*! version 1.0.2 03jul2012 MJC /* History MJC 03jul2012 version 1.0.2 - bug fix, extra touse passed to Mata when addplot() is specified; addplot passed to dummy plots to aid legend(order()) use MJC 21sep2011 version 1.0.1 - bug fix: touse now passed to Mata correctly. Matrices of effect estimates, standard errors and status codes now returned in r(). Label option added to pass to metan. Scheme option added: grayscale and color allowed. legend(order(1 "" ...)) now synched. MJC 21aug2011 version 1.0.0 */ program define extfunnel, rclass version 11.0 syntax varlist(min=2 max=2) [if] [in], [ /// FIXEDI /// -Fixed effect MA using inverse-variance method- RANDOMI /// -Random effect MA using inverse-variance method- /// CPoints(string) /// -No. of points to evaluate contours; Defaults vary between fixed/random- NULL(real 0) /// -Value of null hypothesis- Isquared(numlist asc max=5) /// -I-squared contours- TAUsquared(numlist asc max=5) /// -Tau-squared contours- /// MEASure(string) /// -One of lci/uci/ciwidth- LOE(numlist asc min=2 max=2) /// -Limits of equivalence- LOELine /// -Display limits of clinical equivalence on the plot- /// NEWSTUDYTreatment(string) /// -No. of patients in the treatment group of the new study- NEWSTUDYControl(string) /// -No. of patients in the control group of the new study- OR /// -generate log odds ratios when doing "peacock plot" RR /// -generate log risk ratios when doing "peacock plot" /// XRANGE(numlist asc max=2) /// -Range of effect estimates to evaluate- YRANGE(numlist asc max=2) /// -Range of standard errors to evaluate- SUMD /// -Suppress displaying summary diamond- SUMDPosition(real 0) /// -Vertical position to put summary diamond- PREDiction /// -Show prediction interval on twoway graph- NONULLline /// -Suppress displaying vertical line of no effect- NOPOOLEDline /// -Suppress displaying vertical line of pooled estimate- NOSHADing /// -Suppress displaying shaded regions- NOSCATter /// -Suppress displaying scatter of original meta-analysis- NOMETAN /// -Suppress display of metan output- EFORM /// -Exponentiate effect estimates from metan and in graph- ADDPlot(string) /// -Add an additional twoway graph- Level(cilevel) /// -Statistical significance level- LABel(string) /// -To pass to metan- SCHEME(string) /// -Graphical colour scheme: either grayscale or color- /// * /// -Twoway options- ] tokenize `varlist' marksample touse global touse `touse' /*****************************************************************************************************************************************************/ /* ERROR CHECKS AND DEFAULTS */ capture which metan if _rc >0 { display in yellow "You need to install the command metan. This can be installed using," display in yellow ". {stata ssc install metan}" exit 198 } local nmethods : word count `fixedi' `randomi' if (`nmethods'>1) { di as error "Can only specify one of fixedi or randomi" exit 198 } local nmethods2 : word count `or' `rr' if (`nmethods2'>1) { di as error "Can only specify one of or and rr" exit 198 } if "`or'"=="" & "`rr'"=="" { local stat = "or" } else { local stat = "`or' `rr'" } if (`nmethods'==0) { local method "fixedi" } else { local method = trim("`fixedi' `randomi'") } if "`newstudytreatment'"=="" & "`newstudycontrol'"!="" { di as error "newstudytreatment must be specified as well as newstudycontrol" exit } if "`newstudytreatment'"!="" & "`newstudycontrol'"=="" { di as error "newstudycontrol must be specified as well as newstudytreatment" exit 198 } if "`measure'"!="" & "`loe'"!="" { di as error "Cannot specify both measure and loe" exit 198 } if "`cpoints'" != "" { capture confirm integer number `cpoints' if _rc>0 { display in red "cpoints option must be an integer" exit 198 } } if "`method'"=="fixedi" & "`cpoints'"=="" & "`measure'"=="" & "`loe'"=="" { local cpoints = 3500 } if "`method'"=="randomi" | "`measure'"!="" | "`loe'"!="" { if "`cpoints'"=="" { local cpoints = 100 scalar contpoints = 100 } else { scalar contpoints = `cpoints' } } if "`loe'"!="" | "`measure'"!="" { if `cpoints'>500 { di as error "Maximum number of contour points is 500 when loe/measure is specified" exit 198 } } if "`method'"=="randomi" & "`cpoints'"!="" & `cpoints'>500 { di as error "Maximum number of contour points is 500 when fixed/random/randomi method is specified" exit 198 } qui su `2' if `touse' if (`r(min)'<0) { di as error "Variable `2' is assumed to contain standard errors - a negative value has been found" exit 198 } if ("`isquared'"!="" & "`tausquared'"!="") { di as error "I-squared and Tau-squared contours cannot be specified together" exit 198 } if ("`newstudytreatment'" != "" & "`newstudycontrol'" == "") | ("`newstudytreatment'" == "" & "`newstudycontrol'" != "") { di as error "Both newstudytreatment and newstudycontrol must be specified" exit 198 } if "`newstudytreatment'" != "" & "`newstudycontrol'" != "" { capture confirm integer number `newstudytreatment' if _rc>0 { display as error "newstudytreatment option must be an integer" exit 198 } capture confirm integer number `newstudycontrol' if _rc>0 { display as error "newstudycontrol option must be an integer" exit 198 } } if "`measure'"!="" { if "`measure'"!="lci" & "`measure'"!="uci" & "`measure'"!="ciwidth" { di as error "measure must be one of lci, uci or ciwidth" exit 198 } } if "`loe'"!="" { local lind = 1 foreach limit in `loe' { scalar loe`lind' = `limit' local `++lind' } if `null'loe2 { di as error "loe region must contain null" exit 198 } } if ("`tausquared'"!="" | "`isquared'"!="") & "`eform'"!="" { di as error "Cannot use eform option with tausquared or isquared contours" exit 198 } if "`measure'"=="ciwidth" & `null'<=0 { di as error "Null when measure = ciwidth must be > 0" exit 198 } if "`loeline'"!="" & "`loe'"=="" { di as error "loeline can only be specified with loe" exit 198 } if "`noshading'"!="" & ("`measure'"!="" | "`loe'"!="" | "`newstudytreatment'"!="") { di as error "noshading option cannot be used" exit 198 } /*****************************************************************************************************************************************************/ /* Colour scheme */ if "`scheme'"=="" { local colscheme "grayscale" } else { local l = length("`scheme'") if substr("grayscale",1,max(1,`l')) == "`scheme'" { local colscheme "grayscale" } else if substr("color",1,max(1,`l')) == "`scheme'" { local colscheme "color" } else { di as error "Unknown scheme" exit 198 } } if "`loe'"=="" { if "`colscheme'"=="grayscale" { local col1 "gray*0.3" local col2 "gray*0.8" local col3 "gray*1.2" } else { local col1 "green" local col2 "blue" local col3 "red" } } else { if "`colscheme'"=="grayscale" { local col1 "gray*0.1" local col2 "gray*0.4" local col3 "gray*0.7" local col4 "gray*1" local col5 "gray*1.3" local col6 "gray*1.6" local col7 "gray*1.9" local col8 "gray*2.2" } else { local col1 "green" local col2 "orange" local col3 "blue" local col4 "yellow" local col5 "blue*0.5" local col6 "green*0.5" local col7 "red" local col8 "purple" } } quietly { /*****************************************************************************************************************************************************/ /* PRELIMINARY CALCULATIONS */ preserve local xtitle "Effect Estimate" /* Default x-title */ local ytitle "Standard Error" /* Default y-title */ /* Summarise effect estimates and standard errors and extract: max, min, range and no. of studies */ tokenize `varlist' su `1' if `touse', mean local N_SS = `r(N)' scalar nstudies=`r(N)' local minSS = `r(min)' local maxSS = `r(max)' local rangeSS = (`maxSS') - (`minSS') su `2' if `touse', mean local N_seSS = `r(N)' local minseSS = `r(min)' local maxseSS = `r(max)' local rangeseSS = (`maxseSS') - (`minseSS') local Nobs = _N local N_new = `Nobs'+1 local siglev = abs(invnormal((100-`level')/200)) /* Critical value */ scalar siglev = `siglev' /* Need scalars for mata program */ scalar null = `null' /*** Conduct AD meta-analysis using metan ***/ if "`nometan'"!=""{ local show "quietly" } } if "`label'"!="" { local metanlab "label(`label')" } `show' di as txt "Original meta-analysis results:" `show' metan `1' `2' if `touse', `method' nograph nokeep olevel(`level') `eform' `metanlab' quietly { if "`eform'"=="" { local pooled_est = `r(ES)' /* Pooled estimate */ local tau2= `r(tau2)' /* Estimat of tau2 */ local se_pooled_est = `r(seES)' /* Pooled standard error */ } else{ local pooled_est = log(`r(ES)') /* Pooled estimate */ local tau2= `r(tau2)' /* Estimat of tau2 */ local se_pooled_est = `r(selogES)' /* Pooled standard error */ } /*** Define x/y min/max for area ranges, defaults are scatter plot plus a bit. ***/ local xcount : word count `xrange' local ycount : word count `yrange' tempname xmin xmax ymin ymax if (`xcount' > 0){ tokenize `xrange' scalar `xmin' =`1' scalar `xmax' =`2' } else { scalar `xmin' = `minSS' - (0.5*(`rangeSS')) scalar `xmax' = `maxSS' + (0.5*(`rangeSS')) } if (`ycount' > 0){ tokenize `yrange' scalar `ymin' = 1.0E-06 scalar `ymax' = `2' } else { scalar `ymin' = 1.0E-06 scalar `ymax' = `maxseSS' + 0.5*`rangeseSS' } if (`ymin'<=0){ scalar `ymin' = 1.0E-06 di "Minimum of y must be > 0" } /*****************************************************************************************************************************************************/ /* For mata */ if "`measure'"=="lci" { scalar power = 1 } else if "`measure'"=="uci" { scalar power = 2 } else if "`measure'"=="ciwidth" { scalar power = 3 } else if "`loe'"!="" { scalar power = 4 } else { scalar power = 0 } if "`eform'"!="" { local explabel "xscale(log)" } scalar ymin = `ymin' scalar ymax = `ymax' scalar xmin = `xmin' scalar xmax = `xmax' /*****************************************************************************************************************************************************/ /* Generate basis variables to input into formula */ tempvar tempSS tempseSS range2 `tempSS' `xmin' `xmax' `cpoints' /* range2 creates double formatted variables */ range2 `tempseSS' `ymin' `ymax' `cpoints' local xwidth = `xmax' - `xmin' local ywidth = `ymax' - `ymin' if ("`method'" != "fixedi" | "`measure'"!="" | "`loe'"!="") { /* Return matrices of effect estimates and standard errors */ tempname matSS matseSS mkmat `tempSS', matrix(`matSS') nomissing mkmat `tempseSS' , matrix(`matseSS') nomissing return matrix ESmat = `matSS' return matrix seESmat = `matseSS' } /*****************************************************************************************************************************************************/ /* Dummy plots for legend */ /* Begin twoway label index and label */ local labnum 0 if "`noshading'"=="" & "`newstudytreatment'"=="" { *coords for full graph box shaded area tempvar x2 y2 range `x2' `xmin' `xmax' 2 if "`eform'"!="" { replace `x2' = exp(`x2') } gen `y2'=`ymax' in 1/2 *dummy plots for legend if (power==0) { local labnum `++labnum' local dumplots "(area `y2' `x2', col(`col1') fi(inten100))" local lab `"`lab' `labnum' "Non-sig. effect (`=100-`level''% level)""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col2') fi(inten100))" local lab `"`lab' `labnum' "Sig. effect > NULL (`=100-`level''% level)""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col3') fi(inten100))" local lab `"`lab' `labnum' "Sig. effect < NULL (`=100-`level''% level)""' } else if (power==1) { local labnum `++labnum' local dumplots "(area `y2' `x2', col(`col1') fi(inten100))" local lab `"`lab' `labnum' "Lower CI < `null'""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col3') fi(inten100))" local lab `"`lab' `labnum' "Lower CI >= `null'""' } else if (power==2) { local labnum `++labnum' local dumplots "(area `y2' `x2', col(`col1') fi(inten100))" local lab `"`lab' `labnum' "Upper CI < `null'""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col3') fi(inten100))" local lab `"`lab' `labnum' "Upper CI >= `null'""' } else if (power==3) { local labnum `++labnum' local dumplots "(area `y2' `x2', col(`col1') fi(inten100))" local lab `"`lab' `labnum' "CI width < `null'""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col3') fi(inten100))" local lab `"`lab' `labnum' "CI width >= `null'""' } else if (power==4) { local labnum `++labnum' local dumplots "(area `y2' `x2', col(`col1') fi(inten100))" local lab `"`lab' `labnum' "1. Insufficient evidence to" "confirm or exclude an" "important difference""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col2') fi(inten100))" local lab `"`lab' `labnum' "2. Stat. sig. benefit," "clinical benefit unclear""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col3') fi(inten100))" local lab `"`lab' `labnum' "3. Stat. sig. harm," "clinical harm unclear""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col4') fi(inten100))" local lab `"`lab' `labnum' "4. Stat. sig. benefit," "no clinical benefit""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col5') fi(inten100))" local lab `"`lab' `labnum' "5. Stat. sig. harm," "no clinical harm""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col6') fi(inten100))" local lab `"`lab' `labnum' "6. No evidence of" "an important difference""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col7') fi(inten100))" local lab `"`lab' `labnum' "7. Important benefit""' local labnum `++labnum' local dumplots "`dumplots'(area `y2' `x2', col(`col8') fi(inten100))" local lab `"`lab' `labnum' "8. Important harm""' local extraleg "pos(3) col(1) size(small)" } if (power!=4) { local basearea "(area `y2' `x2', base(0) col(`col1') fi(inten100))" } else { local basearea "(area `y2' `x2', col(`col1') base(0) fi(inten100))" } } /*****************************************************************************************************************************************************/ /* I-SQUARED CONTOURS */ if ("`isquared'"!="" & "`tausquared'"=="") { tokenize `varlist' /*replace `1' = `pooled_est' in `=`N_SS'+1' replace `2' = 0.1 in `=`N_SS'+1' metan `1' `2', `method'i nograph nokeep level(`level') noisily di as txt " Minimum I-squared =" %8.4f `r(i_sq)'*/ *prelim calcs tempvar eval1gen eval2gen eval3gen gen `eval1gen' = ((`1')^2)/((`2')^2) su `eval1gen' if `touse' local eval1 = `r(sum)' gen `eval2gen' = (`1')/((`2')^2) su `eval2gen' if `touse' local eval2 = `r(sum)' gen `eval3gen' = 1/((`2')^2) su `eval3gen' if `touse' local eval3 = `r(sum)' local Ni2 = wordcount("`isquared'") /* Line patterns */ local lpat "shortdash longdash_dot dash_dot shortdash_dot dot" local n 0 foreach i in `lpat' { local lp`++n' `i' } *** Need to generate each contour as a variable, then replace as missing ones outside plot region. local pati 1 tokenize `isquared' forvalues i2 =1/`Ni2' { if (``i2'' > 100) { di as error "I-squared must be between 0% and 100%" exit 198 } tempvar isq`i2'1 isq`i2'2 gen double `isq`i2'1' = (`eval2'/`eval3')+sqrt((((`eval3'+(1/((`tempseSS')^2)))/(`eval3'*(1/(`tempseSS'^2))))*(((`eval2')^2/(`eval3'+(1/((`tempseSS')^2))))+(`N_SS'/(1-(``i2''/100)))-`eval1'))+(`eval2'/`eval3')^2) count if inrange(`isq`i2'1',`xmin',`xmax') local ncont1 = r(N) gen double `isq`i2'2' = (`eval2'/`eval3')-sqrt((((`eval3'+(1/((`tempseSS')^2)))/(`eval3'*(1/(`tempseSS'^2))))*(((`eval2')^2/(`eval3'+(1/((`tempseSS')^2))))+(`N_SS'/(1-(``i2''/100)))-`eval1'))+(`eval2'/`eval3')^2) count if inrange(`isq`i2'2',`xmin',`xmax') local ncont2 = r(N) local percentshown = 100*(`ncont1'+`ncont2')/(`cpoints'*2) *noisily di as txt %8.4f `percentshown' "% of I-squared = ``i2'' is shown" local line "`line' (line `tempseSS' `isq`i2'1' if inrange(`isq`i2'1',`xmin',`xmax'),lpat(`lp`pati'') lcol(black))(line `tempseSS' `isq`i2'2' if inrange(`isq`i2'2',`xmin',`xmax'), lpat(`lp`pati'') lcol(black))" local dumplots "`dumplots' (line `tempseSS' `isq`i2'1' if inrange(`isq`i2'1',`xmin',`xmax'),lpat(`lp`pati'') lcol(black))" local `++labnum' local lab `"`lab' `labnum' "{it:I}-squared = ``i2''%""' local pati = `pati'+1 /*if "`eform'"!="" { replace isq`i2'1 = exp(isq`i2'1) replace isq`i2'2 = exp(isq`i2'2) }*/ } } /*****************************************************************************************************************************************************/ /* TAU-SQUARED CONTOURS */ if ("`isquared'"=="" & "`tausquared'"!="") { tokenize `varlist' /*replace `1' = `pooled_est' in `=`N_SS'+1' replace `2' = 0.1 in `=`N_SS'+1' metan `1' `2', `method'i nograph nokeep level(`level') noisily di as txt " Minimum tau-squared =" %8.4f `r(tau2)'*/ *prelim calcs. tempvar ss2 sess2 evalv evalx evaly evalz gen `ss2' = `1'^2 gen `sess2'=`2'^2 gen `evalv' = (1/`sess2')^2 su `evalv' if `touse', mean local evalvl = `r(sum)' gen `evalx' = `1'*(1/`sess2') su `evalx' if `touse', mean local evalxl = `r(sum)' gen `evaly' = 1/`sess2' su `evaly' if `touse', mean local evalyl = `r(sum)' gen `evalz' = (1/`sess2')*`ss2' su `evalz' if `touse', mean local evalzl = `r(sum)' local sstest1 = `evalxl'/`evalyl' local ntau2 = wordcount("`tausquared'") local lpat "shortdash longdash_dot dash_dot shortdash_dot dot" local tt 0 foreach i in `lpat' { local lp`++tt' `i' } local pati 1 tokenize `tausquared' forvalues tau22 =1/`ntau2' { tempvar tau2`tau22'1 tau2`tau22'2 gen double `tau2`tau22'1' = `sstest1' + sqrt(((`evalyl'+(1/(`tempseSS'^2)))/(`evalyl'*(1/(`tempseSS'^2))))* (`evalyl'+(1/(`tempseSS'^2))-((`evalvl'+(1/(`tempseSS'^2))^2)/(`evalyl'+(1/(`tempseSS'^2)))))*``tau22'' +((`sstest1')^2) - ((`evalyl'+(1/(`tempseSS'^2)))/(`evalyl'*(1/(`tempseSS'^2)))) * (`evalzl'-(((`evalxl')^2)/(`evalyl'+(1/(`tempseSS'^2))))-`N_SS')) count if inrange(`tau2`tau22'1',`xmin',`xmax') local ncont1 = r(N) gen double `tau2`tau22'2' = `sstest1' - sqrt(((`evalyl'+(1/(`tempseSS'^2)))/(`evalyl'*(1/(`tempseSS'^2))))* (`evalyl'+(1/(`tempseSS'^2))-((`evalvl'+(1/(`tempseSS'^2))^2)/(`evalyl'+(1/(`tempseSS'^2)))))*``tau22'' +((`sstest1')^2) - ((`evalyl'+(1/(`tempseSS'^2)))/(`evalyl'*(1/(`tempseSS'^2)))) * (`evalzl'-(((`evalxl')^2)/(`evalyl'+(1/(`tempseSS'^2))))-`N_SS')) count if inrange(`tau2`tau22'2',`xmin',`xmax') local ncont2 = r(N) local percentshown = 100*(`ncont1'+`ncont2')/(`cpoints'*2) *noisily di as txt %8.4f `percentshown' "% of Tau-squared = ``tau22'' is shown" local line "`line' (line `tempseSS' `tau2`tau22'1' if inrange(`tau2`tau22'1',`xmin',`xmax'),lpat(`lp`pati'') lcol(black))(line `tempseSS' `tau2`tau22'2' if inrange(`tau2`tau22'2',`xmin',`xmax'), lpat(`lp`pati'') lcol(black))" local dumplots "`dumplots' (line `tempseSS' `tau2`tau22'1' if inrange(`tau2`tau22'1',`xmin',`xmax'),lpat(`lp`pati'') lcol(black))" local `++labnum' local lab `"`lab' `labnum' "{it:{&tau}} - squared = ``tau22''""' local pati = `pati'+1 /*if "`eform'"!="" { replace tau2`tau22'1 = exp(tau2`tau22'1) replace tau2`tau22'2 = exp(tau2`tau22'2) }*/ } } /*****************************************************************************************************************************************************/ /* SUMMARY PLOTS */ tokenize `varlist' /* SEAGULL PLOT */ if "`newstudytreatment'"!="" & "`newstudycontrol'"!="" { local newobs = `=`newstudytreatment'+1'*`=`newstudycontrol'+1' cap set obs `newobs' tempvar indnew logornew selogornew gen `indnew' = 0 gen `logornew' = . gen `selogornew' = . local templogor "`1'" local tempselogor "`2'" local ii = 1 forvalues i=1/`=`newstudytreatment'+1' { forvalues j=1/`=`newstudycontrol'+1' { local event_t = `i' - 1 local noevent_t = `newstudytreatment' - `event_t' local event_c = `j' - 1 local noevent_c = `newstudycontrol' - `event_c' if `event_t'==0 | `noevent_t'==0 | `event_c'==0 | `noevent_c'==0 { local event_t = `event_t' + 0.5 local noevent_t = `noevent_t' + 0.5 local event_c = `event_c' + 0.5 local noevent_c = `noevent_c' + 0.5 } if "`stat'"!="rr" { local lor = log((`event_t'*`noevent_c')/(`noevent_t'*`event_c')) local selor = sqrt((1/`event_t') + (1/`noevent_t') + (1/`event_c') + (1/`noevent_c')) } else { local p1 = (`event_t'/(`noevent_t'+`event_t')) local p2 = (`event_c'/(`noevent_c'+`event_c')) local lor = log(`p1'/`p2') local selor = sqrt(((1-`p1')/(`p1'*(`noevent_t'+`event_t')))+((1-`p2')/(`p2'*(`noevent_c'+`event_c')))) } replace `logornew' = `lor' in `ii' replace `templogor' = `lor' in `N_new' replace `selogornew' = `selor' in `ii' replace `tempselogor' = `selor' in `N_new' metan `templogor' `tempselogor', `method' nograph notable level(`level') local logor = `r(ES)' local selogor = `r(seES)' local uci = `logor' + `siglev' * `selogor' local lci = `logor' - `siglev' * `selogor' if `lci'<0 & `uci'<0 { replace `indnew' = 1 in `ii' } else if `uci'>0 & `lci'>0 { replace `indnew' = 2 in `ii' } local `++ii' } } if "`eform'"!="" { replace `logornew' = exp(`logornew') } local seagull "(scatter `selogornew' `logornew' if `indnew'==0,ms(o) mcol(`col1'))(scatter `selogornew' `logornew' if `indnew'==1,ms(o) mcol(`col3'))(scatter `selogornew' `logornew' if `indnew'==2,ms(o) mcol(`col2') )" local dumplots "`dumplots' (scatter `selogornew' `logornew' if `indnew'==0,ms(o) mcol(`col1'))(scatter `selogornew' `logornew' if `indnew'==1,ms(o) mcol(`col3'))(scatter `selogornew' `logornew' if `indnew'==2,ms(o) mcol(`col2') )" local labnum `++labnum' local lab `"`lab' `labnum' "Non-sig. effect (`=100-`level''% level)""' local labnum `++labnum' local lab `"`lab' `labnum' "Sig. effect < NULL (`=100-`level''% level)""' local labnum `++labnum' local lab `"`lab' `labnum' "Sig. effect > NULL (`=100-`level''% level)""' qui su `selogornew',mean if `r(max)'>`ymax' { scalar `ymax' = `r(max)' } /* Prediction interval */ if "`prediction'"!="" { local tstat = abs(invttail(`=`N_SS'-2',(100-`level')/200)) tempvar p1 p2 sumdy gen `sumdy' = `sumdposition' gen `p1' = `pooled_est' + `tstat'*sqrt(`se_pooled_est'^2+`tau2') gen `p2' = `pooled_est' - `tstat'*sqrt(`se_pooled_est'^2+`tau2') local predplot "(rcap `p2' `p1' `sumdy' in 1, hor lcol(black))" local dumplots "`dumplots'(rcap `p2' `p1' `sumdy' in 1, hor lcol(black))" if "`eform'"!="" { replace `p1' = exp(`p1') replace `p2' = exp(`p2') } local labnum `++labnum' local lab `"`lab' `labnum' "Prediction interval""' } /* Show summary diamond on plot */ if "`sumd'"!="" { tempvar sumdiam_x sumdiam_y gen `sumdiam_x' = `pooled_est' - `siglev'*`se_pooled_est' in 1 replace `sumdiam_x' = `pooled_est' in 2 replace `sumdiam_x' = `pooled_est'+`siglev'*`se_pooled_est' in 3 replace `sumdiam_x' = `pooled_est' in 4 gen `sumdiam_y' = `sumdposition' in 1 replace `sumdiam_y' = `sumdposition' - 0.04*`ywidth' in 2 replace `sumdiam_y' = `sumdposition' in 3 replace `sumdiam_y' = `sumdposition' + 0.04*`ywidth' in 4 if "`eform'"!="" { replace `sumdiam_x' = exp(`sumdiam_x') } local sumarea "(area `sumdiam_y' `sumdiam_x', nodropbase col(black) fi(inten100))" local dumplots "`dumplots'(area `sumdiam_y' `sumdiam_x', nodropbase col(black) fi(inten100))" local labnum `++labnum' } /* Display y-line at null hypothesis */ if "`nonullline'"=="" & "`newstudytreatment'"!="" { tempvar y_null_line x_null_line range `y_null_line' `ymin' `ymax' 2 gen `x_null_line' = `null' in 1/2 if "`eform'"!="" { replace `x_null_line' = exp(`x_null_line') } local nullline "(line `y_null_line' `x_null_line', lpat(solid) lcol(black))" local dumplots "`dumplots'(line `y_null_line' `x_null_line', lpat(solid) lcol(black))" local labnum `++labnum' local lab `"`lab' `labnum' "Null effect""' } /* Display y-line at pooled estimate of original meta-analysis */ if "`nopooledline'"=="" & "`newstudytreatment'"!="" { tempvar y_pooled_est x_pooled_est range `y_pooled_est' `ymin' `ymax' 2 gen `x_pooled_est' = `pooled_est' in 1/2 if "`eform'"!="" { replace `x_pooled_est' = exp(`x_pooled_est') } local pooledline "(line `y_pooled_est' `x_pooled_est', lpat(dash) lcol(black))" local dumplots "`dumplots'(line `y_pooled_est' `x_pooled_est', lpat(dash) lcol(black))" local labnum `++labnum' local lab `"`lab' `labnum' "Pooled effect""' } } /* Prediction interval */ if "`prediction'"!="" & "`newstudytreatment'"=="" { local tstat = abs(invttail(`=`N_SS'-2',(100-`level')/200)) tempvar p1 p2 sumdy gen `sumdy' = `sumdposition' gen `p1' = `pooled_est' + `tstat'*sqrt(`se_pooled_est'^2+`tau2') gen `p2' = `pooled_est' - `tstat'*sqrt(`se_pooled_est'^2+`tau2') local predplot "(rcap `p2' `p1' `sumdy' in 1, hor lcol(black))" local dumplots "`dumplots'(rcap `p2' `p1' `sumdy' in 1, hor lcol(black))" if "`eform'"!="" { replace `p1' = exp(`p1') replace `p2' = exp(`p2') } local labnum `++labnum' local lab `"`lab' `labnum' "Prediction interval""' } /* Show summary diamond on plot */ if "`sumd'"!="" & "`newstudytreatment'"=="" { tempvar sumdiam_x sumdiam_y gen `sumdiam_x' = `pooled_est' - `siglev'*`se_pooled_est' in 1 replace `sumdiam_x' = `pooled_est' in 2 replace `sumdiam_x' = `pooled_est'+`siglev'*`se_pooled_est' in 3 replace `sumdiam_x' = `pooled_est' in 4 gen `sumdiam_y' = `sumdposition' in 1 replace `sumdiam_y' = `sumdposition' - 0.04*`ywidth' in 2 replace `sumdiam_y' = `sumdposition' in 3 replace `sumdiam_y' = `sumdposition' + 0.04*`ywidth' in 4 if "`eform'"!="" { replace `sumdiam_x' = exp(`sumdiam_x') } local sumarea "(area `sumdiam_y' `sumdiam_x', nodropbase col(black) fi(inten100))" local dumplots "`dumplots'(area `sumdiam_y' `sumdiam_x', nodropbase col(black) fi(inten100))" local labnum `++labnum' } /* Display y-line at null hypothesis */ if "`nonullline'"=="" & "`newstudytreatment'"=="" & "`measure'"!="ciwidth"{ tempvar y_null_line x_null_line range `y_null_line' `ymin' `ymax' 2 gen `x_null_line' = `null' in 1/2 if "`eform'"!="" { replace `x_null_line' = exp(`x_null_line') } local nullline "(line `y_null_line' `x_null_line', lpat(solid) lcol(black))" local dumplots "`dumplots'(line `y_null_line' `x_null_line', lpat(solid) lcol(black))" local labnum `++labnum' local lab `"`lab' `labnum' "Null effect""' } /* Display y-line at pooled estimate of original meta-analysis */ if "`nopooledline'"=="" & "`newstudytreatment'"=="" { tempvar y_pooled_est x_pooled_est range `y_pooled_est' `ymin' `ymax' 2 gen `x_pooled_est' = `pooled_est' in 1/2 if "`eform'"!="" { replace `x_pooled_est' = exp(`x_pooled_est') } local pooledline "(line `y_pooled_est' `x_pooled_est', lpat(dash) lcol(black))" local dumplots "`dumplots'(line `y_pooled_est' `x_pooled_est', lpat(dash) lcol(black))" local labnum `++labnum' local lab `"`lab' `labnum' "Pooled effect""' } /* Show scatterplot of original study results on plot */ if "`noscatter'"=="" { tokenize `varlist' if "`eform'"!="" { tempvar scat1 gen `scat1' = exp(`1') if `touse' local scatterplot "(scatter `2' `scat1' if `touse', ms(o) mcol(black))" local dumplots "`dumplots'(scatter `2' `scat1' if `touse', ms(o) mcol(black))" } else { local scatterplot "(scatter `2' `1' if `touse', ms(o) mcol(black))" local dumplots "`dumplots'(scatter `2' `1' if `touse', ms(o) mcol(black))" } local labnum `++labnum' local lab `"`lab' `labnum' "Study effects""' } if "`loeline'"!="" { tempvar y_loe loe1 loe2 range `y_loe' `ymin' `ymax' 2 local loe11 : word 1 of `loe' gen `loe1' = `loe11' in 1/2 local loe21 : word 2 of `loe' gen `loe2' = `loe21' in 1/2 if "`eform'"!="" { replace `loe1' = exp(`loe1') replace `loe2' = exp(`loe2') } local loeline "(line `y_loe' `loe1', lpat(dot) lcol(black))" local dumplots "`dumplots'(line `y_loe' `loe1', lpat(dot) lcol(black))" local labnum `++labnum' local lab `"`lab' `labnum' "Limits of equivalence""' local loeline "`loeline' (line `y_loe' `loe2', lpat(dot) lcol(black))" local labnum `++labnum' } /*****************************************************************************************************************************************************/ /* CORE CONTOUR CALCULATIONS AND AREA GRAPH BUILDING */ tokenize `varlist' /* re-tokenize varlist */ noisily di noisily di as txt "Building graph:" if "`noshading'"=="" & "`newstudytreatment'"=="" { if ("`method'"=="fixedi" & "`measure'"=="" & "`loe'"==""){ /* Fixed effects can probably be done better */ *calculate weights tempvar weights weight_SS vwt c1SS c2SS gen double `weights' = 1/(`2'^2) if `touse' su `weights' if `touse', meanonly local sum_weights = `r(sum)' *Pre-calcs for formula gen double `weight_SS' = `weights'*`1' if `touse' su `weight_SS' if `touse', meanonly local sum_weights_SS = `r(sum)' gen double `vwt' = 1/(`tempseSS'^2) *Apply formulas gen double `c1SS' = (1/`vwt') * (`null' * (`sum_weights' + `vwt') - `sum_weights_SS' + `siglev'*(`sum_weights'+`vwt')^0.5) gen double `c2SS' = (1/`vwt') * (`null' * (`sum_weights' + `vwt') - `sum_weights_SS' - `siglev'*(`sum_weights'+`vwt')^0.5) ******************************************************************************************************************** *generate dummy repeat vectors for shaded areas (manipulation purposes) tempvar c1SS_2 c2SS_2 tempse1SS_2 tempse2SS_2 gen double `c1SS_2' = `c1SS' gen double `c2SS_2' = `c2SS' gen double `tempse1SS_2' = `tempseSS' gen double `tempse2SS_2' = `tempseSS' ******************************************************************************************************************** ***prepare c1SS_2 and c2SS_2 *append on end of c2SS_2 and tempse2SS_2 up to 4 coords. *(xmin,ymax) and (xmin,ymin) local N_1 = _N local ext1 = _N+1 local ext2 = _N+2 local ext3 = _N+3 cap set obs `ext3' *Locals to find direction of contours local start_c1SS_2 = `c1SS_2'[1] local end_c1SS_2 = `c1SS_2'[`N_1'] local start_c2SS_2 = `c2SS_2'[1] local end_c2SS_2 = `c2SS_2'[`N_1'] *Right hand line swings right if (`end_c1SS_2' > `start_c1SS_2') { replace `c1SS_2' = . if `c1SS_2' < `xmin' replace `c1SS_2' = . if `c1SS_2' > `xmax' replace `tempse1SS_2' = . if `c1SS_2' == . if (`tempse1SS_2'[`N_1']==`ymax') { replace `c1SS_2' = `xmax' in `ext1' replace `tempse1SS_2' = `ymax' in `ext1' } local area1 "(area `tempse1SS_2' `c1SS_2', col(`col2') fi(inten100))" } *Right hand swings left if (`end_c1SS_2' < `start_c1SS_2'){ replace `c1SS_2' =. if `c1SS_2' < `xmin' replace `c1SS_2' =. if `c1SS_2' > `xmax' replace `tempse1SS_2' =. if `c1SS_2' == . if (`tempse1SS_2'[`N_1']==`ymax') { replace `c1SS_2' = `xmax' in `ext1' replace `tempse1SS_2' = `ymax' in `ext1' } else { replace `c1SS_2' = `xmin' in `ext1' replace `tempse1SS_2' = `ymax' in `ext1' replace `c1SS_2' = `xmax' in `ext2' replace `tempse1SS_2' = `ymax' in `ext2' } local area1 "(area `tempse1SS_2' `c1SS_2', col(`col2') fi(inten100))" } *Left hand line swings left if (`end_c2SS_2' < `start_c2SS_2'){ replace `c2SS_2' = . if `c2SS_2' < `xmin' replace `c2SS_2' = . if `c2SS_2' > `xmax' replace `tempse2SS_2' = . if `c2SS_2' == . if (`tempse2SS_2'[`N_1']==`ymax'){ replace `c2SS_2' = `xmin' in `ext1' replace `tempse2SS_2' = `ymax' in `ext1' } local area2 "(area `tempse2SS_2' `c2SS_2', col(`col3') fi(inten100))" } *Left hand swings right if (`end_c2SS_2' > `start_c2SS_2') { replace `c2SS_2' =. if `c2SS_2' < `xmin' replace `c2SS_2' =. if `c2SS_2' > `xmax' replace `tempse2SS_2' =. if `c2SS_2' ==. if (`tempse2SS_2'[`N_1']==`ymax') { replace `c2SS_2' = `xmin' in `ext1' replace `tempse2SS_2' = `ymax' in `ext1' } else { replace `c2SS_2' = `xmax' in `ext1' replace `tempse2SS_2' = `ymax' in `ext1' replace `c2SS_2' = `xmin' in `ext2' replace `tempse2SS_2' = `ymax' in `ext2' } local area2 "(area `tempse2SS_2' `c2SS_2', col(`col3') fi(inten100))" } if "`eform'"!="" { replace `c1SS_2' = exp(`c1SS_2') replace `c2SS_2' = exp(`c2SS_2') } } else { *calculate weights tempvar weights gen `weights' = 1/((`2'^2)+`tau2') if `touse' local obs = (`cpoints'-1)*4 cap set obs `obs' *re code into 4 coord format: tempvar newx newy mat_ind gen double `newx' = . gen double `newy' = . gen double `mat_ind' = 0 *gen id variable, each id of length 4 tempvar id tempind local obs1 = `cpoints'-1 local obs2 = `obs1'-1 gen `tempind' = _n<=`obs' egen double `id'=seq() if `tempind', from(1) to(`obs1') block(4) replace `newx' = `tempSS'[1] in 1/2 replace `newx' = `tempSS'[2] in 3/4 replace `newy' = `tempseSS'[1] in 1/4 replace `newy' = `tempseSS'[2] in 2/3 local add=4 forvalues i = 1/`obs2' { local t1 = `i'+`add' local t2 = `i'+`add'+1 local t3 = `i'+`add'+2 local t4 = `i'+`add'+3 local tt1=`i'+1 local tt2=`i'+2 replace `newx' = `tempSS'[`tt1'] in `t1'/`t4' replace `newx' = `tempSS'[`tt2'] in `t3'/`t4' replace `newy' = `tempseSS'[1] in `t1'/`t4' replace `newy' = `tempseSS'[2] in `t2'/`t3' local add =`add'+3 } local varnames forvalues j = 1/`obs1'{ tempvar newy`j' matind`j' gen double `newy`j''=`tempseSS'[`j'] gen `matind`j'' =. local j2 = `j'+1 replace `newy`j'' = `tempseSS'[`j2'] if `newy' == `tempseSS'[2] local varnames "`varnames' `matind`j''" /* Variables to paste 0/1/... coding into */ } if "`method'"=="randomi" | "`measure'"!="" | "`loe'"!="" { tokenize `varlist' replace `touse' = 0 if `1'==. replace `touse' = 0 if `2'==. mata: dl1mata("`1'","`2'","`method'","`touse'") /* Execute mata program to conduct M-A's */ if "`eform'"!="" { replace `newx' = exp(`newx') } } /* Returns matrix (as variables) of 0,1,2's to identify region */ /*********** REFERENCE CODE **************/ /*else { forvalues p=1/`obs1'{ *first 2 and last 2 local obs_minus1 = `obs'-1 foreach s in 1 2 `obs_minus1' `obs' { qui replace `1' = `newx'[`s'] in `N_new' qui replace `2' = `newy`p''[`s'] in `N_new' metan `1' `2', `method' nograph notable local lci = `r(ES)'-`siglev'*`r(seES)' local uci = `r(ES)'+`siglev'*`r(seES)' if (`uci' < `null' & `lci' < `null') { replace `matind`p'' = 1 in `s' } * if stat. sig. > null else if (`lci' > `null' & `uci' > `null') { replace `matind`p'' = 2 in `s' } } *remaining local tfinal = `obs'-5 forvalues t = 3(4)`tfinal' { local tt2 = `t'+1 local tt3 = `t'+2 local tt4 = `t'+3 qui replace `1' = `newx'[`t'] in `N_new' qui replace `2' = `newy`p''[`t'] in `N_new' metan `1' `2', `method' nograph notable local lci = `r(ES)'-`siglev'*`r(seES)' local uci = `r(ES)'+`siglev'*`r(seES)' if (`uci' < `null' & `lci' < `null') { replace `matind`p'' = 1 in `t' replace `matind`p'' = 1 in `tt4' } * if stat. sig. > null else if (`lci' > `null' & `uci' > `null') { replace `matind`p'' = 2 in `t' replace `matind`p'' = 2 in `tt4' } replace `1' = `newx'[`tt2'] in `N_new' replace `2' = `newy`p''[`tt2'] in `N_new' metan `1' `2', `method' nograph notable local lci = `r(ES)'-`siglev'*`r(seES)' local uci = `r(ES)'+`siglev'*`r(seES)' if (`uci' < `null' & `lci' < `null') { replace `matind`p'' = 1 in `tt2'/`tt3' } * if stat. sig. > null else if (`lci' > `null' & `uci' > `null') { replace `matind`p'' = 2 in `tt2'/`tt3' } } } }*/ /* Return matrices of effect estimates and standard errors */ tempfile temp1 save `temp1', replace tempvar tempind2 tempnew bys `id': gen `tempnew' = `matind`obs1''[2] local varnames "`varnames' `tempnew'" gen `tempind2' = _n==`obs' bys `id': keep if _n==1 | `tempind2'==1 tempname indmat mkmat `varnames', matrix(`indmat') nomissing mat `indmat' = `indmat'' return matrix status = `indmat' use `temp1',clear *make four cells equal in matrix: forvalues i=1/`obs1'{ bys `id': replace `matind`i'' = `matind`i''[1] if (`matind`i''[1]!=`matind`i''[2] | `matind`i''[2]!=`matind`i''[3] | `matind`i''[3]!=`matind`i''[4] | `matind`i''[1]!=`matind`i''[2]) } *Loop over newy's local area "" *replace matind1=1 if _n>4 forvalues w=1/`obs1'{ *w^th row: 1st colour section: local rowN = (`obs1'*4) local rowN2 = (`obs1'*4)-2 local n1 = 1 local n2 = 2 local n3 = 3 local ind1 = `matind`w''[1] local ind2 = `matind`w''[2] local i=3 while `i'<`rowN2'{ while `ind1'==`ind2' { local ind2 = `matind`w''[`i'] local i = `i'+1 if (`i'==`rowN2') { continue, break } } continue, break } *First change in colour detected if (`i'!=`rowN2') { ***need to newy1 now, then do graph *changes code at i-1 local id1 = `i'-1 local id2 = `i' *index for replace with missing in first colour : 1+2 to new_start-3 local miss_replace_end = `i'-4 local end = `i'-2 if (`miss_replace_end'>3) { /*If colour needs to change after first square*/ replace `newy`w'' = . in `n3'/`miss_replace_end' } *need to do graph colours as well if (power!=4) { if (`matind`w''[1]==1) { local area "`area' (area `newy`w'' `newx' in 1/`end', nodropb col(`col3') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==2) { local area "`area' (area `newy`w'' `newx' in 1/`end', nodropb col(`col2') fi(inten100))" /*first area*/ } } else { if (`matind`w''[1]==2) { local area "`area' (area `newy`w'' `newx' in 1/`end', nodropb col(`col2') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==3) { local area "`area' (area `newy`w'' `newx' in 1/`end', nodropb col(`col3') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==4) { local area "`area' (area `newy`w'' `newx' in 1/`end', nodropb col(`col4') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==5) { local area "`area' (area `newy`w'' `newx' in 1/`end', nodropb col(`col5') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==6) { local area "`area' (area `newy`w'' `newx' in 1/`end', nodropb col(`col6') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==7) { local area "`area' (area `newy`w'' `newx' in 1/`end', nodropb col(`col7') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==8) { local area "`area' (area `newy`w'' `newx' in 1/`end', nodropb col(`col8') fi(inten100))" /*first area*/ } } ********************************************** *w^th row: 2nd colour local ind1 = `matind`w''[`id1'] local ind2 = `matind`w''[`id2'] local n1 = `id1' local n2 = `id2' local n3 = `n2'+1 local i = `n3' while `i'<`rowN2' { while `ind1'==`ind2' { local ind2=`matind`w''[`i'] local i=`i'+1 if (`i'==`rowN2') { continue, break } } continue, break } while (`i'<`rowN2') { if (`i'!=`rowN2') { local id1 = `i'-1 local id2 = `i' *2nd colour change detected *index rows to change to missing local missing_replace_end = `i'-4 local end_new_section = `i'-2 *edit newy`' if needed: if (`missing_replace_end'>`n3') { replace `newy`w'' = . in `n3'/`missing_replace_end' } *exit *colours if (power!=4) { if (`matind`w''[`n1']==1) { local area "`area' (area `newy`w'' `newx' in `n1'/`end_new_section', nodropb col(`col3') fi(inten100))" /*second area*/ } else if (`matind`w''[`n1']==2) { local area "`area' (area `newy`w'' `newx' in `n1'/`end_new_section', nodropb col(`col2') fi(inten100))" /*second area*/ } } else { if (`matind`w''[`n1']==2) { local area "`area' (area `newy`w'' `newx' in `n1'/`end_new_section', nodropb col(`col2') fi(inten100))" /*first area*/ } else if (`matind`w''[`n1']==3) { local area "`area' (area `newy`w'' `newx' in `n1'/`end_new_section', nodropb col(`col3') fi(inten100))" /*first area*/ } else if (`matind`w''[`n1']==4) { local area "`area' (area `newy`w'' `newx' in `n1'/`end_new_section', nodropb col(`col4') fi(inten100))" /*first area*/ } else if (`matind`w''[`n1']==5) { local area "`area' (area `newy`w'' `newx' in `n1'/`end_new_section', nodropb col(`col5') fi(inten100))" /*first area*/ } else if (`matind`w''[`n1']==6) { local area "`area' (area `newy`w'' `newx' in `n1'/`end_new_section', nodropb col(`col6') fi(inten100))" /*first area*/ } else if (`matind`w''[`n1']==7) { local area "`area' (area `newy`w'' `newx' in `n1'/`end_new_section', nodropb col(`col7') fi(inten100))" /*first area*/ } else if (`matind`w''[`n1']==8) { local area "`area' (area `newy`w'' `newx' in `n1'/`end_new_section', nodropb col(`col8') fi(inten100))" /*first area*/ } } ********************************************** *w^th row: 3rd colour local ind1 = `matind`w''[`id1'] local ind2 = `matind`w''[`id2'] local n1 = `id1' local n2 = `n1'+1 local n3 = `n2'+1 local i = `n3' while `i'<`rowN2'{ while `ind1'==`ind2'{ local ind2 = `matind`w''[`i'] local i =`i'+1 if (`i'==`rowN2'){ continue, break } } continue, break } } } *colour is same for remainder of row if (`i'==`rowN2') { replace `newy`w'' = . in `n3'/`i' if (power!=4) { if (`matind`w''[`n1']==1) { local area "`area' (area `newy`w'' `newx' in `n1'/`obs', nodropb col(`col3') fi(inten100))" } else if (`matind`w''[`n1']==2) { local area "`area' (area `newy`w'' `newx' in `n1'/`obs', nodropb col(`col2') fi(inten100))" } } else { if (`matind`w''[`n1']==2) { local area "`area' (area `newy`w'' `newx' in `n1'/`obs', nodropb col(`col2') fi(inten100))" } else if (`matind`w''[`n1']==3) { local area "`area' (area `newy`w'' `newx' in `n1'/`obs', nodropb col(`col3') fi(inten100))" } else if (`matind`w''[`n1']==4) { local area "`area' (area `newy`w'' `newx' in `n1'/`obs', nodropb col(`col4') fi(inten100))" } else if (`matind`w''[`n1']==5) { local area "`area' (area `newy`w'' `newx' in `n1'/`obs', nodropb col(`col5') fi(inten100))" } else if (`matind`w''[`n1']==6) { local area "`area' (area `newy`w'' `newx' in `n1'/`obs', nodropb col(`col6') fi(inten100))" } else if (`matind`w''[`n1']==7) { local area "`area' (area `newy`w'' `newx' in `n1'/`obs', nodropb col(`col7') fi(inten100))" } else if (`matind`w''[`n1']==8) { local area "`area' (area `newy`w'' `newx' in `n1'/`obs', nodropb col(`col8') fi(inten100))" } } } } *First colour is same for entire row: else { local miss_replace_end = `obs'-2 replace `newy`w'' = . in 3/`i' if (power!=4) { if (`matind`w''[1]==1) { local area "`area' (area `newy`w'' `newx' in 1/`obs', nodropb col(`col3') fi(inten100))" } else if (`matind`w''[1]==2) { local area "`area' (area `newy`w'' `newx' in 1/`obs', nodropb col(`col2') fi(inten100))" } } else { if (`matind`w''[1]==2) { local area "`area' (area `newy`w'' `newx' in 1/`obs', nodropb col(`col2') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==3) { local area "`area' (area `newy`w'' `newx' in 1/`obs', nodropb col(`col3') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==4) { local area "`area' (area `newy`w'' `newx' in 1/`obs', nodropb col(`col4') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==5) { local area "`area' (area `newy`w'' `newx' in 1/`obs', nodropb col(`col5') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==6) { local area "`area' (area `newy`w'' `newx' in 1/`obs', nodropb col(`col6') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==7) { local area "`area' (area `newy`w'' `newx' in 1/`obs', nodropb col(`col7') fi(inten100))" /*first area*/ } else if (`matind`w''[1]==8) { local area "`area' (area `newy`w'' `newx' in 1/`obs', nodropb col(`col8') fi(inten100))" /*first area*/ } } } /* BUG FIX - if colour change occurs in last square it is ignored above */ local fix1 = `obs' - 4 local fix2 = `obs' - 3 if (`matind`w''[`fix1']!=`matind`w''[`fix2']) { tempvar n2ewy`w' gen `n2ewy`w'' = `newy`w'' replace `n2ewy`w'' = `n2ewy`w''[`obs'] in `fix2' replace `n2ewy`w'' = `n2ewy`w''[`=`obs'-1'] in `=`fix2'+1' if (power!=4) { if (`matind`w''[`obs']==1) { local area "`area' (area `n2ewy`w'' `newx' in `fix2'/`obs', nodropb col(`col3') fi(inten100))" } else if (`matind`w''[`obs']==2) { local area "`area' (area `n2ewy`w'' `newx' in `fix2'/`obs', nodropb col(`col2') fi(inten100))" } else { local area "`area' (area `n2ewy`w'' `newx' in `fix2'/`obs', nodropb col(`col1') fi(inten100))" } } else { if (`matind`w''[`obs']==2) { local area "`area' (area `newy`w'' `newx' in `fix2'/`obs', nodropb col(`col2') fi(inten100))" } else if (`matind`w''[`obs']==3) { local area "`area' (area `newy`w'' `newx' in `fix2'/`obs', nodropb col(`col3') fi(inten100))" } else if (`matind`w''[`obs']==4) { local area "`area' (area `newy`w'' `newx' in `fix2'/`obs', nodropb col(`col4') fi(inten100))" } else if (`matind`w''[`obs']==5) { local area "`area' (area `newy`w'' `newx' in `fix2'/`obs', nodropb col(`col5') fi(inten100))" } else if (`matind`w''[`obs']==6) { local area "`area' (area `newy`w'' `newx' in `fix2'/`obs', nodropb col(`col6') fi(inten100))" } else if (`matind`w''[`obs']==7) { local area "`area' (area `newy`w'' `newx' in `fix2'/`obs', nodropb col(`col7') fi(inten100))" } else if (`matind`w''[`obs']==8) { local area "`area' (area `newy`w'' `newx' in `fix2'/`obs', nodropb col(`col8') fi(inten100))" } else { local area "`area' (area `n2ewy`w'' `newx' in `fix2'/`obs', nodropb col(`col1') fi(inten100))" } } } } } } /*****************************************************************************************************************************************************/ /* FINAL TWOWAY GRAPH */ if ("`method'" == "fixedi" & "`measure'"=="" & "`loe'"==""){ twoway `dumplots' (`addplot') `basearea' `area1' `area2' `line' `seagull' `predplot' `sumarea' /// `nullline' `pooledline' `scatterplot' (`addplot') /// , plotregion(margin(zero)) yscale(reverse) `explabel' xtitle(`xtitle') ytitle(`ytitle') ylabel(,angle(horizontal)) /// legend(order(`lab') symxsize(*0.6)) `xline' title("Contours for impact of a new study") /// `options' } else { twoway `dumplots' (`addplot') `basearea' `area' `line' `seagull' `predplot' `sumarea' /// `nullline' `pooledline' `scatterplot' `loeline' (`addplot') /// , yscale(reverse) `explabel' plotregion(margin(zero)) xtitle(`xtitle') ytitle(`ytitle') ylabel(,angle(horizontal)) /// legend(order(`lab') symxsize(*0.6) `extraleg') `xline' title("Contours for impact of a new study") /// `options' } restore } end /* Stata command `range' + `double' in the generate statement */ program define range2 version 3.1 if "`3'"=="" | "`5'"!="" { error 198 } confirm new var `1' if _N==0 { if "`4'"=="" { error 198 } set obs `4' local o "`4'" } else { if "`4'"!="" { local o "`4'" if `4' > _N { set obs `4' } } else local o=_N } gen double `1'=(_n-1)/(`o'-1)*((`3')-(`2'))+(`2') in 1/`o' end /* Mata program to conduct random effects DerSimonian and Laird meta-analysis */ mata: void dl1mata(string scalar ssvar, string scalar sevar, string scalar method, string scalar touse) { nint = st_numscalar("contpoints") nintminus1 = nint:-1 nint2 = nint:^2 nstudies = st_numscalar("nstudies") siglevel = st_numscalar("siglev") null = st_numscalar("null") studies = st_data(.,ssvar,touse)' sterrors = st_data(.,sevar,touse)' origss_mat = J(nint2[1,1],1,studies) origse_mat = J(nint2[1,1],1,sterrors) setemp = rangen(st_numscalar("ymin"),st_numscalar("ymax"),nint)' //row vectors sstemp = rangen(st_numscalar("xmin"),st_numscalar("xmax"),nint)' serep = J(nint,1,setemp') //each y range repeated and put in column vector ssrep = J(1,cols(sstemp),sstemp[1,1]) for(i=2;i<=cols(setemp);i++){ ssrep = ssrep,J(1,cols(sstemp),sstemp[1,i]) } ssrep=ssrep' newss_mat = origss_mat,ssrep newse_mat = origse_mat,serep N=cols(newss_mat) //no. of studies + new one weights = 1:/(newse_mat:^2) //weights of each study sumweights = quadrowsum(weights) //sum of weights for each M-A weights_times_ss = weights:*newss_mat sum_weights_times_ss = quadrowsum(weights_times_ss) //sum of weights* ss fsum = sum_weights_times_ss:/sumweights ss_minus_fsum_2 = (newss_mat:-J(1,cols(newss_mat),fsum)):^2 q = ss_minus_fsum_2:/(newse_mat:^2) Q = quadrowsum(q) df = N:-1 weights2 = weights:^2 sum_weights2 = quadrowsum(weights2) comp = (Q:-df):/(sumweights:-(sum_weights2:/sumweights)) tau2 = J(nint2,1,0) for (i=1; i<=rows(comp); i++) { if (comp[i,1]>0){ tau2[i,1]=comp[i,1] } } if (method=="fixedi") { wt = 1:/(newse_mat:^2) v1 = wt:*wt:*(newse_mat:^2) } else { wt = 1:/((newse_mat:^2):+tau2) v1 = wt:*wt:*((newse_mat:^2):+tau2) } wt_ss = wt:*newss_mat sum_wt_ss = rowsum(wt_ss) sum_wt = rowsum(wt) summ = sum_wt_ss:/sum_wt //pooled estimate of each M-A sum_v1 = rowsum(v1) varsum = sum_v1:/(sum_wt:^2) //variance of each M-A comb_ss_mat = J(nintminus1:*4,nintminus1,.) comb_var_mat = J(nintminus1:*4,nintminus1,.) //Build matrix of effect estimates index1 = 1 index2 = nintminus1 index3 = 2 index4 = nint //first 2 rows: comb_ss_mat[1,(1::nintminus1)] = summ[(1::nintminus1),1]' comb_ss_mat[2,(1::nintminus1)] = summ[(2::nint),1]' comb_var_mat[1,(1::nintminus1)] = varsum[(1::nintminus1),1]' comb_var_mat[2,(1::nintminus1)] = varsum[(2::nint),1]' row1 = 3 row2 = 4 row3 = 5 row4 = 6 //add 4 sub1start = nint:+2 //12 sub1end = nint:+nint //20 sub2start = nint:+1 //11 sub2end = nint:+nint:-1 //19 for (i=1; i<=(nint:-2); i++) { comb_ss_mat[row1,(1::nintminus1)] = summ[(sub1start::sub1end),1]' comb_ss_mat[(row2::row3),(1::nintminus1)] = J(2,1,summ[(sub2start::sub2end),1]') comb_ss_mat[row4,(1::nintminus1)] = summ[(sub1start::sub1end),1]' comb_var_mat[row1,(1::nintminus1)] = varsum[(sub1start::sub1end),1]' comb_var_mat[(row2::row3),(1::nintminus1)] = J(2,1,varsum[(sub2start::sub2end),1]') comb_var_mat[row4,(1::nintminus1)] = varsum[(sub1start::sub1end),1]' row1 = row1:+4 row2 = row2:+4 row3 = row3:+4 row4 = row4:+4 sub1start = sub1start:+nint sub1end = sub1end:+nint sub2start = sub2start:+nint sub2end = sub2end:+nint } //final 2 rows comb_ss_mat[row1,(1::nintminus1)] = summ[(sub1start::sub1end),1]' comb_ss_mat[row2,(1::nintminus1)] = summ[(sub2start::sub2end),1]' comb_var_mat[row1,(1::nintminus1)] = varsum[(sub1start::sub1end),1]' comb_var_mat[row2,(1::nintminus1)] = varsum[(sub2start::sub2end),1]' //indmat now contains effect estimates of all combinations of the coordinates. //need to do the same for the variance upp_ci_mat = comb_ss_mat:+J(nintminus1:*4,1,siglevel):*sqrt(comb_var_mat) low_ci_mat = comb_ss_mat:-J(nintminus1:*4,1,siglevel):*sqrt(comb_var_mat) mat012 = J(nintminus1:*4,nintminus1,0) pow = st_numscalar("power") if (pow==0) { for (i=1; i<=rows(comb_ss_mat); i++) { for (j=1; j<=cols(comb_ss_mat); j++) { if (upp_ci_mat[i,j]null & upp_ci_mat[i,j]>null) { mat012[i,j] = 2 } } } } else if (pow==1) { for (i=1; i<=rows(comb_ss_mat); i++) { for (j=1; j<=cols(comb_ss_mat); j++) { if (low_ci_mat[i,j]>=null) { mat012[i,j] = 1 } } } } else if (pow==2) { for (i=1; i<=rows(comb_ss_mat); i++) { for (j=1; j<=cols(comb_ss_mat); j++) { if (upp_ci_mat[i,j]>=null) { mat012[i,j] = 1 } } } } else if (pow==3) { ciwidths = upp_ci_mat :- low_ci_mat for (i=1; i<=rows(comb_ss_mat); i++) { for (j=1; j<=cols(comb_ss_mat); j++) { if (ciwidths[i,j]>=null) { mat012[i,j] = 1 } } } } else if (pow==4) { loe1 = st_numscalar("loe1") loe2 = st_numscalar("loe2") for (i=1; i<=rows(comb_ss_mat); i++) { for (j=1; j<=cols(comb_ss_mat); j++) { if (upp_ci_mat[i,j]>=null) { if (upp_ci_mat[i,j]=loe2) { if (low_ci_mat[i,j]=loe2) { if (low_ci_mat[i,j]=loe1) { mat012[i,j] = 1 } } } if (upp_ci_mat[i,j]=loe1) { if (low_ci_mat[i,j]=loe2) { if (low_ci_mat[i,j]>=null) { if (low_ci_mat[i,j]=loe1) { mat012[i,j] = 4 } } if (low_ci_mat[i,j]>=null) { if (upp_ci_mat[i,j]=null) { if (upp_ci_mat[i,j]=loe1) { if (low_ci_mat[i,j]=loe2) { if (upp_ci_mat[i,j]>=loe2) { mat012[i,j] = 8 } } } } } st_view(out=.,.,tokens(st_local("varnames")),st_local("tempind")) out[.,.] = mat012[.,.] } end