/* 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'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=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 (sval10) 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 (temp1){ if (abs(A[temp-1]-sval) `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