/* meta analysis module by Evangelos Kontopantelis
created in STATA v9.2, 26 October 2009
v1.1 10 Jan 2013
-set variables to double for better precision when study variance is extremely small
v1.2 25 Feb 2013
-added bootstrap method (default is 1,000 iterations)
-added sensitivity analysis using pre-set value for I^2
v1.3 16 Oct 2013
-renamed diff scalars since overlap with variables is likely
v1.4 14 Aug 2014
-added option to deal with proportions - Freeman-Tukey transformation
v1.5 05 Sep 2014
-corrected strange STATA error that suppressed heterogeneity outputs: must have been using variable name recently reserved by Stata (tausq)
-changed reported H^2 from H^2M (i.e. H^2 - 1) to H^2 to avoid confusion
-added support for alpha other than 5%
-addded test based CIs for heterogeneity measures
v1.6 18 Sep 2014
-added support for back-transforming effects for binary outcomes to ORs
v1.7 20 Jun 2015
-CIs for I^2 and H^2 added to returned scalars
v2.0 12 Apr 2016
-major update for forest plot to use program _dipgby by Ross Harris & Mike Bradburn
v2.1 12 Apr 2016
-seed number no longer set automatically
v3.0 11 Apr 2017
-by option added for subgroup analyses
v4.0 25 Jul 2018
-added support for binary data (MH, IV or Peto)
v4.1 10 Sep 2018
-added back-transformation for empirical logit
-gave up moving REML code to main ado file (Stata does not allow this for ML algrithms)
-added pscl option for [0,1] or [0,100] scaling of percentages (used with prp or backt options only)
v4.2 17 Sep 2019
-corrected to allow additional twway graph options
v4.3 28 Jan 2020
-four variable syntax now asks for events and non-events, rather than events and totals
*/
/*stata version define*/
version 9
/*metaan calculates the effect and SE for each study (if possible)*/
program define metaan, rclass
/*command syntax*/
syntax namelist(min=2 max=4) [if] [in], /*
*/ [fe dl bdl ml reml pl pe sa exp varc prp label(string) forest /*
*/ reps(integer -127) seed(integer -127) sens(real 80) plplot(string) /*
*/ grpby(varname numeric) backt(string) pscl /*
*/ /*MH and Peto methods*/ mhor mhrr mhrd por poe /*
*/ /*new plot options*/ DP(integer 2) ASTEXT(integer 50) TEXTSize(real 100.0) /*
*/ BOXSCA(real 100.0) noOVERALL NOHET NULL(real 999) NULLOFF NOWARNING XLAbel(passthru) /*
*/ XTick(passthru) FORCE SUMMARYONLY EFFECT(string) FAVOURS(string) /*
*/ DOUBLE BOXOPT(string) CLASSIC DIAMOPT(string) POINTOPT(string) CIOPT(string) /*
*/ OLINEOPT(string) noSTATS noWT lcols(string) rcols(string) /*
*/ /*junk now*/ forestw(real -127.5) * ]
/*temp variables used in all methods*/
tempvar temp1 temp2 temp3 temp4
/*temp scalars*/
tempname btor
/*INITIAL STUFF*/
//get level
local clvl=c(level)
local mltps=(100-`clvl')/200
local mltpl=invnormal(1-(100-`clvl')/200)
//methods
local mlist1 "fe dl bdl ml pl reml pe sa poe"
/*find out if the variables specified are there and include the information they are supposed to have*/
local tswc=wordcount("`namelist'")
forvalues i=1(1)`tswc' {
local var`i' = word("`namelist'",`i')
}
if `tswc'==3 {
di as error "Output variables need to be two or four (see help file)"
error 110
}
/*make sure a user doesn't give the same variable names*/
forvalues i=1(1)`tswc' {
forvalues j=1(1)`tswc' {
if `i'!=`j' & "`var`i''"=="`var`j''" {
di as error "Outcome variables need to have different names!"
error 110
}
}
}
//find out the variable specified is there and includes the information they are supposed to have
forvalues i=1(1)`tswc'{
capture confirm numeric variable `var`i'', exact
if _rc!=0 {
di as error "Variable `var`i'' does not exist or is not numeric"
error 110
}
}
/*check what is happening with model options*/
local mdcnt=0
forvalues i=1(1)`=wordcount("`mlist1'")' {
if "``=word("`mlist1'",`i')''"!="" local mdcnt = `mdcnt' + 1
}
if `mdcnt'>=2 {
di as error "Please select only one model option"
error 110
}
//if provided effect and its variability
if `tswc'==2 {
if `mdcnt' == 0 {
di as error "You must select a model option: fe, dl, bdl, ml, pl, reml, pe ,sa, poe"
error 110
}
//variables that need to be positive
qui sum `var2' `if' `in'
if r(min)<=0 {
di as error "Variable `var2' must only contain positive numbers"
error 110
}
/*prp and varc*/
if "`prp'"!="" {
if "`varc'"!="" {
di as error "Cannot use varc and prp options together: either standard meta-analysis or one of proportions!"
error 110
}
}
//exponentiate or not
scalar `btor'=0
if "`exp'"!="" {
scalar `btor'=1
if "`prp'"!="" {
di as error "Cannot use exp and prp options together: either meta-analysis of dichotomous outcomes or one of proportions!"
error 110
}
}
//additional checks and info if MA of proportions
scalar prpinfo=0
if "`prp'"!="" {
//make sure both variables include integers and positive
qui sum `var1'
scalar minvar1=r(min)
if r(min)<0 {
di as error "Variable `var1' (numerators) must only contain positive or zero numbers"
error 110
}
qui count if mod(`var1', 1)!=0 & `var1'!=.
if r(N)>0 {
di as error "Variable `var1' (numerators) must only contain integers"
error 110
}
qui count if mod(`var2', 1)!=0 & `var2'!=.
if r(N)>0 {
di as error "Variable `var2' (denominators) must only contain integers"
error 110
}
//make sure denominators are not smaller than numerators
qui gen `temp1'=`var2'-`var1'
qui sum `temp1'
if r(min)<0 {
di as error "Numerators appear larger than denominators in some cases"
error 110
}
qui drop `temp1'
//to pass to forest plot for range
scalar prpinfo=1
}
//cannot have Peto or MH
if "`mhor'"!="" | "`mhrr'"!="" | "`mhrd'"!="" | "`por'"!="" {
di as error "Cannot select Peto or Mantel-Haenszel when providing an effect and its variance"
error 110
}
/*Back-transformation options - logit (need percentage) or empirical logit (need percentage and denominator)*/
local btcntr=0
if "`backt'"!="" {
if "`prp'"!="" {
di as error "Options backt and prp cannot be requested together"
error 197
}
if "`exp'"!="" {
di as error "Options backt and exp cannot be requested together"
error 197
}
tokenize "`backt'", parse(" ,")
if "`3'"!="" {
di as error "At most two variable names in backt() for empirical logit: percentages and denominators"
error 197
}
while "`1'"!="" {
cap confirm var `1'
if _rc!=0 {
di as err "Back-transformation variable `1' not defined"
error 109
}
/*store variable names in local variables for future use*/
local btcntr = `btcntr'+1
local btvar`btcntr' = "`1'"
macro shift
}
//make sure first variable is percentages and second denominators
forvalues i=1(1)`btcntr' {
capture confirm numeric variable `btvar`i'', exact
if _rc!=0 {
di as error "Variable `btvar`i'' does not exist or is not numeric"
error 109
}
if `i'==1 {
qui sum `btvar`i''
if r(min)<0 | r(max)>1 {
di as error "Variable `btvar`i'' holds percentages and needs to be in the [0,1] range"
error 109
}
//if only one variable provided (percentages) it cannot have any zeros
if `btcntr'==1 & (r(min)==0 | r(max)==1) {
di as error "Logit back-transformation cannot be applied to variable `btvar`i'' (percentages) since it holds zeros and/or ones"
di as error "Include values in the (0,1) range or provide denominators so that the empirical logit back-transfomration can be applied"
error 109
}
}
if `i'==2 {
qui count if mod(`btvar`i'', 1)!=0
local bttmp=r(N)
qui sum `btvar`i''
if r(min)<=0 | `bttmp'>0 {
di as error "Variable `btvar`i'' holds denominators (positive integers)"
error 109
}
}
}
}
//pscl option can only be used with prp or backt
if "`pscl'"!="" {
if "`backt'"=="" & "`prp'"=="" {
di as error "Option pscl can only be used with backt or prp options"
error 197
}
}
}
//if provided events and non-events variables
if `tswc'==4 {
scalar prpinfo=0
//calculate denominators to comply with changing input
qui gen `temp1'=`var1'+`var2'
qui gen `temp2'=`var3'+`var4'
//exponentiate or not
scalar `btor'=0
if "`exp'"!="" {
scalar `btor'=1
if "`mhrd'"!="" {
di as error "Cannot use exp option with Mantel-Haenszel Risk Difference (option mhrd)"
error 110
}
}
//events
foreach i in 1 3 {
qui sum `var`i''
if r(min)<0 {
di as error "Variable `var`i'' must only contain non-negative numbers"
error 110
}
}
//denominators
foreach i in 1 2 {
qui sum `temp`i''
if r(min)<=0 {
di as error "Events and non-events for group `i' must be positive"
error 110
}
}
//options not allowed
foreach x in prp varc backt {
if "``x''"!="" {
di as error "Cannot use option `x' with events and populations variables"
error 110
}
}
if "`poe'"!="" {
di as error "Cannot use Peto O-E option (poe) with the four variable syntax - needs HR or OR estimates"
error 110
}
//must have one of Peto or MH
local cnt = 0
foreach x in por mhor mhrr mhrd {
if "``x''"!="" local cnt = `cnt' + 1
}
if `cnt' == 0 {
di as error "You must select a calculation option when using events and populations: por, mhor, mhrr, mhrd"
error 110
}
else if `cnt'>=2 {
di as error "Please select only one calculation option (por, mhor, mhrr, mhrd)"
error 110
}
qui drop `temp1' `temp2'
}
//if provided effect and its variability OR a different model has been selected for event data
if "`tswc'"=="2" | ("`tswc'"=="4" & `mdcnt'>0) {
/*if profile likelihood plot is selected*/
if "`plplot'"!="" {
if "`ml'"=="" & "`pl'"=="" & "`reml'"=="" {
di as error "Option plplot can only be used with the maximum-likelihood, profile-likelihood & restricted maximum-likelihood methods"
exit
}
tokenize "`plplot'", parse(" ,")
if "`2'"!="" | ("`1'"!="mu" & "`1'"!="tsq") {
di as error "Please use plplot with either the mu or tsq option i.e. plplot(mu) or plplot(tsq)"
exit
}
local plplotvar = "`1'"
}
/*bootsrap options*/
if `reps'!=-127 & "`bdl'"=="" {
di as error "Bootsrap repetitions option only required for the Bootstrapped DerSimonian-Laird method"
error 110
}
if "`bdl'"!="" {
/*set default to 1000*/
if `reps'==-127 local reps=1000
/*error if few repetitions chosen*/
if `reps'<100 {
di as error "At least 100 repetitions are recommended"
error 110
}
}
/*sensitivity analysis*/
if "`sa'"!="" {
if `sens'<0 | `sens'>=100 {
di as error "Sensitivity analyses performed using an I^2 value and constrained to [0,100)"
error 110
}
}
else {
if `sens'!=80 {
di as error "sens() option can only be used when requesting a sensitivity analysis (sa) model"
error 110
}
}
}
//common in both 2 or 4 input variables
/*to select one of the plot options*/
local cnt=0
foreach x in forest plplot {
if "``x''"!="" local cnt = `cnt' + 1
}
if `forestw'!=-127.5 local cnt = `cnt' + 1
/*if more than one plot options used come back with an error*/
if `cnt'>1 {
di as error "Please select only one of the available plot options"
error 110
}
/*set seed*/
if `seed'!=-127 {
qui set seed `seed'
}
/*by option*/
if "`grpby'" !="" {
qui count if mod(`grpby', 1)!=0
if r(N)>0 {
di as error "by option variable needs to be numeric and include integers only"
exit
}
}
/*label options - title(authors) or year or both*/
local lcntr=0
if "`label'"!="" {
tokenize "`label'", parse(" ,")
if "`3'"!="" {
di as error "At most two variable names in label(): e.g ...label(authors year)"
exit
}
while "`1'"!="" {
cap confirm var `1'
if _rc!=0 {
di as err "Label variable `1' not defined"
exit
}
/*store variable names in local variables for future use*/
local lcntr = `lcntr'+1
local lvar`lcntr' = "`1'"
macro shift
}
}
//NEW FOREST PLOT OPTIONS
global MA_dp = "`dp'"
global MA_AS_TEXT `astext'
global MA_TEXT_SCA `textsize'
global MA_FBSC `boxsca'
global MA_nulloff "`nulloff'"
global MA_nohet "`nohet'"
global MA_nowarning "`nowarning'"
if `null'!=999 {
global MA_NULL=`null'
}
else {
//set default null if none provided
if `btor'==1{
global MA_NULL=1
}
else {
global MA_NULL=0
}
}
global MA_rfdist "" /*future trial option not used*/
global RFL ""
local xcounts "" /*raw counts display not used*/
global MA_G1L ""
global MA_G2L ""
local log "" /*issued with raw counts so not using*/
//if meta-analysis of proportions
if prpinfo==1 {
if "`xlabel'"=="" {
if "`pscl'"=="" {
local xlabel "xlabel(0,0.2,0.4,0.6,0.8,1)"
}
else {
local xlabel "xlabel(0,20,40,60,80,100)"
}
}
local force="force"
global MA_nulloff "nulloff"
}
global MA_rjhby "`grpby'" /*BY option*/
global rjhHetGrp ""
global MA_summaryonly "`summaryonly'"
global MA_params = 0 /*irrelevant - number of parameters as input in metan*/
global IND=`clvl'
global MA_ESLA "`effect'"
global MA_FAVOURS "`favours'"
if "`effect'"=="" {
local sumstat "ES"
if `btor'==1 {
local sumstat "OR"
if "`mhrr'"!="" {
local sumstat "RR"
}
else if "`mhrd'"!="" {
local sumstat "RD"
}
else if "`poe'"!="" {
local sumstat "HR"
}
}
}
else {
local sumstat "`effect'"
}
global MA_firststats ""
global MA_secondstats ""
global MA_userDescM "" /*irrelevant - linked to first() option in metan*/
global MA_userDesc ""
global MA_DOUBLE "`double'"
global MA_efficacy "" /*irrelevant - linked to efficacy() option in metan*/
global MA_OTHEROPTS `"`options'"' /*graph options added - _dispgby can identify them*/
global MA_BOXOPT `"`boxopt'"'
global MA_classic "`classic'"
global MA_DIAMOPT `"`diamopt'"'
global MA_POINTOPT `"`pointopt'"'
global MA_CIOPT `"`ciopt'"'
global MA_OLINEOPT `"`olineopt'"'
/*use temporary file to mess with the data as much as needed (can't use tempvar because need to call in many programs)*/
preserve
/*in and if options - create temp variable that will enable thir usage*/
qui capture drop use
qui gen use = 0
qui replace use = 1 `if' `in'
qui drop if use==0
//if provided use effect and its variance
if `tswc'==2 {
/*rename to avoid using the same variable names*/
capture drop tempeff_007
rename `var1' tempeff_007
capture drop tempeff_007var
rename `var2' tempeff_007var
capture drop eff
capture drop effvar
//different approach for standard MA and one of proportions
if prpinfo==0 {
/*varc variable option and changes*/
qui gen double eff = tempeff_007
if "`varc'"!="" {
qui gen double effvar = tempeff_007var
}
else {
qui gen double effvar = tempeff_007var^2
}
}
else {
//if proportions transform (tempeff_007=num, tempeff_007var=den)
/*effect & SE*/
qui gen double eff = asin(sqrt(tempeff_007/(tempeff_007var+1))) + asin(sqrt((tempeff_007+1)/(tempeff_007var+1)))
qui gen double effvar = (1/(tempeff_007var+1))
}
/*exclude studies where the eff or the SEeff are missing*/
qui replace use=0 if eff==. | effvar==.
qui drop if use==0
}
//if events and denominator calculate
else {
/*rename to avoid using the same variable names*/
capture drop ev1_127
rename `var1' ev1_127
capture drop total1_127
qui gen total1_127=ev1_127+`var2'
capture drop ev2_127
rename `var3' ev2_127
capture drop total2_127
qui gen total2_127=ev2_127+`var4'
/*exclude studies where the eff or the SEeff are missing*/
qui replace use=0 if ev1_127==. | ev2_127==. | total1_127==. | total2_127==.
qui drop if use==0
capture drop eff
capture drop effvar
capture drop weights
//call selected weighting program to return eff and eff var
if "`mhor'"!="" {
calc_MHORf "gen"
}
else if "`mhrr'"!="" {
calc_MHRRf "gen"
}
else if "`mhrd'"!="" {
calc_MHRDf "gen"
}
else if "`por'"!="" {
calc_PORf "gen"
}
*restore, not
*error
}
/*generate and id for the studies*/
qui capture qui drop studyid
qui egen studyid = seq()
/*create label temp variable*/
/*first convert if needed to string variables*/
forvalues i=1(1)`lcntr' {
qui capture confirm string var `lvar`i''
if _rc!=0 {
qui capture drop temp1
qui gen temp1 = string(`lvar`i'')
qui drop `lvar`i''
qui rename temp1 `lvar`i''
}
}
/*then merge*/
qui capture drop labelvar
if `lcntr'==0 {
qui gen str20 labelvar = string(studyid)
}
else if `lcntr'==1 {
qui gen str20 labelvar = trim(`lvar1')
}
else {
qui gen str20 labelvar = trim(`lvar1') + ", " + trim(`lvar2')
}
/*lower and upper CI for studies*/
qui capture drop studyloCI
qui capture drop studyupCI
qui gen studyloCI = eff - `mltpl'*sqrt(effvar)
qui gen studyupCI = eff + `mltpl'*sqrt(effvar)
//loop to accomodate by option
local numloop=1
if "`grpby'"!="" {
qui tab `grpby'
local numloop=r(r)+1
//mess with variable to make sure it's numbered 1(1)max
local lblfoo : val label `grpby'
qui sum `grpby'
local tmpmin=r(min)
local tmpmax=r(max)
local cntr=0
qui gen tgrpby=.
forvalues i=`tmpmin'(1)`tmpmax' {
qui count if `grpby'==`i'
if r(N)>0 {
local cntr=`cntr'+1
qui replace tgrpby=`cntr' if `grpby'==`i'
if "`lblfoo'"=="" {
local tfoo`cntr' "`i'"
}
else {
local tfoo`cntr' : label `lblfoo' `i'
}
label define tgrpbylbl `cntr' "`tfoo`cntr''", add
}
}
label val tgrpby tgrpbylbl
sort tgrpby
}
tempfile master1
qui save `master1', replace
forvalues xlp=1(1)`numloop' {
qui use `master1', clear
if `numloop'>1 {
if `xlp'<`numloop' {
qui keep if tgrpby==`xlp'
di _newline(2) as text "Subgroup analysis: `tfoo`xlp''"
}
else {
di _newline(2) as text "Overall analysis"
}
}
/*CALCULATE ALL METHODS*/
/*number of used studies*/
qui count
scalar k=r(N)
/*********************/
/*Fixed effects model*/
qui gen double `temp1' = eff/effvar
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = 1/effvar
qui sum `temp2'
scalar sum2 = r(sum)
/*mean & var estimates*/
scalar fe_mu = sum1/sum2
scalar fe_var = 1/sum2
/*weights - to be used in display and graphs - different for each method*/
qui capture drop weights
qui gen double weights = `temp2'/sum2
qui drop `temp1' `temp2'
scalar fe_lo = fe_mu - `mltpl'*sqrt(fe_var)
scalar fe_up = fe_mu + `mltpl'*sqrt(fe_var)
/*************************/
/*DerSimonian-Laird model*/
/*calculate Cochran's Qw*/
qui gen double `temp1' = (1/effvar)*(eff-fe_mu)^2
qui sum `temp1'
scalar qw = r(sum)
qui drop `temp1'
/*calculate t^2 estimate*/
qui gen double `temp1' = 1/effvar
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = (1/effvar)^2
qui sum `temp2'
scalar sum2 = r(sum)
/* t^2 can't be negative!*/
scalar tsq_dl = max((qw-(k-1))/(sum1-sum2/sum1),0)
/*calculate the s^2 that will be used in the heterogeneity measures*/
scalar ssq = (k-1)*sum1/(sum1^2-sum2)
qui drop `temp1' `temp2'
/*DL model*/
qui gen double `temp1' = eff/(effvar+tsq_dl)
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = 1/(effvar+tsq_dl)
qui sum `temp2'
scalar sum2 = r(sum)
/*weights - to be used in display and graphs - different for each method*/
if "`dl'"!="" | "`pe'"!="" {
qui capture drop weights
qui gen weights = `temp2'/sum2
}
/*mean & var estimates*/
scalar dl_mu = sum1/sum2
scalar dl_var = 1/sum2
scalar dl_lo = dl_mu - `mltpl'*sqrt(dl_var)
scalar dl_up = dl_mu + `mltpl'*sqrt(dl_var)
qui drop `temp1' `temp2'
/****************/
/*T-test "model"*/
/*only use effect sizes and not effect variances*/
qui sum eff
scalar tt_mu = r(mean)
scalar tt_var = r(Var)
/*calculate CIs*/
scalar tt_lo = tt_mu - invttail(k-1,`mltps')*sqrt(tt_var/k)
scalar tt_up = tt_mu + invttail(k-1,`mltps')*sqrt(tt_var/k)
/*********************/
/*events methods**/
if `tswc'==4 & `mdcnt'==0 {
if "`por'"!="" {
calc_PORf "gen"
scalar peto_mu = r(mu)
scalar peto_var = r(var)
scalar peto_lo = r(lo)
scalar peto_up = r(up)
}
else {
if "`mhor'"!="" {
calc_MHORf "gen"
}
else if "`mhrr'"!="" {
calc_MHRRf "gen"
}
else if "`mhrd'"!="" {
calc_MHRDf "gen"
}
scalar mh_mu = r(mu)
scalar mh_var = r(var)
scalar mh_lo = r(lo)
scalar mh_up = r(up)
}
}
/*********************/
/*O-E Peto method**/
if "`poe'"!="" {
calc_PORf "partgen" "effex"
scalar poe_mu = r(mu)
scalar poe_var = r(var)
scalar poe_lo = r(lo)
scalar poe_up = r(up)
}
/*************************/
/*DerSimonian-Laird bootstrap model*/
/*only if requested*/
if "`bdl'"!="" {
tempvar teffvar teff
forvalues i=1(1)`reps' {
qui gen `teffvar'=.
qui gen `teff'=.
/*NEED A RANDOM SAMPLE WITH REPLACEMENT*/
forvalues j=1(1)`=k' {
local tnum = 1+int((`=k')*runiform())
qui replace `teff'=eff[`tnum'] in `j'
qui replace `teffvar'=effvar[`tnum'] in `j'
}
/*CALCULATE COCHRAN'S Qw */
qui gen double `temp1' = (1/`teffvar')*(`teff'-fe_mu)^2
qui sum `temp1'
scalar _qw = r(sum)
qui drop `temp1'
/*CALCULATE ESTIMATE of t^2*/
qui gen double `temp1' = 1/`teffvar'
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = (1/`teffvar')^2
qui sum `temp2'
scalar sum2 = r(sum)
/* t^2 can't be negative!*/
scalar _tsq_dl = max((_qw-(k-1))/(sum1-sum2/sum1),0)
/*calculate the s^2 that will be used in the heterogeneity measures*/
scalar _ssq = (k-1)*sum1/(sum1^2-sum2)
/*add to MATA matrices*/
if `i'==1 {
qui mata: A = st_numscalar("_tsq_dl")
qui mata: B = st_numscalar("_ssq")
qui mata: C = st_numscalar("_qw")
*qui mata: eff=st_data(.,"`teff'")
*qui mata: effvar=st_data(.,"`teffvar'")
}
else {
qui mata: A = (A \ st_numscalar("_tsq_dl"))
qui mata: B = (B \ st_numscalar("_ssq"))
qui mata: C = (C \ st_numscalar("_qw"))
*qui mata: eff=(eff, st_data(.,"`teff'"))
*qui mata: effvar=(effvar, st_data(.,"`teffvar'"))
}
*list `teff' `teffvar'
qui drop `temp1' `temp2' `teff' `teffvar'
}
/*use MATA matrices to calculate means*/
mata: meanA = mean(A)
mata: meanB = mean(B)
mata: meanC = mean(C)
mata: st_numscalar("tsq_bdl", meanA)
mata: st_numscalar("ssq", meanB)
mata: st_numscalar("qw", meanC)
/*DL model*/
qui gen double `temp1' = eff/(effvar+tsq_bdl)
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = 1/(effvar+tsq_bdl)
qui sum `temp2'
scalar sum2 = r(sum)
/*weights - to be used in display and graphs - different for each method*/
qui capture drop weights
qui gen weights = `temp2'/sum2
/*mean & var estimates*/
scalar bdl_mu = sum1/sum2
scalar bdl_var = 1/sum2
scalar bdl_lo = bdl_mu - `mltpl'*sqrt(bdl_var)
scalar bdl_up = bdl_mu + `mltpl'*sqrt(bdl_var)
qui drop `temp1' `temp2'
}
/*************************/
/*Sensitivity analysis*/
if "`sa'"!="" {
/*get I^2 and H^2 values*/
scalar Isq=`sens'
scalar Hsq = Isq/(100-Isq)+1
/*calculate ssq estimate*/
qui gen double `temp1' = 1/effvar
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = (1/effvar)^2
qui sum `temp2'
scalar sum2 = r(sum)
scalar ssq = (k-1)*sum1/(sum1^2-sum2)
qui drop `temp1' `temp2'
/*calculate the tau^2 estimate from input and ssq*/
scalar tsq_sa = (Hsq)*ssq-ssq
/*calculate Cochran's Qw - not reported under sensitivity analysis*/
scalar qw = .
/*DL model*/
qui gen double `temp1' = eff/(effvar+tsq_sa)
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = 1/(effvar+tsq_sa)
qui sum `temp2'
scalar sum2 = r(sum)
/*weights - to be used in display and graphs - different for each method*/
qui capture drop weights
qui gen weights = `temp2'/sum2
/*mean & var estimates*/
scalar sa_mu = sum1/sum2
scalar sa_var = 1/sum2
scalar sa_lo = sa_mu - `mltpl'*sqrt(sa_var)
scalar sa_up = sa_mu + `mltpl'*sqrt(sa_var)
qui drop `temp1' `temp2'
}
/**************************/
/*Maximum Likelihood model*/
/*calculate ML only if asked for (or if PL is asked) - since it takes time*/
if "`ml'"!="" | "`pl'"!="" {
/*initialise t^2 estimate*/
scalar tsq_ml = max(tsq_dl, 1e-4)
/*calculate initial mu_m (tau estimate)*/
qui gen double `temp1' = eff/(effvar+tsq_ml)
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = 1/(effvar+tsq_ml)
qui sum `temp2'
scalar sum2 = r(sum)
scalar mu_m = sum1/sum2
qui drop `temp1' `temp2'
/*dynamic system: we want it to converge to a fixed point - 500 iterations limit*/
scalar nz=0
scalar diff188_1 = 1
scalar diff188_2 = 1
/*continue with loop till both converge or we reach 500 iterations*/
while (abs(diff188_1)>10^-6 | abs(diff188_2)>10^-6) & nz<500 {
scalar nz = nz + 1
/*calculate new tsq_ml (t^2 estimate)*/
qui gen double `temp1' = ((eff-mu_m)^2-effvar)/(effvar+tsq_ml)^2
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = 1/(effvar+tsq_ml)^2
qui sum `temp2'
scalar sum2 = r(sum)
/*calculate new mu_m (tau estimate)*/
qui gen double `temp3' = eff/(effvar+tsq_ml)
qui sum `temp3'
scalar sum3 = r(sum)
qui gen double `temp4' = 1/(effvar+tsq_ml)
qui sum `temp4'
scalar sum4 = r(sum)
/*calc the differences*/
scalar diff188_1 = tsq_ml - sum1/sum2
scalar diff188_2 = mu_m - sum3/sum4
/*calculate the new values*/
scalar tsq_ml = sum1/sum2
scalar mu_m = sum3/sum4
qui drop `temp1' `temp2' `temp3' `temp4'
}
/*save the results to the appropriate variables - choose accordingly if the estimate is negative or not*/
if tsq_ml>=0 {
scalar ml_mu = mu_m
qui gen double `temp1' = 1/(effvar+tsq_ml)
qui sum `temp1'
scalar ml_var = 1/r(sum)
scalar sum2 = r(sum)
/*weights - to be used in display and graphs - different for each method*/
qui capture drop weights
qui gen double weights = `temp1'/sum2
qui drop `temp1'
scalar ml_info=1
/*if it has not converged update information*/
if abs(diff188_1)>10^-6 | abs(diff188_2)>10^-6 {
scalar ml_info = 0
scalar ml_mu = .
scalar ml_var = .
scalar tsq_ml = .
qui replace weights=.
}
}
else {
scalar ml_mu = fe_mu
scalar ml_var = fe_var
scalar tsq_ml = 0
scalar ml_info = -1
}
/*confidence intervals*/
scalar ml_lo = ml_mu - `mltpl'*sqrt(ml_var)
scalar ml_up = ml_mu + `mltpl'*sqrt(ml_var)
/*calculate the log-likelihood value*/
qui gen double `temp1' = -1/2*ln(2*_pi*(effvar+tsq_ml))-1/2*(eff-ml_mu)^2/(effvar+tsq_ml)
qui sum `temp1'
scalar ml_val = r(sum)
qui drop `temp1'
}
/*calculate PL only if asked for and ML converged - since it takes time*/
if "`pl'"!="" {
if ml_info!=0 {
/********************************************************************/
/*Profile Likelihood method*/
/*ML estimates from previous step - PL estimates are the same only the CI is different*/
scalar pl_mu = ml_mu /*mean effect*/
scalar pl_var = ml_var /*effect variance*/
scalar tsq_pl = tsq_ml
scalar Tm = tsq_ml /*between study variance - only used in CI calculation*/
/*main method calculations: CI for the effect estimate*/
/*Lower bound calculations*/
scalar Mmlo = pl_mu - 20 /*using SE wasn't efficient - many SEs close to 0*/
/*binary search - stop if within 10^-6*/
scalar st_0 = Mmlo
scalar fin_0 = pl_mu
scalar diff188_ = 1
scalar cntr = 0
while (diff188_<0 | diff188_>10^-6) & cntr<500 {
local tempxxx = (st_0+fin_0)/2
scalar cntr = cntr + 1
/*call program that calculates the ML value for fixed tau*/
qui PL_value `tempxxx'
scalar diff188_ = r(mlv) - ml_val + 1.92
if diff188_<0 { /*move upwards*/
scalar st_0 = `tempxxx'
}
if diff188_>10^-6 { /*move downwards*/
scalar fin_0 = `tempxxx'
}
}
scalar pl_lo = `tempxxx'
/*if it hasn't converged update indicator*/
scalar pl_lo_info = 1
if cntr==500 {
scalar pl_lo_info = 0
scalar pl_lo = .
}
/*Upper bound calculations*/
scalar Mmup = pl_mu + 20 /*using SE wasn't efficient - many SEs close to 0*/
/*binary search - stop if within 10^-6*/
scalar st_0 = pl_mu
scalar fin_0 = Mmup
scalar diff188_ = 1
scalar cntr = 0
while (diff188_<0 | diff188_>10^-6) & cntr<500 {
scalar cntr = cntr + 1
local tempxxx = (st_0+fin_0)/2
/*call program that calculates the ML value for fixed tau*/
qui PL_value `tempxxx'
scalar diff188_ = r(mlv) - ml_val + 1.92
if diff188_>10^-6 { /*move upwards*/
scalar st_0 = `tempxxx'
}
if diff188_<0 { /*move downwards*/
scalar fin_0 = `tempxxx'
}
}
scalar pl_up = `tempxxx'
/*if it hasn't converged update indicator*/
scalar pl_up_info = 1
if cntr==500 {
scalar pl_up_info = 0
scalar pl_up = .
}
/*heterogeneity calculations: CI for the t^2 estimate*/
/*Lower bound calculations*/
scalar Tmlo = 0
/*binary search - stop if within 10^-6*/
scalar st_0 = Tmlo
scalar fin_0 = Tm
scalar diff188_ = 1
scalar cntr = 0
while (diff188_<0 | diff188_>10^-6) & cntr<500 {
local tempxxx = (st_0+fin_0)/2
scalar cntr = cntr + 1
/*call program that calculates the ML value for fixed t^2*/
qui PL_value2 `tempxxx'
scalar diff188_ = r(mlv) - ml_val + 1.92
if diff188_<0 { /*move upwards*/
scalar st_0 = `tempxxx'
}
if diff188_>10^-6 { /*move downwards*/
scalar fin_0 = `tempxxx'
}
}
/*take into account lower limit very very close to zero*/
scalar tsq_pl_lo = `tempxxx'
if `tempxxx' < 10^8 scalar tsq_pl_lo=0
/*if it hasn't converged update indicator - but take into account values very close to zero*/
scalar tsq_pl_lo_info = 1
if cntr==500 & tsq_pl_lo!=0{
scalar tsq_pl_lo_info = 0
scalar tsq_pl_lo = .
}
/*Upper bound calculations*/
scalar Tmup = Tm + 10
/*binary search - stop if within 10^-6*/
scalar st_0 = Tm
scalar fin_0 = Tmup
scalar diff188_ = 1
scalar cntr = 0
while (diff188_<0 | diff188_>10^-6) & cntr<500{
scalar cntr = cntr + 1
local tempxxx = (st_0+fin_0)/2
/*call program that calculates the ML value for fixed t^2*/
qui PL_value2 `tempxxx'
scalar diff188_ = r(mlv) - ml_val + 1.92
if diff188_>10^-6 { /*move upwards*/
scalar st_0 = `tempxxx'
}
if diff188_<0 { /*move downwards*/
scalar fin_0 = `tempxxx'
}
}
scalar tsq_pl_up = `tempxxx'
scalar tsq_pl_up_info = 1
/*if it hasn't converged update indicator*/
if cntr==500 {
scalar tsq_pl_up_info = 0
scalar tsq_pl_up = .
}
}
/*if ML has not converged then PL is not executed*/
else {
scalar pl_mu = .
scalar pl_var = .
scalar pl_lo_info = .
scalar pl_lo = .
scalar pl_up_info = .
scalar pl_up = .
scalar tsq_pl_lo_info = .
scalar tsq_pl_lo = .
scalar tsq_pl_up_info = .
scalar tsq_pl_up = .
scalar tsq_pl = .
}
}
/**************************/
/*Restricted Maximum Likelihood model*/
/*calculate REML only if asked for - since it may take time*/
if "`reml'"!="" {
scalar tsq_reml = max(tsq_dl, 1e-4)
/* estimate tau2 using ml */
qui ml model d0 metaan_reml (eff effvar = ), maximize init(_cons=`=tsq_reml') search(off) nopreserve noscvars /*
*/ log iterate(100)
scalar remlconvinfo = e(converged)
scalar tsq_reml = max(_b[_cons],0)
/*save the results to the appropriate variables - choose accordingly if the estimate is negative or not*/
if tsq_reml>=0 {
qui gen double `temp1' = eff/(effvar+tsq_reml)
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = 1/(effvar+tsq_reml)
qui sum `temp2'
scalar sum2 = r(sum)
/*weights - to be used in display and graphs - different for each method*/
qui capture drop weights
qui gen double weights = `temp2'/sum2
/*mean & var estimates*/
scalar reml_mu = sum1/sum2
scalar reml_var = 1/sum2
qui drop `temp1' `temp2'
scalar reml_info=1
/*if it has not converged update information*/
if remlconvinfo!=1 {
scalar reml_info = 0
scalar reml_mu = .
scalar reml_var = .
scalar tsq_reml = .
/*qui replace weights=.*/
}
}
else {
scalar reml_mu = fe_mu
scalar reml_var = fe_var
scalar tsq_reml = 0
scalar reml_info = -1
}
/*confidence intervals*/
scalar reml_lo = reml_mu - `mltpl'*sqrt(reml_var)
scalar reml_up = reml_mu + `mltpl'*sqrt(reml_var)
}
/*calculate PE only if asked for - since it takes time*/
if "`pe'"!="" {
/*Permutations method (Follmann, 1999)*/
/*drop if there and re-generate*/
scalar pe_pval=.
scalar pe_mu = .
scalar pe_var = .
scalar pe_lo=.
scalar pe_up=.
/*the method does not work for k<6*/
if k>=6 {
/*create locals from scalars to manage calls...*/
local tk = k
local dlmu = dl_mu
/*create appropriate permutation (or random) matrix in memory. for k<=10 use permutations; above that random (2^k rows)*/
if k<=10 {
mata: permtable = permmat(`tk')
}
else {
mata: permtable = randmat(`tk',11)
}
/*get the number of rows of permtable to speed up code*/
mata: rnum = rows(permtable)
/*create an empty m table that will be fed into the main mata function*/
mata: matempty = emptymat(rnum)
/*call function that will return the variances and means in a matrix*/
mata: varmu = createmat(`tk')
/*call the mata function that does all the work and returns a CI matrix*/
mata: pe_p(`tk',varmu,permtable,matempty,0,`dlmu',`mltps')
/*convert the scalar returned to values in case i of pe_prob*/
scalar pe_pval = r(pe_p)
scalar pe_mu = dl_mu
scalar pe_lo = r(pe_lo)
scalar pe_up = r(pe_up)
/*weights taken care of in DL*/
}
}
/*select the method to be displayed and create appropriate strings*/
if "`fe'"!="" {
local headoutstr = "Fixed-effect method selected"
local modsel = "fe"
}
else if "`dl'"!="" {
local headoutstr = "DerSimonian-Laird random-effects method selected"
local modsel = "dl"
}
else if "`bdl'"!="" {
local headoutstr = "Bootstrapped DerSimonian-Laird random-effects method selected"
local modsel = "bdl"
}
else if "`ml'"!="" {
local headoutstr = "Maximum Likelihood method selected"
local modsel = "ml"
}
else if "`pl'"!="" {
local headoutstr = "Profile Likelihood method selected"
local modsel = "pl"
}
else if "`reml'"!="" {
local headoutstr = "Restricted Maximum Likelihood (REML) method selected"
local modsel = "reml"
}
else if "`pe'"!="" {
local headoutstr = "Permutations method selected (Follmann & Proschan)"
local modsel = "pe"
}
else if "`sa'"!="" {
local headoutstr = "Sensitivity analysis selected with preset heterogeneity (I^2=`sens'%)"
local modsel = "sa"
}
else if "`poe'"!="" {
local headoutstr = "Peto fixed-effect for O-E data selected (log(HR) or log(OR))"
local modsel = "poe"
}
/*if MH/Peto used in conjuction with one of the above*/
if "`mhor'"!="" {
if `mdcnt'==0 {
local headoutstr = "Mantel-Haenszel fixed-effect method (based on Odds-Ratio)"
local modsel = "mh"
}
else {
local headoutstr = "`headoutstr' (MH Odds-Ratio)"
}
}
else if "`mhrr'"!="" {
if `mdcnt'==0 {
local headoutstr = "Mantel-Haenszel fixed-effect method (based on Risk-Ratio)"
local modsel = "mh"
}
else {
local headoutstr = "`headoutstr' (MH Risk-Ratio)"
}
}
else if "`mhrd'"!="" {
if `mdcnt'==0 {
local headoutstr = "Mantel-Haenszel fixed-effect method (based on Risk-Difference)"
local modsel = "mh"
}
else {
local headoutstr = "`headoutstr' (MH Risk-Difference)"
}
}
else if "`por'"!="" {
if `mdcnt'==0 {
local headoutstr = "Peto Odds-Ratio (fixed-effect method)"
local modsel = "peto"
}
else {
local headoutstr = "`headoutstr' (Peto Odds-Ratio)"
}
}
/*information on methods - post texts*/
/*Bootstrap DL on number of repetitions*/
if "`modsel'"=="bdl" {
local posttext1 = "bootstrap of `reps' repetitions"
local ptxttype1 = "text"
}
/*ML information on conergence or not*/
if "`modsel'"=="ml" | "`modsel'"=="pl" {
if ml_info==1 {
local posttext1 = "ML method succesfully converged"
local ptxttype1 = "text"
}
else if ml_info==0 {
local posttext1 = "Warning: ML convergence was not achieved after 500 iterations"
local ptxttype1 = "error"
}
else if ml_info==-1 {
local posttext1 = "Warning: tau^2 estimate using ML was negative and fixed-effect method estimates are used instead"
local ptxttype1 = "error"
}
}
/*PL information on conergence or not*/
if "`modsel'"=="pl" {
if ml_info==0 {
local posttext2 = "PL method not executed since prerequisite ML method did not converge"
local ptxttype2 = "error"
}
else {
if pl_lo_info==1 & pl_up_info==1 {
local posttext2 = "PL method succesfully converged for both upper and lower CI limits"
local ptxttype2 = "text"
}
else {
local posttext2 = "Warning: PL method did not converge for either the lower or upper CI limit of the overall effect, after 500 iterations"
local ptxttype2 = "error"
}
}
}
/*REML information on conergence or not*/
if "`modsel'"=="reml" {
if reml_info==1 {
local posttext3 = "REML method succesfully converged"
local ptxttype3 = "text"
}
else if reml_info==0 {
local posttext3 = "Warning: REML convergence was not achieved after 100 iterations"
local ptxttype3 = "error"
}
else if reml_info==-1 {
local posttext3 = "Warning: tau^2 estimate using REML was negative and fixed-effect method estimates are used instead"
local ptxttype3 = "error"
}
}
/*PE information: below or above 6*/
if "`modsel'"=="pe" {
if k>=6 {
scalar execpe = 1
local posttext4 = "PE method succesfully computed since the number of studies is 6 or above"
local ptxttype4 = "text"
}
else {
scalar execpe = 0
local posttext4 = "Warning: PE method could not be computed since the number of studies is below 6"
local ptxttype4 = "error"
}
}
/*SA information - for Qw and below heterogeneity measures*/
if "`modsel'"=="sa" {
local posttext5 = "Cochrane's Q not reported under a sensitivity analysis"
local ptxttype5 = "text"
}
//MH and Peto information
if "`modsel'"=="mh" | "`modsel'"=="peto" {
local posttext7 = "Heterogeneity measures reported under an inverse-variance weighting assumption"
local ptxttype7 = "text"
}
//if exponentiated
if "`exp'"!="" {
local posttext6 = "Effects have been exponentiated"
local ptxttype6 = "text"
}
//restore, not
//error 110
//back-transform here if meta-analysis of proportions
if prpinfo==1 {
/*overall effect*/
scalar `modsel'_mu = (sin(`modsel'_mu/2))^2
scalar `modsel'_lo = (sin(`modsel'_lo/2))^2
scalar `modsel'_up = (sin(`modsel'_up/2))^2
/*variable back-transform*/
qui replace eff = (sin(eff/2))^2
qui replace studyloCI = (sin(studyloCI/2))^2
qui replace studyupCI = (sin(studyupCI/2))^2
}
//back-transform here if exponentiation requested with `exp' option
if `btor'==1 {
/*overall effect*/
scalar `modsel'_mu = exp(`modsel'_mu)
scalar `modsel'_lo = exp(`modsel'_lo)
scalar `modsel'_up = exp(`modsel'_up)
/*variable back-transform*/
qui replace eff = exp(eff)
qui replace studyloCI = exp(studyloCI)
qui replace studyupCI = exp(studyupCI)
}
//back-transform from logit or empirical logit using
if "`backt'"!="" {
//unweighted or weighted "anchor" for overall effect
if `btcntr'==1 {
qui sum `btvar1'
}
else {
qui sum `btvar1' [fweight=`btvar2']
}
/*overall effect - use the mean percentage to anchor*/
scalar `modsel'_mu = 1/(exp(-`modsel'_mu)*(1-r(mean))/r(mean)+1)-r(mean)
scalar `modsel'_lo = 1/(exp(-`modsel'_lo)*(1-r(mean))/r(mean)+1)-r(mean)
scalar `modsel'_up = 1/(exp(-`modsel'_up)*(1-r(mean))/r(mean)+1)-r(mean)
//if only percentage provided => simple logit for study estimates (no zeros or ones allowed, checked earlier)
if `btcntr'==1 {
/*variable back-transform*/
qui replace eff = 1/(exp(-eff)*(1-`btvar1')/`btvar1'+1)-`btvar1'
qui replace studyloCI = 1/(exp(-studyloCI)*(1-`btvar1')/`btvar1'+1)-`btvar1'
qui replace studyupCI = 1/(exp(-studyupCI)*(1-`btvar1')/`btvar1'+1)-`btvar1'
}
//if denominators are provided => empirical logit for study estimates where zeros or ones available
else {
/*variable back-transform*/
qui replace eff = 1/(exp(-eff)*(1-`btvar1')/`btvar1'+1)-`btvar1' if `btvar1'>0 & `btvar1'<1
qui replace studyloCI = 1/(exp(-studyloCI)*(1-`btvar1')/`btvar1'+1)-`btvar1' if `btvar1'>0 & `btvar1'<1
qui replace studyupCI = 1/(exp(-studyupCI)*(1-`btvar1')/`btvar1'+1)-`btvar1' if `btvar1'>0 & `btvar1'<1
qui replace eff = 1/(exp(-eff)*(1-`btvar1'+0.5/`btvar2')/(`btvar1'+0.5/`btvar2')+1)-`btvar1' if `btvar1'==0 | `btvar1'==1
qui replace studyloCI = 1/(exp(-studyloCI)*(1-`btvar1'+0.5/`btvar2')/(`btvar1'+0.5/`btvar2')+1)-`btvar1' if `btvar1'==0 | `btvar1'==1
qui replace studyupCI = 1/(exp(-studyupCI)*(1-`btvar1'+0.5/`btvar2')/(`btvar1'+0.5/`btvar2')+1)-`btvar1' if `btvar1'==0 | `btvar1'==1
}
}
//if meta-analysis of proportions set to % for forest plot
if "`pscl'"!="" {
foreach x of varlist eff studyloCI studyupCI {
qui replace `x'=`x'*100
}
foreach x in mu lo up {
scalar `modsel'_`x'=`modsel'_`x'*100
}
}
if "`grpby'"=="" {
di as text _newline(2) "`headoutstr'"
}
else {
di as text "`headoutstr'"
}
di as text "{hline 21}{c TT}{hline 45}
di as text "{col 9}Study{col 22}{c |}{col 26}Effect{col 35}[`clvl'% Conf. Interval]{col 57} % Weight"
di as text "{hline 21}{c +}{hline 45}
local stnum = k
forvalues i=1(1)`stnum' {
di as text %-20s labelvar[`i'] "{col 22}{c |}" as result _col(25) %7.3f eff[`i'] _col(36) %7.3f studyloCI[`i'] /*
*/_col(46) %7.3f studyupCI[`i'] _col(57) %7.2f weights[`i']*100
}
di as text "{hline 21}{c +}{hline 45}
qui sum weights
scalar sumweights = 100*r(sum)
di as text %-20s "Overall effect (`modsel'){col 22}{c |}" as result _col(25) %7.3f `modsel'_mu _col(36) %7.3f `modsel'_lo /*
*/_col(46) %7.3f `modsel'_up _col(57) %7.2f sumweights
di as text "{hline 21}{c BT}{hline 45}
/*post texts*/
if "`exp'"!="" {
di as `ptxttype6' "`posttext6'"
}
if "`modsel'"=="ml" | "`modsel'"=="pl" | "`modsel'"=="bdl" {
di as `ptxttype1' "`posttext1'"
}
if "`modsel'"=="pl" {
di as `ptxttype2' "`posttext2'"
}
if "`modsel'"=="reml" {
di as `ptxttype3' "`posttext3'"
}
if "`modsel'"=="pe" {
di as `ptxttype4' "`posttext4'"
}
/*calculate and display heterogeneity measures*/
/*cochran's Q*/
scalar qpval = 1-chi2(k-1,qw)
/*t^2 estimate selection*/
/*the default is DL but will change if ML, PL or boot DL is selected*/
scalar tausqx = tsq_dl
local tsqmdl = "dl"
if "`modsel'"=="bdl" {
scalar tausqx = tsq_bdl
local tsqmdl = "bdl"
}
else if "`modsel'"=="ml" {
scalar tausqx = tsq_ml
local tsqmdl = "ml"
}
else if "`modsel'"=="pl" {
scalar tausqx = tsq_ml
local tsqmdl = "ml"
/*if pl is selected we will use t^2 CIs and convergence information*/
if tsq_pl_lo_info==1 & tsq_pl_up_info==1 {
local posttext5 = "PL method succesfully converged for both upper and lower CI limits of the tau^2 estimate"
local ptxttype5 = "text"
}
else {
local posttext5 = "Warning: PL method did not converge for either the lower or upper CI limit of the tau^2 estimate, after 500 iterations"
local ptxttype5 = "error"
}
}
else if "`modsel'"=="reml" {
scalar tausqx = tsq_reml
local tsqmdl = "reml"
}
else if "`modsel'"=="sa" {
scalar tausqx = tsq_sa
local tsqmdl = "sa"
}
/*I^2 and H^2 - from Higgins 2002 paper */
/*scalar Hsq=(qw-(k-1))/(k-1) relies on DL - not used*/
scalar Hsq = ((tausqx + ssq)/ssq) /*this is H^2 as described by Mittlboeck*/
if Hsq<1 scalar Hsq=1
/*scalar Isq=100*(qw-(k-1))/qw relies on DL method - not used*/
scalar Isq = 100*(Hsq-1)/(Hsq)
if Isq<0 scalar Isq=0
//calculate CIs for H^2 and I^2 using the test-based confidence interval
tempname Qart selnh1 selnh2 setouse Hsqlo Hsqup Isqlo Isqup
scalar `Qart'=(Hsq)*(k-1)
scalar `selnh1' = 0.5*(ln(`Qart')-ln(k-1))/(sqrt(2*`Qart')-sqrt(2*k-3))
scalar `selnh2' = sqrt(1/(2*(k-2))*(1-1/(3*(k-2)^2)))
scalar `setouse'=`selnh1'
if `Qart'<k scalar `setouse'=`selnh2'
scalar `Hsqlo' = max(exp((ln(sqrt(Hsq))-`mltpl'*`setouse'))^2,0)
scalar `Hsqup' = exp((ln(sqrt(Hsq))+`mltpl'*`setouse'))^2
scalar `Isqlo' = max(100*(`Hsqlo'-1)/(`Hsqlo'),0)
scalar `Isqup' = min(100*(`Hsqup'-1)/(`Hsqup'),100)
/*if forest selected and overall analysis - collect some needed data*/
if ("`forest'"!="" | `forestw'!=-127.5) {
local `modsel'_mu`xlp'=`modsel'_mu
local `modsel'_lo`xlp'=`modsel'_lo
local `modsel'_up`xlp'=`modsel'_up
local tausqx`xlp'=tausqx
local Isq`xlp'=Isq
local k`xlp' = k
}
/*if forest selected and overall analysis - call forest plot*/
if ("`forest'"!="" | `forestw'!=-127.5) & `xlp'==`numloop' {
if `modsel'_mu==. | `modsel'_lo==. | `modsel'_up==. {
di as error "Forest plot not produced since the model did not converge"
}
/*NEW FOREST PLOT*/
//edits and options
qui replace weights=weights*100
tempvar rawdata tau2 df
//linked to unused option COUNTS
qui gen str1 `rawdata'=""
//linked to unused option RFDIST
qui gen `tau2'=.
qui gen `df'=.
//additional globals - post analyses
global MA_method1 `=upper("`modsel'")'
global MA_method2 "IV"
/*add the overall effects as an extra observations if requested*/
if "`overall'"=="" {
qui set obs `=k+`numloop'+2*(`numloop'-1)'
forvalues sg=1(1)`numloop' {
qui replace eff = ``modsel'_mu`sg'' in `=k+`sg''
qui replace studyloCI = ``modsel'_lo`sg'' in `=k+`sg''
qui replace studyupCI = ``modsel'_up`sg'' in `=k+`sg''
qui replace studyid = 0 in `=k+`sg''
if `sg'<`numloop' {
//results
qui replace use=3 in `=k+`sg''
qui sum weights if tgrpby==`sg'
qui replace weights=r(sum) in `=k+`sg''
if "$MA_nohet" != "" {
qui replace labelvar="Subtotal" in `=k+`sg''
}
else {
qui replace labelvar="Subtotal (I{superscript:2} = " + string(`Isq`sg'',"%5.1f") + "%)" in `=k+`sg''
}
qui replace tgrpby=`sg' in `=k+`sg''
}
else {
qui replace use=5 in `=k+`sg''
qui replace weights=100 in `=k+`sg''
if "$MA_nohet" != "" {
qui replace labelvar="Overall" in `=k+`sg''
}
else {
qui replace labelvar="Overall (I{superscript:2} = " + string(`Isq`sg'',"%5.1f") + "%)" in `=k+`sg''
}
}
qui replace `tau2' = `tausqx`sg'' in `=k+`sg''
qui replace `df' = `k`sg'' in `=k+`sg''
}
//add labels and space
forvalues sg=1(1)`=`numloop'-1' {
qui replace use=0 in `=k+`numloop'+2*`sg'-1'
qui replace labelvar="`tfoo`sg''" in `=k+`numloop'+2*`sg'-1'
qui replace tgrpby=`sg' in `=k+`numloop'+2*`sg'-1'
qui replace use=36 in `=k+`numloop'+2*`sg''
qui replace tgrpby=`sg' in `=k+`numloop'+2*`sg''
}
capture confirm variable tgrpby
if _rc==0 {
sort tgrpby use studyid
}
else {
sort use studyid
}
qui replace use=0 if use==36
}
*restore, not
*drop author outcome
*list
*error
//if meta-analysis exponentiated tweak label if not set - not needed
/*
if `btor'==1 & "`xlabel'"=="" {
qui sum studyupCI
local tempup = ceil(r(max))
if `tempup'<=2 {
local xlabel "xlabel(0,1,2)"
}
else {
local xlabel "xlabel(0,1,2,`tempup')"
}
local force="force"
}*/
_dispgby eff studyloCI studyupCI weights use labelvar `rawdata' `tau2' `df', `log' /*
*/ `xlabel' `xtick' `force' sumstat(`sumstat') `saving' `box' t1("`t1'") /*
*/ t2("`t2'") b1("`b1'") b2("`b2'") lcols("`lcols'") rcols("`rcols'") `overall' `wt' `stats' `xcounts' `eform' /*
*/ `groupla' `cornfield'
}
/*if plplot is selected and overall analysis call*/
if "`plplot'"!="" & `xlp'==`numloop' {
if `modsel'_mu==. | tsq_`modsel'==. {
di as error "Likelihood plot not produced since the model did not converge"
}
plplot `modsel'_mu tsq_`modsel' `plplotvar' `modsel' `btor'
}
/*more displays*/
di as text _newline(2) "Heterogeneity Measures"
/*cochran's Q*/
di as text "{hline 15}{c TT}{hline 35}
di as text "{col 16}{c |}{col 22}value{col 32}df{col 39}p-value"
di as text "{hline 15}{c +}{hline 35}
di as text %15s "Cochrane Q {col 16}{c |}" as result _col(20) %8.2f qw _col(29) %6.0f k-1 _col(38) %7.3f qpval
di as text "{hline 15}{c BT}{hline 35}
di _newline(1) as text "{hline 15}{c TT}{hline 35}
di as text "{col 16}{c |}{col 22}value{col 32}[`clvl'% Conf. Interval]"
di as text "{hline 15}{c +}{hline 35}
di as text %15s "I^2(%) {col 16}{c |}" as result _col(20) %8.2f Isq _col(29) %8.2f `Isqlo' _col(38) %8.2f `Isqup'
di as text %15s "H^2 {col 16}{c |}" as result _col(20) %8.2f Hsq _col(29) %8.2f `Hsqlo' _col(38) %8.2f `Hsqup'
if "`modsel'"!="pl" {
di as text %15s "tau^2(`tsqmdl') {col 16}{c |}" as result _col(20) %8.3f tausqx
}
di as text "{hline 15}{c BT}{hline 35}
/*if PL model was selected display in different box since it has different fields (CI interval)*/
if "`modsel'"=="pl" {
di _newline(1) as text "{hline 15}{c TT}{hline 35}
di as text "{col 16}{c |}{col 22}value{col 30}[`clvl'% Conf. Interval]"
di as text "{hline 15}{c +}{hline 35}
di as text %15s "tau^2 {col 16}{c |}" as result _col(20) %8.3f tausqx /*
*/ _col(29) %8.3f tsq_pl_lo _col(38) %8.3f tsq_pl_up
di as text "{hline 15}{c BT}{hline 35}
if ml_info != 0 {
di as text "Estimate obtained with Maximum likelihood - Profile likelihood provides the CI"
di as `ptxttype5' "`posttext5'"
}
}
else if "`modsel'"=="sa" {
di as `ptxttype5' "`posttext5'"
}
if "`modsel'"=="mh" | "`modsel'"=="peto" {
di as `ptxttype7' "`posttext7'"
}
//if analysis by group report comparison across subgroups
if "`grpby'"!="" & `xlp'==1 {
scalar sumqw=0
}
if "`grpby'"!="" & `xlp'==`numloop' {
scalar qdiff=qw-sumqw
scalar qpdiff = 1-chi2(`numloop'-2,qdiff)
scalar Hdiff=(qdiff-(`numloop'-2))/(`numloop'-2) /*relies on DL*/
if Hdiff<1 scalar Hdiff=1
scalar Idiff = 100*(Hdiff-1)/(Hdiff)
if Idiff<0 scalar Idiff=0
di as text "Test for subgroup differences: Chi^2=" %4.2f as res qdiff as text ", df=" %1.0f as res `=`numloop'-2' _continue
di as text ", p-value=" as res %5.3f qpdiff as text ", H^2=" as res %4.2f Hdiff as text ", I^2=" as res %4.2f Idiff
}
//save qw in locals to be used for subgroup differences evaluation, if grpby() option used
if "`grpby'"!="" {
scalar qw`xlp'=qw
scalar sumqw=sumqw+qw`xlp'
}
}
/*return list - irrelevant to grpby option, will only return overall results*/
/*universal*/
return scalar eff = `modsel'_mu
return scalar efflo = `modsel'_lo
return scalar effup = `modsel'_up
return scalar effvar = `modsel'_var
return scalar Q = qw
return scalar df = k-1
return scalar Qpval = qpval
return scalar Isq = Isq
return scalar Isq_lo = `Isqlo'
return scalar Isq_up = `Isqup'
return scalar Hsq = Hsq
return scalar Hsq_lo = `Hsqlo'
return scalar Hsq_up = `Hsqup'
return scalar tausq_dl = tsq_dl
/*specific*/
if "`bdl'"!="" {
return scalar tausq_bdl = tsq_bdl
}
else if "`sa'"!="" {
return scalar tausq_sa = tsq_sa
}
else if "`ml'"!="" {
return scalar tausq_ml = tsq_ml
return scalar conv_ml = ml_info
}
else if "`pl'"!="" {
return scalar tausq_pl = tsq_ml
return scalar tausqlo_pl = tsq_pl_lo
return scalar tausqup_pl = tsq_pl_up
return scalar conv_ml = ml_info
return scalar cloeff_pl = pl_lo_info
return scalar cupeff_pl = pl_up_info
return scalar ctausqlo_pl = tsq_pl_lo_info
return scalar ctausqup_pl = tsq_pl_up_info
}
else if "`reml'"!="" {
return scalar tausq_reml = tsq_reml
return scalar conv_reml = reml_info
}
else if "`pe'"!="" {
return scalar exec_pe = execpe
}
/*un-comment to get lots of created variables you can toy with*/
*restore, not
end
/*profile likelihood plot*/
program plplot
preserve
/*arguments*/
scalar mu = `1'
scalar tsq = `2'
local plottype = "`3'" /*mu or tsq*/
local modsel = "`4'"
if "`modsel'" == "reml" {
local tstr "REML"
}
else {
local tstr "ML/PL"
}
local expinf=`5'
if `expinf'==1 {
scalar mu=ln(mu)
local txtstr1=" (exponentiated)"
local txtstr2=" (not exponentiated)"
qui replace eff=ln(eff)
}
tempvar temp1
/*loop to get 100 points around the parameter of interest: mu-1 to mu+1 or 0 to 2*tsq*/
if "`plottype'"=="tsq" {
/*bounds*/
local lowerb = round(`=mu-1', 0.1)
local upperb = round(`=mu+1', 0.1)
local cntr=0
forvalues i=`lowerb'(0.01)`upperb' {
/*generate a counter so that we won't lose track!*/
local cntr = `cntr'+1
/*calculate the likelihood*/
qui gen double `temp1' = -1/2*ln(2*_pi*(effvar+tsq))-1/2*(eff-`i')^2/(effvar+tsq)
qui sum `temp1'
scalar yval`cntr' = r(sum)
scalar xval`cntr' = `i'
qui drop `temp1'
}
local sttle "for tau^2 fixed to the `tstr' estimate"
local xttle = "mu values`txtstr1'"
}
else {
/*different number of decimals depending on tsq value*/
scalar dlev = 0.1
if tsq < 0.1 scalar dlev = 0.01
/*bounds*/
local lowerb = 0
local upperb = round(max(`=2*tsq',0.1), dlev)
local cntr=0
forvalues i=`lowerb'(`=dlev/10')`upperb' {
/*generate a counter so that we won't lose track!*/
local cntr = `cntr'+1
/*calculate the likelihood*/
qui gen double `temp1' = -1/2*ln(2*_pi*(effvar+`i'))-1/2*(eff-mu)^2/(effvar+`i')
qui sum `temp1'
scalar yval`cntr' = r(sum)
scalar xval`cntr' = `i'
qui drop `temp1'
}
local sttle "for mu fixed to the `tstr' estimate`txtstr2'"
local xttle = "tau^2 values"
}
local ttle = "Likelihood plot"
local yttle = "log-likelihood"
/*now generate X and Y variables to create the graph*/
qui clear
qui set obs `cntr'
qui gen xvar=.
qui gen yvar=.
forvalues i=1(1)`cntr' {
qui replace xvar=xval`i' in `i'
qui replace yvar=yval`i' in `i'
}
//exponentiate only for the mu
if `expinf'==1 & "`plottype'"=="tsq"{
qui replace xvar=exp(xvar)
}
/*the graph*/
twoway mspline yvar xvar, lcolor(maroon) lwidth(medthick) /*
*/ ylabel(, labsize(2) angle(0)) ytitle("`yttle'", size(small)) /*
*/ xlabel(, labsize(2.5) angle(0)) xtitle("`xttle'", size(small)) /*
*/ title("`ttle'", size("medsmall"))/*
*/ subtitle("`sttle'", size("small"))
end
/*Profile Likelihood extras*/
/*program that estimates t^2 with mu fixed for Profile Likelihood method (used in tau CI estimation)*/
program PL_value, rclass
/*temporary variables*/
tempvar temp1 temp2
/*grab argument (the fixed mean) and assign*/
local tempx = `1'
scalar mu = `tempx'
/*starting value for the t^2 estimate*/
scalar t_sq_m = 0.001
/*dynamic system: we want it to converge to a fixed point - 500 iterations limit*/
scalar nx = 0
scalar diff188_1 = 1
/*continue with loop till t^2 estimate converges or we reach 500 iterations*/
while abs(diff188_1)>10^-6 & nx<500 {
scalar nx = nx + 1
/*calculate the new t_sq_m (t^2 estimate)*/
qui gen double `temp1' = ((eff-mu)^2-effvar)/(effvar+t_sq_m)^2
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = 1/(effvar+t_sq_m)^2
qui sum `temp2'
scalar sum2 = r(sum)
qui drop `temp1' `temp2'
/*calc the difference*/
scalar diff188_1 = t_sq_m - sum1/sum2
/*calculate the new value*/
scalar t_sq_m = sum1/sum2
}
/*calculate the log-likelihood value*/
qui gen double `temp1' = -1/2*ln(2*_pi*(effvar+t_sq_m))-1/2*(eff-mu)^2/(effvar+t_sq_m)
qui sum `temp1'
scalar mlval = r(sum)
qui drop `temp1'
/*return values to main program*/
return scalar mlv = mlval
end
/*program that estimates mu with t^2 fixed for Profile Likelihood method (used in t^2 CI estimation)*/
program PL_value2, rclass
/*temporary variables*/
tempvar temp1 temp2
/*grab argument (the fixed tau^2) and assign*/
local tempx = `1'
scalar t_sq_m = `tempx'
/*starting value for the mu estimate*/
scalar mu = 0
/*dynamic system: we want it to converge to a fixed point - 500 iterations limit*/
scalar nx = 0
scalar diff188_1 = 1
/*continue with loop till mu estimate converges or we reach 500 iterations*/
while abs(diff188_1)>10^-6 & nx<500 {
scalar nx = nx + 1
/*calculate the new mu estimate*/
qui gen double `temp1' = eff/(effvar+t_sq_m)
qui sum `temp1'
scalar sum1 = r(sum)
qui gen double `temp2' = 1/(effvar+t_sq_m)
qui sum `temp2'
scalar sum2 = r(sum)
qui drop `temp1' `temp2'
/*calc the difference*/
scalar diff188_1 = mu - sum1/sum2
/*calculate the new value*/
scalar mu = sum1/sum2
}
/*calculate the log-likelihood value*/
qui gen double `temp1' = -1/2*ln(2*_pi*(effvar+t_sq_m))-1/2*(eff-mu)^2/(effvar+t_sq_m)
qui sum `temp1'
scalar mlval = r(sum)
qui drop `temp1'
/*return values to main program*/
return scalar mlv = mlval
end
/*Permutation method extras*/
mata:
mata clear
/*calculates all the permutations for n columns*/
real matrix permmat(real scalar n)
{
real matrix A
real matrix tempA
real matrix permtable
real matrix pos
string matrix fin_0
string matrix found
/*first row of the matrix - all 1's*/
tempA = 1
A = 1
for (i=1; i< n; i++) {
A = (A,tempA)
}
/*up to this point created a line of n elements in A*/
permtable = A
pos = n
fin_0 = "false"
while (fin_0=="false") {
A[1,pos] = -1
/*make sure that when reversed at pos, all to the right (if any) are reset to 1*/
for (temp=pos+1; temp<=n; temp++) {
A[1,temp] = 1
}
/*reset the position to the end*/
pos = n
/*append the current permutation to the final table*/
permtable = (permtable\A)
/*check last element and change if it's 1*/
if (A[1,n] == 1) {
A[1,n] = -1
permtable = (permtable\A)
}
/*find first 1 on the left and alter pos accordingly*/
tpos = pos - 1
found ="false"
while (found=="false" & tpos>0) {
if (A[1,tpos]==1) {
found="true"
}
else {
tpos = tpos - 1
}
}
/*set new position*/
pos = tpos
if (pos==0) {
fin_0 = "true"
}
}
return(permtable)
}
/*creates random permutation table*/
real matrix randmat(real scalar n, real scalar k)
{
real scalar r
real matrix permtable
/*create random table - perhaps in the future get time to change seed number based on time*/
uniformseed(1976)
r = 2^k
permtable = uniform(r, n)
/*round numbers since they are in the (0,1) range - get them to be either 1 or -1*/
permtable = round(permtable)
permtable = permtable:*2
permtable = permtable:-1
return(permtable)
}
/*creates and returns an 1XN empty matrix*/
real matrix emptymat(real scalar N)
{
real matrix m
m = 0
for (i=2; i<= N; i++) {
m = (m\0)
}
return(m)
}
/*creates and returns a 2XN matrix with values for the meta-analysis*/
real matrix createmat(real scalar N)
{
real matrix A
real scalar temp, cntr, i
/*create first row*/
A = .
for (i=2; i<= N; i++) {
A = (A,.)
}
/*double up, for second row*/
A = (A\A)
/*fill up with effect sizes and variances data*/
i=0
cntr = 0
while (i<N) {
cntr = cntr + 1
temp = st_data(cntr,"use")
/*if study is supposed to be used then grab data and increment counter*/
if (temp==1) {
i = i + 1
A[1,i] = st_data(cntr,"effvar") /*variances in 1st row*/
A[2,i] = st_data(cntr,"eff") /*means in 2nd row*/
}
}
return(A)
}
/*returns the probability, upperCI and lowerCI in a 1X2 matrix*/
void pe_p(real scalar k, real matrix varmu, real matrix permtable, real matrix m,real scalar mu, real scalar re_mu, real scalar mltps)
{
real matrix means
real matrix Xperm
real matrix var
real matrix Xtemp
real matrix utemp
real matrix btemp
real matrix permtemp
real matrix b
real matrix cp
real matrix cptemp
real scalar numrows
real scalar posx
real scalar probab
real scalar pe_lo
real scalar pe_up
real matrix prcpu
real matrix prcpl
real matrix prcpt1
real matrix prcpt2
/*create the Xp table from tables permtable & the second line of varmu*/
means = abs(varmu[2,.]:-mu)
re_mu = re_mu-mu
Xperm = permtable:*means
/*grab the variances too*/
var = varmu[1,.]
/*find number of rows*/
numrows = rows(Xperm)
/*initialise the matrices that we will edit cell by cell*/
b = m:*0
cp = m:*0
prcpu = m:*0
prcpl = m:*0
/*grab each row and calculate the tau prediction for it - same dimensions as matrix m*/
for (i=1; i<= numrows; i++) {
Xtemp = Xperm[i,.]
/*calculate Sum(u(i))*/
utemp = 1:/var
su1 = sum(utemp)
/*calculate Sum(u(i)^2)*/
utemp = 1:/(var:*var)
su2 = sum(utemp)
/*calculate Sum(u(i)*X(i))*/
utemp = (Xtemp:*(1:/var)):/su1
m0 = sum(utemp)
/*calculate Sum{u(i)*(X(i)-m0)^2}*/
utemp = (1:/var):*((Xtemp:-m0):^2)
anum = sum(utemp)
/*calculate t_square estimate*/
t_sq = (anum-(k-1))/(su1-su2/su1)
if (t_sq<0) {
t_sq = 0
}
/*calculate m prediction for case i*/
/*calculate SUM(1/(t^2+varj))*/
utemp = 1:/(var:+t_sq)
tempS = sum(utemp)
utemp = ((1:/(var:+t_sq)):/tempS):*Xtemp
m[i] = sum(utemp)
/*needed for CI calculation*/
permtemp = permtable[i,.]
btemp = ((1:/(var:+t_sq)):/tempS):*permtemp
tempS = sum(btemp)
b[i] = tempS
tempS = (re_mu-m[i])/(1-b[i])
cp[i] = tempS
}
//Wolfgang
//permtable
//cp
/*more needed for the CI calculation*/
for (i=1; i<= numrows; i++) {
/*trying to count how many permutations satisfy mu-mu(p)>=c(1-b(p)). sqrt will set negatives to missing*/
cptemp = sqrt(cp:-cp[i]):+1
cptemp = round(cptemp:/cptemp)
tempS = sum(cptemp)/numrows
prcpu[i] = tempS
/*trying to count how many permutations satisfy mu-mu(p)<=c(1-b(p)).*/
cptemp = sqrt(-cp:+cp[i]):+1
cptemp = round(cptemp:/cptemp)
tempS = sum(cptemp)/numrows
prcpl[i] = tempS
}
//prcpu
//prcpl
/*find the smallest cp(i) for which prcpu(i)<=0.025 - upper CI*/
prcpt1 = sqrt(mltps:-prcpu):+1
prcpt2 = round(prcpt1:/prcpt1)
cptemp = prcpt2:*cp
pe_up = min(cptemp)
/*find the largest cp(i) for which prcpl(i)<=0.025 - lower CI*/
prcpt1 = sqrt(mltps:-prcpl):+1
prcpt2 = round(prcpt1:/prcpt1)
cptemp = prcpt2:*cp
pe_lo = max(cptemp)
//prcpl
//prcpt1
//prcpt2
//cptemp
/*BACK TO P value - sort the table*/
_sort(m,1)
/*call binary search to find the position*/
posx = bsearch(re_mu, m)
/*now from the position calculate probability*/
if (posx>numrows/2) { /*if the position in the lower half of the table...*/
probab = 2*(numrows-posx+1)/numrows
}
else { /*if the position in the upper half of the table...*/
probab = 2*posx/numrows
}
st_numscalar("r(pe_p)",probab)
st_numscalar("r(pe_lo)",pe_lo)
st_numscalar("r(pe_up)",pe_up)
}
/*binary search*/
real scalar bsearch(real scalar sval, real matrix A)
{
real scalar stpos
real scalar endpos
real scalar k
real scalar found
real scalar minsofar
real scalar temp
/*iniitalise start & end position*/
stpos = 1
endpos = rows(A)
found = 0 /*not found yet*/
/*main binary search loop*/
while (stpos<=endpos & found==0) {
k = trunc((stpos+endpos)/2)
if (abs(A[k]-sval)<10^-8) {
found=1
}
else {
if (sval<A[k]) {
endpos = k - 1
}
else {
stpos = k + 1
}
}
}
/*if it is not found, i.e. if re_mu was not included in the random matrix (for k>10)
we need to look for the closest value*/
if (found==0) {
/*look if A[k+1]is closer to re_mu(sval) than A[k]*/
minsofar = abs(A[k]-sval)
temp = k
if (temp<rows(A)){
if (abs(A[temp+1]-sval)<minsofar) {
minsofar = abs(A[temp+1]-sval)
k = temp + 1
}
}
if (temp>1){
if (abs(A[temp-1]-sval)<minsofar) {
minsofar = abs(A[temp-1]-sval)
k = temp - 1
}
}
}
return(k)
}
end
/*Mantel-Haenszel Odds-ratio fixed-effect method*/
program calc_MHORf, rclass
//not needed but remnant from older code
local i=127
//CI level
local clvl=c(level)
local mltpl=invnormal(1-(100-`clvl')/200)
//temp variables
tempvar a`i' b`i' c`i' d`i' N`i' OR`i' varOR`i' w`i'
/*only use OR, effmeas= 1 or 2(reversed)*/
qui gen `a`i'' = .
qui gen `b`i'' = .
qui gen `c`i'' = .
qui gen `d`i'' = .
/*events and non events*/
qui replace `a`i'' = ev1_`i'
qui replace `b`i'' = total1_`i'-ev1_`i'
qui replace `c`i'' = ev2_`i'
qui replace `d`i'' = total2_`i'-ev2_`i'
/*correction for zero cells*/
foreach x in a b c d {
qui replace ``x'`i''=``x'`i''+0.5 if `a`i''<1 | `c`i''<1 | `b`i''<1 | `d`i''<1
}
qui gen `N`i'' = `a`i''+`b`i''+`c`i''+`d`i''
/*but if there are too many that are zero, set all to zero*/
foreach x in a b c d {
qui replace ``x'`i''=0 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
}
/*more calculations*/
qui gen double `OR`i'' = (`a`i''*`d`i'')/(`b`i''*`c`i'')
qui gen double `varOR`i'' = (1/`a`i''+1/`d`i''+1/`b`i''+1/`c`i'')
qui gen double `w`i'' = (`b`i''*`c`i'')/`N`i''
/*if one of the cells is still zero set the numerators in the sums to zero to avoid problems*/
qui replace `w`i''=1 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
qui replace `N`i''=1 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
qui replace `OR`i''=0 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
qui replace `varOR`i''=1 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
//intermediate variables
tempvar temp1 temp2 vE vF vG vH vR vS
qui gen double `temp1' = `w`i''*`OR`i''
qui gen double `temp2' = `w`i''
qui gen double `vR' = (`a`i''*`d`i'')/`N`i''
qui gen double `vS' = (`b`i''*`c`i'')/`N`i''
qui gen double `vE' = ((`a`i''+`d`i'')*`a`i''*`d`i'')/(`N`i''^2)
qui gen double `vF' = ((`a`i''+`d`i'')*`b`i''*`c`i'')/(`N`i''^2)
qui gen double `vG' = ((`b`i''+`c`i'')*`a`i''*`d`i'')/(`N`i''^2)
qui gen double `vH' = ((`b`i''+`c`i'')*`b`i''*`c`i'')/(`N`i''^2)
foreach x in temp1 temp2 vE vF vG vH vR vS {
qui sum ``x''
local l`x'=r(sum)
}
/*mean estimate (coefficient)*/
local MHf_mu = ln(`ltemp1'/`ltemp2')
/*variance estimate*/
local MHf_var =0.5*(`lvE'/(`lvR'^2) + (`lvF'+`lvG')/(`lvR'*`lvS') + `lvH'/(`lvS'^2))
local MHf_lo = `MHf_mu'-`mltpl'*sqrt(`MHf_var')
local MHf_up = `MHf_mu'+`mltpl'*sqrt(`MHf_var')
/*scalars*/
foreach x in mu var lo up {
scalar `x'=`MHf_`x''
}
//create study variables only when requested to do so
if "`1'"=="gen" {
capture drop eff effvar weights
gen double eff=ln(`OR`i'')
gen double effvar=`varOR`i''
qui sum `w`i''
gen double weights=`w`i''/r(sum)
}
/*return overall to main program*/
foreach x in mu var lo up {
return scalar `x'=`x'
}
end
/*Mantel-Haenszel Risk-ratio fixed-effect method*/
program calc_MHRRf, rclass
//not needed but remnant from older code
local i=127
//CI level
local clvl=c(level)
local mltpl=invnormal(1-(100-`clvl')/200)
//temp variables
tempvar a`i' b`i' c`i' d`i' N`i' RR`i' varRR`i' w`i'
/*only use OR, effmeas= 1 or 2(reversed)*/
qui gen `a`i'' = .
qui gen `b`i'' = .
qui gen `c`i'' = .
qui gen `d`i'' = .
/*events and non events*/
qui replace `a`i'' = ev1_`i'
qui replace `b`i'' = total1_`i'-ev1_`i'
qui replace `c`i'' = ev2_`i'
qui replace `d`i'' = total2_`i'-ev2_`i'
/*correction for zero cells*/
foreach x in a b c d {
qui replace ``x'`i''=``x'`i''+0.5 if `a`i''<1 | `c`i''<1 | `b`i''<1 | `d`i''<1
}
qui gen `N`i'' = `a`i''+`b`i''+`c`i''+`d`i''
/*but if there are too many that are zero, set all to zero*/
foreach x in a b c d {
qui replace ``x'`i''=0 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
}
/*more calculations*/
qui gen double `RR`i'' = (`a`i''/(`a`i''+`b`i''))/(`c`i''/(`c`i''+`d`i''))
qui gen double `varRR`i'' = (1/`a`i'' + 1/`c`i'' - 1/(`a`i''+`b`i'') - 1/(`c`i''+`d`i''))
qui gen double `w`i'' = `c`i''*(`a`i''+`b`i'')/`N`i''
/*if one of the cells is still zero set the numerators in the sums to zero to avoid problems*/
qui replace `w`i''=1 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
qui replace `N`i''=1 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
qui replace `RR`i''=0 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
qui replace `varRR`i''=1 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
//intermediate variables
tempvar temp1 temp2 vP vR vS
qui gen double `temp1' = `w`i''*`RR`i''
qui gen double `temp2' = `w`i''
qui gen double `vP' = ((`a`i''+`b`i'')*(`c`i''+`d`i'')*(`a`i''+`c`i'') - `a`i''*`c`i''*`N`i'')/(`N`i''^2)
qui gen double `vR' = (`a`i''*(`c`i''+`d`i''))/`N`i''
qui gen double `vS' = (`c`i''*(`a`i''+`b`i''))/`N`i''
foreach x in temp1 temp2 vP vR vS {
qui sum ``x''
local l`x'=r(sum)
}
/*mean estimate (coefficient)*/
local MHf_mu = ln(`ltemp1'/`ltemp2')
/*variance estimate*/
local MHf_var =(`lvP'/(`lvR'*`lvS'))
local MHf_lo = `MHf_mu'-`mltpl'*sqrt(`MHf_var')
local MHf_up = `MHf_mu'+`mltpl'*sqrt(`MHf_var')
/*scalars*/
foreach x in mu var lo up {
scalar `x'=`MHf_`x''
}
//create study variables only when requested to do so
if "`1'"=="gen" {
capture drop eff effvar weights
gen double eff=ln(`RR`i'')
gen double effvar=`varRR`i''
qui sum `w`i''
gen double weights=`w`i''/r(sum)
}
/*return overall to main program*/
foreach x in mu var lo up {
return scalar `x'=`x'
}
end
//METHODS FOR 4 VARIABLE SYNTAX: EVENTS AND POPULATIONS
/*Mantel-Haenszel Risk difference fixed-effect method*/
program calc_MHRDf, rclass
//not needed but remnant from older code
local i=127
//CI level
local clvl=c(level)
local mltpl=invnormal(1-(100-`clvl')/200)
//temp variables
tempvar a`i' b`i' c`i' d`i' N`i' RD`i' varRD`i' w`i'
/*only use OR, effmeas= 1 or 2(reversed)*/
qui gen `a`i'' = .
qui gen `b`i'' = .
qui gen `c`i'' = .
qui gen `d`i'' = .
/*events and non events*/
qui replace `a`i'' = ev1_`i'
qui replace `b`i'' = total1_`i'-ev1_`i'
qui replace `c`i'' = ev2_`i'
qui replace `d`i'' = total2_`i'-ev2_`i'
/*correction for zero cells*/
foreach x in a b c d {
/*`a`i''<1 | `c`i''<1 | `b`i''<1 | `d`i''<1*/
qui replace ``x'`i''=``x'`i''+0.5 if (`a`i''<1 & `c`i''<1) | (`a`i''<1 & `d`i''<1) | (`b`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
}
qui gen `N`i'' = `a`i''+`b`i''+`c`i''+`d`i''
/*more calculations*/
qui gen double `RD`i'' = (`a`i''/(`a`i''+`b`i'')) - (`c`i''/(`c`i''+`d`i''))
qui gen double `varRD`i'' = (`a`i''*`b`i'')/((`a`i''+`b`i'')^3) + (`c`i''*`d`i'')/((`c`i''+`d`i'')^3)
qui gen double `w`i'' = ((`a`i''+`b`i'')*(`c`i''+`d`i''))/`N`i''
//intermediate variables
tempvar temp1 temp2 vJ vK
qui gen double `temp1' = `w`i''*`RD`i''
qui gen double `temp2' = `w`i''
qui gen double `vJ' = (`a`i''*`b`i''*((`c`i''+`d`i'')^3)+`c`i''*`d`i''*((`a`i''+`b`i'')^3))/((`a`i''+`b`i'')*(`c`i''+`d`i'')*(`N`i''^2))
qui gen double `vK' = ((`a`i''+`b`i'')*(`c`i''+`d`i''))/`N`i''
foreach x in temp1 temp2 vJ vK {
qui sum ``x''
local l`x'=r(sum)
}
/*mean estimate*/
local MHf_mu = (`ltemp1'/`ltemp2')
/*variance estimate*/
local MHf_var =(`lvJ'/(`lvK'^2))
local MHf_lo = `MHf_mu'-`mltpl'*sqrt(`MHf_var')
local MHf_up = `MHf_mu'+`mltpl'*sqrt(`MHf_var')
/*scalars*/
foreach x in mu var lo up {
scalar `x'=`MHf_`x''
}
//create study variables only when requested to do so
if "`1'"=="gen" {
capture drop eff effvar weights
gen double eff=(`RD`i'')
gen double effvar=`varRD`i''
qui sum `w`i''
gen double weights=`w`i''/r(sum)
}
/*return overall to main program*/
foreach x in mu var lo up {
return scalar `x'=`x'
}
end
/*Peto Odds-ratio fixed-effect method*/
program calc_PORf, rclass
//not needed but remnant from older code
local i=127
//CI level
local clvl=c(level)
local mltpl=invnormal(1-(100-`clvl')/200)
//temp variables
tempvar a`i' b`i' c`i' d`i' N`i' OR`i' varOR`i' Z`i' V`i'
//if effect already there don't calculate
if "`2'"!="effex" {
/*only use OR, effmeas= 1 or 2(reversed)*/
qui gen `a`i'' = .
qui gen `b`i'' = .
qui gen `c`i'' = .
qui gen `d`i'' = .
/*events and non events*/
qui replace `a`i'' = ev1_`i'
qui replace `b`i'' = total1_`i'-ev1_`i'
qui replace `c`i'' = ev2_`i'
qui replace `d`i'' = total2_`i'-ev2_`i'
/*correction for zero cells*/
foreach x in a b c d {
*qui replace ``x'`i''=``x'`i''+0.5 if `a`i''<1 | `c`i''<1 | `b`i''<1 | `d`i''<1
}
qui gen `N`i'' = `a`i''+`b`i''+`c`i''+`d`i''
/*but if there are too many that are zero, set all to zero*/
foreach x in a b c d {
*qui replace ``x'`i''=0 if (`a`i''<1 & `c`i''<1) | (`b`i''<1 & `d`i''<1)
}
/*more calculations*/
qui gen double `Z`i'' = `a`i'' - (`a`i''+`b`i'')*(`a`i''+`c`i'')/`N`i''
qui gen double `V`i'' = (`a`i''+`b`i'')*(`c`i''+`d`i'')*(`a`i''+`c`i'')*(`b`i''+`d`i'')/(`N`i''^2*(`N`i''-1))
qui gen double `OR`i'' = exp(`Z`i''/`V`i'')
qui gen double `varOR`i'' = 1/`V`i''
/*can't have any problems with zero cells: only when a=c=b=d=0 which is impossible*/
}
else {
qui gen double `V`i'' = 1/effvar
qui gen double `OR`i'' = exp(eff)
qui gen double `varOR`i'' = 1/`V`i''
}
//intermediate variables
tempvar temp1 temp2
qui gen double `temp1' = `V`i''*ln(`OR`i'')
qui gen double `temp2' = `V`i''
foreach x in temp1 temp2 {
qui sum ``x''
local l`x'=r(sum)
}
/*mean estimate (coefficients)*/
local Pf_mu = (`ltemp1'/`ltemp2')
/*variance estimate*/
local Pf_var = 1/`ltemp2'
local Pf_lo = `Pf_mu' - `mltpl'*sqrt(`Pf_var')
local Pf_up = `Pf_mu' + `mltpl'*sqrt(`Pf_var')
/*scalars*/
foreach x in mu var lo up {
scalar `x'=`Pf_`x''
}
//create study variables only when requested to do so
if "`1'"=="gen" {
capture drop eff effvar
gen double eff=ln(`OR`i'')
gen double effvar=`varOR`i''
}
if "`1'"=="gen" | "`1'"=="partgen" {
capture drop weights
qui sum `V`i''
gen double weights=`V`i''/r(sum)
}
/*return overall to main program*/
foreach x in mu var lo up {
return scalar `x'=`x'
}
end
**********************************************************
*** ***
*** NEW ***
*** _DISPGBY PROGRAM ***
*** ROSS HARRIS ***
*** JULY 2006 ***
*** * * * ***
*** ***
**********************************************************
//"appropriated" with many thanks
program define _dispgby
version 9.0
// AXmin AXmax ARE THE OVERALL LEFT AND RIGHT COORDS
// DXmin dxMAX ARE THE LEFT AND RIGHT COORDS OF THE GRAPH PART
#delimit ;
syntax varlist(min=6 max=10 default=none ) [if] [in] [,
LOG XLAbel(string) XTICK(string) FORCE SAVING(string) noBOX SUMSTAT(string)
T1(string) T2(string) B1(string) B2(string) LCOLS(string) /* JUNK NOW */
RCOLS(string) noOVERALL noWT noSTATS COUNTS EFORM
noGROUPLA CORNFIELD];
#delimit cr
tempvar effect lci uci weight wtdisp use label tlabel id yrange xrange Ghsqrwt rawdata i2 mylabel
tokenize "`varlist'", parse(" ")
qui{
gen `effect'=`1'
gen `lci' =`2'
gen `uci' =`3'
gen `weight'=`4' // was 4
gen byte `use'=`5'
gen str `label'=`6'
gen str `mylabel'=`6'
if "`lcols'" == ""{
local lcols "`mylabel'"
label var `mylabel' "Study ID"
}
gen str80 `rawdata' = `7'
compress `rawdata'
if "`8'"!="" & "$MA_rjhby" != ""{
*gen `wtdisp'=`8'
gen `wtdisp'=`weight'
}
else {
gen `wtdisp'=`weight'
}
*if "`10'" != "" & "$MA_rjhby" != ""{
if "$MA_rjhby" != ""{
tempvar tau2 df
gen `tau2' = `8'
gen `df' = `9'
}
if "`9'" != "" & "$MA_rjhby" == ""{ // DIFFERENT IF FROM metan OR metanby
tempvar tau2 df
gen `tau2' = `8'
gen `df' = `9'
}
replace `weight' = `wtdisp' // bodge solu for SG weights
if "$MA_summaryonly" != ""{
drop if `use' == 1
}
// SET UP EXTENDED CIs FOR RANDOM EFFECTS DISTRIBUTION
// THIS CODE IS A BIT NASTY AS I SET THIS UP BADLY INITIALLY
// REQUIRES MAJOR REWORK IDEALLY...
tempvar tauLCI tauUCI SE tauLCIinf tauUCIinf
*replace `tau2' = .a if `tau2' == 0 // no heterogeneity
replace `tau2' = .b if `df'-1 == 0 // inestimable predictive distribution
//replace `tau2' = . if (`use' == 5 | `use' == 3) & "$MA_method1" != "D+L"
//replace `tau2' = . if (`use' == 17 | `use' == 19) & "$MA_method2" != "D+L"
replace `tau2' = . if (`use' == 5 | `use' == 3) & "$MA_method1" == "FE"
replace `tau2' = . if (`use' == 17 | `use' == 19) & "$MA_method2" == "FE"
gen `tauLCI' = .
gen `tauUCI' = .
gen `tauLCIinf' = .
gen `tauUCIinf' = .
gen `SE' = .
// modified so rf CI (rflevel) used
if "$MA_rfdist" != ""{
if ( ("`sumstat'"=="OR" | "`sumstat'"=="RR" | "`sumstat'"=="HR") & ("`log'"=="") ) | ("`eform'"!="") {
replace `SE' = (ln(`uci')-ln(`lci')) / (invnorm($RFL/200+0.5)*2)
replace `tauLCI' = exp( ln(`effect') - invttail((`df'-1), 0.5-$RFL/200)*sqrt( `tau2' +`SE'^2 ) )
replace `tauUCI' = exp( ln(`effect') + invttail((`df'-1), 0.5-$RFL/200)*sqrt( `tau2' +`SE'^2 ) )
replace `tauLCI' = 1e-9 if `tau2' == .b
replace `tauUCI' = 1e9 if `tau2' == .b
}
else{
replace `SE' = (`uci'-`lci') / (invnorm($RFL/200+0.5)*2)
replace `tauLCI' = `effect'-invttail((`df'-1), 0.5-$RFL/200)*sqrt(`tau2'+`SE'^2)
replace `tauUCI' = `effect'+invttail((`df'-1), 0.5-$RFL/200)*sqrt(`tau2'+`SE'^2)
replace `tauLCI' = -1e9 if `tau2' == .b
replace `tauUCI' = 1e9 if `tau2' == .b
}
}
if "$MA_rfdist" != ""{
qui count
local prevN = r(N)
tempvar expTau orderTau
gen `orderTau' = _n
gen `expTau' = 1
replace `expTau' = 2 if `tau2' != . // but expand if .a or .b
expand `expTau'
replace `use' = 4 if _n > `prevN'
replace `orderTau' = `orderTau' + 0.5 if _n > `prevN'
sort `orderTau'
}
tempvar estText weightText RFdistText RFdistLabel
local dp = $MA_dp
gen str `estText' = string(`effect', "%10.`dp'f") + " (" + string(`lci', "%10.`dp'f") + ", " +string(`uci', "%10.`dp'f") + ")"
replace `estText' = "(Excluded)" if `use' == 2
// don't show effect size again, just CI
gen `RFdistLabel' = "with estimated predictive interval" if `use' == 4 & `tau2' < .
gen `RFdistText' = /* string(`effect', "%10.`dp'f") + */ ". (" + string(`tauLCI', "%10.`dp'f") + ", " +string(`tauUCI', "%10.`dp'f") ///
+ ")" if `use' == 4 & `tau2' < .
/* not used
replace `RFdistLabel' = "No observed heterogeneity" if `use' == 4 & `tau2' == .a
replace `RFdistText' = string(`effect', "%10.`dp'f") + " (" + string(`lci', "%10.`dp'f") + ", " +string(`uci', "%10.`dp'f") ///
+ ")" if `use' == 4 & `tau2' == .a
*/
// don't show effect size again, just CI
replace `RFdistLabel' = "Inestimable predictive distribution with <3 studies" if `use' == 4 & `tau2' == .b
replace `RFdistText' = /* string(`effect', "%4.2f") + */ ". ( - , - )" if `use' == 4 & `tau2' == .b
qui replace `estText' = " " + `estText' if `effect' >= 0 & `use' != 4
gen str `weightText' = string(`weight', "%4.2f")
replace `weightText' = "" if `use' == 17 | `use' == 19 // can cause confusion and not necessary
replace `rawdata' = "" if `use' == 17 | `use' == 19
if "`counts'" != ""{
if $MA_params == 6{
local type "N, mean (SD);"
}
else{
local type "Events,"
}
tempvar raw1 raw2
gen str `raw1' = substr(`rawdata',1,(strpos(`rawdata',";")-1) )
gen str `raw2' = substr(`rawdata',(strpos(`rawdata',";")+1), (length(`rawdata')-strpos(`rawdata',";")) )
label var `raw1' "`type' $MA_G1L"
label var `raw2' "`type' $MA_G2L"
}
/* RJH - probably a better way to get this but I've just used globals from earlier */
if "`overall'" == "" & "$MA_nohet" == ""{
if "$MA_method1" == "USER"{
if "$MA_firststats" != ""{
replace `label' = "Overall ($MA_firststats)" if `use'==5
}
else{
replace `label' = "Overall" if `use'==5
}
}
replace `label' = "Overall ($MA_secondstats)" if `use' == 17 & "$MA_method2" == "USER" & "$MA_secondstats" != ""
replace `label' = "Overall" if `use' == 17 & "$MA_method2" == "USER" & "$MA_secondstats" == ""
}
if "`overall'" == "" & "$MA_nohet" != ""{
replace `label' = "Overall" if `use' == 5 | `use' == 17
}
tempvar hetGroupLabel expandOverall orderOverall
if "$MA_rjhby" != "" & "$MA_nohet" == "" & "$MA_method1" == "IV"{
* replace `label' = `label' + ";" if `use' == 5
qui count
local prevMax = r(N)
gen `orderOverall' = _n
gen `expandOverall' = 1
replace `expandOverall' = 2 if `use' == 5
expand `expandOverall'
replace `orderOverall' = `orderOverall' -0.5 if _n > `prevMax'
gen `hetGroupLabel' = "Heterogeneity between groups: p = " + ///
string($rjhHetGrp, "%5.3f") if _n > `prevMax'
replace `use' = 4 if _n > `prevMax'
sort `orderOverall'
}
else{
gen `hetGroupLabel' = .
}
replace `label' = "Overall" if `use' == 17 & "$MA_method2" != "USER"
replace `label' = "Subtotal" if `use' == 19
qui count if (`use'==1 | `use'==2)
local ntrials=r(N)
qui count if (`use'>=0 & `use'<=5)
local ymax=r(N)
gen `id'=`ymax'-_n+1 if `use'<9 | `use' == 17 | `use' == 19
if "$MA_method2" != "" | "$MA_method1" == "USER" {
local dispM1 = "$MA_method1"
local dispM2 = "$MA_method2"
if "$MA_method1" == "USER"{
local dispM1 "$MA_userDescM"
}
if "$MA_method2" == "USER"{
local dispM2 "$MA_userDesc"
}
replace `label' = "`dispM1'" + " " + `label' if (`use' == 3 | `use' == 5) & substr(`label',1,3) != "het"
replace `label' = "`dispM2'" + " " + `label' if `use' == 17 | `use' == 19
}
// GET MIN AND MAX DISPLAY
// SORT OUT TICKS- CODE PINCHED FROM MIKE AND FIDDLED. TURNS OUT I'VE BEEN USING SIMILAR NAMES...
// AS SUGGESTED BY JS JUST ACCEPT ANYTHING AS TICKS AND RESPONSIBILITY IS TO USER!
qui summ `lci', detail
local DXmin = r(min)
qui summ `uci', detail
local DXmax = r(max)
local h0 = 0
// MIKE MAKES A MAX VALUE IF SOMETHING EXTREME OCCURS...
if (( ("`sumstat'"=="OR" | "`sumstat'"=="RR" | "`sumstat'"=="HR") & ("`log'"=="") ) | ("`eform'"!="")) {
local h0=1
local Glog "xlog"
local xlog "log"
local xexp "exp"
replace `lci'=1e-9 if `lci'<1e-8
replace `lci'=1e9 if `lci'>1e8 & `lci'!=.
replace `uci'=1e-9 if `uci'<1e-8
replace `uci'=1e9 if `uci'>1e8 & `uci'!=.
if `DXmin'<1e-8 {
local DXmin=1e-8
}
if `DXmax'>1e8 {
local DXmax=1e8
}
}
if "$MA_NULL" != ""{
local h0 = $MA_NULL
}
if `h0' != 0 & `h0' != 1{
noi di "Null specified as `h0' in graph- for most effect measures null is 0 or 1"
}
if "`cornfield'"!="" {
replace `lci'=`log'(1e-9) if ( (`lci'==. | `lci'==0) & (`effect'!=. & `use'==1) )
replace `uci'=`log'(1e9) if ( (`uci'==.) & (`effect'!=. & `use'==1) )
}
// THIS BIT CHANGED- THE USER CAN PUT ANYTHING IN
local flag1=0
if ("`xlabel'"=="" | "`xtick'" == "") & "$MA_nulloff" == ""{ // if no xlabel or tick
local xtick "`h0'"
}
if "`xlabel'"==""{
local Gmodxhi=max( abs(`xlog'(`DXmin')),abs(`xlog'(`DXmax')))
if `Gmodxhi'==. {
local Gmodxhi=2
}
local DXmin=`xexp'(-`Gmodxhi')
local DXmax=`xexp'( `Gmodxhi')
if "$MA_nulloff" == ""{
local xlabel "`DXmin',`h0',`DXmax'"
}
else{
local xlabel "`DXmin',`DXmax'"
}
}
local DXmin2 = min(`xlabel',`DXmin')
local DXmax2 = max(`xlabel',`DXmax')
if "`force'" == ""{
local Gmodxhi=max( abs(`xlog'(`DXmin')), abs(`xlog'(`DXmax')), ///
abs(`xlog'(`DXmin2')), abs(`xlog'(`DXmax2')) )
if `Gmodxhi'==. {
local Gmodxhi=2
}
local DXmin=`xexp'(-`Gmodxhi')
local DXmax=`xexp'( `Gmodxhi')
if "`xlabel'" != "" & "$MA_nulloff" == ""{
local xlabel "`h0',`xlabel'"
}
}
if "`force'" != ""{
local DXmin = min(`xlabel')
local DXmax = max(`xlabel')
if "$MA_nulloff" == ""{
local xlabel "`h0',`xlabel'"
}
}
// LABELS- DON'T ALLOW SILLY NO. OF DECIMAL PLACES
local lblcmd ""
tokenize "`xlabel'", parse(",")
while "`1'" != ""{
if "`1'" != ","{
local lbl = string(`1',"%7.3g")
local val = `1'
local lblcmd `lblcmd' `val' "`lbl'"
}
mac shift
}
if "`xtick'" == ""{
local xtick = "`xlabel'"
}
local xtick2 = ""
tokenize "`xtick'", parse(",")
while "`1'" != ""{
if "`1'" != ","{
local xtick2 = "`xtick2' " + string(`1')
}
if "`1'" == ","{
local xtick2 = "`xtick2'`1'"
}
mac shift
}
local xtick = "`xtick2'"
local DXmin=`xlog'(min(`xlabel',`xtick',`DXmin'))
local DXmax=`xlog'(max(`xlabel',`xtick',`DXmax'))
if ("`eform'" != "" | "`xlog'" != "") {
local lblcmd ""
tokenize "`xlabel'", parse(",")
while "`1'" != ""{
if "`1'" != ","{
local lbl = string(`1',"%7.3g")
local val = ln(`1')
local lblcmd `lblcmd' `val' "`lbl'"
}
mac shift
}
replace `effect' = ln(`effect')
replace `lci' = ln(`lci')
replace `uci' = ln(`uci')
replace `tauLCI' = ln(`tauLCI')
replace `tauUCI' = ln(`tauUCI')
local xtick2 ""
tokenize "`xtick'", parse(",")
while "`1'" != ""{
if "`1'" != ","{
local ln = ln(`1')
local xtick2 "`xtick2' `ln'"
}
if "`1'" == ","{
local xtick2 "`xtick2'`1'"
}
mac shift
}
local xtick "`xtick2'"
local h0 = 0
}
// JUNK
*noi di "min: `DXmin', `DXminLab'; h0: `h0', `h0Lab'; max: `DXmax', `DXmaxLab'"
local DXwidth = `DXmax'-`DXmin'
if `DXmin' > 0{
local h0 = 1
}
} // END QUI
// END OF TICKS AND LABLES
// MAKE OFF-SCALE ARROWS
qui{
tempvar offLeftX offLeftX2 offRightX offRightX2 offYlo offYhi
local arrowWidth = 0.02 // FRACTION OF GRAPH WIDTH
local arrowHeight = 0.5/2 // Y SCALE IS JUST ORDERED NUMBER- 2x0.25 IS 0.5 OF AVAILABLE SPACE
gen `offLeftX' = `DXmin' if `lci' < `DXmin' | `tauLCI' < `DXmin'
gen `offLeftX2' = `DXmin' + `DXwidth'*`arrowWidth' if `lci' < `DXmin' | `tauLCI' < `DXmin'
gen `offRightX' = `DXmax' if `uci' > `DXmax' | (`tauUCI' > `DXmax' & `tauLCI' < .)
gen `offRightX2' = `DXmax' - `DXwidth'*`arrowWidth' if `uci' > `DXmax' | (`tauUCI' > `DXmax' & `tauLCI' < .)
gen `offYlo' = `id' - `arrowHeight'
gen `offYhi' = `id' + `arrowHeight'
replace `lci' = `DXmin' if `lci' < `DXmin' & (`use' == 1 | `use' == 2)
replace `uci' = `DXmax' if `uci' > `DXmax' & (`use' == 1 | `use' == 2)
replace `lci' = . if `uci' < `DXmin' & (`use' == 1 | `use' == 2)
replace `uci' = . if `lci' > `DXmax' & (`use' == 1 | `use' == 2)
replace `effect' = . if `effect' < `DXmin' & (`use' == 1 | `use' == 2)
replace `effect' = . if `effect' > `DXmax' & (`use' == 1 | `use' == 2)
} // end qui
************************
** COLUMNS **
************************
// OPTIONS FOR L-R JUSTIFY?
// HAVE ONE MORE COL POSITION THAN NECESSARY, COULD THEN R-JUSTIFY
// BY ADDING 1 TO LOOP, ALSO HAVE MAX DIST FOR OUTER EDGE
// HAVE USER SPECIFY % OF GRAPH USED FOR TEXT?
qui{ // KEEP QUIET UNTIL AFTER DIAMONDS
local titleOff = 0
if "`lcols'" == ""{
local lcols = "`label'"
local titleOff = 1
}
// DOUBLE LINE OPTION
if "$MA_DOUBLE" != "" & ("`lcols'" != "" | "`rcols'" != ""){
tempvar expand orig
gen `orig' = _n
gen `expand' = 1
replace `expand' = 2 if `use' == 1
expand `expand'
sort `orig'
replace `id' = `id' - 0.45 if `id' == `id'[_n-1]
replace `use' = 2 if mod(`id',1) != 0 & `use' != 5
replace `effect' = . if mod(`id',1) != 0
replace `lci' = . if mod(`id',1) != 0
replace `uci' = . if mod(`id',1) != 0
replace `estText' = "" if mod(`id',1) != 0
cap replace `raw1' = "" if mod(`id',1) != 0
cap replace `raw2' = "" if mod(`id',1) != 0
replace `weightText' = "" if mod(`id',1) != 0
foreach var of varlist `lcols' `rcols'{
cap confirm string var `var'
if _rc == 0{
tempvar length words tosplit splitwhere best
gen `splitwhere' = 0
gen `best' = .
gen `length' = length(`var')
summ `length', det
gen `words' = wordcount(`var')
gen `tosplit' = 1 if `length' > r(max)/2+1 & `words' >= 2
summ `words', det
local max = r(max)
forvalues i = 1/`max'{
replace `splitwhere' = strpos(`var',word(`var',`i')) ///
if abs( strpos(`var',word(`var',`i')) - length(`var')/2 ) < `best' ///
& `tosplit' == 1
replace `best' = abs(strpos(`var',word(`var',`i')) - length(`var')/2) ///
if abs(strpos(`var',word(`var',`i')) - length(`var')/2) < `best'
}
replace `var' = substr(`var',1,(`splitwhere'-1)) if `tosplit' == 1 & mod(`id',1) == 0
replace `var' = substr(`var',`splitwhere',length(`var')) if `tosplit' == 1 & mod(`id',1) != 0
replace `var' = "" if `tosplit' != 1 & mod(`id',1) != 0 & `use' != 5
drop `length' `words' `tosplit' `splitwhere' `best'
}
if _rc != 0{
replace `var' = . if mod(`id',1) != 0 & `use' != 5
}
}
}
summ `id' if `use' != 9
local max = r(max)
local new = r(N)+4
if `new' > _N {
set obs `new'
}
forvalues i = 1/4{ // up to four lines for titles
local multip = 1
local add = 0
if "$MA_DOUBLE" != ""{ // DOUBLE OPTION- CLOSER TOGETHER, GAP BENEATH
local multip = 0.45
local add = 0.5
}
local idNew`i' = `max' + `i'*`multip' + `add'
local Nnew`i'=r(N)+`i'
local tmp = `Nnew`i''
replace `id' = `idNew`i'' + 1 in `tmp'
replace `use' = 1 in `tmp'
if `i' == 1{
global borderline = `idNew`i''-0.25
}
}
local maxline = 1
if "`lcols'" != ""{
tokenize "`lcols'"
local lcolsN = 0
while "`1'" != ""{
cap confirm var `1'
if _rc!=0 {
di in re "Variable `1' not defined"
exit _rc
}
local lcolsN = `lcolsN' + 1
tempvar left`lcolsN' leftLB`lcolsN' leftWD`lcolsN'
cap confirm string var `1'
if _rc == 0{
gen str `leftLB`lcolsN'' = `1'
}
if _rc != 0{
cap decode `1', gen(`leftLB`lcolsN'')
if _rc != 0{
local f: format `1'
gen str `leftLB`lcolsN'' = string(`1', "`f'")
replace `leftLB`lcolsN'' = "" if `leftLB`lcolsN'' == "."
}
}
replace `leftLB`lcolsN'' = "" if (`use' != 1 & `use' != 2)
local colName: variable label `1'
if "`colName'"==""{
local colName = "`1'"
}
// WORK OUT IF TITLE IS BIGGER THAN THE VARIABLE
// SPREAD OVER UP TO FOUR LINES IF NECESSARY
local titleln = length("`colName'")
tempvar tmpln
gen `tmpln' = length(`leftLB`lcolsN'')
qui summ `tmpln' if `use' != 0
local otherln = r(max)
drop `tmpln'
// NOW HAVE LENGTH OF TITLE AND MAX LENGTH OF VARIABLE
local spread = int(`titleln'/`otherln')+1
if `spread'>4{
local spread = 4
}
local line = 1
local end = 0
local count = -1
local c2 = -2
local first = word("`colName'",1)
local last = word("`colName'",`count')
local nextlast = word("`colName'",`c2')
while `end' == 0{
replace `leftLB`lcolsN'' = "`last'" + " " + `leftLB`lcolsN'' in `Nnew`line''
local check = `leftLB`lcolsN''[`Nnew`line''] + " `nextlast'" // what next will be
local count = `count'-1
local last = word("`colName'",`count')
if "`last'" == ""{
local end = 1
}
if length(`leftLB`lcolsN''[`Nnew`line'']) > `titleln'/`spread' | ///
length("`check'") > `titleln'/`spread' & "`first'" == "`nextlast'"{
if `end' == 0{
local line = `line'+1
}
}
}
if `line' > `maxline'{
local maxline = `line'
}
mac shift
}
}
if `titleOff' == 1 {
forvalues i = 1/4{
replace `leftLB1' = "" in `Nnew`i'' // get rid of horrible __var name
}
}
replace `leftLB1' = `label' if `use' != 1 & `use' != 2 // put titles back in (overall, sub est etc.)
// STUFF ADDED FOR JS TO INCLUDE EFFICACY AS COLUMN WITH OVERALL
*effect lci uci tempvars
if "$MA_efficacy" != ""{
tempvar vetemp ucivetemp lcivetemp vaccine_efficacy
qui {
gen `vetemp'=100*(1-exp(`effect'))
tostring `vetemp', replace force format(%4.0f)
gen `ucivetemp'=100*(1-exp(`lci'))
tostring `ucivetemp', replace force format(%4.0f)
gen `lcivetemp'=100*(1-exp(`uci'))
tostring `lcivetemp', replace force format(%4.0f)
gen str30 `vaccine_efficacy'=`vetemp'+" ("+`lcivetemp'+", "+`ucivetemp'+")" if `effect' != .
label var `vaccine_efficacy' "Vaccine efficacy (%)"
local rcols = "`vaccine_efficacy' " + "`rcols' "
}
}
if "`wt'" == ""{
local rcols = "`weightText' " + "`rcols'"
if "$MA_method2" != ""{
label var `weightText' "% Weight ($MA_method1)"
}
else{
label var `weightText' "% Weight"
}
}
if "`counts'" != ""{
local rcols = "`raw1' " + "`raw2' " + "`rcols'"
}
if "`stats'" == ""{
local rcols = "`estText' " + "`rcols'"
if "$MA_ESLA" == ""{
global MA_ESLA = "`sumstat'"
}
label var `estText' "$MA_ESLA ($IND% CI)"
}
tempvar extra
gen `extra' = ""
label var `extra' " "
local rcols = "`rcols' `extra'"
local rcolsN = 0
if "`rcols'" != ""{
tokenize "`rcols'"
local rcolsN = 0
while "`1'" != ""{
cap confirm var `1'
if _rc!=0 {
di in re "Variable `1' not defined"
exit _rc
}
local rcolsN = `rcolsN' + 1
tempvar right`rcolsN' rightLB`rcolsN' rightWD`rcolsN'
cap confirm string var `1'
if _rc == 0{
gen str `rightLB`rcolsN'' = `1'
}
if _rc != 0{
local f: format `1'
gen str `rightLB`rcolsN'' = string(`1', "`f'")
replace `rightLB`rcolsN'' = "" if `rightLB`rcolsN'' == "."
}
local colName: variable label `1'
if "`colName'"==""{
local colName = "`1'"
}
// WORK OUT IF TITLE IS BIGGER THAN THE VARIABLE
// SPREAD OVER UP TO FOUR LINES IF NECESSARY
local titleln = length("`colName'")
tempvar tmpln
gen `tmpln' = length(`rightLB`rcolsN'')
qui summ `tmpln' if `use' != 0
local otherln = r(max)
drop `tmpln'
// NOW HAVE LENGTH OF TITLE AND MAX LENGTH OF VARIABLE
local spread = int(`titleln'/`otherln')+1
if `spread'>4{
local spread = 4
}
local line = 1
local end = 0
local count = -1
local c2 = -2
local first = word("`colName'",1)
local last = word("`colName'",`count')
local nextlast = word("`colName'",`c2')
while `end' == 0{
replace `rightLB`rcolsN'' = "`last'" + " " + `rightLB`rcolsN'' in `Nnew`line''
local check = `rightLB`rcolsN''[`Nnew`line''] + " `nextlast'" // what next will be
local count = `count'-1
local last = word("`colName'",`count')
if "`last'" == ""{
local end = 1
}
if length(`rightLB`rcolsN''[`Nnew`line'']) > `titleln'/`spread' | ///
length("`check'") > `titleln'/`spread' & "`first'" == "`nextlast'"{
if `end' == 0{
local line = `line'+1
}
}
}
if `line' > `maxline'{
local maxline = `line'
}
mac shift
}
}
// now get rid of extra title rows if they weren't used
if `maxline'==3{
drop in `Nnew4'
}
if `maxline'==2{
drop in `Nnew3'/`Nnew4'
}
if `maxline'==1{
drop in `Nnew2'/`Nnew4'
}
/* BODGE SOLU- EXTRA COLS */
while `rcolsN' < 2{
local rcolsN = `rcolsN' + 1
tempvar right`rcolsN' rightLB`rcolsN' rightWD`rcolsN'
gen str `rightLB`rcolsN'' = " "
}
local skip = 1
if "`stats'" == "" & "`wt'" == ""{ // sort out titles for stats and weight, if there
local skip = 3
}
if "`stats'" != "" & "`wt'" == ""{
local skip = 2
}
if "`stats'" == "" & "`wt'" != ""{
local skip = 2
}
if "`counts'" != ""{
local skip = `skip' + 2
}
if "$MA_efficacy" != ""{
local skip = `skip' + 1
}
/* SET TWO DUMMY RCOLS IF NOSTATS NOWEIGHT */
forvalues i = `skip'/`rcolsN'{ // get rid of junk if not weight, stats or counts
replace `rightLB`i'' = "" if (`use' != 1 & `use' != 2)
}
forvalues i = 1/`rcolsN'{
replace `rightLB`i'' = "" if (`use' ==0)
}
local leftWDtot = 0
local rightWDtot = 0
local leftWDtotNoTi = 0
forvalues i = 1/`lcolsN'{
getWidth `leftLB`i'' `leftWD`i''
qui summ `leftWD`i'' if `use' != 0 & `use' != 4 & `use' != 3 & `use' != 5 & ///
`use' != 17 & `use' != 19 // DON'T INCLUDE OVERALL STATS AT THIS POINT
local maxL = r(max)
local leftWDtotNoTi = `leftWDtotNoTi' + `maxL'
replace `leftWD`i'' = `maxL'
}
tempvar titleLN // CHECK IF OVERALL LENGTH BIGGER THAN REST OF LCOLS
getWidth `leftLB1' `titleLN'
qui summ `titleLN' if `use' != 0 & `use' != 4
local leftWDtot = max(`leftWDtotNoTi', r(max))
forvalues i = 1/`rcolsN'{
getWidth `rightLB`i'' `rightWD`i''
qui summ `rightWD`i'' if `use' != 0 & `use' != 4
replace `rightWD`i'' = r(max)
local rightWDtot = `rightWDtot' + r(max)
}
// CHECK IF NOT WIDE ENOUGH (I.E., OVERALL INFO TOO WIDE)
// LOOK FOR EDGE OF DIAMOND summ `lci' if `use' == ...
tempvar maxLeft
getWidth `leftLB1' `maxLeft'
qui count if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
if r(N) > 0{
summ `maxLeft' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19 // NOT TITLES THOUGH!
local max = r(max)
if `max' > `leftWDtotNoTi'{
// WORK OUT HOW FAR INTO PLOT CAN EXTEND
// WIDTH OF LEFT COLUMNS AS FRACTION OF WHOLE GRAPH
local x = `leftWDtot'*($MA_AS_TEXT/100)/(`leftWDtot'+`rightWDtot')
tempvar y
// SPACE TO LEFT OF DIAMOND WITHIN PLOT (FRAC OF GRAPH)
gen `y' = ((100-$MA_AS_TEXT)/100)*(`lci'-`DXmin') / (`DXmax'-`DXmin')
qui summ `y' if `use' == 3 | `use' == 5
local extend = 1*(r(min)+`x')/`x'
local leftWDtot = max(`leftWDtot'/`extend',`leftWDtotNoTi') // TRIM TO KEEP ON SAFE SIDE
// ALSO MAKE SURE NOT LESS THAN BEFORE!
}
}
global LEFT_WD = `leftWDtot'
global RIGHT_WD = `rightWDtot'
local ratio = $MA_AS_TEXT // USER SPECIFIED- % OF GRAPH TAKEN BY TEXT (ELSE NUM COLS CALC?)
local textWD = (`DXwidth'/(1-`ratio'/100)-`DXwidth') /(`leftWDtot'+`rightWDtot')
forvalues i = 1/`lcolsN'{
gen `left`i'' = `DXmin' - `leftWDtot'*`textWD'
local leftWDtot = `leftWDtot'-`leftWD`i''
}
gen `right1' = `DXmax'
forvalues i = 2/`rcolsN'{
local r2 = `i'-1
gen `right`i'' = `right`r2'' + `rightWD`r2''*`textWD'
}
local AXmin = `left1'
local AXmax = `DXmax' + `rightWDtot'*`textWD'
foreach type in "" "inf"{
replace `tauLCI`inf'' = `DXmin' if `tauLCI' < `DXmin' & `tauLCI`inf'' != .
replace `tauLCI`inf'' = . if `lci' < `DXmin'
replace `tauLCI`inf'' = . if `tauLCI`inf'' > `lci'
replace `tauUCI`inf'' = `DXmax' if `tauUCI`inf'' > `DXmax' & `tauUCI`inf'' != .
replace `tauUCI`inf'' = . if `uci' > `DXmax'
replace `tauUCI`inf'' = . if `tauUCI`inf'' < `uci'
//replace `tauLCI`inf'' = . if (`use' == 3 | `use' == 5) & "$MA_method1" != "D+L"
//replace `tauUCI`inf'' = . if (`use' == 3 | `use' == 5) & "$MA_method1" != "D+L"
//replace `tauLCI`inf'' = . if (`use' == 17 | `use' == 19) & "$MA_method2" != "D+L"
//replace `tauUCI`inf'' = . if (`use' == 17 | `use' == 19) & "$MA_method2" != "D+L"
replace `tauLCI`inf'' = . if (`use' == 3 | `use' == 5) & "$MA_method1" == "FE"
replace `tauUCI`inf'' = . if (`use' == 3 | `use' == 5) & "$MA_method1" == "FE"
replace `tauLCI`inf'' = . if (`use' == 17 | `use' == 19) & "$MA_method2" == "FE"
replace `tauUCI`inf'' = . if (`use' == 17 | `use' == 19) & "$MA_method2" == "FE"
}
// DIAMONDS TAKE FOREVER...I DON'T THINK THIS IS WHAT MIKE DID
tempvar DIAMleftX DIAMrightX DIAMbottomX DIAMtopX DIAMleftY1 DIAMrightY1 DIAMleftY2 DIAMrightY2 DIAMbottomY DIAMtopY
gen `DIAMleftX' = `lci' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
replace `DIAMleftX' = `DXmin' if `lci' < `DXmin' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMleftX' = . if `effect' < `DXmin' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
gen `DIAMleftY1' = `id' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
replace `DIAMleftY1' = `id' + 0.4*( abs((`DXmin'-`lci')/(`effect'-`lci')) ) if `lci' < `DXmin' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMleftY1' = . if `effect' < `DXmin' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
gen `DIAMleftY2' = `id' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
replace `DIAMleftY2' = `id' - 0.4*( abs((`DXmin'-`lci')/(`effect'-`lci')) ) if `lci' < `DXmin' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMleftY2' = . if `effect' < `DXmin' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
gen `DIAMrightX' = `uci' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
replace `DIAMrightX' = `DXmax' if `uci' > `DXmax' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMrightX' = . if `effect' > `DXmax' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
gen `DIAMrightY1' = `id' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
replace `DIAMrightY1' = `id' + 0.4*( abs((`uci'-`DXmax')/(`uci'-`effect')) ) if `uci' > `DXmax' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMrightY1' = . if `effect' > `DXmax' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
gen `DIAMrightY2' = `id' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
replace `DIAMrightY2' = `id' - 0.4*( abs((`uci'-`DXmax')/(`uci'-`effect')) ) if `uci' > `DXmax' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMrightY2' = . if `effect' > `DXmax' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
gen `DIAMbottomY' = `id' - 0.4 if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
replace `DIAMbottomY' = `id' - 0.4*( abs((`uci'-`DXmin')/(`uci'-`effect')) ) if `effect' < `DXmin' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMbottomY' = `id' - 0.4*( abs((`DXmax'-`lci')/(`effect'-`lci')) ) if `effect' > `DXmax' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
gen `DIAMtopY' = `id' + 0.4 if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
replace `DIAMtopY' = `id' + 0.4*( abs((`uci'-`DXmin')/(`uci'-`effect')) ) if `effect' < `DXmin' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMtopY' = `id' + 0.4*( abs((`DXmax'-`lci')/(`effect'-`lci')) ) if `effect' > `DXmax' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
gen `DIAMtopX' = `effect' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19
replace `DIAMtopX' = `DXmin' if `effect' < `DXmin' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMtopX' = `DXmax' if `effect' > `DXmax' & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
replace `DIAMtopX' = . if (`uci' < `DXmin' | `lci' > `DXmax') & (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19)
gen `DIAMbottomX' = `DIAMtopX'
} // END QUI
// v1.11 TEXT SIZE SOLU
// v1.16 TRYING AGAIN!
// IF aspect IS USED IN "$MA_OTHEROPTS" (OTHER GRAPH OPTS) THEN THIS HELPS TO CALCULATE TEXT SIZE
// IF NO ASPECT, BUT xsize AND ysize USED THEN FIND RATIO MANUALLY
// STATA ALWAYS TRIES TO PRODUCE A GRAPH WITH ASPECT ABOUT 0.77 - TRY TO FIND "NATURAL ASPECT"
local aspect = .
if strpos(`"$MA_OTHEROPTS"',"aspect") > 0{
local aspectTXT = substr( `"$MA_OTHEROPTS"', (strpos(`"$MA_OTHEROPTS"',"aspect")), (length(`"$MA_OTHEROPTS"')) )
local aspectTXT = substr( "`aspectTXT'", 1, ( strpos("`aspectTXT'",")")) )
local aspect = real( substr( "`aspectTXT'", ( strpos("`aspectTXT'","(") +1 ), ///
( strpos("`aspectTXT'",")") - strpos("`aspectTXT'","(") -1 ) ))
}
if strpos(`"$MA_OTHEROPTS"',"xsize") > 0 ///
& strpos(`"$MA_OTHEROPTS"',"ysize") > 0 ///
& strpos(`"$MA_OTHEROPTS"',"aspect") == 0{
local xsizeTXT = substr( `"$MA_OTHEROPTS"', (strpos(`"$MA_OTHEROPTS"',"xsize")), (length(`"$MA_OTHEROPTS"')) )
// Ian White's bug fix!
local xsizeTXT = substr( `"`xsizeTXT'"', 1, ( strpos(`"`xsizeTXT'"',")")) )
local xsize = real( substr( `"`xsizeTXT'"', ( strpos(`"`xsizeTXT'"',"(") +1 ), ///
( strpos(`"`xsizeTXT'"',")") - strpos(`"`xsizeTXT'"',"(") -1 ) ))
local ysizeTXT = substr( `"$MA_OTHEROPTS"', (strpos(`"$MA_OTHEROPTS"',"ysize")), (length(`"$MA_OTHEROPTS"')) )
local ysizeTXT = substr( `"`ysizeTXT'"', 1, ( strpos(`"`ysizeTXT'"',")")) )
local ysize = real( substr( `"`ysizeTXT'"', ( strpos(`"`ysizeTXT'"',"(") +1 ), ///
( strpos(`"`ysizeTXT'"',")") - strpos(`"`ysizeTXT'"',"(") -1 ) ))
local aspect = `ysize'/`xsize'
}
local approx_chars = ($LEFT_WD + $RIGHT_WD)/($MA_AS_TEXT/100)
qui count if `use' != 9
local height = r(N)
local natu_aspect = 1.3*`height'/`approx_chars'
if `aspect' == .{
// sort out relative to text, but not to ridiculous degree
local new_asp = 0.5*`natu_aspect' + 0.5*1
global MA_OTHEROPTS `"$MA_OTHEROPTS aspect(`new_asp')"'
local aspectRat = max( `new_asp'/`natu_aspect' , `natu_aspect'/`new_asp' )
}
if `aspect' != .{
local aspectRat = max( `aspect'/`natu_aspect' , `natu_aspect'/`aspect' )
}
local adj = 1.25
if `natu_aspect' > 0.7{
local adj = 1/(`natu_aspect'^1.3+0.2)
}
local textSize = `adj' * $MA_TEXT_SCA / (`approx_chars' * sqrt(`aspectRat') )
local textSize2 = `adj' * $MA_TEXT_SCA / (`approx_chars' * sqrt(`aspectRat') )
forvalues i = 1/`lcolsN'{
local lcolCommands`i' "(scatter `id' `left`i'' if `use' != 4, msymbol(none) mlabel(`leftLB`i'') mlabcolor(black) mlabpos(3) mlabsize(`textSize')) "
}
forvalues i = 1/`rcolsN'{
local rcolCommands`i' "(scatter `id' `right`i'' if `use' != 4, msymbol(none) mlabel(`rightLB`i'') mlabcolor(black) mlabpos(3) mlabsize(`textSize')) "
}
if "$MA_rfdist" != ""{
if "`stats'" == ""{
local predIntCmd "(scatter `id' `right1' if `use' == 4, msymbol(none) mlabel(`RFdistText') mlabcolor(black) mlabpos(3) mlabsize(`textSize')) "
}
if "$MA_nohet" == ""{
local predIntCmd2 "(scatter `id' `left1' if `use' == 4, msymbol(none) mlabel(`RFdistLabel') mlabcolor(black) mlabpos(3) mlabsize(`textSize')) "
}
}
if "$MA_nohet" == "" & "$MA_rjhby" != ""{
local hetGroupCmd "(scatter `id' `left1' if `use' == 4, msymbol(none) mlabel(`hetGroupLabel') mlabcolor(black) mlabpos(3) mlabsize(`textSize')) "
}
// OTHER BITS AND BOBS
local dispBox "none"
if "`nobox'" == ""{
local dispBox "square "
}
local boxsize = $MA_FBSC/150
if "$MA_FAVOURS" != ""{
local pos = strpos("$MA_FAVOURS", "#")
local leftfav = substr("$MA_FAVOURS",1,(`pos'-1))
local rightfav = substr("$MA_FAVOURS",(`pos'+1),(length("$MA_FAVOURS")-`pos'+1) )
}
if `h0' != . & "$MA_nulloff" == ""{
local leftfp = `DXmin' + (`h0'-`DXmin')/2
local rightfp = `h0' + (`DXmax'-`h0')/2
}
else{
local leftfp = `DXmin'
local rightfp = `DXmax'
}
// GRAPH APPEARANCE OPTIONS- ADDED v1.15
/*
if `"$MA_OPT"' != "" & strpos(`"$MA_OPT"',"m") == 0{(
global MA_OPT = `"$MA_OPT m()"'
}
*/
if `"$MA_BOXOPT"' != "" & strpos(`"$MA_BOXOPT"',"msymbol") == 0{ // make defaults if unspecified
global MA_BOXOPT = `"$MA_BOXOPT msymbol(square)"'
}
if `"$MA_BOXOPT"' != "" & strpos(`"$MA_BOXOPT"',"mcolor") == 0{ // make defaults if unspecified
global MA_BOXOPT = `"$MA_BOXOPT mcolor("180 180 180")"'
}
if `"$MA_BOXOPT"' == ""{
local boxopt "msymbol(`dispBox') msize(`boxsize') mcolor("180 180 180")"
}
else{
if strpos(`"$MA_BOXOPT"',"mla") != 0{
di as error "Option mlabel() not allowed in boxopt()"
exit
}
if strpos(`"$MA_BOXOPT"',"msi") != 0{
di as error "Option msize() not allowed in boxopt()"
exit
}
local boxopt `"$MA_BOXOPT msize(`boxsize')"'
}
if "$MA_classic" != ""{
local boxopt "mcolor(black) msymbol(square) msize(`boxsize')"
}
if "`box'" != ""{
local boxopt "msymbol(none)"
}
if `"$MA_DIAMOPT"' == ""{
local diamopt "lcolor("0 0 100")"
}
else{
if strpos(`"$MA_DIAMOPT"',"hor") != 0 | strpos(`"$MA_DIAMOPT"',"vert") != 0{
di as error "Options horizontal/vertical not allowed in diamopt()"
exit
}
if strpos(`"$MA_DIAMOPT"',"con") != 0{
di as error "Option connect() not allowed in diamopt()"
exit
}
if strpos(`"$MA_DIAMOPT"',"lp") != 0{
di as error "Option lpattern() not allowed in diamopt()"
exit
}
local diamopt `"$MA_DIAMOPT"'
}
if `"$MA_POINTOPT"' != "" & strpos(`"$MA_POINTOPT"',"msymbol") == 0{(
global MA_POINTOPT = `"$MA_POINTOPT msymbol(diamond)"'
}
if `"$MA_POINTOPT"' != "" & strpos(`"$MA_POINTOPT"',"msize") == 0{(
global MA_POINTOPT = `"$MA_POINTOPT msize(vsmall)"'
}
if `"$MA_POINTOPT"' != "" & strpos(`"$MA_POINTOPT"',"mcolor") == 0{(
global MA_POINTOPT = `"$MA_POINTOPT mcolor(black)"'
}
if `"$MA_POINTOPT"' == ""{
local pointopt "msymbol(diamond) msize(vsmall) mcolor("0 0 0")"
}
else{
local pointopt `"$MA_POINTOPT"'
}
if "$MA_classic" != "" & "`box'" == ""{
local pointopt "msymbol(none)"
}
if `"$MA_CIOPT"' != "" & strpos(`"$MA_CIOPT"',"lcolor") == 0{(
global MA_CIOPT = `"$MA_CIOPT lcolor(black)"'
}
if `"$MA_CIOPT"' == ""{
local ciopt "lcolor("0 0 0")"
}
else{
if strpos(`"$MA_CIOPT"',"hor") != 0 | strpos(`"$MA_CIOPT"',"vert") != 0{
di as error "Options horizontal/vertical not allowed in ciopt()"
exit
}
if strpos(`"$MA_CIOPT"',"con") != 0{
di as error "Option connect() not allowed in ciopt()"
exit
}
if strpos(`"$MA_CIOPT"',"lp") != 0{
di as error "Option lpattern() not allowed in ciopt()"
exit
}
local ciopt `"$MA_CIOPT"'
}
// END GRAPH OPTS
//if "$MA_method1" == "D+L"{
if "$MA_method1" != "FE"{
tempvar noteposx noteposy notelab
qui{
summ `id'
gen `noteposy' = r(min) -1.5 in 1
summ `left1'
gen `noteposx' = r(mean) in 1
gen `notelab' = "NOTE: Weights are from random effects analysis" in 1
local notecmd "(scatter `noteposy' `noteposx', msymbol(none) mlabel(`notelab') mlabcolor(black) mlabpos(3) mlabsize(`textSize')) "
if "$MA_method1" == "MH" {
replace `notelab' = "NOTE: Mantel-Haenszel fixed-effect weights" in 1
}
if "$MA_method1" == "PETO" {
replace `notelab' = "NOTE: Peto fixed-effect weights" in 1
}
}
if "$MA_nowarning" != ""{
local notecmd
}
}
if "`overall'" != ""{
local overallCommand ""
qui drop if `use' == 5
qui summ `id'
local DYmin = r(min)
cap replace `noteposy' = r(min) -.5 in 1
}
// quick bodge to get overall- can't find log version!
tempvar tempOv ovLine ovMin ovMax h0Line
qui gen `tempOv' = `effect' if `use' == 5
sort `tempOv'
qui summ `id'
local DYmin = r(min)-2
local DYmax = r(max)+1
qui gen `ovLine' = `tempOv' in 1
qui gen `ovMin' = r(min)-2 in 1
qui gen `ovMax' = $borderline in 1
qui gen `h0Line' = `h0' in 1
if `"$MA_OLINEOPT"' == ""{
local overallCommand " (pcspike `ovMin' `ovLine' `ovMax' `ovLine', lwidth(thin) lcolor(maroon) lpattern(shortdash)) "
}
else{
local overallCommand `" (pcspike `ovMin' `ovLine' `ovMax' `ovLine', $MA_OLINEOPT) "'
}
if `ovLine' > `DXmax' | `ovLine' < `DXmin' | "`overall'" != ""{ // ditch if not on graph
local overallCommand ""
}
local nullCommand " (pcspike `ovMin' `h0Line' `ovMax' `h0Line', lwidth(thin) lcolor(black) ) "
// gap if "favours" used
if "`leftfav'" != "" | "`rightfav'" != ""{
local gap = "labgap(5)"
}
// if summary only must not have weights
local awweight "[aw= `weight']"
if "$MA_summaryonly" != ""{
local awweight ""
}
qui summ `weight'
if r(N) == 0{
local awweight ""
}
// rfdist off scale arrows only used when appropriate
qui{
tempvar rfarrow
gen `rfarrow' = 0
if "$MA_rfdist" != ""{
//if "$MA_method1" == "D+L"{
if "$MA_method1" != "FE"{
replace `rfarrow' = 1 if `use' == 3 | `use' == 5
}
//if "$MA_method2" == "D+L"{
if "$MA_method2" != "FE"{
replace `rfarrow' = 1 if `use' == 17 | `use' == 19
}
}
} // end qui
// final addition- if aspect() given but not xsize() ysize(), put these in to get rid of gaps
// need to fiddle to allow space for bottom title
// should this just replace the aspect option?
// suppose good to keep- most people hopefully using xsize and ysize and can always change themselves if using aspect
if strpos(`"$MA_OTHEROPTS"',"xsize") == 0 & strpos(`"$MA_OTHEROPTS"',"ysize") == 0 ///
& strpos(`"$MA_OTHEROPTS"',"aspect") > 0 {
local aspct = substr(`"$MA_OTHEROPTS"', (strpos(`"$MA_OTHEROPTS"',"aspect(")+7 ) , length(`"$MA_OTHEROPTS"') )
local aspct = substr(`"`aspct'"', 1, (strpos(`"`aspct'"',")")-1) )
if `aspct' > 1{
local xx = (11.5+(2-2*1/`aspct'))/`aspct'
local yy = 12
}
if `aspct' <= 1{
local yy = 12*`aspct'
local xx = 11.5-(2-2*`aspct')
}
global MA_OTHEROPTS = `"$MA_OTHEROPTS"' + " xsize(`xx') ysize(`yy')"
}
// switch off null if wanted
if "$MA_nulloff" != ""{
local nullCommand ""
}
***************************
*** GRAPH ***
***************************
#delimit ;
twoway
/* NOTE FOR RF, OVERALL AND NULL LINES FIRST */
`notecmd' `overallCommand' `nullCommand' `predIntCmd' `predIntCmd2' `hetGroupCmd'
/* PLOT BOXES AND PUT ALL THE GRAPH OPTIONS IN THERE */
(scatter `id' `effect' `awweight' if `use' == 1,
`boxopt'
yscale(range(`DYmin' `DYmax') noline )
ylabel(none) ytitle("")
xscale(range(`AXmin' `AXmax'))
xlabel(`lblcmd', labsize(`textSize2') )
yline($borderline, lwidth(thin) lcolor(gs12))
/* THIS BIT DOES favours. NOTE SPACES TO SUPPRESS IF THIS IS NOT USED */
xmlabel(`leftfp' "`leftfav' " `rightfp' "`rightfav' ", noticks labels labsize(`textSize')
`gap' /* PUT LABELS UNDER xticks? Yes as labels now extended */ )
xtitle("") legend(off) xtick("`xtick'") )
/* END OF FIRST SCATTER */
/* HERE ARE THE CONFIDENCE INTERVALS */
(pcspike `id' `lci' `id' `uci' if `use' == 1, `ciopt')
/* ADD ARROWS IF OFFSCALE USING offLeftX offLeftX2 offRightX offRightX2 offYlo offYhi */
(pcspike `id' `offLeftX' `offYlo' `offLeftX2' if `use' == 1, `ciopt')
(pcspike `id' `offLeftX' `offYhi' `offLeftX2' if `use' == 1, `ciopt')
(pcspike `id' `offRightX' `offYlo' `offRightX2' if `use' == 1, `ciopt')
(pcspike `id' `offRightX' `offYhi' `offRightX2' if `use' == 1, `ciopt')
/* DIAMONDS FOR SUMMARY ESTIMATES -START FROM 9 O'CLOCK */
(pcspike `DIAMleftY1' `DIAMleftX' `DIAMtopY' `DIAMtopX' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19, `diamopt')
(pcspike `DIAMtopY' `DIAMtopX' `DIAMrightY1' `DIAMrightX' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19, `diamopt')
(pcspike `DIAMrightY2' `DIAMrightX' `DIAMbottomY' `DIAMbottomX' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19, `diamopt')
(pcspike `DIAMbottomY' `DIAMbottomX' `DIAMleftY2' `DIAMleftX' if `use' == 3 | `use' == 5 | `use' == 17 | `use' == 19, `diamopt')
/* EXTENDED CI FOR RANDOM EFFECTS, SHOW DISTRIBUTION AS RECOMMENDED BY JULIAN HIGGINS
DOTTED LINES FOR INESTIMABLE DISTRIBUTION */
(pcspike `DIAMleftY1' `DIAMleftX' `DIAMleftY1' `tauLCI' if (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19) & `tau2' < ., `diamopt')
(pcspike `DIAMrightY1' `DIAMrightX' `DIAMrightY1' `tauUCI' if (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19) & `tau2' < ., `diamopt')
(pcspike `DIAMleftY1' `DIAMleftX' `DIAMleftY1' `tauLCI' if (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19) & `tau2' ==.b, `diamopt' lpattern(shortdash))
(pcspike `DIAMrightY1' `DIAMrightX' `DIAMrightY1' `tauUCI' if (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19) & `tau2' ==.b, `diamopt' lpattern(shortdash))
/* DIAMOND EXTENSION FOR RF DIST ALSO HAS ARROWS... */
(pcspike `id' `offLeftX' `offYlo' `offLeftX2' if (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19) & `rfarrow' == 1, `diamopt')
(pcspike `id' `offLeftX' `offYhi' `offLeftX2' if (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19) & `rfarrow' == 1, `diamopt')
(pcspike `id' `offRightX' `offYlo' `offRightX2' if (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19) & `rfarrow' == 1, `diamopt')
(pcspike `id' `offRightX' `offYhi' `offRightX2' if (`use' == 3 | `use' == 5 | `use' == 17 | `use' == 19) & `rfarrow' == 1, `diamopt')
/* COLUMN VARIBLES */
`lcolCommands1' `lcolCommands2' `lcolCommands3' `lcolCommands4' `lcolCommands5' `lcolCommands6'
`lcolCommands7' `lcolCommands8' `lcolCommands9' `lcolCommands10' `lcolCommands11' `lcolCommands12'
`rcolCommands1' `rcolCommands2' `rcolCommands3' `rcolCommands4' `rcolCommands5' `rcolCommands6'
`rcolCommands7' `rcolCommands8' `rcolCommands9' `rcolCommands10' `rcolCommands11' `rcolCommands12'
(scatter `id' `right1' if `use' != 4 & `use' != 0,
msymbol(none) mlabel(`rightLB1') mlabcolor("0 0 0") mlabpos(3) mlabsize(`textSize'))
(scatter `id' `right2' if `use' != 4 & `use' != 0,
msymbol(none) mlabel(`rightLB2') mlabcolor("0 0 0") mlabpos(3) mlabsize(`textSize'))
/* (scatter `id' `right2', mlabel(`use')) JUNK, TO SEE WHAT'S WHERE */
/* LAST OF ALL PLOT EFFECT MARKERS TO CLARIFY AND OVERALL EFFECT LINE */
(scatter `id' `effect' if `use' == 1, `pointopt' )
, $MA_OTHEROPTS /* RMH added */ plotregion(margin(zero));
#delimit cr
end
program define getWidth
version 9.0
// ROSS HARRIS, 13TH JULY 2006
// TEXT SIZES VARY DEPENDING ON CHARACTER
// THIS PROGRAM GENERATES APPROXIMATE DISPLAY WIDTH OF A STRING
// FIRST ARG IS STRING TO MEASURE, SECOND THE NEW VARIABLE
// PREVIOUS CODE DROPPED COMPLETELY AND REPLACED WITH SUGGESTION
// FROM Jeff Pitblado
qui{
gen `2' = 0
count
local N = r(N)
forvalues i = 1/`N'{
local this = `1'[`i']
local width: _length "`this'"
replace `2' = `width' +1 in `i'
}
} // end qui
end