*! version 1.4.6 15dec2022 JM. Domenech, R. Sesma /* META-ANALYSIS: OR,RR,RD,IR,ID,B,MD,R COMBINED */ program define mar, rclass version 12 syntax anything(id="risk measure and type of data") [if] [in], [TABles INFlu METhod(string) CUM(string) /* */ BEGG EGGer MACask(string) NUS(numlist max=1 >0 <1) INTerval(string) FORest(string) HISTogram RADial /* */ cumulative abbe FUNnel FUNNELH GALBraith Level(numlist max=1 >50 <100) COHen HEDges GLAss keepmissing /* */ keep zero(string) exp date(string) nst(string)] *Tokenize anything parameter tokenize `anything' local type `1' local stat `2' *Check parameters if ("`stat'"!="or" & "`stat'"!="rr" & "`stat'"!="ir" & "`stat'"!="rd" & /* */ "`stat'"!="id" & "`stat'"!="b" & "`stat'"!="md" & "`stat'"!="r") print_error "statistic `stat' invalid" if ("`type'"!="raw" & "`type'"!="sst" & "`type'"!="mix") print_error "type `type' invalid" if ("`method'"=="") local method "fem" if ("`method'"!="fem" & "`method'"!="rem" & "`method'"!="mh" & "`method'"!="peto") print_error "method `method' invalid" if ("`type'"!="raw" & ("`method'"=="mh" | "`method'"=="peto")) print_error "method `method' is only available with raw data" if ("`method'"=="peto" & "`stat'"!="or") print_error "method `method' is only compatible with or statistic" if (("`method'"=="peto" | "`method'"=="mh") & ("`stat'"=="rd" | "`stat'"=="id")) print_error "method `method' is not compatible with `stat' statistic" if ("`level'"=="") local level 95 if ("`zero'"=="") local zero "n" if ("`zero'"!="n" & "`zero'"!="c" & "`zero'"!="p" & "`zero'"!="e") print_error "zero invalid" if ("`exp'"!="" & "`stat'"!="b") print_error "exp option incompatible with `stat' statistic" if ("`date'"=="") local date "YMD" if ("`date'"!="YMD" & "`date'"!="YM" & "`date'"!="Y") print_error "date invalid" if ("`cumulative'"!="" & "`cum'"=="") print_error "cumulative analisys (cum option) is needed for the cumulative graphic" if ("`type'"!="raw" & "`abbe'"!="") print_error "raw data type is needed for l'Abbe plot" if ("`type'"=="raw" & "`stat'"=="md" & "`abbe'"!="") print_error "md statistic is not compatible with l'Abbe plot" if ("`type'"=="raw" & "`stat'"=="b") print_error "b statistic is not compatible with raw data" if ("`type'"=="raw" & "`stat'"=="r") print_error "r statistic is not compatible with raw data" if ("`type'"=="mix" & "`stat'"=="r") print_error "r statistic is not compatible with mixed data" if ("`cohen'"!="" & ("`hedges'"!="" | "`glass'"!="")) print_error "only one of cohen, hedges, glass options is allowed" if ("`hedges'"!="" & "`glass'"!="") print_error "only one of cohen, hedges, glass options is allowed" if ("`cohen'"=="" & "`hedges'"=="" & "`glass'"=="") local md_std 0 //No standardization if ("`cohen'"!="") local md_std 1 //Mean difference standardization: Cohen's d if ("`hedges'"!="") local md_std 2 //Mean difference standardization: Hedges' g if ("`glass'"!="") local md_std 3 //Mean difference standardization: Glass's delta tempname vzeros *Study is mandatory checkvars study, errmsg local study_type : type study if (substr("`study_type'",1,3)!="str") qui tostring(study), replace force //Transform numeric study variable to string marksample touse //ifin marksample preserve //Preserve original data quietly keep if `touse' //Select cases if (_N == 0) print_error "no data available" //Check there's data on the dataset *Check for auxiliary variables: date var and interval, forest plot vars checkvars date, local(ldate) local ldate = cond(`ldate'==1,"date","") if ("`interval'"!="" & "`interval'"!="se") checkvars `interval', errmsg if ("`forest'"!="" & "`forest'"!="se") checkvars `forest', errmsg local lzeros 0 if ("`type'"=="raw") { data_raw `stat', zero(`zero') vzeros(`vzeros') md_std(`md_std') level(`level') local lzeros = r(zeros) local ormh_zero = r(ormh_zero) local md_exists = r(md_exists) } if ("`type'"=="sst") { if ("`stat'"=="r") { tempname rho capture g `rho' = r } data_sst `stat',level(`level') } local keep_mix = "" if ("`type'"=="mix") { restore quietly { *Look for variable names in the dataset tempname E matrix `E' = J(1,21,.) local i 1 foreach v in `stat' a1 a0 b1 b0 n1 n0 n i1 i0 t1 t0 se lo up cl pexp pexp_nc pres pres_ue _mar_ { capture confirm variable `v' matrix `E'[1,`i'] = (!_rc) local i = `i'+1 } /*E gets existence of variables: 1: statistic, 2: a1, 3: a0, 4: b1, 5: b0, 6: n1, 7: n0, 8: n, 9: i1, 10: i0 11: t1, 12: t0, 13: se, 14: lo, 15: up, 16: cl, 17: pexp, 18: pexp_nc, 19: pres, 20: pres_ue 21: _mar_*/ *Check for data errors if the dataset is not keeped if (!`E'[1,21]) { quietly { capture drop __err g __err = 0 if (`E'[1,1] & `E'[1,2] & `E'[1,3]) replace __err = 1 if (`stat'<. & (a1<. | a0<.)) if (`E'[1,2] & `E'[1,3] & `E'[1,13]) replace __err = 2 if (se<. & (a1<. | a0<.)) if (`E'[1,2] & `E'[1,3] & `E'[1,14] & `E'[1,15]) replace __err = 3 if ((lo<. | up<.) & (a1<. | a0<.)) if (`E'[1,13] & `E'[1,14] & `E'[1,15]) replace __err = 4 if ((lo<. | up<.) & se<.) if (`E'[1,2] & `E'[1,3] & `E'[1,8]) replace __err = 5 if (n<. & (a1<. | a0<.)) if (`E'[1,4] & `E'[1,5] & `E'[1,6] & `E'[1,7]) replace __err = 6 if ((b1<. | b0<.) & (n1<. | n0<.)) if (`E'[1,9] & `E'[1,10] & `E'[1,11] & `E'[1,12]) replace __err = 7 if ((i1<. | i0<.) & (t1<. | t0<.)) if (`E'[1,2] & `E'[1,3] & `E'[1,17]) replace __err = 8 if (pexp<. & (a1<. | a0<.)) if (`E'[1,2] & `E'[1,3] & `E'[1,18]) replace __err = 8 if (pexp_nc<. & (a1<. | a0<.)) if (`E'[1,2] & `E'[1,3] & `E'[1,19]) replace __err = 8 if (pres<. & (a1<. | a0<.)) if (`E'[1,2] & `E'[1,3] & `E'[1,20]) replace __err = 8 if (pres_ue<. & (a1<. | a0<.)) label variable __err "Data error" label define derr 0 "Correct" 1 "Error: a1,a0 incompatible with `stat'" 2 "Error: a1,a0 incompatible with se" /// 3 "Error: a1,a0 incompatible with lo,up" 4 "Error: se incompatible with lo,up" /// 5 "Error: n incompatible with a1,a0" 6 "Error: b1,b0 incompatible with n1,n0" /// 7 "Error: i1,i0 incompatible with t1,t0" 8 "Error: pexp,pexp_nc,pres,pres_ue incompatible with a1,a0", replace label values __err derr summarize __err if (r(sum)>0) print_error "some observations have data errors -- see __err variable" drop __err label drop derr } } preserve *Mixed type: stat variable is needed if (!`E'[1,1]) print_error "missing `stat' variable for mixed type of data" if (!`E'[1,13]) g se=. //If se doesn't exists, it's created missing if ((`E'[1,2] & !`E'[1,3]) | (!`E'[1,2] & `E'[1,3])) print_error "for raw data observations variables a1 and a0 are needed" if (`E'[1,2] & `E'[1,3]) { if ("`stat'"=="b" | "`stat'"=="md") print_error "raw data observations are incompatible with `stat' statistic" *Raw data observations: complete the missing data and compute statistic if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="rd") { if ((`E'[1,4] & !`E'[1,5]) | (!`E'[1,4] & `E'[1,5])) print_error "missing variable b1 or b0 for raw data observations" if (!`E'[1,4] & !`E'[1,5]) { g b1 = . //Create b1, b0 variables missing g b0 = . } if ((`E'[1,6] & !`E'[1,7]) | (!`E'[1,6] & `E'[1,7])) print_error "missing variable n1 or n0 for raw data observations" if (!`E'[1,6] & !`E'[1,7]) { g n1 = . //Create n1, n0 variables missing g n0 = . } replace n1 = a1 + b1 if (n1==.) replace n0 = a0 + b0 if (n0==.) replace b1 = n1 - a1 if (b1==.) replace b0 = n0 - a0 if (b0==.) tempname vraw g `vraw' = !missing(a1) //raw observations must have a valid value on a1 zero_correction `zero', vzeros(`vzeros') mixed vraw(`vraw') local lzeros = r(zeros) local ormh_zero = r(ormh) if (`lzeros') local keep_mix = "a1 a0 b1 b0 n1 n0 n" if ("`stat'" == "or") { replace or = (a1 / b1) / (a0 / b0) if or==. replace se = sqrt(1/a1 + 1/b1 + 1/a0 + 1/b0) if (se==. & a1!=.) } if ("`stat'" == "rr") { replace rr = (a1 / n1) / (a0 / n0) if rr==. replace se = sqrt(1/a1 - 1/n1 + 1/a0 - 1/n0) if (se==. & a1!=.) } if ("`stat'" == "rd") { replace rd = (a1 / n1) - (a0 / n0) if rd==. replace se = sqrt(((a1/n1)*(1-(a1/n1)))/n1 + ((a0/n0)*(1-(a0/n0)))/n0) if (se==. & a1!=.) } if (`E'[1,8]) replace n = n1 + n0 if (n==.) if (`E'[1,17]) replace pexp = n1/(n1+n0) if (pexp==.) if (`E'[1,18]) replace pexp_nc = b1/(b0+b1) if (pexp_nc==.) if (`E'[1,19]) replace pres = (a0+a1)/(n1+n0) if (pres==.) if (`E'[1,20]) replace pres_ue = a0/n0 if (pres_ue==.) if (`E'[1,17] & `lzeros') local keep_mix = "`keep_mix' pexp" if (`E'[1,18] & `lzeros') local keep_mix = "`keep_mix' pexp_nc" if (`E'[1,19] & `lzeros') local keep_mix = "`keep_mix' pres" if (`E'[1,20] & `lzeros') local keep_mix = "`keep_mix' pres_ue" } if ("`stat'"=="ir" | "`stat'"=="id") { if ((`E'[1,11] & !`E'[1,12]) | (!`E'[1,11] & `E'[1,12])) print_error "missing variable b1 or b0 for raw data observations" if (!`E'[1,11] & !`E'[1,12]) { g t1 = . //Create b1, b0 variables missing g t0 = . } if ((`E'[1,9] & !`E'[1,10]) | (!`E'[1,9] & `E'[1,10])) print_error "missing variable n1 or n0 for raw data observations" if (!`E'[1,9] & !`E'[1,10]) { g i1 = . //Create n1, n0 variables missing g i0 = . } replace i1 = a1/t1 if (i1==.) replace i0 = a0/t0 if (i0==.) replace t1 = a1/i1 if (t1==.) replace t0 = a0/i0 if (t0==.) if ("`stat'" == "ir") { replace ir = (a1/t1) / (a0/t0) if ir==. replace se = sqrt(1/a1 + 1/a0) if (se==. & a1!=.) } if ("`stat'" == "id") { replace id = (a1/t1) - (a0/t0) if id==. replace se = sqrt(a1/t1^2 + a0/t0^2) if (se==. & a1!=.) } if (`E'[1,17]) replace pexp = t1/(t1+t0) if (pexp==.) if (`E'[1,19]) replace pres = (a0+a1)/(t1+t0) if (pres==.) if (`E'[1,20]) replace pres_ue = a0/t0 if (pres_ue==.) } } if ((`E'[1,14] & !`E'[1,15]) | (!`E'[1,14] & `E'[1,15])) print_error "missing variable lo or up for sst data observations" if (!`E'[1,14] & !`E'[1,15]) { qui gen lo = . //Create lo, up variables missing qui gen up = . } if (!`E'[1,16]) g cl=. replace cl=95 if cl==. if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir") { replace se = (ln(up)-ln(lo))/(2*invnormal((cl+100)/200)) if se==. replace lo = `stat'*exp(-invnormal((cl+100)/200)*se) if lo==. replace up = `stat'*exp(invnormal((cl+100)/200)*se) if up==. } else { replace se = (up-lo)/(2*invnormal((cl+100)/200)) if se==. replace lo = `stat' - invnormal((`level'+100)/200)*se if lo==. replace up = `stat' + invnormal((`level'+100)/200)*se if up==. } } *From this point on the command must behave as if type = sst local type = "sst" } *Cumulative sort variable qui gen __cumsort = . local ldesc 1 //Descendent (greater -> lower) order by default if ("`cum'"!="") { local lreverse = (substr("`cum'",1,1)=="-") //Reverse order? if (`lreverse') local cum = substr("`cum'",2,.) //Variable name without the reverse sign if ("`cum'"=="date" | "`cum'"=="n") local ldesc 0 //date | n variable: sort ascedending (lower -> greater) if (`lreverse') local ldesc = !`ldesc' //Reverse order qui replace __cumsort = `cum' } *Check for missing variables local ntotal = _N if ("`keepmissing'"=="") { *By default, studies with missing data are deleted (keepmissing option prevents the default) if ("`type'"=="raw") { if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="rd") local lvars = "a1,b1,n1,a0,b0,n0" if ("`stat'"=="ir" | "`stat'"=="id") local lvars = "a1,i1,t1,a0,i0,t0" if ("`stat'"=="md") local lvars = "m1,sd1,n1,m0,sd0,n0" } if ("`type'"=="sst") { if ("`stat'"!="r") local lvars = "`stat',se,lo,up" else local lvars = "r,n" } qui keep if !missing(`lvars') } local nmiss = _N local ldeleted = `ntotal'!=`nmiss' local ndeleted = `ntotal'-`nmiss' local len = _N local len_1 = _N+1 local raw_md = ("`type'"=="raw" & "`stat'"=="md") //md studies with raw data if (`raw_md') local type = "sst" //md studies with raw data must behave like sst data //Results using mata mata: init("`type'","`stat'",`len',`level',("`influ'"!=""),"`method'",("`cum'"!=""),`ldesc') //Read mata results into stata matrix tempname d st stp rw sw peto inv invlo invup mh iv ivr inf cu matrix `d' = __mar_data matrix `st' = __mar_st matrix `stp' = __mar_stp matrix `rw' = __mar_rw matrix `sw' = __mar_sw matrix `peto' = __mar_peto matrix `mh' = __mar_mh matrix `iv' = __mar_iv matrix `ivr' = __mar_ivr matrix `inf' = __mar_influ matrix `cu' = __mar_cum ****PRINT results qui g str dates = "" if ("`ldate'"!="") { *if date not daily, it's the numeric year local fmt : format date if (substr("`fmt'",1,3)=="%td" & "`date'"=="YMD") local fmt = "%tdCCYY-NN-DD" if (substr("`fmt'",1,3)=="%td" & "`date'"=="YM") local fmt = "%tdCCYY-NN" if (substr("`fmt'",1,3)=="%td" & "`date'"=="Y") local fmt = "%tdCCYY" if (substr("`fmt'",1,3)!="%td") local fmt = "%4.0f" qui replace dates = string(date,"`fmt'") } *Print caption di as res "META-ANALYSIS: OR,RR,RD,IR,ID,B,MD,R COMBINED" if ("`nst'"!="") di as txt "{bf:STUDY:} `nst'" *Get study names & zero notes forvalues i = 1/`len_1'{ local n`i' = cond(`i'<`len_1',abbrev(study[`i'],12),"GLOBAL") local z`i' " " if ("`type'"=="raw") { if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="rd") { if ("`zero'"!="n" & `i'!=`len_1' & `vzeros'[`i']==1) local z`i' "*" } } } local title = upper("`stat'") if ("`type'"=="raw" & "`tables'"!="") { //Print 2x2 tables (if available and asked) di if ("`stat'"=="or") di as txt _col(25)"{c |}" _col(37)"{c |}" _col(49)"{c |}" _col(61)"{c |}{center 11:OR}{c |}{center 11:OR_peto}" di as txt "{lalign 24:STUDY}{c |} UnExposed {c |} Exposed {c |}{center 11:TOTAL}{c |}" _c if ("`stat'"=="or") di as txt "{center 11:(SE_lnOR)}{c |}{center 11:(SE_lnORp)}" if ("`stat'"=="rr" | "`stat'"=="ir") di as txt " `title' (SE_ln`title')" if ("`stat'"=="rd" | "`stat'"=="id") di as txt " `title' (SE)" di as txt "{hline 24}{c +}{hline 11}{c +}{hline 11}{c +}{hline 11}{c +}" _c if ("`stat'"=="or") di as txt "{hline 11}{c +}{hline 11}" if ("`stat'"=="rr" | "`stat'"=="ir") di as txt "{hline 13}" if ("`stat'"=="rd" | "`stat'"=="id") di as txt "{hline 11}" forvalues i = 1/`len_1'{ di as txt "{lalign 14:`n`i''}" _col(19) "Cases {c |} " %9.0g `d'[`i',2] " {c |} " %9.0g `d'[`i',1] " {c |} " /* */ %9.0g `d'[`i',7] " {c |} " as res %9.0g `st'[`i',1] as txt "`z`i''" _c if ("`stat'"=="or") di as res "{c |} " %9.0g `stp'[`i',1] as txt "`z`i''" else di as txt " " di as txt %-10s dates[`i'] _c if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="rd") di _col(16) "NonCases {c |} " _c if ("`stat'"=="ir" | "`stat'"=="id") di _col(13) "Person-time {c |} " _c di as txt %9.0g `d'[`i',4] " {c |} " %9.0g `d'[`i',3] " {c |} " %9.0g `d'[`i',8] as txt " {c |}(" as res %9.0g `st'[`i',2] as txt ")" _c if ("`stat'"=="or") di as txt "{c |}(" as res %9.0g `stp'[`i',2] as txt ")" else di as txt " " di as txt _col(25) "{c LT}{hline 11}{c +}{hline 11}{c +}{hline 11}{c RT}" _c if ("`stat'"=="or") di as txt _col(73) "{c |}" else di as txt " " if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="rd") di as txt _col(19) "TOTAL {c |} " _c if ("`stat'"=="ir" | "`stat'"=="id") di _col(10) "Incidence Rate {c |} " _c di as txt %9.0g `d'[`i',6] " {c |} " %9.0g `d'[`i',5] " {c |} " %9.0g `d'[`i',9] /* */ as txt " {c |} p= " as res %6.4f `st'[`i',5] _c if ("`stat'"=="or") di as txt " {c |} p= " as res %6.4f `stp'[`i',5] else di as txt " " local sep "+" if (`i'==`len_1') local sep "BT" di as txt "{hline 24}{c `sep'}{hline 11}{c `sep'}{hline 11}{c `sep'}{hline 11}{c `sep'}{hline 11}" _c if ("`stat'"=="or") di as txt "{c `sep'}{hline 11}" else di as txt " " } if (("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="rd") & (`lzeros' & "`zero'"!="n")) { if ("`zero'"=="c") di as txt "(*) Computed with a constant continuity correction (k=0.5)" if ("`zero'"=="p") di as txt "(*) Computed with a proportional continuity correction" if ("`zero'"=="e") di as txt "(*) Computed with an empirical continuity correction (ORmh= " %9.0g `ormh_zero' ")" di as txt " for studies with zero events." } } //Print statistics table /*RAW WEIGHTS: -or: Peto MH IV_FEM IV_REM -rr, ir: MH IV_FEM IV_REM -rd, id: IV_FEM IV_REM if there are uncorrected zeros: IV_FEM & IV_REM weights not printed for or, rr SST WEIGHTS: IV_FEM IV_REM*/ local zunc = ("`type'"=="raw" & ("`stat'"=="or" | "`stat'"=="rr")) & (`lzeros' & "`zero'"=="n") di if ("`type'"=="raw" & "`stat'"=="or") { if ("`method'"=="peto") local title = "OR_peto" local pos = cond(!`zunc',66,60) di as txt "{hline 24}{c TT}{hline 11}{c TT}{hline 21}{c TT}{hline 17}" cond(!`zunc',"{hline 14}"," ") di as txt _col(25) "{c |}" _col(37) "{c |}" %4.0g `level' "% Conf. Interval {c |}" _col(`pos') "Relative WEIGHTS" di as txt "STUDY" _col(14) upper("`ldate'") _col(25) "{c |}{center 11:`title'}{c |} {ralign 9:Lower} {ralign 9:Upper} {c |}" /* */ cond(!`zunc'," Peto MH IV_FEM IV_REM"," Peto MH") di as txt "{hline 24}{c +}{hline 11}{c +}{hline 21}{c +}" /* */ cond(!`zunc',"{hline 7}{c TT}{hline 7}{c TT}{hline 7}{c TT}{hline 7}","{hline 7}{c TT}{hline 9}") forvalues i = 1/`len_1'{ if (`i'==`len_1') di as txt "{hline 24}{c +}{hline 11}{c +}{hline 21}{c +}" /* */ cond(!`zunc',"{hline 7}{c +}{hline 7}{c +}{hline 7}{c +}{hline 7}","{hline 7}{c +}{hline 9}") di as txt "{lalign 12:`n`i''} " %-10s dates[`i'] _col(25) "{c |} " /* */ as res %9.0g cond("`method'"!="peto",`st'[`i',1],`stp'[`i',1]) "`z`i''{c |} " /* */ %9.0g cond("`method'"!="peto",`st'[`i',3],`stp'[`i',3]) " " /* */ %9.0g cond("`method'"!="peto",`st'[`i',4],`stp'[`i',4]) " {c |}" _c *Print relative weights local lnum 3 2 if (!`zunc') local lnum `lnum' 1 4 foreach j of numlist `lnum' { if (`i'==`len_1') local p = 100 else local p = `rw'[`i',`j'] print_pct `p' if (`j'==3 | (`j'==2 & !`zunc') | `j'==1) di as txt " {c |}" _c } di " " } di as txt "{hline 24}{c BT}{hline 11}{c BT}{hline 21}{c +}" /* */ cond(!`zunc',"{hline 7}{c +}{hline 7}{c +}{hline 7}{c +}{hline 7}","{hline 7}{c +}{hline 9}") di as txt "{ralign 58:TOTAL Weights }{c |}" as res %6.0g `sw'[1,3] " {c |}" %6.0g `sw'[1,2] _c if (!`zunc') di as res " {c |}" %6.0g `sw'[1,1] " {c |}" %6.0g `sw'[1,4] else di " " di as txt "{hline 58}{c BT}" _c di as txt cond(!`zunc',"{hline 7}{c BT}{hline 7}{c BT}{hline 7}{c BT}{hline 7}","{hline 7}{c BT}{hline 9}") } if ("`type'"=="raw" & ("`stat'"=="rr" | "`stat'"=="ir")) { di as txt "{hline 24}{c TT}{hline 11}{c TT}{hline 21}{c TT}{hline 12}" cond(!`zunc',"{hline 11}"," ") di as txt _col(25) "{c |}" _col(37) "{c |}" %4.0g `level' "% Conf. Interval {c |} " _c if (!`zunc') di _col(61) "Relative WEIGHTS" else di as txt "MH Relative" di as txt "STUDY" _col(14) upper("`ldate'") _col(25) "{c |}{center 11:`title'}{c |} {ralign 9:Lower} {ralign 9:Upper} {c |}" /* */ cond(!`zunc'," MH IV_FEM IV_REM"," Weights") di as txt "{hline 24}{c +}{hline 11}{c +}{hline 21}{c +}" cond(!`zunc',"{hline 7}{c TT}{hline 7}{c TT}{hline 7}","{hline 12}") forvalues i = 1/`len_1'{ if (`i'==`len_1') di as txt "{hline 24}{c +}{hline 11}{c +}{hline 21}{c +}" cond(!`zunc',"{hline 7}{c +}{hline 7}{c +}{hline 7}","{hline 12}") di as txt "{lalign 12:`n`i''} " %-10s dates[`i'] _col(25) "{c |} " as res %9.0g `st'[`i',1] "`z`i''{c |} " /* */ %9.0g `st'[`i',3] " " %9.0g `st'[`i',4] " {c |}" _c *Print relative weights local lnum 2 if (!`zunc') local lnum `lnum' 1 4 foreach j of numlist `lnum' { if (`i'==`len_1') local p = 100 else local p = `rw'[`i',`j'] if (`zunc') di " " _c print_pct `p' if ((`j'==2 & !`zunc') | `j'==1) di as txt " {c |}" _c } di " " } di as txt "{hline 24}{c BT}{hline 11}{c BT}{hline 21}{c +}" /* */ cond(!`zunc',"{hline 7}{c +}{hline 7}{c +}{hline 7}","{hline 12}") local fmt = cond(!`zunc',"%6.0g","%8.0g") di as txt "{ralign 58:TOTAL Weights }{c |}" as res `fmt' `sw'[1,2] _c if (!`zunc') di as res " {c |}" `fmt' `sw'[1,1] " {c |}" `fmt' `sw'[1,4] else di " " di as txt "{hline 58}{c BT}" _c di as txt cond(!`zunc',"{hline 7}{c BT}{hline 7}{c BT}{hline 7}","{hline 12}") } if ("`type'"=="raw" & ("`stat'"=="rd" | "`stat'"=="id")) { di as txt "{hline 24}{c TT}{hline 11}{c TT}{hline 21}{c TT}{hline 23}" di as txt _col(25) "{c |}" _col(37) "{c |}" %4.0g `level' "% Conf. Interval {c |} Relative IoV Weights" di as txt "STUDY" _col(14) upper("`ldate'") _col(25) "{c |}{center 11:`title'}{c |} {ralign 9:Lower} {ralign 9:Upper} {c |}" /* */ " Fixed Eff Random Eff" di as txt "{hline 24}{c +}{hline 11}{c +}{hline 21}{c +}{hline 11}{c +}{hline 11}" forvalues i = 1/`len_1'{ if (`i'==`len_1') di as txt "{hline 24}{c +}{hline 11}{c +}{hline 21}{c +}{hline 11}{c +}{hline 11}" di as txt "{lalign 12:`n`i''} " %-10s dates[`i'] _col(25) "{c |} " as res %9.0g `st'[`i',1] "`z`i''{c |} " /* */ %9.0g `st'[`i',3] " " %9.0g `st'[`i',4] " {c |}" _c *Print relative weights foreach j of numlist 1 4 { if (`i'==`len_1') local p = 100 else local p = `rw'[`i',`j'] di " " _c print_pct `p' if (`j'==1) di as txt " {c |}" _c } di " " } di as txt "{hline 24}{c BT}{hline 11}{c BT}{hline 21}{c +}{hline 11}{c +}{hline 11}" di as txt "{ralign 58:TOTAL Weights }{c |}" as res %10.0g `sw'[1,1] " {c |}" %10.0g `sw'[1,4] di as txt "{hline 58}{c BT}{hline 11}{c BT}{hline 11}" } if (("`type'"=="raw") & ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="rd") & (`lzeros' & "`zero'"!="n")) { if ("`zero'"=="c") di as txt "(*) Computed with a constant continuity correction (k=0.5)" if ("`zero'"=="p") di as txt "(*) Computed with a proportional continuity correction" if ("`zero'"=="e") di as txt "(*) Computed with an empirical continuity correction (ORmh= " %9.0g `ormh_zero' ")" di as txt " for studies with zero events." } if ("`type'"=="sst") { if ("`stat'"=="b") local title = cond("`exp'"==""," b","Exp(b)") if ("`stat'"=="r") local title = " r" if ("`stat'"=="md") local title = cond(`md_std'==1,"d",cond(`md_std'==2,"g",cond(`md_std'==3,"Delta","MD"))) local md_type = "" if ("`stat'"=="md") local md_type = cond(`md_std'==1,"Cohen's",cond(`md_std'==2,"Hedges'",cond(`md_std'==3,"Glass's",""))) di as txt "{hline 24}{c TT}{hline 11}{c TT}{hline 21}{c TT}{hline 5}{c TT}{hline 23}" di as txt _col(25) "{c |}{center 11:`md_type'}{c |} Confidence Interval {c |}" _col(65) "{c |} Relative IoV Weights" di as txt "STUDY" _col(14) upper("`ldate'") _col(25) "{c |}{center 11:`title'}{c |} {ralign 9:Lower} {ralign 9:Upper} {c |}" /* */ " CL {c |} Fixed Eff Random Eff" di as txt "{hline 24}{c +}{hline 11}{c +}{hline 21}{c +}{hline 5}{c +}{hline 11}{c TT}{hline 11}" forvalues i = 1/`len'{ di as txt "{lalign 12:`n`i''} " %-10s dates[`i'] _col(25) "{c |} " /* */ as res %9.0g cond("`stat'"=="b" & "`exp'"!="",exp(`st'[`i',1]),`st'[`i',1]) " {c |} " /* */ %9.0g cond("`stat'"=="b" & "`exp'"!="",exp(`st'[`i',3]),`st'[`i',3]) " " /* */ %9.0g cond("`stat'"=="b" & "`exp'"!="",exp(`st'[`i',4]),`st'[`i',4]) " " /* */ "{c |}" %4.0g cl[`i'] " {c |}" _c *Print relative weights foreach j of numlist 1 4 { local p = `rw'[`i',`j'] di " " _c print_pct `p' if (`j'==1) di as txt " {c |}" _c } di " " } di as txt "{hline 24}{c BT}{hline 11}{c BT}{hline 21}{c BT}{hline 5}{c +}{hline 11}{c +}{hline 11}" di as txt "{ralign 64:TOTAL}{c |}" as res %9.0f 100 as txt "% {c |}" as res %9.0f 100 as txt "%" di as txt _col(65) "{c LT}{hline 11}{c +}{hline 11}" di as txt "{ralign 64:TOTAL Weights }{c |}" as res %10.0g `sw'[1,1] " {c |}" %10.0g `sw'[1,4] di as txt "{hline 64}{c BT}{hline 11}{c BT}{hline 11}" } if (`ldeleted') { di as txt "WARNING: `ndeleted' studies were excluded due to missing data" } ***Print estimations table if !("`type'"=="raw" & "`stat'"=="rr" & `lzeros' & "`zero'"=="n") { *For raw data & rr statistic with zero events, there's no results available *For any other combination, print estimations results di di as txt "{hline 24}{c TT}{hline 12}{c TT}{hline 9}{c TT}{hline 22}" di as txt _col(25)"{c |}" _col(38)"{c |}" _col(48)"{c |} " %4.0g `level' "% Conf. Interval" di as txt _col(25)"{c |} Estimation {c |} Signif. {c |} {ralign 9:Lower} {ralign 9:Upper}" di as txt "{hline 24}{c +}{hline 12}{c +}{hline 9}{c +}{hline 22}" } if ("`type'"=="raw" & "`stat'"=="or") local title = "OR" local z " " if ("`type'"=="raw" & `lzeros' & "`zero'"!="n") local z "*" if ("`type'"=="raw" & "`stat'"=="or") { *raw data + stat=or: print Peto combined results di as txt "PETO COMBINED: OR_peto`z'{c |} " as res %9.0g `peto'[1,1] " {c |} " /* */ %6.4f `peto'[1,5] " {c |} " %9.0g `peto'[1,3] " " %9.0g `peto'[1,4] if (`peto'[1,1] < 1) di as txt "{ralign 23:1/OR} {c |} " as res %9.0g 1/`peto'[1,1] " {c |} " /* */ _col(48)"{c |} " %9.0g 1/`peto'[1,4] " " %9.0g 1/`peto'[1,3] di as txt "{ralign 24:SE(lnOR_peto)}{c |} " as res %9.0g `peto'[1,2] " {c |} " _col(48) "{c |}" di as txt " Homogeneity Chi-Square {c |} " as res %9.0g `peto'[1,6] " {c |} " %6.4f `peto'[1,7] " {c |}" return scalar stat_peto = `peto'[1,1] //RETURNED RESULTS return scalar se_peto = `peto'[1,2] return scalar stat_peto_lo = `peto'[1,3] return scalar stat_peto_up = `peto'[1,4] return scalar w_peto = `sw'[1,3] } if ("`type'"=="raw" & ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir") & !(`lzeros' & "`zero'"=="n")) { if ("`stat'"=="or") di as txt "{hline 24}{c +}{hline 12}{c +}{hline 9}{c +}{hline 22}" di as txt "M-H COMBINED:{ralign 10:`title'} {c |} " as res %9.0g `mh'[1,1] " {c |} " /* */ %6.4f `mh'[1,5] " {c |} " %9.0g `mh'[1,3] " " %9.0g `mh'[1,4] if (`mh'[1,1] < 1) di as txt "{ralign 23:1/`title'} {c |} " as res %9.0g 1/`mh'[1,1] " {c |} " /* */ _col(48) "{c |} " %9.0g 1/`mh'[1,4] " " %9.0g 1/`mh'[1,3] di as txt "{ralign 24:SE(ln`title')}{c |} " as res %9.0g `mh'[1,2] " {c |}" _col(48) "{c |}" di as txt " Homogeneity Chi-Square {c |} " as res %9.0g `mh'[1,6] " {c |} " %6.4f `mh'[1,7] " {c |}" di as txt "{hline 24}{c +}{hline 12}{c +}{hline 9}{c +}{hline 22}" return scalar stat_mh = `mh'[1,1] //RETURNED RESULTS return scalar se_mh = `mh'[1,2] return scalar stat_mh_lo = `mh'[1,3] return scalar stat_mh_up = `mh'[1,4] return scalar w_mh = `sw'[1,2] } if (`lzeros' & "`zero'"=="n" & ("`stat'"=="or" | "`stat'"=="rr")) { if ("`type'"=="raw" & "`stat'"=="or") di as txt "{hline 24}{c BT}{hline 12}{c BT}{hline 9}{c BT}{hline 22}" display "" display as txt "{bf:NOTE:} M-H combined, Fixed and Random Effects models can't be computed" _n "with studies with zero events." } else { if (`iv'[1,15]>0) di as txt "{lalign 24:FIXED EFFECTS MODEL:}{c |}" _col(38) "{c |}" _col(48) "{c |}" if (`iv'[1,15]==0) di as txt "{lalign 24:FIXED & RANDOM EFFECTS:}{c |}" _col(38) "{c |}" _col(48) "{c |}" //If Tau2=0, there are no random effects local c = cond("`stat'"=="b" & "`exp'"!="","Exp(b)","`title'") di as txt "{ralign 23:IoV Weighted `c'} {c |} " as res /* */ %9.0g cond("`stat'"=="b" & "`exp'"!="",exp(`iv'[1,1]),`iv'[1,1]) " {c |} " %6.4f `iv'[1,5] " {c |} "/* */ %9.0g cond("`stat'"=="b" & "`exp'"!="",exp(`iv'[1,3]),`iv'[1,3]) " " /* */ %9.0g cond("`stat'"=="b" & "`exp'"!="",exp(`iv'[1,4]),`iv'[1,4]) if (`iv'[1,1] < 1 & ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir")) { di as txt "{ralign 23:1/`title'} {c |} " as res %9.0g 1/`iv'[1,1] " {c |} " /* */ _col(48) "{c |} " %9.0g 1/`iv'[1,4] " " %9.0g 1/`iv'[1,3] } local c "SE" if ("`stat'"=="r") local c "SE(Zr)" if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir") local c "SE(ln`title')" di as txt "{ralign 23: `c'} {c |} " as res %9.0g `iv'[1,2] " {c |} " _col(48) "{c |}" di as txt "Heterogeneity Measures: {c |}" _col(38) "{c |}" _col(48) "{c |}" di as txt " Relative excess H {c |} " as res %9.0g `iv'[1,8] " {c |}" _col(48) "{c |} " _c if !(`len'<3 & `len'>=`iv'[1,6]) di as res %9.0g `iv'[1,10] " " %9.0g `iv'[1,11] else di " " if !(`len'<3 & `len'>=`iv'[1,6]) di as txt "{ralign 24:SE(lnH)}{c |} " as res %9.0g `iv'[1,9] " {c |}" _col(48) "{c |}" di as txt " Pct. of variation I^2 {c |} " as res %9.0g `iv'[1,12] " {c |}" _col(48) "{c |} " _c if (`iv'[1,12]>0 & !(`len'<3 & `len'>=`iv'[1,6])) di as res %9.0g `iv'[1,13] " " %9.0g `iv'[1,14] else di "" di as txt " due to heterogeneity {c |}" _col(38) "{c |}" _col(48) "{c |}" di as txt " Homogeneity Chi-Square {c |} " as res %9.0g `iv'[1,6] " {c |} " %6.4f `iv'[1,7] " {c |}" di as txt "VAR between stud. Tau^2 {c |}" as res %11.0g `iv'[1,15] " {c |}" _col(48) "{c |}" return scalar stat_fem = `iv'[1,1] //RETURNED RESULTS return scalar se_fem = `iv'[1,2] return scalar stat_fem_lo = `iv'[1,3] return scalar stat_fem_up = `iv'[1,4] return scalar tau2 = `iv'[1,15] return scalar w_fem = `sw'[1,1] ***Print IOV model - Random effects results (only if Tau2>0) if (`iv'[1,15]>0) { di as txt "{hline 24}{c +}{hline 12}{c +}{hline 9}{c +}{hline 22}" di as txt "{lalign 24:RANDOM EFFECTS MODEL:}{c |}" _col(38) "{c |}" _col(48) "{c |}" local c = cond("`stat'"=="b" & "`exp'"!="","Exp(b)","`title'") di as txt "{ralign 23:IoV Weighted `c'} {c |} " as res /* */ %9.0g cond("`stat'"=="b" & "`exp'"!="",exp(`ivr'[1,1]),`ivr'[1,1]) " {c |} " %6.4f `ivr'[1,5] " {c |} "/* */ %9.0g cond("`stat'"=="b" & "`exp'"!="",exp(`ivr'[1,3]),`ivr'[1,3]) " " /* */ %9.0g cond("`stat'"=="b" & "`exp'"!="",exp(`ivr'[1,4]),`ivr'[1,4]) if (`ivr'[1,1] < 1 & ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir")) { di as txt "{ralign 23:1/`title'} {c |} " as res %9.0g 1/`ivr'[1,1] " {c |} " /* */ _col(48) "{c |} " %9.0g 1/`ivr'[1,4] " " %9.0g 1/`ivr'[1,3] } local c "SE" if ("`stat'"=="r") local c "SE(Zr)" if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir") local c "SE(ln`title')" di as txt "{ralign 23: `c'} {c |} " as res %9.0g `ivr'[1,2] " {c |} " _col(48) "{c |}" di as txt "Heterogeneity Measures: {c |}" _col(38) "{c |}" _col(48) "{c |}" di as txt " Relative excess H {c |} " as res %9.0g `ivr'[1,8] " {c |}" _col(48) "{c |} " _c if !(`len'<3 & `len'>=`ivr'[1,6]) di as res %9.0g `ivr'[1,10] " " %9.0g `ivr'[1,11] else di " " if !(`len'<3 & `len'>=`ivr'[1,6]) di as txt "{ralign 24:SE(lnH)}{c |} " as res %9.0g `ivr'[1,9] " {c |}" _col(48) "{c |}" di as txt " Pct. of variation I^2 {c |} " as res %9.0g `ivr'[1,12] " {c |}" _col(48) "{c |} " _c if (`ivr'[1,12]>0 & !(`len'<3 & `len'>=`ivr'[1,6])) di as res %9.0g `ivr'[1,13] " " %9.0g `ivr'[1,14] else di "" di as txt " due to heterogeneity {c |}" _col(38) "{c |}" _col(48) "{c |}" di as txt " Homogeneity Chi-Square {c |} " as res %9.0g `ivr'[1,6] " {c |} " %6.4f `ivr'[1,7] " {c |}" return scalar stat_rem = `ivr'[1,1] //RETURNED RESULTS return scalar se_rem = `ivr'[1,2] return scalar stat_rem_lo = `ivr'[1,3] return scalar stat_rem_up = `ivr'[1,4] return scalar w_rem = `sw'[1,4] } di as txt "{hline 24}{c BT}{hline 12}{c BT}{hline 9}{c BT}{hline 22}" if ("`type'"=="raw" & "`stat'"=="or" & `lzeros' & "`zero'"!="n") { di as txt "(*) OR of PETO COMBINED does not need continuity correction," di as txt " but this estimation is computed with corrected data" } } ***INFLUENCE ANALYSIS if ("`influ'"!="" & !(`lzeros' & "`zero'"=="n" & ("`stat'"=="or" | "`stat'"=="rr") & ("`method'"=="fem" | "`method'"=="rem"))) { tempname m if ("`method'"=="peto" | "`method'"=="mh") local i 44 if ("`method'"=="fem" | "`method'"=="rem") local i 59 if ("`method'"=="peto") { matrix `m' = `peto' local c1 "PETO" local c2 "PETO" } if ("`method'"=="mh") { matrix `m' = `mh' local c1 "MH" local c2 "MH" } if ("`method'"=="fem") { matrix `m' = `iv' local c1 "IoV Fixed Effect Model" local c2 "IoV FEM" } if ("`method'"=="rem") { matrix `m' = `ivr' local c1 "IoV Random Effect Model" local c2 "IoV REM" } if ("`type'"=="raw" & "`stat'"=="or") local title = "OR" di di as txt "INFLUENCE ANALYSIS" if !("`stat'"=="b" & "`exp'"!="") { di as txt "{hline 12}{c TT}{hline 44}" _c if ("`method'"=="fem" | "`method'"=="rem") di as txt "{hline 15}" _c di as txt "{c TT}{hline 9}" di as txt _col(13) "{c |}{center `i':`c1' estimation if the study is deleted}{c |} Relative" di as txt "{lalign 12:STUDY}{c |}" _col(23) %4.0g `level' "% Confidence Interval" _col(50) "{center 8:`title'w}" _c if ("`method'"=="fem" | "`method'"=="rem") di as txt _col(59) "Heterogeneity " _c di as txt "{c |} Weights" di as txt "{lalign 12:deleted}{c |} {ralign 9:`title'w} {ralign 9:Lower} {ralign 9:Upper} {c |}{center 8:change}" _c if ("`method'"=="fem" | "`method'"=="rem") di as txt "{c |}{ralign 5:I^2} {ralign 6:p} " _c di as txt "{c |}{center 9:`c2'}" di as txt "{hline 12}{c +}{hline 35}{c +}{hline 8}" _c if ("`method'"=="fem" | "`method'"=="rem") di as txt "{c +}{hline 14}" _c di as txt "{c +}{hline 9}" forvalues i = 1/`len'{ di as txt "{lalign 12:`n`i''}{c |} " as res %9.0g `inf'[`i',1] " " %9.0g `inf'[`i',2] as txt " to " /* */ as res %9.0g `inf'[`i',3] " {c |}" _c local p = `inf'[`i',4] print_pct `p', a(6) di as txt " {c |}" _c if ("`method'"=="fem" | "`method'"=="rem") { local p = `inf'[`i',5] print_pct `p' di as txt " " _c di as res %6.4f `inf'[`i',6] " {c |}" _c } local j = cond("`method'"=="peto",3,cond("`method'"=="mh",2,cond("`method'"=="fem",1,4))) local p = `rw'[`i',`j'] di " " _c print_pct `p' di " " } di as txt "{hline 12}{c +}{hline 35}{c +}{hline 8}" _c if ("`method'"=="fem" | "`method'"=="rem") di as txt "{c +}{hline 14}" _c di as txt "{c +}{hline 9}" di as txt "{ralign 11:NONE} {c |} " as res %9.0g `m'[1,1] " " %9.0g `m'[1,3] as txt " to " /* */ as res %9.0g `m'[1,4] " {c |}" %6.0f 0 as txt "% {c |}" _c if ("`method'"=="fem" | "`method'"=="rem") di as res %5.1f `m'[1,12] as txt "% " /* */ as res %6.4f `m'[1,7] " {c |}" _c di as res %6.0f 100 as txt "%" di as txt "{hline 12}{c BT}{hline 35}{c BT}{hline 8}" _c if ("`method'"=="fem" | "`method'"=="rem") di as txt "{c BT}{hline 14}" _c di as txt "{c BT}{hline 9}" } else { di as txt "{hline 12}{c TT}{hline 78}{c TT}{hline 9}" di as txt _col(13) "{c |}{center 78:`c1' Model estimation if the study is deleted}{c |} Relative" di as txt "{lalign 12:STUDY}{c |}" _col(34) %4.0g `level' "% Confidence Interval" _col(61) "{ralign 5:bw}" /* */ _col(70) "Exp(bw) Heterogeneity {c |} Weights" di as txt "{lalign 12:deleted}{c |} {ralign 9:bw} {ralign 9:Exp(bw)} {ralign 9:Lower} {ralign 9:Upper} {c |}"/* */ "{ralign 7:change}{ralign 8:change} {c |}{ralign 5:I^2} {ralign 6:p} {c |}{center 9:`c2'}" di as txt "{hline 12}{c +}{hline 46}{c +}{hline 16}{c +}{hline 14}{c +}{hline 9}" forvalues i = 1/`len'{ di as txt "{lalign 12:`n`i''}{c |} " as res %9.0g `inf'[`i',1] " " %9.0g exp(`inf'[`i',1]) /* */ " " %9.0g exp(`inf'[`i',2]) as txt " to " as res %9.0g exp(`inf'[`i',3]) " {c |} " _c local p = `inf'[`i',4] print_pct `p' di " " _c local p = 100*(exp(`inf'[`i',1]) - exp(`m'[1,1]))/(exp(`m'[1,1])) print_pct `p' di " {c |}" _c local p = `inf'[`i',5] print_pct `p' di " " as res %6.4f `inf'[`i',6] " {c |} " _c local j = cond("`method'"=="fem",1,4) local p = `rw'[`i',`j'] print_pct `p' di " " } di as txt "{hline 12}{c +}{hline 46}{c +}{hline 16}{c +}{hline 14}{c +}{hline 9}" di as txt "{ralign 11:NONE} {c |} " as res %9.0g `m'[1,1] " " %9.0g exp(`m'[1,1]) /* */ " " %9.0g exp(`m'[1,3]) as txt " to " as res %9.0g exp(`m'[1,4]) " {c |}" /* */ %6.0f 0 as txt "% " as res %6.0f 0 as txt "% {c |}" /* */ as res %5.1f `m'[1,12] as txt "% " as res %6.4f `m'[1,7] " {c |}" as res %6.0f 100 as txt "%" di as txt "{hline 12}{c BT}{hline 46}{c BT}{hline 16}{c BT}{hline 14}{c BT}{hline 9}" } } ***CUMULATIVE ANALYSIS if ("`cum'"!="" & !(`lzeros' & "`zero'"=="n" & ("`stat'"=="or" | "`stat'"=="rr") & ("`method'"=="fem" | "`method'"=="rem"))) { quietly { tempvar cs s if ("`cum'"=="date") g __cums = dates if ("`cum'"!="date") g __cums = string(__cumsort,"%6.0g") g `s' = _n if (`ldesc') gsort -__cumsort else gsort +__cumsort forvalues i = 1/`len'{ local nc`i' = abbrev(study[`i'],12) } } if ("`method'"=="fem") local c "IoV estimation: Fixed effects model" if ("`method'"=="rem") local c "IoV estimation: Random effects model" if ("`method'"=="mh") local c "Mantel-Haenszel estimation" if ("`method'"=="peto") local c "Peto estimation" di di as txt "{hline 36}{c TT}{hline 9} CUMULATIVE META-ANALYSIS {hline 9}" di as txt _col(37) "{c |}{center 44:`c'}" di as txt _col(37) "{c |}Number of{c |}" _col(59) "{c |} " %4.0g `level' "% Conf. Interval" di as txt "{lalign 12:STUDY} " cond("`cum'"!="date",upper("`cum'"),"DATE") _col(25) /* */ "{c |}{rcenter 11:`title'}{c |} studies {c |}{rcenter 11:`title'}{c |} {ralign 9:Lower} {ralign 9:Upper}" di as txt "{hline 24}{c +}{hline 11}{c +}{hline 9}{c +}{hline 11}{c +}{hline 22}" forvalues i = 1/`len'{ forvalues j = 1/4 { if ("`stat'"=="b" & "`exp'"!="") matrix `cu'[`i',`j'] = exp(`cu'[`i',`j']) } di as txt "{lalign 12:`nc`i''} " __cums[`i'] _col(25) "{c |} " as res %9.0g `cu'[`i',4] /* */ " {c |} " %7.0f `i' " {c |} " %9.0g `cu'[`i',1] " {c |} " %9.0g `cu'[`i',2] " " %9.0g `cu'[`i',3] } di as txt "{hline 24}{c BT}{hline 11}{c BT}{hline 9}{c BT}{hline 11}{c BT}{hline 22}" qui sort `s' } *Influence / Cumulative don't shown note if (("`influ'"!="" | "`cum'"!="") & (`lzeros' & "`zero'"=="n" & ("`stat'"=="or" | "`stat'"=="rr") & ("`method'"=="fem" | "`method'"=="rem"))) { di as txt _col(7) cond("`influ'"!="" & "`cum'"!="","Influence and Cumulative Meta-",cond("`influ'"!="","Influence ","Cumulative Meta-")) _c di as txt "Analysis can't be computed with studies with zero events." di as txt _col(7) "To obtain these results, you may select only studies without zero events" di as txt _col(7) "or apply a continuity correction to include all studies (see zero option)." } *Publication Bias if (("`begg'"!="" | "`egger'"!="" | "`macask'"!="" | "`nus'"!="")) { local estitle = upper("`stat'") quietly { tempname esi sei pi gen `esi' = . gen `sei' = . gen `pi' = . forvalues i = 1/`len'{ replace `esi' = __mar_st[`i',1] in `i' replace `sei' = __mar_st[`i',2] in `i' replace `pi' = __mar_st[`i',5] in `i' } if ("`stat'"=="r") replace `esi' = ln(sqrt((1+`esi')/(1-`esi'))) } if ("`begg'"!="") { tempname rho_vt rho_n n_vt n_n p_vt p_n vi ti es_begg se_begg n_begg if ("`method'"=="peto") { scalar `es_begg' = `peto'[1,1] scalar `se_begg' = `peto'[1,2] } if ("`method'"=="mh") { scalar `es_begg' = `mh'[1,1] scalar `se_begg' = `mh'[1,2] } if ("`method'"=="fem") { scalar `es_begg' = `iv'[1,1] scalar `se_begg' = `iv'[1,2] } if ("`method'"=="rem") { if (`iv'[1,15]>0) { scalar `es_begg' = `ivr'[1,1] scalar `se_begg' = `ivr'[1,2] } else { scalar `es_begg' = `iv'[1,1] scalar `se_begg' = `iv'[1,2] } } if ("`stat'"=="r") scalar `es_begg' = ln(sqrt((1+`es_begg')/(1-`es_begg'))) qui gen `vi' = `sei'^2 if (`raw_md') { qui gen `n_begg' = n1+n0 } else { qui gen `n_begg' = n } if ("`stat'"=="rd" | "`stat'"=="id" | "`stat'"=="md" | "`stat'"=="r") { qui gen `ti' = (`esi' - `es_begg')/sqrt(`vi' - `se_begg'^2) } else { qui gen `ti' = (ln(`esi') - ln(`es_begg'))/sqrt(`vi' - `se_begg'^2) } di di as res "Begg Publication Bias" qui ktau `ti' `vi' scalar `rho_vt' = r(tau_b) scalar `p_vt' = r(p) scalar `n_vt' = r(N) qui ktau `ti' `n_begg' scalar `rho_n' = r(tau_b) scalar `p_n' = r(p) scalar `n_n' = r(N) di as txt _col(53) "{c |} Vi {c |} N" di as txt "{hline 52}{c +}{hline 9}{c +}{hline 9}" di as txt "Kendall's tau_b Ti Correlation Coefficient {c |} " as res %7.4f `rho_vt' " {c |} " %7.4f `rho_n' di as txt _col(23) "Prob > |z| (continuity corr.) {c |} " as res %7.4f `p_vt' " {c |} " %7.4f `p_n' di as txt _col(51) "N {c |} " as res %7.0f `n_vt' " {c |} " %7.0f `n_n' di as txt "{hline 52}{c BT}{hline 9}{c BT}{hline 9}" } if ("`egger'"!="") { di di as res "Egger Publication Bias" qui gen Precision = 1 / `sei' if ("`stat'"=="rd" | "`stat'"=="id" | "`stat'"=="md" | "`stat'"=="r") { if ("`stat'"=="r") { qui gen Zr_SE = `esi' / `sei' regress Zr_SE Precision, noheader } else { qui gen `estitle'_SE = `esi' / `sei' regress `estitle'_SE Precision, noheader } } else { qui gen ln`estitle'_SE = ln(`esi') / `sei' regress ln`estitle'_SE Precision, noheader } } if ("`macask'"!="" & "`stat'"!="rd" & "`stat'"!="id" & "`stat'"!="md" & "`stat'"!="r") { *Macask weight variable is needed if ("`type'"!="raw") { checkvars `macask', errmsg checkvars n, errmsg } else { if ("`macask'"!="pres" & "`macask'"!="pres_ue" & "`macask'"!="pexp" & "`macask'"=="pexp_nc") { checkvars `macask', errmsg } } qui gen w = n*`macask'*(1-`macask') qui gen inv_n = 1 / n qui gen ln`estitle' = ln(`esi') di di as res "Macaskill Publication Bias" di as txt "regress ln`estitle' n [aw=w], noheader" regress ln`estitle' n [aw=w], noheader di regress ln`estitle' inv_n [aw=w], noheader } if ("`nus'"!="") { if ("`type'"=="raw") { qui gen p_nus = `pi' if ("`method'"=="peto") { matname __mar_stp es_p se_p lo_p up_p p_peto chi2_p, columns(.) explicit svmat __mar_stp, names(col) qui gen p_nus = p_peto } if ("`stat'"=="or" | "`stat'"=="rr") qui gen z_nus = sqrt((n-1)*((a1*n-n1*m1)^2)/(m1*m0*n1*n0)) if ("`stat'"=="ir") qui gen z_nus = sqrt(((a1*t-t1*m1)^2)/(m1*t1*t0)) } if ("`type'"=="sst") { if ("`stat'"=="rd" | "`stat'"=="id" | "`stat'"=="md" | "`stat'"=="r") { qui gen p_nus = 2*(1-normal(abs(`esi'/`sei'))) qui gen z_nus = `esi' / `sei' } else { qui gen p_nus = 2*(1-normal(abs(ln(`esi')/`sei'))) qui gen z_nus = ln(`esi') / `sei' } } quietly { gen __sort_nus = _n sort p_nus gen _rp = _n gen _flt = p_nus < `nus' if p_nus < . egen _mnus = mean(_flt) egen _pmax = max(p_nus) egen _pmin = min(p_nus) egen _sumz = total(z_nus) replace _mnus = _mnus*_N gen _f = ceil((_sumz^2)/2.706 - _N) gen _pcut0 = . replace _pcut0 = p_nus if (_rp == _mnus) egen _pcut = max(_pcut0) } di di as res "Publication Bias" *Print Rosenthal method di di as txt "RESISTANCE FOR FUTURE NULL RESULTS (ROSENTHAL model)" di as txt "Number of additional non-significant studies that is necessary" di as txt "to add a significant meta-analysis for becoming non-significant ......." as res %7.0f _f[1] local num = _f[1] local criteria = _N*5+10 if (_f[1]<=`criteria') { di as txt "WARNING! The number of unpublished studies `num' is <= `criteria' (k*5+10)" } else { di as txt "The number of unpublished studies `num' is > `criteria' (k*5+10)" } di di as txt "GLESER-OLKIN Models" di as txt "{hline 12}{c TT}{hline 4}{c TT}{hline 8}" di as txt "Study {c |} RP {c |} p" di as txt "{hline 12}{c +}{hline 4}{c +}{hline 8}" forvalues i = 1/`len'{ di as txt %-11s study[`i'] " {c |} " as res %2.0f _rp[`i'] " {c |} " %5.4f p_nus[`i'] } di as txt "{hline 12}{c BT}{hline 4}{c BT}{hline 8}" di as txt "Number of reported studies..............: " as res %2.0f _N di as txt "Largest reported p-value................: " as res %5.4f _pmax[1] di as txt "Smallest reported p-value...............: " as res %5.4f _pmin[1] di as txt "P cut for select smallest p-values......: " as res %3.2f `nus' di as txt "Number of studies with smallest p-values: " as res %2.0f _mnus[1] di as txt "Maximum of the smallest p-values........: " as res %5.4f _pcut[1] tempname knus pmax mnus pcut ngo ngog q qg f idf nl nlg scalar `knus' = _N scalar `pmax' = _pmax[1] scalar `mnus' = _mnus[1] scalar `pcut' = _pcut[1] if (`mnus' > 0) { *Gleser-Olkin model scalar `ngo' = trunc((`knus'-1)/`pmax' - `knus') scalar `ngog' = trunc((`mnus'-1)/`pcut' - `knus') *Gleser-Olkin model: lower bound of the CI scalar `q' = -999 scalar `qg' = -999 forvalues i = 0/999 { *Simple model scalar `f' = (`i'+1)*`pmax'/(`knus'*(1-`pmax')) scalar `idf' = invF(2*`knus',2*(`i'+1),(100-`level')/100) if (`q'<0 & `idf'<`f') scalar `q' = `i' *Generalized simple model scalar `f' = (`i'+1)*`pcut'/(`mnus'*(1-`pcut')) scalar `idf' = invF(2*`mnus',2*(`i'+1),(100-`level')/100) if (`qg'<0 & `idf'<`f') scalar `qg' = `i' if (`qg'>0 & `q'>0) continue, break } scalar `nl' = `q' scalar `nlg' = `qg' - `knus' + `mnus' di as txt "{hline 21}{c TT}{hline 29}" di as txt _col(22) "{c |}NUMBER OF UNPUBLISHED STUDIES" di as txt _col(22) "{c |} {c |} Generalized" di as txt _col(22) "{c |} Simple model {c |} Simple model" di as txt "{hline 21}{c +}{hline 14}{c +}{hline 14}" di as txt " Unbiased Estimation {c |}" as res %9.0f `ngo' _col(37) "{c |}" %9.0f `ngog' if (`ngog'>=0) di as txt " `level'% CI Lower bound {c |}" as res %9.0f `nl' _col(37) "{c |}" %9.0f `nlg' else di as txt " `level'% CI Lower bound {c |}" as res %9.0f `nl' _col(37) "{c |} unestimable" di as txt "{hline 21}{c BT}{hline 14}{c BT}{hline 14}" } else { di as txt "NUMBER OF UNPUBLISHED STUDIES = Not estimable" } } } tempfile dirsave if ("`keep'"!="") qui save `dirsave', replace ***GRAPHICS if ("`interval'"!="" | "`forest'"!="" | "`histogram'"!="" | "`radial'"!="" | "`cumulative'"!="" | /* */ "`abbe'"!="" | "`funnel'"!="" | "`funnelh'"!="" | "`galbraith'"!="") { local int_var "" if ("`interval'"!="" & "`interval'"!="se") local int_var "`interval'" local for_var "" if ("`forest'"!="" & "`forest'"!="se") local for_var "`forest'" local cum_var "" if ("`cumulative'"!="") local cum_var "__cums __cumsort" local abbe_vars "" if ("`abbe'"!="") local abbe_vars "a0 n0 a1 n1" keep study `cum_var' `abbe_vars' `int_var' `for_var' quietly { *Build the dataset for the graphs gen `stat' = . gen se = . gen lo = . gen up = . gen rw = . gen _comb = 0 forvalues i = 1/`len'{ replace `stat' = __mar_st[`i',1] in `i' replace se = __mar_st[`i',2] in `i' replace lo = __mar_st[`i',3] in `i' replace up = __mar_st[`i',4] in `i' if ("`method'"=="fem") replace rw = __mar_rw[`i',1] in `i' if ("`method'"=="mh") replace rw = __mar_rw[`i',2] in `i' if ("`method'"=="peto") replace rw = __mar_rw[`i',3] in `i' if ("`method'"=="rem") replace rw = __mar_rw[`i',4] in `i' } if ("`method'"=="fem") gen sw = __mar_sw[1,1] if ("`method'"=="mh") gen sw = __mar_sw[1,2] if ("`method'"=="peto") gen sw = __mar_sw[1,3] if ("`method'"=="rem") gen sw = __mar_sw[1,4] gen wght = round(rw*sw/100) *Add the COMBINE results (estimated estatistic) if ("`method'"=="fem") matrix __mar_temp = __mar_iv if ("`method'"=="mh") matrix __mar_temp = __mar_mh if ("`method'"=="peto") matrix __mar_temp = __mar_peto if ("`method'"=="rem") matrix __mar_temp = __mar_ivr set obs `len_1' replace study = "Overall" in `len_1' replace `stat' = __mar_temp[1,1] in `len_1' replace se = __mar_temp[1,2] in `len_1' replace lo = __mar_temp[1,3] in `len_1' replace up = __mar_temp[1,4] in `len_1' replace _comb = 1 in `len_1' matrix drop __mar_temp if ("`method'"=="fem") local ref_es = __mar_iv[1,1] if ("`method'"=="mh") local ref_es = __mar_mh[1,1] if ("`method'"=="peto") local ref_es = __mar_peto[1,1] if ("`method'"=="rem") local ref_es = __mar_ivr[1,1] *Correlation graphics: use Fisher's Z wiht SE if ("`stat'"=="r") { replace `stat' = ln(sqrt((1+`stat')/(1-`stat'))) replace lo = ln(sqrt((1+lo)/(1-lo))) replace up = ln(sqrt((1+up)/(1-up))) local ref_es = ln(sqrt((1+`ref_es')/(1-`ref_es'))) } } if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir" | ("`stat'"=="b" & "`exp'"!="")) { if ("`stat'"=="b" & "`exp'"!="") local xtitle = "Exp(b)" else local xtitle = upper("`stat'") local xtitle = "`xtitle' (logarithmic scale)" local ref_line = 1 local xscale = "xscale(log)" } else { if ("`stat'"=="md") { if (`md_std'==0) local xtitle "Mean difference" if (`md_std'==1) local xtitle "Cohen's d" if (`md_std'==2) local xtitle "Hedges' g" if (`md_std'==3) local xtitle "Glass's Delta" } else if ("`stat'"=="r") { local xtitle "Zr" } else { local xtitle = upper("`stat'") } local ref_line = 0 local xscale = "" } *Create graphics if ("`interval'"!="") { quietly { gen __sort_g = _n gsort -_comb +`interval' gen id_int = _n forvalues i=1/`len_1' { local val = study[`i'] label define d_int `i' "`val'", add } label values id_int d_int if ("`stat'"=="b" & "`exp'"!="") { replace `stat' = exp(`stat') replace lo = exp(lo) replace up = exp(up) local ref_es = exp(`ref_es') } /*gen xpos_st_i = minlo - 0.1 (scatter id xpos_st_i, msymbol(none) mlabel(study) mlabcolor(black) mlabpos(3) mlabsize(small)), */ noisily twoway (rcap lo up id_int, horizontal lwidth(medthick) lcolor(black)) /// (scatter id_int `stat', mfcolor(gs10) msize(medlarge) mcolor(black)), /// xline(`ref_line', lcolor(black)) xline(`ref_es', lcolor(black) lpattern(dash)) /// xtitle(`xtitle', size(medlarge)) ylabel(#`len_1', nogrid valuelabels angle(horizontal) labsize(medsmall)) /// legend(off) ytitle("") `xscale' name(interval,replace) title(Interval Plot) if ("`stat'"=="b" & "`exp'"!="") { replace `stat' = ln(`stat') replace lo = ln(lo) replace up = ln(up) local ref_es = ln(`ref_es') } sort __sort_g drop __sort_g } } if ("`forest'"!="") { quietly { gen __sort_g = _n gsort -_comb +`forest' gen id_for = _n gen ytop = id_for + 0.4 if _comb==1 gen ydown = id_for - 0.4 if _comb==1 forvalues i=1/`len_1' { local val = study[`i'] label define d_for `i' "`val'", add } label values id_for d_for if ("`stat'"=="b" & "`exp'"!="") { replace `stat' = exp(`stat') replace lo = exp(lo) replace up = exp(up) local ref_es = exp(`ref_es') } /*gen xpos_st = minlo - 0.1 (scatter id xpos_st, msymbol(none) mlabel(study) mlabcolor(black) mlabpos(3) mlabsize(small))*/ twoway (scatter id_for `stat' if _comb==0 [aweight=rw], msymbol(square) mfcolor("180 180 180") mlcolor(black)) /// (pcspike id_for lo id_for up if _comb==0, lcolor(black) lwidth(medthick)) /// (pcspike id_for lo ytop `stat' if _comb==1, lcolor(black) lwidth(medthick)) /// (pcspike ytop `stat' id_for up if _comb==1, lcolor(black) lwidth(medthick)) /// (pcspike id_for up ydown `stat' if _comb==1, lcolor(black) lwidth(medthick)) /// (pcspike ydown `stat' id_for lo if _comb==1, lcolor(black) lwidth(medthick)) /// (scatter id_for `stat' if _comb==0, msymbol(diamond) msize(vsmall) mcolor(black)), /// `xscale' legend(off) ylabel(#`len_1', nogrid valuelabels angle(horizontal) labsize(medsmall)) /// xtitle(`xtitle') xline(`ref_line', lcolor(black)) xline(`ref_es', lcolor(black) lpattern(dash)) /// ytitle("") title(Forest Plot) name(forest,replace) if ("`stat'"=="b" & "`exp'"!="") { replace `stat' = ln(`stat') replace lo = ln(lo) replace up = ln(up) local ref_es = ln(`ref_es') } sort __sort_g drop __sort_g } } if ("`radial'"!="") { if ("`stat'"=="rd" | "`stat'"=="id" | "`stat'"=="md" | "`stat'"=="r") { local xtitler = "`xtitle' / SE" qui gen x_rad = `stat'/se } else { local xtitler = "ln" + upper("`stat'") + " / SE" qui gen x_rad = ln(`stat')/se } qui gen y_rad = 1/se twoway (scatter y_rad x_rad), xtitle(`xtitler', size(medlarge)) /// ytitle(1 / SE, size(medlarge)) title(Radial Plot) name(radial,replace) } if ("`histogram'"!="") { qui drop if _comb==1 /*local xtitleh = upper("`stat'") if ("`stat'"=="md") { if (`md_std'==1) local xtitleh "Cohen's d" if (`md_std'==2) local xtitleh "Hedges' g" if (`md_std'==3) local xtitleh "Glass's Delta" }*/ qui histogram `stat' [w=wght], frequency xtitle(`xtitle', size(medlarge)) /// ytitle("Study weight sum", size(medlarge)) ylabel(,angle(horizontal)) /// title("Histogram (weighted)") name(histogram,replace) } if ("`cumulative'"!="") { quietly { gen __sort_g = _n drop if _comb==1 if (`ldesc') gsort -__cumsort else gsort +__cumsort gen id_cum = _n gen es_cum = . gen lo_cum = . gen up_cum = . forvalues i = 1/`len'{ replace es_cum = __mar_cum[`i',1] in `i' replace lo_cum = __mar_cum[`i',2] in `i' replace up_cum = __mar_cum[`i',3] in `i' } gen cum_title = study + " (" + __cums + ")" forvalues i=1/`len' { local val = cum_title[`i'] label define d_cum `i' "`val'", add } label values id_cum d_cum } if ("`stat'"=="b" & "`exp'"!="") { qui replace es_cum = exp(es_cum) qui replace lo_cum = exp(lo_cum) qui replace up_cum = exp(up_cum) } /*(scatter id xpos_title, msymbol(none) mlabel(cum_title) mlabcolor(black) mlabpos(3) mlabsize(small))*/ twoway (scatter id_cum es_cum, msymbol(circle) msize(small) mcolor(black)) /// (pcspike id_cum lo_cum id_cum up_cum, lcolor(navy)), /// `xscale' legend(off) yscale(reverse) ylabel(#`len', nogrid valuelabels angle(horizontal) labsize(small)) ytitle("") /// xtitle(`xtitle') xline(`ref_line', lcolor(black)) name(cumulative,replace) legend(off) title(Cumulative Plot) sort __sort_g drop __sort_g } if ("`abbe'"!="") { quietly { drop if _comb==1 gen Px = (a0/n0)*100 gen Py = (a1/n1)*100 set obs `len_1' replace Px = 0 in `len_1' } qui twoway (line Px Px) (scatter Py Px [aweight=rw], msymbol(square) mlcolor(black) mfcolor(none)) , /// legend(off) xtitle(Px) ytitle(Py) name(abbe,replace) title(Abbe Plot) } if ("`galbraith'"!="") { quietly { drop if _comb==1 if ("`stat'"=="rd" | "`stat'"=="id" | "`stat'"=="md" | "`stat'"=="r") { gen y_galb = `stat'/se local ln_se = `ref_es' local ytitle = "`xtitle' / SE" } else { gen y_galb = ln(`stat')/se local ln_se = ln(`ref_es') local ytitle = "ln" + upper("`stat'") + " / SE" } gen x_galb = 1/se egen x_galb_max = max(x_galb) } local max_galb = x_galb_max[1] twoway (scatter y_galb x_galb, msize(medium) mlabel(study) mlabposition(9)) /// (function y=`ln_se'*x, range(0 `max_galb') lcolor(black) lpattern(solid)) /// (function y=`ln_se'*x+2, range(0 `max_galb') lcolor(black) lpattern(dash)) /// (function y=`ln_se'*x-2, range(0 `max_galb') lcolor(black) lpattern(dash)), /// legend(off) xtitle(1 / SE) ytitle(`ytitle') name(galbraith,replace) title(Galbraith (radial) Plot) } if (("`funnel'"!="" | "`funnelh'"!="") & ("`stat'"!="rd" & "`stat'"!="id" & "`stat'"!="md")) { quietly { drop if _comb==1 if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir") { gen fun_es = ln(`stat') local fun_ref = ln(`ref_es') local xtitle = "ln" + upper("`stat'") local ytitle = "SE(ln" + upper("`stat')") } else if ("`stat'"=="rd" | "`stat'"=="id" | "`stat'"=="md" | "`stat'"=="b" | "`stat'"=="r") { gen fun_es = `stat' local fun_ref = `ref_es' local xtitle = upper("`stat'") local ytitle = upper("SE(`stat')") if ("`stat'"=="r") { local xtitle "Zr" local ytitle "SE(Zr)" } } gen fun_se = se forval i = 1/6 { local obs`i' = _N + `i' } set obs `obs6' egen maxse = max(se) if ("`funnel'"!="") { *Estimated statistic line gen y_es = 0 if _n==`obs1' replace y_es = maxse if _n==`obs2' gen x_fun = `fun_ref' if _n==`obs1' | _n==`obs2' *Pseudo 95% confidence limits - lo gen y_lo = 0 if _n==`obs3' replace y_lo = maxse if _n==`obs4' replace x_fun = `fun_ref' if _n==`obs3' replace x_fun = `fun_ref' - invnorm(0.975) * maxse if _n==`obs4' *Pseudo 95% confidence limits - up gen y_up = 0 if _n==`obs5' replace y_up = maxse if _n==`obs6' replace x_fun = `fun_ref' if _n==`obs5' replace x_fun = `fun_ref' + invnorm(0.975) * maxse if _n==`obs6' } if ("`funnelh'"!="") { *Estimated statistic line gen y_esh = `fun_ref' if _n==`obs1' | _n==`obs2' gen x_funh = 0 if _n==`obs1' replace x_funh = maxse if _n==`obs2' *Pseudo 95% confidence limits - lo gen y_loh = `fun_ref' if _n==`obs3' replace y_loh = `fun_ref' - invnorm(0.975) * maxse if _n==`obs4' replace x_funh = 0 if _n==`obs3' replace x_funh = maxse if _n==`obs4' *Pseudo 95% confidence limits - up gen y_uph = `fun_ref' if _n==`obs5' replace y_uph = `fun_ref' + invnorm(0.975) * maxse if _n==`obs6' replace x_funh = 0 if _n==`obs5' replace x_funh = maxse if _n==`obs6' } } if ("`funnel'"!="") { twoway (scatter fun_se fun_es, yscale(reverse)) /// (line y_lo y_up y_es x_fun, msymbol(none none none) /// clcolor(black black black) clpat(dash dash solid) /// clwidth(medium medium medium)), /// legend(off) xtitle(`xtitle') ytitle(`ytitle') /// title(Funnel plot with pseudo 95% confidence limits) name(funnel,replace) } if ("`funnelh'"!="") { twoway (scatter fun_es fun_se) /// (line y_loh y_uph y_esh x_funh, msymbol(none none none) /// clcolor(black black black) clpat(dash dash solid) /// clwidth(medium medium medium)), /// legend(off) ytitle(`xtitle') xtitle(`ytitle') /// title(Funnel plot with pseudo 95% confidence limits) name(funnelh,replace) } } } *Save variables to keep if ("`keep'"!="") { qui use `dirsave', clear local keep_vars study if ("`ldate'"!="") local keep_vars = "`keep_vars' date" if (`raw_md') { local keep_vars = "`keep_vars' m1 sd1 n1 m0 sd0 n0" if (`md_exists') local keep_vars = "`keep_vars' _md0_orig _se0_orig _lo0_orig _up0_orig" } if ("`type'"=="raw") { if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="rd") local keep_vars = "`keep_vars' a1 a0 b1 b0 m1 m0 n1 n0 n pres pres_ue pres_ex pexp pexp_nc pexp_ca" if ("`stat'"=="ir" | "`stat'"=="id") local keep_vars = "`keep_vars' a1 a0 t1 t0 i1 i0 m1 t pres pres_ue pres_ex pexp pexp_ca" erasevars `stat' se lo up quietly { *Add matrix variables gen `stat' = . gen se = . gen lo = . gen up = . forvalues i = 1/`len'{ replace `stat' = __mar_st[`i',1] in `i' replace se = __mar_st[`i',2] in `i' replace lo = __mar_st[`i',3] in `i' replace up = __mar_st[`i',4] in `i' } if ("`stat'"=="or") { erasevars or_peto lnor_peto se_or_peto lo_or_peto up_or_peto wp gen or_peto = . gen lnor_peto = . gen se_or_peto = . gen lo_or_peto = . gen up_or_peto = . gen wp = . forvalues i = 1/`len'{ replace or_peto = __mar_stp[`i',1] in `i' replace se_or_peto = __mar_stp[`i',2] in `i' replace lo_or_peto = __mar_stp[`i',3] in `i' replace up_or_peto = __mar_stp[`i',4] in `i' replace wp = __mar_rw[`i',3] * __mar_sw[1,3]/100 in `i' } replace lnor_peto = ln(or_peto) local keep_vars = "`keep_vars' or_peto lnor_peto se_or_peto lo_or_peto up_or_peto wp" } if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir") { erasevars wmh gen wmh = . forvalues i = 1/`len'{ replace wmh = __mar_rw[`i',2] * __mar_sw[1,2]/100 in `i' } local keep_vars = "`keep_vars' wmh" } } } if ("`type'"=="sst") local keep_vars = "`keep_vars' cl" if ("`stat'"=="r") { qui { erasevars z rename r z rename `rho' r order r, before(z) replace lo = (exp(2*lo)-1)/(exp(2*lo)+1) replace up = (exp(2*up)-1)/(exp(2*up)+1) } } quietly { erasevars wivf wivr gen wivf = . gen wivr = . forvalues i = 1/`len'{ replace wivf = __mar_rw[`i',1] * __mar_sw[1,1]/100 in `i' replace wivr = __mar_rw[`i',4] * __mar_sw[1,4]/100 in `i' } } local keep_vars = "`keep_vars' `stat' se lo up wivf wivr" if ("`stat'"=="r") local keep_vars = "`keep_vars' z" local keep_vars = "`keep_vars' `keep_mix'" keep `keep_vars' if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir") { *Add ln(stat) for or, rr, ir quietly generate ln`stat' = ln(`stat') order ln`stat', after(`stat') } qui save `dirsave', replace } restore //Restore original dataset *Add variables to the original dataset if ("`keep'"!="") { tempname sort_orig qui gen `sort_orig' = _n qui merge 1:1 study using `dirsave', update replace nogenerate sort `sort_orig' capture confirm variable _mar_ if (_rc!=0) qui gen _mar_ = 1 //Variable to identify the dataset as 'keeped' if (`raw_md') { if (`md_std'>0) { if (`md_std'==1) local vars_md "d se_d lo_d up_d" if (`md_std'==2) local vars_md "g se_g lo_g up_g" if (`md_std'==3) local vars_md "Delta se_Delta lo_Delta up_Delta" foreach v in `vars_md' { capture confirm variable `v', exact if (_rc==0) drop `v' } rename (`stat' se lo up) (`vars_md') if (`md_exists') rename (_md0_orig _se0_orig _lo0_orig _up0_orig) (md se lo up) } else { if (`md_exists') drop _md0_orig _se0_orig _lo0_orig _up0_orig } } *Label variables quietly{ label variable study "Study" if ("`ldate'"!="") label variable date "Date of Study" if ("`type'"=="raw") { label variable a1 "Cases Exposed" label variable a0 "Cases Unexposed" label variable pres "Proportion of response in Total sample" label variable pres_ue "Proportion of response in Unexposed" label variable pres_ex "Proportion of response in Exposed" label variable pexp "Proportion of exposure in Total sample" label variable pexp_ca "Proportion of exposure in Cases" if ("`stat'"=="ir" | "`stat'"=="id") { label variable t1 "Person-time Exposed" label variable t0 "Person-time Unexposed" label variable i1 "Incidence Exposed" label variable i0 "Incidence Unexposed" label variable m1 "Total Cases" label variable t "Total Person-time" } else { label variable pexp_nc "Proportion of exposure in Non Cases" label variable b1 "Non Cases Exposed" label variable b0 "Non Cases Unexposed" label variable n1 "Total Exposed" label variable n0 "Total Unexposed" label variable m1 "Total Cases" label variable m0 "Total Non Cases" label variable n "Total" } if ("`stat'"=="or" | "`stat'"=="rr") { *v1.3.9: for raw data & or | rr stat: always keep or, or_peto, rr if ("`stat'"=="or") { erasevars se_or lo_or up_or rr lnrr se_rr lo_rr up_rr rename (se lo up) (se_or lo_or up_or) qui g rr = (a1 / n1) / (a0 / n0) qui g lnrr = ln(rr) qui g se_rr = sqrt(1/a1 - 1/n1 + 1/a0 - 1/n0) qui g lo_rr = rr * exp(-invnormal((`level'+100)/200) * se_rr) qui g up_rr = rr * exp(invnormal((`level'+100)/200) * se_rr) order or lnor se_or lo_or up_or or_peto lnor_peto se_or_peto lo_or_peto up_or_peto rr lnrr se_rr lo_rr up_rr, after(pres_ex) } if ("`stat'"=="rr") { erasevars se_rr lo_rr up_rr or lnor se_or lo_or up_or rename (se lo up) (se_rr lo_rr up_rr) qui g or = (a1 / b1) / (a0 / b0) qui g lnor = ln(or) qui g se_or = sqrt(1/a1 + 1/b1 + 1/a0 + 1/b0) qui g lo_or = or * exp(-invnormal((`level'+100)/200) * se_or) qui g up_or = or * exp(invnormal((`level'+100)/200) * se_or) erasevars or_peto lnor_peto se_or_peto lo_or_peto up_or_peto qui g se_or_peto = (m1 * (n1/n)) * ((n0 * m0) / (n * (n-1))) qui g or_peto = exp((a1 - (m1 * (n1/n))) / se_or_peto) qui g lnor_peto = ln(or_peto) qui replace se_or_peto = sqrt(1/se_or_peto) qui g lo_or_peto = or_peto * exp(-invnormal((`level'+100)/200) * se_or_peto) qui g up_or_peto = or_peto * exp(invnormal((`level'+100)/200) * se_or_peto) order rr lnrr se_rr lo_rr up_rr or lnor se_or lo_or up_or or_peto lnor_peto se_or_peto lo_or_peto up_or_peto, after(pres_ex) } label variable or "Odds Ratio" label variable lnor "ln(Odds Ratio)" label variable se_or "Standard Error of ln(OR)" label variable lo_or "OR Lower Bound (`level'% CI)" label variable up_or "OR Upper Bound (`level'% CI)" label variable or_peto "Odds Ratio (Peto)" label variable lnor_peto "ln(Odds Ratio) (Peto)" label variable se_or_peto "Standard Error of ln(OR) (Peto)" label variable lo_or_peto "OR (Peto) Lower Bound (`level'% CI)" label variable up_or_peto "OR (Peto) Upper Bound (`level'% CI)" label variable rr "Risk Ratio" label variable lnrr "ln(Risk Ratio)" label variable se_rr "Standard Error of ln(RR)" label variable lo_rr "RR Lower Bound (`level'% CI)" label variable up_rr "RR Upper Bound (`level'% CI)" } } if ("`type'"!="raw" & "`stat'"=="or") { label variable or "Odds Ratio" label variable lnor "ln(Odds Ratio)" label variable se "Standard Error of ln(OR)" label variable lo "OR Lower Bound (`level'% CI)" label variable up "OR Upper Bound (`level'% CI)" } if ("`type'"!="raw" & "`stat'"=="rr") { label variable rr "Risk Ratio" label variable lnrr "ln(Risk Ratio)" label variable se "Standard Error of ln(RR)" label variable lo "RR Lower Bound (`level'% CI)" label variable up "RR Upper Bound (`level'% CI)" } if ("`stat'"=="ir") { label variable ir "Rate Ratio" label variable lnir "ln(Rate Ratio)" label variable se "Standard Error of ln(IR)" label variable lo "IR Lower Bound (`level'% CI)" label variable up "IR Upper Bound (`level'% CI)" } if ("`stat'"=="rd") { label variable rd "Risk Difference" label variable se "Standard Error of RD" label variable lo "RD Lower Bound (`level'% CI)" label variable up "RD Upper Bound (`level'% CI)" } if ("`stat'"=="id") { label variable id "Rate Difference" label variable se "Standard Error of ID" label variable lo "ID Lower Bound (`level'% CI)" label variable up "ID Upper Bound (`level'% CI)" } if ("`stat'"=="b") { label variable b "Slope" label variable se "Standard Error of B" label variable lo "B Lower Bound (`level'% CI)" label variable up "B Upper Bound (`level'% CI)" } if ("`stat'"=="md") { if (`raw_md') { if (`md_std'==0) { label variable md "Mean Difference" label variable se "Standard Error of MD" label variable lo "MD Lower Bound (`level'% CI)" label variable up "MD Upper Bound (`level'% CI)" } if (`md_std'==1) { label variable d "Cohen's d standardised mean difference" label variable se_d "Standard Error of d" label variable lo_d "d Lower Bound (`level'% CI)" label variable up_d "d Upper Bound (`level'% CI)" } if (`md_std'==2) { label variable g "Hedges' adjusted g standardised mean difference" label variable se_g "Standard Error of g" label variable lo_g "g Lower Bound (`level'% CI)" label variable up_g "g Upper Bound (`level'% CI)" } if (`md_std'==3) { label variable Delta "Glass's Delta standardised mean difference" label variable se_Delta "Standard Error of Delta" label variable lo_Delta "Delta Lower Bound (`level'% CI)" label variable up_Delta "Delta Upper Bound (`level'% CI)" } } else { if (`md_std'==0) label variable md "Mean Difference" if (`md_std'==1) label variable md "Mean Difference (Cohen's d)" if (`md_std'==2) label variable md "Mean Difference (Hedges' adjusted g)" if (`md_std'==3) label variable md "Mean Difference (Glass's Delta)" label variable se "Standard Error of MD" label variable lo "MD Lower Bound (`level'% CI)" label variable up "MD Upper Bound (`level'% CI)" } } if ("`stat'"=="r") { label variable r "Correlation coefficient" label variable n "Sample Size" label variable z "Fisher's z Transformed Correlation coefficient (Zr)" label variable se "Standard Error of Zr" label variable lo "Correlation Coefficient Lower Bound (`level'% CI)" label variable up "Correlation Coefficient Upper Bound (`level'% CI)" } if ("`type'"=="sst") label variable cl "Confidence Level" if ("`type'"=="raw" & "`stat'"=="or") label variable wp "Peto Weight" if ("`type'"=="raw" & ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir")) label variable wmh "Mantel-Haenszel Weight" label variable wivf "Inverse of Variance Fixed Effects Model Weight" label variable wivr "Inverse of Variance Random Effects Model Weight" } } *Drop temporary matrix matrix drop __mar_data __mar_st __mar_stp __mar_rw __mar_sw __mar_peto __mar_mh __mar_iv __mar_ivr __mar_influ __mar_cum end program data_raw, rclass syntax anything(name=stat), zero(string) vzeros(name) md_std(integer) level(real) *Prepare raw data for analysis local lzeros 0 //default values for zero correction flags local ormh_zero 0 local md_exists 0 quietly { if ("`stat'"!="md") checkvars a1 a0, errmsg if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="rd") { checkvars b1 b0, local(found_b) checkvars n1 n0, local(found_n) if (!`found_b' & !`found_n') print_error "missing b1,b0 or n1,n0 variable(s)" if (`found_b') { erasevars n1 n0 g n1 = a1 + b1 g n0 = a0 + b0 } if (`found_n') { erasevars b1 b0 g b1 = n1 - a1 g b0 = n0 - a0 } erasevars m1 m0 n g m1 = a0 + a1 g m0 = b0 + b1 g n = n1 + n0 zero_correction `zero', vzeros(`vzeros') local lzeros = r(zeros) local ormh_zero = r(ormh) erasevars pexp pexp_nc pexp_ca pres pres_ue pres_ex g pexp = n1/n g pexp_nc = b1/m0 g pexp_ca = a1/m1 //v1.3.9: add pexp_ca, pres_ex variables g pres = m1/n g pres_ue = a0/n0 g pres_ex = a1/n1 //v1.3.9: add pexp_ca, pres_ex variables } if ("`stat'"=="ir" | "`stat'"=="id") { checkvars t1 t0, local(found_t) checkvars i1 i0, local(found_i) if (!`found_t' & !`found_i') print_error "missing t1,t0 or i1,i0 variable(s)" if (`found_t') { erasevars i1 i0 g i1 = a1 / t1 g i0 = a0 / t0 } if (`found_i') { erasevars t1 t0 g t1 = a1 / i1 g t0 = a0 / i0 } erasevars m1 t i pexp pexp_ca pres pres_ue pres_ex g m1 = a0 + a1 g t = t1 + t0 g i = m1 / t g pexp = t1/t g pexp_ca = a1/m1 //v1.3.9: add pexp_ca, pres_ex variables g pres = m1/t g pres_ue = a0/t0 g pres_ex = a1/t1 //v1.3.9: add pexp_ca, pres_ex variables } if ("`stat'"=="md") { //check for variable names and erase existent variables checkvars m1 sd1 n1 m0 sd0 n0, errmsg checkvars md se lo up, local(md_exists) if (`md_exists') { g _md0_orig = md g _se0_orig = se g _lo0_orig = lo g _up0_orig = up } erasevars md se lo up cl tempname N s g `N' = n1+n0 //Total number of participants g `s' = sqrt(((n1-1)*sd1^2+(n0-1)*sd0^2)/(`N'-2)) //Pooled standard deviation of the two groups if (`md_std'==0) { g md = m1 - m0 //No standardization g se = sqrt(sd1^2/n1 + sd0^2/n0) } else if (`md_std'==1) { g md = (m1 - m0)/`s' //Cohen's d g se = sqrt(`N'/(n1*n0) + md^2/(2*(`N'-2))) } else if (`md_std'==2) { g md = ((m1 - m0)/`s')*(1-3/(4*`N'-9)) //Hedge's adjusted g g se = sqrt(`N'/(n1*n0) + md^2/(2*(`N'-3.94))) } else if (`md_std'==3) { g md = (m1 - m0)/sd0 //Glass's delta g se = sqrt(`N'/(n1*n0) + md^2/(2*(n0-1))) } g cl = `level' g lo = md - invnormal((cl+100)/200)*se g up = md + invnormal((cl+100)/200)*se } } return scalar zeros = `lzeros' return scalar ormh_zero = `ormh_zero' return scalar md_exists = `md_exists' end program data_sst syntax anything(name=stat), level(real) checkvars `stat', errmsg quietly{ if ("`stat'"!="r") { checkvars lo up, local(found_ci) if (!`found_ci') checkvars se, local(found_se) else local found_se 0 if (!`found_ci' & !`found_se') print_error "missing lo,up or se variable(s)" } else { *Correlation coefficient: r, n variables checkvars n, errmsg capture drop se replace r = ln(sqrt((1+r)/(1-r))) //Fisher's Z g se = 1/sqrt(n-3) //Fisher's Z Standard Error local found_ci 0 local found_se 1 } *Compute Confidece Interval if (`found_ci') { capture drop se checkvars cl if (!r(ok)) g cl = `level' if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir") g se = (ln(up)-ln(lo))/(2*invnormal((cl+100)/200)) else g se = (up-lo)/(2*invnormal((cl+100)/200)) } if (`found_se') { erasevars lo up if ("`stat'"=="or" | "`stat'"=="rr" | "`stat'"=="ir") { g lo = `stat'*exp(-invnormal((`level'+100)/200)*se) g up = `stat'*exp(invnormal((`level'+100)/200)*se) } else { g lo = `stat' - invnormal((`level'+100)/200)*se g up = `stat' + invnormal((`level'+100)/200)*se } capture drop cl g cl = `level' } } end program checkvars, rclass syntax anything(name=vars), [errmsg local(str)] foreach v in `vars' { capture confirm variable `v' local ok = (_rc==0) if (!`ok') { if ("`error_msg'"!="") print_error "missing `v' variable" break } } if "`local'" != "" { c_local `local' `"`ok'"' } return scalar ok = `ok' end program erasevars syntax anything(name=vars) foreach v in `vars' { capture confirm variable `v', exact if !_rc { qui drop `v' } } end program zero_correction, rclass syntax anything(name=type), vzeros(string) [mixed vraw(string)] tempname k1 k0 nmhi dmhi nmh dmh ormh local lzeros 0 if ("`mixed'"=="") { //raw data qui gen `vzeros' = (a1==0) | (a0==0) | (b1==0) | (b0==0) } else { //mixed data qui gen `vzeros' = (`vraw'==1) & ((a1==0) | (a0==0) | (b1==0) | (b0==0)) } qui sum `vzeros' local lzeros = r(max) scalar `ormh' = 0 if ("`type'"!="n" & `lzeros'==1) { quietly { if ("`type'"=="c") { gen `k1' = 0.5 if (`vzeros'==1) gen `k0' = 0.5 if (`vzeros'==1) } if ("`type'"=="p") { gen `k1' = n1/(n1+n0) if (`vzeros'==1) gen `k0' = n0/(n1+n0) if (`vzeros'==1) } if ("`type'"=="e") { //Compute the ORmh with the studies without zero events gen `nmhi' = a1*b0/(n1+n0) if (`vzeros'==0) gen `dmhi' = a0*b1/(n1+n0) if (`vzeros'==0) sum `nmhi' scalar `nmh' = r(sum) sum `dmhi' scalar `dmh' = r(sum) scalar `ormh' = `nmh'/`dmh' gen `k1' = (`ormh' * n1)/(n0 + `ormh'*n1) if (`vzeros'==1) gen `k0' = n0/(n0 + `ormh'*n1) if (`vzeros'==1) } replace a1 = a1 + `k1' if (`vzeros'==1) replace a0 = a0 + `k0' if (`vzeros'==1) replace b1 = b1 + `k1' if (`vzeros'==1) replace b0 = b0 + `k0' if (`vzeros'==1) replace n1 = a1 + b1 if (`vzeros'==1) replace n0 = a0 + b0 if (`vzeros'==1) if ("`mixed'"=="") { replace m1 = a0 + a1 if (`vzeros'==1) replace m0 = b0 + b1 if (`vzeros'==1) } replace n = n1 + n0 if (`vzeros'==1) } } return scalar zeros = `lzeros' return scalar ormh = `ormh' end program define print_pct syntax anything, [col(numlist max=1) a(real 5) nopercent] local p = `anything' local v = abs(`p') if (`v'<0.01 & `v'>0) local fmt = "%`a'.3f" else local fmt = cond(abs(`p')<10,cond(abs(`p')<1,"%`a'.0g","%`a'.2f"),"%`a'.1f") if ("`col'"=="") di as res `fmt' `p' _c else di as res _col(`col') `fmt' `p' _c if ("`percent'"=="") di as txt "%" _c end program define print_error args message noisily display in red "`message'" exit 198 end version 12 mata: //# struct mar_results struct mar_results { real matrix data //2x2 data: 1 a1; 2 a0; 3 b1; 4 b0; 5 n1; 6 n0; 7 m1; 8 m0; 9 n real matrix st //stat: 1 stat; 2 se; 3 loCI; 4 upCI; 5 p; 6 chi2 real matrix st_peto //PETO stat: 1 or_peto; 2 se_peto; 3 loCI_peto; 4 upCI_peto; 5 p_peto real matrix w //weights: 1 iv; 2 mh; 3 peto; 4 ivr real rowvector sw //Sum of total weights real matrix rw //relative weights: 1 iv; 2 mh; 3 peto; 4 ivr real matrix influ //influence analysis: 1 es; 2 lo; 3 up; 4 change; 5 I2; 6 p real matrix cum //cumulative analysis: 1 es_cum; 2 lo; 3 up; 4; es real rowvector est_peto //Peto combined estimations real rowvector est_mh //MH combined estimations real rowvector est_iv //IoV - Fixed effects estimations real rowvector est_ivr //IoV - Random effects estimations real rowvector a1 //2x2 tables: rowvector with terms real rowvector a0 real rowvector b1 real rowvector b0 real rowvector n1 real rowvector n0 real rowvector m1 real rowvector m0 real rowvector n real rowvector t1 real rowvector t0 real rowvector i1 real rowvector i0 real rowvector t real rowvector i string scalar type string scalar stat real scalar len real scalar level real scalar z real scalar get_peto real scalar get_mh real scalar get_influ string scalar method real scalar get_cum real scalar cum_descent } end version 12 mata: //# init void init(string scalar type, string scalar stat, real scalar nstudy, real scalar level,real scalar influ,string scalar method,real scalar cum,real scalar cumdescent) { struct mar_results scalar r real matrix temp r.type = type r.stat = stat r.level = level r.z = invnormal((r.level+100)/200) r.len = nstudy if (r.type == "raw") r.st = J(r.len+1,6,.) if (r.type == "sst") r.st = J(r.len,6,.) r.st_peto = J(r.len+1,5,.) r.w = J(r.len,4,.) r.get_peto = 0 r.get_mh = 0 r.get_influ = influ r.method = method r.influ = J(r.len,6,.) r.get_cum = cum r.cum = J(r.len,4,.) r.cum_descent = cumdescent if (r.type == "raw") { if (r.stat=="or" | r.stat=="rr" | r.stat=="rd") r.data = st_data(.,("a1","a0","b1","b0","n1","n0","m1","m0","n","__cumsort")) if (r.stat=="ir" | r.stat=="id") r.data = st_data(.,("a1","a0","t1","t0","i1","i0","m1","t","i","__cumsort")) r.data = r.data \ colsum(r.data) //Add TOTAL row if (r.stat=="ir" | r.stat=="id") { r.data[r.len+1,5] = r.data[r.len+1,1]:/r.data[r.len+1,3] r.data[r.len+1,6] = r.data[r.len+1,2]:/r.data[r.len+1,4] r.data[r.len+1,9] = r.data[r.len+1,7]:/r.data[r.len+1,8] } if (r.stat=="or" | r.stat=="rr" | r.stat=="ir") r.get_mh = 1 if (r.stat=="or") r.get_peto = 1 getterms(r) //Get individual terms getstdata(r) //Build st matrix with statistic data st_matrix("__mar_data",r.data) r.data = r.data[1::r.len,.] //Delete TOTAL row getterms(r) //Get individual terms - without the TOTAL row } else { r.data = st_data(.,(r.stat,"se","lo","up","cl","__cumsort")) getstdata(r) //Build st matrix with statistic data st_matrix("__mar_data",r.data) } getweights(r) //Compute weights if (r.type == "raw" & (r.stat == "or" | r.stat == "rr" | r.stat == "ir")) { if (r.stat == "or") { r.est_peto = peto(r,r.w[.,3],r.st_peto[1::r.len,1],r.a1,r.m1,r.n1,r.n) //Peto estimations } st_matrix("__mar_peto",r.est_peto) r.est_mh = mh(r,r.st[1::r.len,1],r.st[1::r.len,2],r.a1,r.a0,r.b1,r.b0,r.t1,r.t0,r.n,r.t,0) //MH estimations st_matrix("__mar_mh",r.est_mh) } else { st_matrix("__mar_peto",J(1,1,.)) st_matrix("__mar_mh",J(1,1,.)) } r.est_iv = iov(r,r.w[.,1],r.st[1::r.len,1]) //IoV - fixed effects estimations if (r.stat == "r") { r.est_iv[1,1] = ZtoRho(r.est_iv[1,1]) r.est_iv[1,3] = ZtoRho(r.est_iv[1,3]) r.est_iv[1,4] = ZtoRho(r.est_iv[1,4]) } st_matrix("__mar_iv",r.est_iv) getwivr(r) //Get IoV - random effects weights r.est_ivr = iov(r,r.w[.,4],r.st[1::r.len,1]) //IoV - random effects estimations if (r.stat == "r") { r.est_ivr[1,1] = ZtoRho(r.est_ivr[1,1]) r.est_ivr[1,3] = ZtoRho(r.est_ivr[1,3]) r.est_ivr[1,4] = ZtoRho(r.est_ivr[1,4]) } st_matrix("__mar_ivr",r.est_ivr) if (r.get_influ) influ(r) //Influence analysis if (r.get_cum) cum(r) //Cumulative analysis if (r.stat == "r") { r.st[.,1] = ZtoRhoM(r.st[.,1]) r.st[.,3] = ZtoRhoM(r.st[.,3]) r.st[.,4] = ZtoRhoM(r.st[.,4]) } st_matrix("__mar_st",r.st) st_matrix("__mar_stp",r.st_peto) st_matrix("__mar_rw",r.rw) st_matrix("__mar_sw",r.sw) st_matrix("__mar_influ",r.influ) st_matrix("__mar_cum",r.cum) } end version 12 mata: //# ZtoRho real scalar ZtoRho(real scalar z) { //Fisher's Z to rho correlation coefficient return( (exp(2*z)-1)/(exp(2*z)+1) ) } end version 12 mata: //# ZtoRhoM real matrix ZtoRhoM(real matrix Z) { //Fisher's Z to rho correlation coefficient return((exp(2:*Z) :- 1) :/ (exp(2:*Z) :+ 1)) } end version 12 mata: //# cum void cum(struct mar_results scalar r) { real matrix D, T, R, temp, S real scalar order temp = J(r.len,1,.) //Build the data matrix to compute the cumulative analysis and sort if (r.method=="fem") D = (r.data[.,cols(r.data)],r.w[.,1],r.st[1::r.len,1]) if (r.method=="rem") D = (r.data[.,cols(r.data)],r.w[.,1],r.st[1::r.len,1],r.st[1::r.len,2]) if (r.method=="mh") { if (r.stat=="or" | r.stat=="rr" | r.stat=="rd") { D = (r.data[.,cols(r.data)],r.st[1::r.len,1],r.st[1::r.len,2],r.a1,r.a0,r.b1,r.b0,temp,temp,r.n,temp) } if (r.stat=="ir" | r.stat=="id") { D = (r.data[.,cols(r.data)],r.st[1::r.len,1],r.st[1::r.len,2],r.a1,r.a0,temp,temp,r.t1,r.t0,temp,r.t) } } if (r.method=="peto") D = (r.data[.,cols(r.data)],r.w[.,3],r.st[1::r.len,3],r.a1,r.m1,r.n1,r.n) S = (r.data[.,cols(r.data)],r.st[1::r.len,1],r.st[1::r.len,3],r.st[1::r.len,4]) if (r.method=="peto") S = (r.data[.,cols(r.data)],r.st_peto[1::r.len,1],r.st_peto[1::r.len,3],r.st_peto[1::r.len,4]) //Sort matrix order = -1 if (r.cum_descent==0) order = 1 D = sort(D,order) S = sort(S,order) r.cum[1,1] = S[1,2] //The first study needs no calculations r.cum[1,2] = S[1,3] r.cum[1,3] = S[1,4] r.cum[1,4] = S[1,2] for (i=2; i<=r.len; i++) { T = D[|1,1 \ i, cols(D)|] //Get the accumulated matrix //Get the results for each method if (r.method=="fem") R = iov(r, T[.,2], T[.,3]) if (r.method=="mh") R = mh(r, T[.,2], T[.,3], T[.,4], T[.,5], T[.,6], T[.,7], T[.,8], T[.,9], T[.,10], T[.,11], 0) if (r.method=="peto") R = peto(r, T[.,2], T[.,3], T[.,4], T[.,5], T[.,6], T[.,7]) if (r.method=="rem") { R = iov(r, T[.,2], T[.,3]) temp = 1 :/ ((T[.,4]):^2 :+ R[1,15]) //IoV weight - Random effects R = iov(r, temp, T[.,3]) } r.cum[i,1] = R[1,1] r.cum[i,2] = R[1,3] r.cum[i,3] = R[1,4] r.cum[i,4] = S[i,2] } if (r.stat=="r") { r.cum[.,1] = ZtoRhoM(r.cum[.,1]) r.cum[.,2] = ZtoRhoM(r.cum[.,2]) r.cum[.,3] = ZtoRhoM(r.cum[.,3]) r.cum[.,4] = ZtoRhoM(r.cum[.,4]) } } end version 12 mata: //# delrow real matrix delrow(real matrix M, real scalar i, real scalar n) { //Return matrix M without the i row. The matrix is nxn real matrix R if (i==1) { R = M[|2,1 \ n, cols(M)|] //Delete the first row } else if (i==n) { R = M[|1,1 \ n-1, cols(M)|] //Delete the last row } else { R = M[|1,1 \ i-1, cols(M)|] \ M[|i+1,1 \ n, cols(M)|] //Delete the ith row } return(R) } end version 12 mata: //# influ void influ(struct mar_results scalar r) { real rowvector R real matrix w, st, a1, a0, b1, b0, n, t1, t0, t //INFLUENCE ANALYSIS for (i=1; i<=r.len; i++) { if (r.method=="fem") { w = delrow(r.w,i,r.len) st = delrow(r.st,i,r.len) R = iov(r, w[.,1], st[.,1]) if (r.stat == "r") { R[1,1] = ZtoRho(R[1,1]) R[1,3] = ZtoRho(R[1,3]) R[1,4] = ZtoRho(R[1,4]) } r.influ[i,1] = R[1,1] r.influ[i,2] = R[1,3] r.influ[i,3] = R[1,4] r.influ[i,4] = 100*(R[1,1]-r.est_iv[1,1])/r.est_iv[1,1] //Percentage of change of the statistic r.influ[i,5] = R[1,12] r.influ[i,6] = R[1,7] } if (r.method=="rem") { w = delrow(r.w,i,r.len) st = delrow(r.st,i,r.len) R = iov(r, w[.,1], st[.,1]) w = 1 :/ (st[.,2]:^2 :+ R[1,15]) R = iov(r, w, st[.,1]) if (r.stat == "r") { R[1,1] = ZtoRho(R[1,1]) R[1,3] = ZtoRho(R[1,3]) R[1,4] = ZtoRho(R[1,4]) } r.influ[i,1] = R[1,1] r.influ[i,2] = R[1,3] r.influ[i,3] = R[1,4] r.influ[i,4] = 100*(R[1,1]-r.est_ivr[1,1])/r.est_ivr[1,1] //Percentage of change of the statistic r.influ[i,5] = R[1,12] r.influ[i,6] = R[1,7] } if (r.method=="mh") { st = delrow(r.st,i,r.len) a1 = delrow(r.a1,i,r.len) a0 = delrow(r.a0,i,r.len) if (r.stat!="ir") { b1 = delrow(r.b1,i,r.len) b0 = delrow(r.b0,i,r.len) n = delrow(r.n,i,r.len) t1 = J(r.len-1,1,.) t0 = J(r.len-1,1,.) t = J(r.len-1,1,.) } else { t1 = delrow(r.t1,i,r.len) t0 = delrow(r.t0,i,r.len) t = delrow(r.t,i,r.len) b1 = J(r.len-1,1,.) b0 = J(r.len-1,1,.) n = J(r.len-1,1,.) } R = mh(r, st[.,1], st[.,2], a1, a0, b1, b0, t1, t0, n, t, 0) r.influ[i,1] = R[1,1] r.influ[i,2] = R[1,3] r.influ[i,3] = R[1,4] r.influ[i,4] = 100*(R[1,1]-r.est_mh[1,1])/r.est_mh[1,1] //Percentage of change of the statistic } if (r.method=="peto") { w = delrow(r.w,i,r.len) st = delrow(r.st,i,r.len) a1 = delrow(r.a1,i,r.len) m1 = delrow(r.m1,i,r.len) n1 = delrow(r.n1,i,r.len) n = delrow(r.n,i,r.len) R = peto(r, w[.,3], st[.,3], a1, m1, n1, n) r.influ[i,1] = R[1,1] r.influ[i,2] = R[1,3] r.influ[i,3] = R[1,4] r.influ[i,4] = 100*(R[1,1]-r.est_peto[1,1])/r.est_peto[1,1] //Percentage of change of the statistic } } } end version 12 mata: //# getterms void getterms(struct mar_results scalar r) { if (r.stat=="or" | r.stat=="rr" | r.stat=="rd") { r.a1 = r.data[.,1] r.a0 = r.data[.,2] r.b1 = r.data[.,3] r.b0 = r.data[.,4] r.n1 = r.data[.,5] r.n0 = r.data[.,6] r.m1 = r.data[.,7] r.m0 = r.data[.,8] r.n = r.data[.,9] } if (r.stat=="ir" | r.stat=="id") { r.a1 = r.data[.,1] r.a0 = r.data[.,2] r.t1 = r.data[.,3] r.t0 = r.data[.,4] r.i1 = r.data[.,5] r.i0 = r.data[.,6] r.m1 = r.data[.,7] r.t = r.data[.,8] r.i = r.data[.,9] } } end version 12 mata: //# mh real rowvector mh(struct mar_results scalar r, real matrix st, real matrix s_e, real matrix a1, /* */ real matrix a0, real matrix b1, real matrix b0, real matrix t1, real matrix t0, real matrix n, real matrix t, real scalar zerocor) { real matrix R, num, den, rmh, s, pr, psqr, qs, numse, n1, n0, m1 real scalar se, es, lo, up, p, Q, pQ if (r.stat == "or" | zerocor) { num = a1:*b0:/n den = a0:*b1:/n rmh = (a1:*b0):/n s = (a0:*b1):/n pr = ((a1:+b0):/n):*rmh psqr = ((a1:+b0):/n):*s :+ (1:-((a1:+b0):/n)):*rmh qs = (1:-((a1:+b0):/n)):*s //Standard error (or) se = sqrt(colsum(pr)/(2*colsum(rmh)^2) + colsum(psqr)/(2*colsum(rmh)*colsum(s)) + /* */ colsum(qs)/(2*colsum(s)^2)) } else { if (r.stat == "rr") { n1 = a1 :+ b1 n0 = a0 :+ b0 m1 = a0 :+ a1 num = a1:*n0:/n den = a0:*n1:/n numse = (m1:*n1:*n0:/n:^2) :- a1:*a0:/n } if (r.stat == "ir") { m1 = a0 :+ a1 num = a1:*t0:/t den = a0:*t1:/t numse = (m1:*t1:*t0:/t:^2) :- a1:*a0:/t } if (r.stat == "rr" | r.stat == "ir") { se = sqrt(colsum(numse)/(colsum(num)*colsum(den))) //Standard Error (rr, ir) } } es = colsum(num)/colsum(den) //Estimation of the statistic lo = es*exp(-r.z*se) //CI up = es*exp(r.z*se) p = 2*(1-normal(abs(ln(es)/se))) //Significance if (missing(st)==0) { Q = colsum((ln(st):-ln(es)):^2:/(s_e:^2)) pQ = 1-chi2(r.len-1, abs(Q)) } R = (es,se,lo,up,p,Q,pQ) return(R) } end version 12 mata: //# peto real rowvector peto(struct mar_results scalar r, real matrix w, real matrix st, real matrix a1, real matrix m1, real matrix n1, real matrix n) { real rowvector R, ae, Q real scalar sw, es, se, lo, up, p, Qpeto, pQpeto ae = (a1 :- (m1:*n1:/n)) //Peto combined: terms for the summatory sw = colsum(w) es = exp(colsum(ae)/sw) //Estimation of the OR Peto se = sqrt(1/sw) //SE lo = es*exp(-r.z*se) //CI up = es*exp(r.z*se) p = 2*(1-normal(abs(ln(es)/se))) //Significance Q = ((ln(st) :- ln(es)):^2) :* w //Homogeneity Q Test (Peto) Qpeto = colsum(Q) pQpeto = (1-chi2(r.len-1,abs(Qpeto))) R = (es,se,lo,up,p,Qpeto,pQpeto) return(R) } end version 12 mata: //# iov real rowvector iov(struct mar_results scalar r, real matrix w, real matrix st) { real rowvector R, num, Q real scalar se, sw, es, lo, up, p, Qiv, pQiv real scalar k, H, swH, loH, upH, I2, upI2, loI2, Tau2 //IoV - fixed & random effects k = rows(w) sw = colsum(w) se = sqrt(1/sw) //SE if (r.stat=="or" | r.stat=="rr" | r.stat=="ir") { num = w :* ln(st) es = exp(colsum(num)/sw) //Weigthed statistic lo = es * exp(-r.z*se) //CI up = es * exp(r.z*se) p = 2*(1-normal(abs(ln(es)/se))) //Significance Q = w :* (ln(st) :- ln(es)):^2 //Homogeneity Q Test } else { num = w :* st es = colsum(num)/sw //Weigthed statistic lo = es - r.z*se //CI up = es + r.z*se p = 2*(1-normal(abs(es/se))) //Significance Q = w :* (st :- es):^2 //Homogeneity Q Test } Qiv = colsum(Q) pQiv = (1- chi2(k-1,abs(Qiv))) //Heterogeneity Measures //Relative excess H = sqrt(Qiv/(k-1)) seH = 0 if (Qiv>k) { seH = (ln(Qiv)-ln(k-1))/(2*(sqrt(2*Qiv) - sqrt(2*k-3))) } else { if (k>2) seH = sqrt((1/(2*(k-2))) * (1-1/(3*(k-2)^2))) } loH = H*exp(-r.z*seH) upH = H*exp(r.z*seH) //% of Variance that is due heterogeneity between studies: I2 Statistic I2 = 100*(H^2-1)/H^2 if (I2<0) I2 = 0 loI2 = 100*(loH^2-1)/loH^2 if (loI2<1) loI2 = 0 upI2 = 100*(upH^2-1)/upH^2 if (upI2<1) upI2 = 0 //Variance between studies (Tau2) for IOV model Tau2 = (Qiv-k+1)/(sw-colsum(w:^2)/sw) if (Tau2<0) Tau2 = 0 R = (es,se,lo,up,p,Qiv,pQiv,H,seH,loH,upH,I2,loI2,upI2,Tau2) return(R) } end version 12 mata: //# getstdata void getstdata(struct mar_results scalar r) { real rowvector t if (r.type == "raw") { if (r.stat == "or") { r.st[.,1] = (r.a1 :/ r.b1) :/ (r.a0 :/ r.b0) //stat=or,Odds Ratio r.st[.,2] = sqrt(1:/r.a1 + 1:/r.b1 + 1:/r.a0 + 1:/r.b0) t = (r.m1 :* (r.n1:/r.n)) :* ((r.n0 :* r.m0) :/ (r.n :* (r.n:-1))) //Peto Weight - temporary r.st_peto[.,1] = exp((r.a1 - (r.m1 :* (r.n1:/r.n))) :/ t) //Peto OR r.st_peto[.,2] = sqrt(1:/t) //Peto Standard Error } if (r.stat == "rr") { r.st[.,1] = (r.a1 :/ r.n1) :/ (r.a0 :/ r.n0) //stat=rr,Relative Risk r.st[.,2] = sqrt(1:/r.a1 - 1:/r.n1 + 1:/r.a0 - 1:/r.n0) } if (r.stat == "ir") { r.st[.,1] = (r.a1 :/ r.t1) :/ (r.a0 :/ r.t0) //stat=ir,Rate Ratio r.st[.,2] = sqrt(1:/r.a1 :+ 1:/r.a0) } if (r.stat == "rd") { r.st[.,1] = (r.a1 :/ r.n1) - (r.a0 :/ r.n0) //stat=rd,Risk Difference r.st[.,2] = sqrt(((r.a1:/r.n1):*(1:-(r.a1:/r.n1))):/r.n1 + ((r.a0:/r.n0):*(1:-(r.a0:/r.n0))):/r.n0) } if (r.stat == "id") { r.st[.,1] = (r.a1 :/ r.t1) :- (r.a0 :/ r.t0) //stat=id,Rate Difference r.st[.,2] = sqrt(r.a1:/r.t1:^2 :+ r.a0:/r.t0:^2) } //Compute bounds for the statistic if (r.stat=="or" | r.stat=="rr" | r.stat=="ir") { r.st[.,3] = r.st[.,1] :* exp(-r.z * r.st[.,2]) r.st[.,4] = r.st[.,1] :* exp(r.z * r.st[.,2]) if (r.get_peto) { r.st_peto[.,3] = r.st_peto[.,1] :* exp(-r.z * r.st_peto[.,2]) r.st_peto[.,4] = r.st_peto[.,1] :* exp(r.z * r.st_peto[.,2]) } } else { r.st[.,3] = r.st[.,1] - (r.z * r.st[.,2]) r.st[.,4] = r.st[.,1] + (r.z * r.st[.,2]) } //MH Chi-Square and signification (p) if (r.stat=="or" | r.stat=="rr" | r.stat=="rd") { r.st[.,6] = (r.n:-1) :* ((r.a1:*r.n :- r.n1:*r.m1):^2) :/ (r.m1:*r.m0:*r.n1:*r.n0) } else { r.st[.,6] = ((r.a1:*r.t - r.t1:*r.m1):^2) :/ (r.m1:*r.t1:*r.t0) } r.st[.,5] = 1 :- chi2(1,abs(r.st[.,6])) if (r.get_peto) r.st_peto[.,5] = 2 :* (1 :- normal(abs(ln(r.st_peto[.,1]):/r.st_peto[.,2]))) } if (r.type == "sst") { r.st[.,1] = r.data[.,1] r.st[.,2] = r.data[.,2] r.st[.,3] = r.data[.,3] r.st[.,4] = r.data[.,4] } } end version 12 mata: //# getweights void getweights(struct mar_results scalar r) { real rowvector sw r.w[.,1] = 1 :/ (r.st[1::r.len,2]):^2 //IoV weight - fixed effects if (r.get_mh) { //MH weight if (r.stat=="or") r.w[.,2] = r.a0 :* r.b1 :/ r.n if (r.stat=="rr") r.w[.,2] = r.a0 :* r.n1 :/ r.n if (r.stat=="ir") r.w[.,2] = r.a0 :* r.t1 :/ r.t } if (r.get_peto) r.w[.,3] = (r.m1 :* (r.n1:/r.n)) :* ((r.n0 :* r.m0) :/ (r.n :* (r.n:-1))) //Peto Weight //Relative weights r.sw = colsum(r.w) r.rw = J(r.len,4,.) r.rw[.,1] = (r.w[.,1] :/ r.sw[1])*100 r.rw[.,2] = (r.w[.,2] :/ r.sw[2])*100 r.rw[.,3] = (r.w[.,3] :/ r.sw[3])*100 //IoV weights - Random effects must be computed after IoV - Fixed: see getwivr } end version 12 mata: //# getwivr void getwivr(struct mar_results scalar r) { real scalar Tau2 Tau2 = r.est_iv[15] r.w[.,4] = 1 :/ ((r.st[1::r.len,2]):^2 :+ Tau2) //IoV weight - Random effects //Relative weights r.sw[4] = colsum(r.w[.,4]) r.rw[.,4] = (r.w[.,4] :/ r.sw[4])*100 } end