// Stand-alone version of subroutine metan.PerformPoolingIV // contains in-line Mata code for on-the-fly compiling // for use with Stata versions prior to 16.1 // If version 16.1+, -metan- can use the pre-compiled Mata routines in lmetan.mlib // Note: // - Content of main routine metan_pooling is identical to that of subroutine metan.PerformPoolingIV // - Content of subroutine metan_pooling.Heterogi is identical to that of subroutine metan.Heterogi // - Content of Mata subroutines is identical to that of compiled Mata library lmetan.mlib *! version 4.05 29nov2021 *! Current version by David Fisher *! Previous versions by Ross Harris and Michael Bradburn program define metan_pooling, rclass syntax varlist(numeric min=2 max=2) [if] [in], MODEL(name) /// [SUMMSTAT(name) TESTSTAT(name) QSTAT(passthru) TESTBased ISQParam /// OEVLIST(varlist numeric min=2 max=2) INVLIST(varlist numeric min=2 max=6) /// NPTS(varname numeric) WGT(varname numeric) WTVAR(varname numeric) /// HKsj KRoger BArtlett SKovgaard RObust LOGRank PRoportion TN(string) /*noTRUNCate*/ TRUNCate(string) EIM OIM /// ISQSA(real 80) TSQSA(real -99) QWT(varname numeric) INIT(name) OLevel(cilevel) HLevel(cilevel) RFLevel(cilevel) /// ITOL(real 1.0x-1a) MAXTausq(real -9) REPS(real 1000) MAXITer(real 1000) QUADPTS(real 100) DIFficult TECHnique(string) * ] // N.B. extra options should just be those allowed for PerformPoolingMH marksample touse // note: *NO* "novarlist" option here local pvlist `varlist' // for clarity // if no wtvar, gen as tempvar if `"`wtvar'"'==`""' { local wtvar tempvar wtvar qui gen `wtvar' = . } // Firstly, check whether only one study // if so, cancel random-effects and set to defaults: iv, cochranq, testbased // t-critval ==> es, se, lci, uci returned but nothing else tempname k qui count if `touse' scalar `k' = r(N) if `k' == 1 { if "`model'"!="peto" { local model iv local isqparam } if "`teststat'"!="z" local teststat z local hksj c_local nsg = 1 } ** Chi-squared test (OR only; includes Peto OR) or logrank HR if "`oevlist'"!="" { tokenize `oevlist' args oe va tempname OE VA chi2 summ `oe' if `touse', meanonly scalar `OE' = cond(r(N), r(sum), .) summ `va' if `touse', meanonly scalar `VA' = cond(r(N), r(sum), .) scalar `chi2' = (`OE'^2 )/`VA' if "`model'"=="peto" | "`logrank'"!="" { return scalar OE = `OE' return scalar V = `VA' } } ************************************* * Standard inverse-variance methods * ************************************* tokenize `pvlist' args _ES _seES tempname eff se_eff crit pvalue qui replace `wtvar' = 1/`_seES'^2 if `touse' summ `_ES' [aw=`wtvar'] if `touse', meanonly scalar `eff' = r(mean) scalar `se_eff' = 1/sqrt(r(sum_w)) // I-V common-effect SE // Derive Cochran's Q tempvar qhet qui gen double `qhet' = `wtvar'*((`_ES' - `eff')^2) summ `qhet' if `touse', meanonly tempname Q Qdf scalar `Q' = cond(r(N), r(sum), .) scalar `Qdf' = cond(r(N), r(N)-1, .) // Derive sigmasq and tausq tempname c sigmasq tausq summ `wtvar' [aw=`wtvar'] if `touse', meanonly scalar `c' = r(sum_w) - r(mean) scalar `sigmasq' = `Qdf'/`c' // [general note: can this be generalised to other (non-IV) methods?] scalar `tausq' = max(0, (`Q' - `Qdf')/`c') // default: D+L estimator ********************************** * Non-iterative tausq estimators * ********************************** // (other than D+L, already derived above) ** Setup two-stage estimators sj2s and dk2s // consider *initial* estimate of tsq if inlist("`model'", "sj2s", "dk2s") { local final `model' local model `"`init'"' if substr(trim(`"`model'"'), 1, 2)==`"sa"' { tempname tausq0 scalar `tausq0' = `tausq' _parse comma model 0 : model syntax [, ISQ(string) TAUSQ(string)] if `"`tausq'"'!=`""' & `"`isq'"'==`""' { nois disp as err `"Only one of {bf:isq()} or {bf:tausq()} may be supplied as suboptions to {bf:sa()}"' exit 184 } else if `"`tausq'"'!=`""' { cap confirm number `tausq' if _rc { nois disp as err `"Error in {bf:tausq()} suboption to {bf:sa()}; a single number was expected"' exit _rc } if `tausq' < 0 { nois disp as err `"tau{c 178} value for sensitivity analysis cannot be negative"' exit 198 } local tsqsa = `tausq' local isqsa } else { if `"`isq'"'==`""' local isq = 80 else { cap confirm number `isq' if _rc { nois disp as err `"Error in {bf:isq()} suboption to {bf:sa()}; a single number was expected"' exit _rc } if `isq'<0 | `isq'>=100 { nois disp as err `"I{c 178} value for sensitivity analysis must be at least 0% and less than 100%"' exit 198 } } local isqsa = `isq' local tsqsa = -99 } scalar `tausq' = `tausq0' } } ** Hartung-Makambi estimator (>0) if "`model'"=="hm" { scalar `tausq' = `Q'^2 / (`c'*(`Q' + 2*`Qdf')) } ** Non-iterative, making use of the sampling variance of _ES else if inlist("`model'", "ev", "he", "b0", "bp") { tempname var_eff meanv qui summ `_ES' if `touse' scalar `var_eff' = r(Var) tempvar v qui gen double `v' = `_seES'^2 summ `v' if `touse', meanonly scalar `meanv' = r(mean) // empirical variance (>0) if "`model'"=="ev" { tempvar residsq qui gen double `residsq' = (`_ES' - r(mean))^2 summ `residsq', meanonly scalar `tausq' = r(mean) } // Hedges aka "variance component" aka Cochran ANOVA-type estimator else if "`model'"=="he" { scalar `tausq' = `var_eff' - `meanv' } // Rukhin Bayes estimators else if inlist("`model'", "b0", "bp") { scalar `tausq' = `var_eff'*(`k' - 1)/(`k' + 1) if "`model'"=="b0" { summ `npts' if `touse', meanonly scalar `tausq' = `tausq' - ( (`r(sum)' - `k')*`Qdf'*`meanv'/((`k' + 1)*(`r(sum)' - `k' + 2)) ) } } scalar `tausq' = max(0, `tausq') // truncate at zero } // Sensitivity analysis: use given Isq/tausq and sigmasq to generate tausq/Isq else if "`model'"=="sa" { if `tsqsa'==-99 scalar `tausq' = `isqsa'*`sigmasq'/(100 - `isqsa') else scalar `tausq' = `tsqsa' } ****************************** * Iterative tausq estimators * ****************************** // Check validity of iteropts cap assert (`maxtausq'>=0 & !missing(`maxtausq')) | `maxtausq'==-9 if _rc { nois disp as err "maxtausq() cannot be negative" exit 198 } cap assert `itol'>=0 & !missing(`itol') if _rc { nois disp as err "itol() cannot be negative" exit 198 } cap { assert (`maxiter'>0 & !missing(`maxiter')) assert round(`maxiter')==`maxiter' } if _rc { nois disp as err "maxiter() must be an integer greater than zero" exit 198 } // maxtausq: use 10*`tausq' if not specified // (and 10 times that for uci -- done in Mata) local maxtausq = cond(`maxtausq'==-9, max(10*`tausq', 100), `maxtausq') // Iterative, using Mata if inlist("`model'", "dlb", "mp", "ml", "pl", "reml") { // Bootstrap D+L // (Kontopantelis PLoS ONE 2013) if "`model'"=="dlb" { cap { assert (`reps'>0 & !missing(`reps')) assert round(`reps')==`reps' } if _rc { nois disp as err "reps() must be an integer greater than zero" exit 198 } cap nois mata: DLb("`_ES' `_seES'", "`touse'", `olevel', `reps') } // Mandel-Paule aka empirical Bayes // (DerSimonian and Kacker CCT 2007) // N.B. Mata routine also performs the Viechtbauer Q-profiling routine for tausq CI // (Viechtbauer Stat Med 2007; 26: 37-52) else if "`model'"=="mp" { cap nois mata: GenQ("`_ES' `_seES'", "`touse'", `hlevel', (`maxtausq', `itol', `maxiter')) } // REML // N.B. Mata routine also performs likelihood profiling to give tausq CI else if "`model'"=="reml" { local hmethod = cond("`difficult'"!="", "hybrid", "m-marquardt") // default = m-marquardt if "`technique'"=="" local technique nr // default = nr cap nois mata: REML("`_ES' `_seES'", "`touse'", `hlevel', (`maxtausq', `itol', `maxiter'), "`hmethod'", "`technique'") if `"`r(ll_negtsq)'"'!=`""' { nois disp `"{error}tau-squared value from last iteration was negative, so has been set to zero"' } return scalar converged = r(converged) return scalar tsq_var = r(tsq_var) return scalar ll = r(ll) } // ML, including Profile Likelihood // with optional Bartlett's (Huizenga Br J Math Stat Psychol 2011) or Skovgaard's (Guolo Stat Med 2012) correction to the likelihood // N.B. Mata routine also performs likelihood profiling to give tausq CI else if inlist("`model'", "ml", "pl") { local mlpl `model' if "`bartlett'"!="" local mlpl plbart else if "`skovgaard'"!="" local mlpl plskov local hmethod = cond("`difficult'"!="", "hybrid", "m-marquardt") // default = m-marquardt if "`technique'"=="" local technique nr // default = nr cap nois mata: MLPL("`_ES' `_seES'", "`touse'", (`olevel', `hlevel'), (`maxtausq', `itol', `maxiter'), "`hmethod'", "`technique'", "`mlpl'") if `"`r(ll_negtsq)'"'!=`""' { nois disp `"{error}tau-squared value from last iteration was negative, so has been set to zero"' } return scalar converged = r(converged) return scalar tsq_var = r(tsq_var) return scalar ll = r(ll) if "`model'"=="pl" { return scalar eff_lci = r(eff_lci) return scalar eff_uci = r(eff_uci) return scalar rc_eff_lci = r(rc_eff_lci) return scalar rc_eff_uci = r(rc_eff_uci) // Need to store these as scalars, in order to calculate critical values tempname chi2 z scalar `chi2' = r(lr) // Likelihood ratio test statistic scalar `z' = r(sll) // Signed log-likelihood test statistic if "`teststat'"=="chi2" return scalar chi2 = r(lr) // Bartlett's correction to the likelihood else return scalar z = r(sll) // Skovgaard's correction to the likelihood } } if _rc { if _rc==1 exit _rc // User break else if _rc==2000 exit _rc // No studies found with sufficient data to be analysed else if _rc>=3000 { nois disp as err "Mata compile-time or run-time error" exit _rc } else if _rc disp `"{error}Error(s) detected during running of Mata code; please check output"' } scalar `tausq' = r(tausq) // check tausq limits and set to missing if necessary tempname tsq_lci tsq_uci scalar `tsq_lci' = r(tsq_lci) scalar `tsq_uci' = r(tsq_uci) if "`model'"!="dlb" { scalar `tsq_lci' = cond(r(rc_tsq_lci)>1 & r(tsq_lci)!=0, ., r(tsq_lci)) scalar `tsq_uci' = cond(r(rc_tsq_uci)>1, ., r(tsq_uci)) } // return extra scalars return scalar maxtausq = `maxtausq' return scalar tsq_lci = `tsq_lci' return scalar tsq_uci = `tsq_uci' return scalar rc_tausq = r(rc_tausq) return scalar rc_tsq_lci = r(rc_tsq_lci) return scalar rc_tsq_uci = r(rc_tsq_uci) } // end if inlist("`model'", "dlb", "mp", "ml", "pl", "reml") // [i.e. iterative tausq estimators] // SEP 2020: Remove this option. If MP CI is desired, the entire MP model can be run. // (Doesn't affect forestplot, so why is it needed? just makes options more complicated) /* // Viechtbauer Q-profiling routine for tausq CI, if *not* Mandel-Paule tsq estimator // (Viechtbauer Stat Med 2007; 26: 37-52) if "`hetci'"=="qprofile" & "`model'"!="mp" { cap nois mata: GenQ("`_ES' `_seES'", "`touse'", `hlevel', (`maxtausq', `itol', `maxiter')) if _rc { if _rc==1 exit _rc // User break else if _rc==2000 exit _rc // No studies found with sufficient data to be analysed else if _rc>=3000 { nois disp as err "Mata compile-time or run-time error" exit _rc } else if _rc disp `"{error}Error(s) detected during running of Mata code; please check output"' } tempname tsq_lci tsq_uci scalar `tsq_lci' = cond(r(rc_tsq_lci)>1 & r(tsq_lci)!=0, ., r(tsq_lci)) scalar `tsq_uci' = cond(r(rc_tsq_uci)>1, ., r(tsq_uci)) // return extra scalars return scalar tsq_lci = `tsq_lci' return scalar tsq_uci = `tsq_uci' return scalar rc_tsq_lci = r(rc_tsq_lci) return scalar rc_tsq_uci = r(rc_tsq_uci) } */ // end of "Iterative, using Mata" section ****************************************************** * User-defined weights; finalise two-step estimators * ****************************************************** if `"`wgt'"'!=`""' { qui replace `wtvar' = `wgt' if `touse' } tempname Qr // will also be used for post-hoc variance correction if "`final'"!="" { qui replace `qhet' = ((`_ES' - `eff')^2)/((`_seES'^2) + `tausq') summ `qhet' if `touse', meanonly scalar `Qr' = r(sum) if "`final'"=="sj2s" { // two-step Sidik-Jonkman // scalar `tausq' = cond(`tausq'==0, `sigmasq'/99, `tausq') * `Qr'/`Qdf' // March 2018: if tsq=0, use Isq=1% scalar `tausq' = `tausq' * `Qr'/`Qdf' /* // Sidik-Jonkman's suggested confidence interval for tausq; not recommended for use tempname tsq_lci tsq_uci scalar `tsq_lci' = `tausq' * `Qdf' / invchi2(`Qdf', .5 - `hlevel'/200) scalar `tsq_uci' = `tausq' * `Qdf' / invchi2(`Qdf', .5 + `hlevel'/200) */ } else if "`final'"=="dk2s" { // two-step DerSimonian-Kacker (MM only) tempname wi1 wi2 wis1 wis2 summ `wtvar' if `touse', meanonly scalar `wi1' = r(sum) // sum of weights summ `wtvar' [aw=`wtvar'] if `touse', meanonly scalar `wi2' = r(sum) // sum of squared weights summ `wtvar' [aw=`_seES'^2] if `touse', meanonly scalar `wis1' = r(sum) // sum of weight * variance summ `wtvar' [aw=`wtvar' * (`_seES'^2)] if `touse', meanonly scalar `wis2' = r(sum) // sum of squared weight * variance scalar `tausq' = (`Qr' - (`wis1' - `wis2'/`wi1')) / (`wi1' - `wi2'/`wi1') scalar `tausq' = max(0, `tausq') // truncate at zero } local model `final' // switch back, so that `model' contains dk2s or sj2s again } ********************************* * Alternative weighting schemes * ********************************* // (not user-defined) // Quality effects (QE) model (extension of IVhet to incorporate quality scores) // (Doi et al, Contemporary Clinical Trials 2015; 45: 123-9) if "`model'"=="qe" { tempvar newqe tauqe // re-scale scores relative to highest value confirm numeric variable `qwt' summ `qwt' if `touse', meanonly qui gen double `newqe' = `qwt'/r(max) // taui and tauhati qui gen double `tauqe' = (1 - `newqe')/(`_seES'*`_seES'*`Qdf') summ `tauqe' if `touse', meanonly local sumtauqe = r(sum) summ `newqe' if `touse', meanonly if r(min) < 1 { // additional correction if any `newqe' are < 1, to avoid neg. weights tempvar newqe_adj qui gen double `newqe_adj' = `newqe' + r(sum)*`tauqe'/(`sumtauqe'*`Qdf') summ `newqe_adj' if `touse', meanonly qui replace `tauqe' = (`sumtauqe'*`k'*`newqe_adj'/r(sum)) - `tauqe' } else qui replace `tauqe' = (`sumtauqe'*`k'*`newqe'/r(sum)) - `tauqe' // Point estimate uses weights = qi/vi + tauhati qui replace `wtvar' = (`newqe'/(`_seES'^2)) + `tauqe' if `touse' } // Biggerstaff and Tweedie approximate Gamma-based weighting // (also derives a variance and confidence interval for tausq_DL) else if "`model'"=="bt" { cap nois mata: BTGamma("`_ES' `_seES'", "`touse'", "`wtvar'", `hlevel', (`maxtausq', `itol', `maxiter', `quadpts')) if _rc { if _rc==1 exit _rc else if _rc>=3000 { nois disp as err "Mata compile-time or run-time error" exit _rc } else if _rc disp `"{error}Error(s) detected during running of Mata code; please check output"' } // check tausq limits and set to missing if necessary tempname tsq_lci tsq_uci scalar `tsq_lci' = r(tsq_lci) scalar `tsq_uci' = r(tsq_uci) scalar `tsq_lci' = cond(r(rc_tsq_lci)>1 & `tsq_lci'!=0, ., `tsq_lci') scalar `tsq_uci' = cond(r(rc_tsq_uci)>1, ., `tsq_uci') // return extra scalars return scalar maxtausq = `maxtausq' return scalar rc_tausq = r(rc_tausq) return scalar tsq_var = r(tsq_var) return scalar tsq_lci = `tsq_lci' return scalar tsq_uci = `tsq_uci' return scalar rc_tsq_lci = r(rc_tsq_lci) return scalar rc_tsq_uci = r(rc_tsq_uci) } // Henmi and Copas method also belongs here // (Henmi and Copas, Stat Med 2010; DOI: 10.1002/sim.4029) // Begins along the same lines as IVhet; that is, a RE model with inv-variance weighting // but goes on to estimate the distribution of pivotal quantity U using a Gamma distribution (c.f. Biggerstaff & Tweedie). // `se_eff' is the same as IVhet, but conf. interval around `eff' is different. else if "`model'"=="hc" { cap nois mata: HC("`_ES' `_seES'", "`touse'", `olevel', (`itol', `maxiter', `quadpts')) if _rc { if _rc==1 exit _rc else if _rc>=3000 { nois disp as err "Mata compile-time or run-time error" exit _rc } else if _rc disp `"{error}Error(s) detected during running of Mata code; please check output"' } return scalar u = r(u) scalar `crit' = r(crit) scalar `pvalue' = r(p) } // end of "Alternative weighting schemes" section ********************************** * Generate pooled eff and se_eff * ********************************** // Alternative or user-defined weighting if `"`wgt'"'!=`""' | inlist("`model'", "ivhet", "qe", "bt", "hc") { // Apply weighting summ `_ES' [aw=`wtvar'] if `touse', meanonly scalar `eff' = r(mean) // Specify underlying model: IV common-effect, or random-effects with additive heterogeneity // (N.B. if *multiplicative* heterogeneity, factor simply multiplies the final pooled variance) local vi = cond("`model'"=="iv", "`_seES'^2", "`_seES'^2 + `tausq'") tempvar wtvce summ `wtvar' if `touse', meanonly qui gen double `wtvce' = (`vi') * `wtvar'^2 / r(sum)^2 summ `wtvce' if `touse', meanonly scalar `se_eff' = sqrt(r(sum)) // May 2020: // Similarly to M-H and Peto methods, re-calculate Q based on standard variance weights // but with respect to the *weighted* pooled effect size if `"`wgt'"'!=`""' { qui replace `qhet' = ((`_ES' - `eff') / `_seES')^2 summ `qhet' if `touse', meanonly scalar `Q' = cond(r(N), r(sum), .) } } // Standard weighting based on additive tau-squared // (N.B. if iv or mu, eff and se_eff have already been calculated) else if !inlist("`model'", "iv", "peto", "mu") { qui replace `wtvar' = 1/(`_seES'^2 + `tausq') if `touse' summ `_ES' [aw=`wtvar'] if `touse', meanonly scalar `eff' = r(mean) scalar `se_eff' = 1/sqrt(r(sum_w)) } // Return weights for CumInfLoop summ `wtvar' if `touse', meanonly return scalar totwt = cond(r(N), r(sum), .) // sum of (non-normalised) weights ********************************* * Post-hoc variance corrections * ********************************* // First, calculate "generalised" (i.e. random-effects) version of Cochran's Q. // Note that the multiplier sqrt(`Q'/`Qdf') is equal to Higgins & Thompson's (Stat Med 2002) `H' statistic // and that van Aert & Jackson (2019) use H* to refer to a "generalised/random-effects" H-statistic, similar to Qr. tempname Hstar scalar `Qr' = `Q' if !inlist("`model'", "iv", "peto", "mu") | "`wgt'"!="" { // Note: if I-V common-effect (e.g. for "mu"), Qr = Q and Hstar = H qui replace `qhet' = `wtvar'*((`_ES' - `eff')^2) summ `qhet' if `touse', meanonly scalar `Qr' = cond(r(N), r(sum), .) return scalar Qr = `Qr' } scalar `Hstar' = sqrt(`Qr'/`Qdf') // Multiplicative heterogeneity (e.g. Thompson and Sharp, Stat Med 1999) // (equivalent to the "full variance" estimator suggested by Sandercock // (https://metasurv.wordpress.com/2013/04/26/ // fixed-or-random-effects-how-about-the-full-variance-model-resolving-a-decades-old-bunfight) // Hartung-Knapp-Sidik-Jonkman variance estimator // (Roever et al, BMC Med Res Methodol 2015; Jackson et al, Stat Med 2017; van Aert & Jackson, Stat Med 2019) if "`model'"=="mu" | "`hksj'"!="" { tempname tcrit zcrit scalar `zcrit' = invnormal(.5 + `olevel'/200) scalar `tcrit' = invttail(`Qdf', .5 - `olevel'/200) // van Aert & Jackson 2019: truncate at z/t if `"`truncate'"'==`"zovert"' scalar `Hstar' = max(`zcrit'/`tcrit', `Hstar') else { if "`hksj'"!="" & `Hstar' < `zcrit'/`tcrit' c_local nzt = 1 // setup error display for later // (e.g.) Roever 2015: truncate at 1 // i.e. don't use if *under* dispersion present if inlist(`"`truncate'"', `"one"', `"1"') scalar `Hstar' = max(1, `Hstar') else if `"`truncate'"'!=`""' { nois disp as err `"invalid use of {bf:truncate()} option"' exit 184 } } scalar `se_eff' = `se_eff' * `Hstar' } return scalar Hstar = `Hstar' // Sidik-Jonkman robust ("sandwich-like") variance estimator // (Sidik and Jonkman, Comp Stat Data Analysis 2006) // (N.B. HKSJ estimator also described in the same paper) else if "`robust'"!="" { tempname sumwi tempvar vr_part summ `wtvar' if `touse', meanonly scalar `sumwi' = r(sum) qui gen double `vr_part' = `wtvar' * `wtvar' * ((`_ES' - `eff')^2) / (1 - (`wtvar'/`sumwi')) summ `vr_part' if `touse', meanonly scalar `se_eff' = sqrt(r(sum))/`sumwi' } // Kenward-Roger variance inflation method // (Morris et al, Stat Med 2018) else if "`kroger'"!="" { tempname wi1 wi2 wi3 nwi2 nwi3 summ `wtvar' if `touse', meanonly scalar `wi1' = r(sum) // sum of weights summ `wtvar' [aw=`wtvar'] if `touse', meanonly scalar `wi2' = r(sum) // sum of squared weights summ `wtvar' [aw=`wtvar'^2] if `touse', meanonly scalar `wi3' = r(sum) // sum of cubed weights scalar `nwi2' = `wi2'/`wi1' // "normalised" sum of squared weights [i.e. sum(wi:^2)/sum(wi)] scalar `nwi3' = `wi3'/`wi1' // "normalised" sum of cubed weights [i.e. sum(wi:^3)/sum(wi)] // expected information tempname I scalar `I' = `wi2'/2 - `nwi3' + (`nwi2'^2)/2 // observed information if "`oim'"!="" { tempvar resid resid2 tempname q2 q3 qui gen double `resid' = `_ES' - `eff' summ `resid' [aw=`wtvar'^2] if `touse', meanonly scalar `q2' = r(sum) // quadratic involving squared weights and residual qui gen double `resid2' = `resid'^2 summ `resid2' [aw=`wtvar'^3] if `touse', meanonly scalar `q3' = r(sum) // quadratic involving cubed weights and squared residual scalar `I' = max(0, (`q2'^2)/`wi1' + `q3' - `I') } // corrected se_eff [sqrt(Phi_A) in Kenward-Roger papers] tempname W V scalar `W' = 1/`I' // approximation of var(tausq) scalar `V' = (1/`wi1') + 2*`W'*(`wi3' - (`wi2'^2)/`wi1')/(`wi1'^2) scalar `se_eff' = sqrt(`V') // denominator degrees of freedom tempname A df_kr scalar `A' = `W' * (`V'*`wi2')^2 scalar `df_kr' = 2 / `A' // return scalar df_kr = `df_kr' } // check for successful pooling if missing(`eff', `se_eff') exit 2002 ********************************************** * Critical values, test statistics, p-values * ********************************************** // Predictive intervals // (uses k-2 df, c.f. Higgins & Thompson 2009; but also see e.g. http://www.metafor-project.org/doku.php/faq#for_random-effects_models_fitt) if `k' < 3 c_local nrfd = 1 // setup error display for later else { tempname rfcritval rflci rfuci scalar `rfcritval' = invttail(`k'-2, .5 - `rflevel'/200) scalar `rflci' = `eff' - `rfcritval' * sqrt(`tausq' + `se_eff'^2) scalar `rfuci' = `eff' + `rfcritval' * sqrt(`tausq' + `se_eff'^2) return scalar rflci = `rflci' return scalar rfuci = `rfuci' } // Proportions if "`proportion'"!="" { tempname z eff_lci eff_uci scalar `z' = `eff'/`se_eff' scalar `crit' = invnormal(.5 + `olevel'/200) scalar `pvalue' = 2*normal(-abs(`z')) scalar `eff_lci' = `eff' - `crit' * `se_eff' scalar `eff_uci' = `eff' + `crit' * `se_eff' if "`summstat'"=="ftukey" { tokenize `invlist' args succ _NN tempname hmean mintes maxtes prop_eff prop_lci prop_uci qui ameans `_NN' if `touse' if "`tn'"=="arithmetic" scalar `hmean' = r(mean) // Arithmetic mean else if "`tn'"=="geometric" scalar `hmean' = r(mean_g) // Geometric mean else if inlist("`tn'", "", "harmonic") scalar `hmean' = r(mean_h) // Harmonic mean (Miller 1978; default) else if "`tn'"=="ivariance" scalar `hmean' = 1/`se_eff'^2 // Barendregt & Doi's suggestion: inverse of pooled variance else { confirm number `tn' scalar `hmean' = `tn' } // recall: transform is = asin(sqrt(`succ' / (`_NN' + 1 ))) + asin(sqrt((`succ' + 1 ) / (`_NN' + 1 ))) // so to get our limits `mintes' and `maxtes', we subsitute `hmean' for `_NN', and let `succ' vary from 0 to `hmean'. scalar `mintes' = /*asin(sqrt(0 /(`hmean' + 1))) + */ asin(sqrt((0 + 1)/(`hmean' + 1 ))) scalar `maxtes' = asin(sqrt(`hmean'/(`hmean' + 1))) + asin(1) /*asin(sqrt((`hmean' + 1)/(`hmean' + 1 )))*/ // Barendregt & Doi use s/v < 2 or (1-s)/v < 2 // where s = sin(eff/2)^2 ~= d/n // and where v = se_eff ~= 1/n // ==> s/v ~= d; (1-s)/v ~= n-d // ... in order to avoid "blow up" when sin(eff) is near zero // Personal communication: "blow up" may result in confidence limits which do not include the point estimate // Therefore, test for this; and use simplified formula sin(eff)^2 in those cases if `eff' < `mintes' scalar `prop_eff' = 0 else if `eff' > `maxtes' scalar `prop_eff' = 1 else scalar `prop_eff' = 0.5 * (1 - sign(cos(`eff')) * sqrt(1 - (sin(`eff') + (sin(`eff') - 1/sin(`eff')) / `hmean')^2 ) ) if `eff_lci' < `mintes' scalar `prop_lci' = 0 else if `eff_lci' > `maxtes' scalar `prop_lci' = 1 else scalar `prop_lci' = 0.5 * (1 - sign(cos(`eff_lci')) * sqrt(1 - (sin(`eff_lci') + (sin(`eff_lci') - 1/sin(`eff_lci')) / `hmean')^2 ) ) if `eff_uci' < `mintes' scalar `prop_uci' = 0 else if `eff_uci' > `maxtes' scalar `prop_uci' = 1 else scalar `prop_uci' = 0.5 * (1 - sign(cos(`eff_uci')) * sqrt(1 - (sin(`eff_uci') + (sin(`eff_uci') - 1/sin(`eff_uci')) / `hmean')^2 ) ) cap { assert `prop_lci' <= `prop_eff' & `prop_eff' <= `prop_uci' assert !missing(`prop_eff', `prop_lci', `prop_uci') } if _rc { scalar `prop_eff' = sin(`eff' /2)^2 scalar `prop_lci' = sin(`eff_lci'/2)^2 scalar `prop_uci' = sin(`eff_uci'/2)^2 } scalar `z' = 0 if `eff' > `mintes' scalar `z' = abs(`eff' - `mintes') / `se_eff' return scalar prop_eff = `prop_eff' return scalar prop_lci = `prop_lci' return scalar prop_uci = `prop_uci' // Predictive intervals if `k' >= 3 { tempname prop_rflci prop_rfuci if `rflci' < `mintes' scalar `prop_rflci' = 0 else if `rflci' > `maxtes' scalar `prop_rflci' = 1 else scalar `prop_rflci' = 0.5 * (1 - sign(cos(`rflci')) * sqrt(1 - (sin(`rflci') + (sin(`rflci') - 1/sin(`rflci')) / `hmean')^2 ) ) if `rfuci' < `mintes' scalar `prop_rfuci' = 0 else if `rfuci' > `maxtes' scalar `prop_rfuci' = 1 else scalar `prop_rfuci' = 0.5 * (1 - sign(cos(`rfuci')) * sqrt(1 - (sin(`rfuci') + (sin(`rfuci') - 1/sin(`rfuci')) / `hmean')^2 ) ) if _rc { // from earlier "cap assert" scalar `prop_rflci' = sin(`rflci'/2)^2 scalar `prop_rfuci' = sin(`rfuci'/2)^2 } return scalar prop_rflci = `prop_rflci' return scalar prop_rfuci = `prop_rfuci' } } else if "`summstat'"=="arcsine" { return scalar prop_eff = sin(`eff')^2 return scalar prop_lci = sin(`eff_lci')^2 return scalar prop_uci = sin(`eff_uci')^2 // Predictive intervals if `k' >= 3 { return scalar prop_rflci = sin(`rflci')^2 return scalar prop_rfuci = sin(`rfuci')^2 } } else if "`summstat'"=="logit" { return scalar prop_eff = invlogit(`eff') return scalar prop_lci = invlogit(`eff_lci') return scalar prop_uci = invlogit(`eff_uci') // Predictive intervals if `k' >= 3 { return scalar prop_rflci = invlogit(`rflci') return scalar prop_rfuci = invlogit(`rfuci') } } else { return scalar prop_eff = `eff' return scalar prop_lci = `eff_lci' return scalar prop_uci = `eff_uci' // Predictive intervals if `k' >= 3 { return scalar prop_rflci = `rflci' return scalar prop_rfuci = `rfuci' } } return scalar z = `z' return scalar eff_lci = `eff_lci' return scalar eff_uci = `eff_uci' } // All other data types else { if "`model'"=="pl" { // N.B. PL confidence limits have already been calculated if "`teststat'"=="chi2" { scalar `crit' = invchi2(1, `olevel'/100) scalar `pvalue' = chi2tail(1, `chi2') } else { scalar `crit' = invnormal(.5 + `olevel'/200) scalar `pvalue' = 2*normal(-abs(`z')) } } else { if "`teststat'"=="chi2" { scalar `crit' = invchi2(1, `olevel'/100) scalar `pvalue' = chi2tail(1, `chi2') return scalar chi2 = `chi2' } else if "`kroger'"!="" { tempname t scalar `crit' = invttail(`df_kr', .5 - `olevel'/200) scalar `t' = `eff'/`se_eff' scalar `pvalue' = 2*ttail(`df_kr', abs(`t')) return scalar t = `t' } else if "`teststat'"=="t" { tempname t scalar `crit' = invttail(`Qdf', .5 - `olevel'/200) scalar `t' = `eff'/`se_eff' scalar `pvalue' = 2*ttail(`Qdf', abs(`eff'/`se_eff')) return scalar t = `t' } else if "`model'"!="hc" { // N.B. HC crit + p-value have already been calculated tempname z scalar `crit' = invnormal(.5 + `olevel'/200) scalar `z' = `eff'/`se_eff' scalar `pvalue' = 2*normal(-abs(`z')) return scalar z = `z' } // Confidence intervals if "`oevlist'"!="" { // crit.value is chi2, but CI is based on z return scalar eff_lci = `eff' - invnormal(.5 + `olevel'/200) * `se_eff' return scalar eff_uci = `eff' + invnormal(.5 + `olevel'/200) * `se_eff' } else { // else we can use crit.value (z or t, or u if HC) return scalar eff_lci = `eff' - `crit' * `se_eff' return scalar eff_uci = `eff' + `crit' * `se_eff' } } } ***************************************** * Derive other heterogeneity statistics * ***************************************** // e.g. H, I-squared and (modified) H-squared; plus Q-based confidence intervals // Sensitivity analysis // (Note: tausq has already been established, whether `tsqsa' or `Isqsa') if "`model'"=="sa" { tempname H Isqval HsqM if `tsqsa' == -99 { scalar `H' = sqrt(100 / (100 - `isqsa')) scalar `Isqval' = `isqsa' scalar `HsqM' = `isqsa'/(100 - `isqsa') return scalar H = `H' return scalar Isq = `Isqval' return scalar HsqM = float(`HsqM') // If user-defined I^2 is a round(ish) number, so should H^2 be } else { scalar `H' = sqrt((`tsqsa' + `sigmasq') / `sigmasq') scalar `Isqval' = 100*`tsqsa'/(`tsqsa' + `sigmasq') scalar `HsqM' = `tsqsa'/`sigmasq' return scalar H = `H' return scalar Isq = `Isqval' return scalar HsqM = `HsqM' } // [Sep 2020] Also save values in matrix `hetstats', same as if `isqparam' (see subroutine -heterogi- ) local t2rownames tausq tsq_lci tsq_uci H H_lci H_uci Isq Isq_lci Isq_uci HsqM HsqM_lci HsqM_uci tempname hetstats local r : word count `t2rownames' matrix define `hetstats' = J(`r', 1, .) matrix rownames `hetstats' = `t2rownames' matrix `hetstats'[rownumb(`hetstats', "tausq"), 1] = `tausq' matrix `hetstats'[rownumb(`hetstats', "H"), 1] = `H' matrix `hetstats'[rownumb(`hetstats', "Isq"), 1] = `Isqval' matrix `hetstats'[rownumb(`hetstats', "HsqM"), 1] = `HsqM' return matrix hetstats = `hetstats' } else { if !inlist("`model'", "iv", "peto", "mu") { local tausqlist `tausq' `tsq_lci' `tsq_uci' cap assert "`tausqlist'"!="" if "`isqparam'"!="" if _rc { nois disp as err "Heterogeneity confidence interval not valid" exit 198 } } cap nois Heterogi `Q' `Qdf' if `touse', `testbased' `isqparam' /// stderr(`_seES') tausqlist(`tausqlist') level(`hlevel') if _rc { if _rc==1 nois disp as err `"User break in {bf:metan.Heterogi}"' else nois disp as err `"Error in {bf:metan.Heterogi}"' c_local err noerr // tell -metan- not to also report an "error in metan.PerformPoolingIV" exit _rc } return add } // Return scalars return scalar eff = `eff' return scalar se_eff = `se_eff' return scalar crit = `crit' return scalar pvalue = `pvalue' if "`kroger'"!="" return scalar df = `df_kr' else if "`teststat'"=="t" return scalar df = `Qdf' return scalar k = `k' // k = number of studies (= count if `touse') return scalar Q = `Q' // Cochran's Q heterogeneity statistic return scalar Qdf = `Qdf' // Q degrees of freedom (= `k' - 1) return scalar sigmasq = `sigmasq' // "typical" within-study variance (Higgins & Thompson 2002) return scalar tausq = `tausq' // between-study heterogeneity variance return scalar c = `c' // scaling factor end // Based on heterogi.ado from SSC, with release notes: // version 2.0 N.Orsini, I. Buchan, 25 Jan 06 // version 1.0 N.Orsini, J.Higgins, M.Bottai, 16 Feb 2005 // (c.f. Higgins & Thompson Stat Med 2002, "Quantifying heterogeneity") program define Heterogi, rclass syntax anything [if] [in], [ TESTBased ISQParam /// STDERR(varname numeric) TAUSQLIST(namelist min=1 max=3) LEVEL(cilevel) ] marksample touse tokenize `anything' assert `"`3'"'==`""' args Q Qdf // setup W1, W2 for tausq if stderr available (e.g. not for M-H) if "`stderr'"!="" { tempvar wtvar qui gen double `wtvar' = 1/`stderr'^2 tempname W1 W2 summ `wtvar' if `touse', meanonly scalar `W1' = r(sum) // sum of weights summ `wtvar' [aw=`wtvar'] if `touse', meanonly scalar `W2' = r(sum) // sum of squared weights tempname sigmasq scalar `sigmasq' = (r(N) - 1) / (`W1' - `W2'/`W1') } ******************** * Standard Q-based * ******************** tempname Q_lci Q_uci scalar `Q_lci' = . scalar `Q_uci' = . ** Confidence intervals: // Test-based interval for ln(Q) [ or, equivalently, ln(H) ] // (Higgins & Thompson, Stats in Medicine 2002) if "`testbased'"!="" { tempname k selogQ scalar `k' = `Qdf' + 1 // Formula 26.4.13 of Abramowitz and Stegun (1965): // Z = sqrt(2Q) - sqrt(2k - 3) is standard normal // Now, expected value of Q is k-1, so form a standard normal variate as follows (taking logs to reduce skew): // Z = [ ln(Q) - ln(k-1) ] / se[ ln(Q) ] // ==> se[ ln(Q) ] = [ ln(Q) - ln(k-1) ] / [ sqrt(2Q) - sqrt(2k - 3) ] scalar `selogQ' = (ln(`Q') - ln(`Qdf')) / ( sqrt(2*`Q') - sqrt(2*`k' - 3) ) // Formula 26.4.36 of Abramowitz and Stegun (1965): // Var[ ln(Q/k-1) ] = [ 2/(k-2) ] * [ 1 - (1/ {3(k-2)^2} ) ] // (use if Q <= k) if `Q' <= `k' { scalar `selogQ' = sqrt( ( 2/(`k'-2)) * (1 - 1/(3*(`k'-2)^2)) ) } tempname Q_lci Q_uci scalar `Q_lci' = max(0, exp( ln(`Q') - invnormal(.5 + `level'/200) * `selogQ' )) scalar `Q_uci' = exp( ln(`Q') + invnormal(.5 + `level'/200) * `selogQ' ) /* // Original code from heterogi.ado // used confidence intervals for lnH rather than for lnQ, but these differ only by a constant: // If, as above, Var[ ln(Q/k-1) ] = [ 2/(k-2) ] * [ 1 - (1/ {3(k-2)^2} ) ] // then if ln(H) = .5 * ln(Q/k-1), then Var[ ln(H) ] = .25 * Var[ ln(Q/k-1) ] = [ 1/ 2(k-2) ] * [ 1 - (1/ {3(k-2)^2} ) ] scalar `selogH' = cond(`Q' > `k', /// .5*( (ln(`Q') - ln(`Qdf')) / ( sqrt(2*`Q') - sqrt(2*`k' - 3) ) ), /// sqrt( ( 1/(2*(`k'-2)) * (1 - 1/(3*(`k'-2)^2)) ) )) scalar `H_lci' = max(1, exp( ln(`H') - invnormal(.5 + `level'/200) * `selogH' )) scalar `H_uci' = exp( ln(`H') + invnormal(.5 + `level'/200) * `selogH' ) */ } // Q-based confidence intervals // using ncchi2 if fixed-effect; Gamma-based if random-effects (ref: Hedges & Pigott, 2001) // ncchi2 previously recommended by JPTH based on personal communications else { if "`tausqlist'"=="" { // fixed (common) effect tempname nc scalar `nc' = max(0, `Q' - `Qdf') // If Q < df, no need to seek the lower bound tempname Q_lci Q_uci scalar `Q_lci' = cond(`nc'==0, 0, invnchi2(`Qdf', `nc', .5 - `level'/200)) scalar `Q_uci' = invnchi2(`Qdf', `nc', .5 + `level'/200) } else { // random-effects cap assert "`stderr'"!="" if _rc { nois disp as err "Heterogeneity confidence interval not valid" exit 198 } tempname W3 tsq_dl summ `wtvar' [aw=`wtvar'^2] if `touse', meanonly scalar `W3' = r(sum) // sum of cubed weights scalar `tsq_dl' = (`Q' - `Qdf') / (`W1' - `W2'/`W1') // non-truncated tsq_DL tempname btVarQ scalar `btVarQ' = 2*`Qdf' + 4*`tsq_dl'*(`W1' - `W2'/`W1') + 2*(`tsq_dl'^2)*(`W2' - 2*`W3'/`W1' + (`W2'/`W1')^2) // If Q < df, no need to seek the lower bound tempname Q_lci Q_uci scalar `Q_lci' = cond(`Q' < `Qdf', 0, invgammap(`Q'^2 / `btVarQ', .5 - `level'/200) * `btVarQ' / `Q') scalar `Q_uci' = invgammap(`Q'^2 / `btVarQ', .5 + `level'/200) * `btVarQ' / `Q' } } // standard, transformed CIs for Isq, as outputted by heterogi.ado // Taken from heterogi.ado by N.Orsini, J.Higgins, M.Bottai, N.Buchan (2005-2006) return scalar Q_lci = `Q_lci' return scalar Q_uci = `Q_uci' return scalar H = max(1, sqrt(`Q' / `Qdf')) return scalar H_lci = max(1, sqrt(`Q_lci' / `Qdf')) return scalar H_uci = sqrt(`Q_uci' / `Qdf') return scalar Isq = 100* max(0, (`Q' - `Qdf') / `Q') return scalar Isq_lci = 100* max(0, (`Q_lci' - `Qdf') / `Q_lci') return scalar Isq_uci = 100* min(1, (`Q_uci' - `Qdf') / `Q_uci') return scalar HsqM = max(0, (`Q' - `Qdf') / `Qdf') return scalar HsqM_lci = max(0, (`Q_lci' - `Qdf') / `Qdf') return scalar HsqM_uci = max(0, (`Q_uci' - `Qdf') / `Qdf') ********************* * Tau-squared based * ********************* if "`isqparam'"!="" { tokenize `tausqlist' args tausq tsq_lci tsq_uci // Save values in matrix `hetstats' local t2rownames tausq tsq_lci tsq_uci H H_lci H_uci Isq Isq_lci Isq_uci HsqM HsqM_lci HsqM_uci tempname hetstats local r : word count `t2rownames' matrix define `hetstats' = J(`r', 1, .) matrix rownames `hetstats' = `t2rownames' if "`tausqlist'"!="" { matrix `hetstats'[rownumb(`hetstats', "tausq"), 1] = `tausq' matrix `hetstats'[rownumb(`hetstats', "H"), 1] = sqrt((`tausq' + `sigmasq') / `sigmasq') matrix `hetstats'[rownumb(`hetstats', "Isq"), 1] = 100* `tausq' / (`tausq' + `sigmasq') matrix `hetstats'[rownumb(`hetstats', "HsqM"), 1] = `tausq' / `sigmasq' } // If `tausq' not defined for this model, store H, Isq and HsqM (& CIs) based on Q instead else { matrix `hetstats'[rownumb(`hetstats', "H"), 1] = max(1, sqrt(`Q' / `Qdf')) matrix `hetstats'[rownumb(`hetstats', "Isq"), 1] = 100* max(0, (`Q' - `Qdf') / `Q') matrix `hetstats'[rownumb(`hetstats', "HsqM"), 1] = max(0, (`Q' - `Qdf') / `Qdf') matrix `hetstats'[rownumb(`hetstats', "H_lci"), 1] = max(1, sqrt(`Q_lci' / `Qdf')) matrix `hetstats'[rownumb(`hetstats', "Isq_lci"), 1] = 100* max(0, (`Q_lci' - `Qdf') / `Q_lci') matrix `hetstats'[rownumb(`hetstats', "HsqM_lci"), 1] = max(0, (`Q_lci' - `Qdf') / `Qdf') matrix `hetstats'[rownumb(`hetstats', "H_uci"), 1] = max(1, sqrt(`Q_uci' / `Qdf')) matrix `hetstats'[rownumb(`hetstats', "Isq_uci"), 1] = 100* max(0, (`Q_uci' - `Qdf') / `Q_uci') matrix `hetstats'[rownumb(`hetstats', "HsqM_uci"), 1] = max(0, (`Q_uci' - `Qdf') / `Qdf') } // Confidence intervals, if appropriate if `"`tsq_lci'"'!=`""' { matrix `hetstats'[rownumb(`hetstats', "tsq_lci"), 1] = `tsq_lci' matrix `hetstats'[rownumb(`hetstats', "H_lci"), 1] = sqrt((`tsq_lci' + `sigmasq') / `sigmasq') matrix `hetstats'[rownumb(`hetstats', "Isq_lci"), 1] = 100* `tsq_lci' / (`tsq_lci' + `sigmasq') matrix `hetstats'[rownumb(`hetstats', "HsqM_lci"), 1] = `tsq_lci' / `sigmasq' matrix `hetstats'[rownumb(`hetstats', "tsq_uci"), 1] = `tsq_uci' matrix `hetstats'[rownumb(`hetstats', "H_uci"), 1] = sqrt((`tsq_uci' + `sigmasq') / `sigmasq') matrix `hetstats'[rownumb(`hetstats', "Isq_uci"), 1] = 100* `tsq_uci' / (`tsq_uci' + `sigmasq') matrix `hetstats'[rownumb(`hetstats', "HsqM_uci"), 1] = `tsq_uci' / `sigmasq' } return matrix hetstats = `hetstats' } end ******************************************************************************** ******************** * Mata subroutines * (for iterative methods) ******************** version 11.0 mata: /* Kontopantelis's bootstrap DerSimonian-Laird estimator */ // (PLoS ONE 2013; 8(7): e69930, and also implemented in metaan) // N.B. using originally estimated ES within the re-samples, as in Kontopantelis's paper */ void DLb(string scalar varlist, string scalar touse, real scalar level, real scalar reps) { // setup real colvector yi, se, vi, wi varlist = tokens(varlist) st_view(yi=., ., varlist[1], touse) if(length(yi)==0) exit(error(2000)) st_view(se=., ., varlist[2], touse) vi = se:^2 wi = 1:/vi // calculate I-V Common eff real scalar eff eff = mean(yi, wi) // carry out bootstrap procedure transmorphic B, J real colvector report B = mm_bs(&ftausq(), (yi, vi), 1, reps, 0, 1, ., ., ., eff) J = mm_jk(&ftausq(), (yi, vi), 1, 1, ., ., ., ., ., eff) report = mm_bs_report(B, ("mean", "bca"), level, 0, J) // truncate at zero report = report:*(report:>0) // return tausq and confidence limits real scalar tausq tausq = report[1] st_numscalar("r(tausq)", tausq) st_numscalar("r(tsq_lci)", report[2]) st_numscalar("r(tsq_uci)", report[3]) } real scalar ftausq(real matrix coeffs, real colvector weight, real scalar eff) { real colvector yi, vi, wi real scalar k, Q, c, tausq yi = select(coeffs[,1], weight) vi = select(coeffs[,2], weight) k = length(yi) wi = 1:/vi Q = crossdev(yi, eff, wi, yi, eff) c = sum(wi) - mean(wi, wi) tausq = max((0, (Q-(k-1))/c)) return(tausq) } /* "Generalised Q" methods */ void GenQ(string scalar varlist, string scalar touse, real scalar hlevel, real rowvector iteropts) { // setup real colvector yi, se, vi, wi varlist = tokens(varlist) st_view(yi=., ., varlist[1], touse) if(length(yi)==0) exit(error(2000)) st_view(se=., ., varlist[2], touse) vi = se:^2 wi = 1:/vi real scalar maxtausq, itol, maxiter maxtausq = iteropts[1] itol = iteropts[2] maxiter = iteropts[3] real scalar k k = length(yi) /* Mandel-Paule estimator of tausq (J Res Natl Bur Stand 1982; 87: 377-85) */ // (also DerSimonian & Kacker, Contemporary Clinical Trials 2007; 28: 105-114) // ... can be shown to be equivalent to the "empirical Bayes" estimator // (e.g. Sidik & Jonkman Stat Med 2007; 26: 1964-81) // and converges more quickly real scalar rc_tausq, tausq rc_tausq = mm_root(tausq=., &Q_crit(), 0, maxtausq, itol, maxiter, yi, vi, k, k-1) st_numscalar("r(tausq)", tausq) st_numscalar("r(rc_tausq)", rc_tausq) /* Confidence interval for tausq by generalised Q-profiling */ // Viechtbauer Stat Med 2007; 26: 37-52 // (N.B. most natural point estimate is Mandel-Paule, but any estimate will do) real scalar eff, Qmin, Qmax eff = mean(yi, wi) // I-V common-effect estimate Qmin = crossdev(yi, eff, wi, yi, eff) // Q(0) = standard Cochran's Q heterogeneity statistic (when tausq=0) wi = 1:/(vi:+maxtausq) eff = mean(yi, wi) Qmax = crossdev(yi, eff, wi, yi, eff) // estimate tausq confidence limits real scalar Q_crit_hi, Q_crit_lo, tsq_lci, rc_tsq_lci, tsq_uci, rc_tsq_uci Q_crit_hi = invchi2(k-1, .5 + hlevel/200) // higher critical value (0.975) to compare GenQ against (for *lower* bound of tausq) Q_crit_lo = invchi2(k-1, .5 - hlevel/200) // lower critical value (0.025) to compare GenQ against (for *upper* bound of tausq) if (Qmin < Q_crit_lo) { // if Q(0) is less the lower critical value, interval is set to null rc_tsq_lci = 2 rc_tsq_uci = 2 tsq_lci = 0 tsq_uci = 0 } else { if (Qmax > Q_crit_lo) { // If Q(maxtausq) is larger than the lower critical value... rc_tsq_uci = 2 tsq_uci = maxtausq // ...upper bound for tausq is tausqmax } else { rc_tsq_uci = mm_root(tsq_uci=., &Q_crit(), 0, maxtausq, itol, maxiter, yi, vi, k, Q_crit_lo) } } if (Qmax > Q_crit_hi) { // If Q(maxtausq) is larger than the higher critical value, interval is set to null rc_tsq_lci = 2 rc_tsq_uci = 2 tsq_lci = maxtausq tsq_uci = maxtausq } else { if (Qmin < Q_crit_hi) { // If Q(0) is less than the higher critical value... rc_tsq_lci = 2 tsq_lci = 0 // ...lower bound for tausq is 0 } else { rc_tsq_lci = mm_root(tsq_lci=., &Q_crit(), 0, maxtausq, itol, maxiter, yi, vi, k, Q_crit_hi) } } // return confidence limits and rc codes st_numscalar("r(tsq_lci)", tsq_lci) st_numscalar("r(tsq_uci)", tsq_uci) st_numscalar("r(rc_tsq_lci)", rc_tsq_lci) st_numscalar("r(rc_tsq_uci)", rc_tsq_uci) } real scalar Q_crit(real scalar tausq, real colvector yi, real colvector vi, real scalar k, real scalar crit) { real colvector wi real scalar eff, newtausq wi = 1:/(vi:+tausq) eff = mean(yi, wi) newtausq = (k/crit)*crossdev(yi, eff, wi, yi, eff)/sum(wi) - mean(vi, wi) // corrected June 2015 return(tausq - newtausq) } /* ML + optional PL (for likelihood profiling for ES CI) */ // (N.B. pass wi back-and-forth as it needs to be calculated anyway for tausq likelihood profiling) void MLPL(string scalar varlist, string scalar touse, real rowvector levels, real rowvector iteropts, string scalar hmethod, string scalar technique, string scalar model) { // setup real colvector yi, se, vi, wi varlist = tokens(varlist) st_view(yi=., ., varlist[1], touse) if(length(yi)==0) exit(error(111)) st_view(se=., ., varlist[2], touse) vi = se:^2 wi = 1:/vi // Initialize real scalar eff0, tausq, eff eff0 = mean(yi, wi) // Effect size with zero tausq tausq = max((0, quadvariance(yi) - mean(vi))) // Initialize tausq using Hedges estimator wi = 1:/(vi:+tausq) eff = mean(yi, wi) // Initialize eff using Hedges estimator // Maximize log-likelihood for ML transmorphic S real rowvector p S = optimize_init() optimize_init_evaluator(S, &ML_est()) optimize_init_evaluatortype(S, "d2") optimize_init_params(S, (eff, tausq)) optimize_init_argument(S, 1, yi) optimize_init_argument(S, 2, vi) optimize_init_argument(S, 3, .) // ML_est() can also estimate tausq with eff held constant; not relevant here optimize_init_technique(S, technique) optimize_init_singularHmethod(S, hmethod) optimize_init_tracelevel(S, "none") p = optimize(S) real scalar rc rc = optimize_result_returncode(S) if(rc) exit(error(rc)) real scalar ll ll = optimize_result_value(S) eff = p[1] tausq = p[2] if(tausq < 0) { tausq = 0 eff = eff0 st_numscalar("r(ll_negtsq)", ll) ll = sum(lnnormalden(yi, eff, sqrt(vi))) } wi = 1:/(vi:+tausq) st_numscalar("r(tausq)", tausq) st_numscalar("r(converged)", optimize_result_converged(S)) st_numscalar("r(ll)", ll) // Variance of tausq (using inverse Fisher information) real scalar tsq_var tsq_var = optimize_result_V(S)[2,2] st_numscalar("r(tsq_var)", tsq_var) // Confidence interval for tausq using likelihood profiling real scalar maxtausq, itol, maxiter maxtausq = iteropts[1] itol = iteropts[2] maxiter = iteropts[3] real scalar level, hlevel, crit level = levels[1] hlevel = levels[2] crit = ll - invchi2(1, hlevel/100)/2 real scalar tsq_lci, rc_tsq_lci, tsq_uci, rc_tsq_uci rc_tsq_lci = mm_root(tsq_lci=., &ML_profile_tausq(), 0, tausq - itol, itol, maxiter, yi, vi, crit) st_numscalar("r(tsq_lci)", tsq_lci) st_numscalar("r(rc_tsq_lci)", rc_tsq_lci) rc_tsq_uci = mm_root(tsq_uci=., &ML_profile_tausq(), tausq + itol, 10*maxtausq, itol, maxiter, yi, vi, crit) st_numscalar("r(tsq_uci)", tsq_uci) st_numscalar("r(rc_tsq_uci)", rc_tsq_uci) // Profile likelihood if (model!="ml") { // Bartlett's correction // (see e.g. Huizenga et al, Br J Math Stat Psychol 2011) real scalar BCFinv BCFinv = 1 if (model=="plbart") { BCFinv = 1 + 2*mean(wi, wi:^2)/sum(wi) - 0.5*mean(wi, wi)/sum(wi) st_numscalar("r(BCF)", 1/BCFinv) } // Log-likelihood based test statistic // (evaluated at b = 0) crit = ll - invchi2(1, level/100)*BCFinv/2 S = optimize_init() optimize_init_evaluator(S, &ML_est()) optimize_init_evaluatortype(S, "d2") optimize_init_params(S, tausq) optimize_init_argument(S, 1, yi) optimize_init_argument(S, 2, vi) optimize_init_argument(S, 3, 0) // estimate tausq with b held constant at zero optimize_init_technique(S, technique) optimize_init_singularHmethod(S, hmethod) optimize_init_tracelevel(S, "none") real scalar tausq0, rc_ll0, ll0, lr tausq0 = optimize(S) rc_ll0 = optimize_result_returncode(S) if(rc_ll0) exit(error(rc_ll0)) ll0 = optimize_result_value(S) if(tausq0 < 0) { tausq0 = 0 ll0 = sum(lnnormalden(yi, 0, sqrt(vi))) } if (abs(ll0 - ll) <= itol) lr = 0 // in case ll, ll_b are very close (within itol) and/or rounding error results in a negative value else lr = 2*(ll - ll0) / BCFinv // Signed log-likelihood statistic // (evaluated at b = 0) real scalar sll if (lr==0) sll = 0 else sll = sign(eff)*sqrt(lr) // Confidence interval for ES using likelihood profiling // (use ten times the ML lci and uci for search limits) real scalar llim, ulim, eff_lci, eff_uci, rc_eff_lci, rc_eff_uci llim = eff - 19.6/sqrt(sum(wi)) ulim = eff + 19.6/sqrt(sum(wi)) // Skovgaard's correction to the signed likelihood statistic if (model=="plskov") { // Collect ML values of eff, tausq, ll to send to ML_skov() real rowvector params params = (eff, tausq, ll) // can't directly correct the critical value, due to the square root (i.e. expression is non-linear) // so instead need to pass the critical value to the iteration procedure, and correct afterwards crit = invnormal(.5 + level/200) sll = ML_skov(0, yi, vi, wi, params, crit, iteropts, hmethod, technique) // find SLL for b fixed at zero sll = sll + crit // ML_skov() returns sll-crit, so add crit back on rc_eff_lci = mm_root(eff_lci=., &ML_skov(), llim, eff-itol, itol, maxiter, yi, vi, wi, params, crit, iteropts, hmethod, technique) rc_eff_uci = mm_root(eff_uci=., &ML_skov(), eff+itol, ulim, itol, maxiter, yi, vi, wi, params, -crit, iteropts, hmethod, technique) st_numscalar("r(eff_lci)", eff_lci) st_numscalar("r(eff_uci)", eff_uci) st_numscalar("r(rc_eff_lci)", rc_eff_lci) st_numscalar("r(rc_eff_uci)", rc_eff_uci) } // Otherwise, use the (squared) likelihood statistic LR = SLL^2 else { rc_eff_lci = mm_root(eff_lci=., &ML_profile_eff(), llim, eff, itol, maxiter, yi, vi, crit, tausq, iteropts, hmethod, technique) rc_eff_uci = mm_root(eff_uci=., &ML_profile_eff(), eff, ulim, itol, maxiter, yi, vi, crit, tausq, iteropts, hmethod, technique) st_numscalar("r(eff_lci)", eff_lci) st_numscalar("r(eff_uci)", eff_uci) st_numscalar("r(rc_eff_lci)", rc_eff_lci) st_numscalar("r(rc_eff_uci)", rc_eff_uci) } st_numscalar("r(ll)", ll) st_numscalar("r(lr)", lr) st_numscalar("r(sll)", sll) } } void ML_est(todo, p, yi, vi, eff_cons, lnf, S, H) { real scalar eff, tausq real colvector wi if(cols(p)==2) { eff = p[1] tausq = max((0, p[2])) lnf = sum(lnnormalden(yi, eff, sqrt(vi:+tausq))) if(todo>=1) { wi = 1:/(vi:+tausq) S = J(1, 2, .) S[1] = quadcross(wi, yi:-eff) S[2] = -0.5*sum(wi) + 0.5*quadcrossdev(yi, eff, wi:^2, yi, eff) if(todo>=2) { H = J(2, 2, .) H[1, 1] = -sum(wi) H[2, 1] = -quadcross(wi:^2, yi:-eff) H[2, 2] = 0.5*quadcross(wi, wi) - quadcrossdev(yi, eff, wi:^3, yi, eff) _makesymmetric(H) } } } else { eff = eff_cons // estimate tausq with eff held constant tausq = max((0, p)) lnf = sum(lnnormalden(yi, eff, sqrt(vi:+tausq))) if(todo>=1) { wi = 1:/(vi:+tausq) S = -0.5*sum(wi) + 0.5*quadcrossdev(yi, eff, wi:^2, yi, eff) if(todo>=2) { H = 0.5*quadcross(wi, wi) - quadcrossdev(yi, eff, wi:^3, yi, eff) } } } } real scalar ML_profile_tausq(real scalar tausq, real colvector yi, real colvector vi, real scalar crit) { real colvector wi real scalar eff, ll wi = 1:/(vi:+tausq) eff = mean(yi, wi) ll = sum(lnnormalden(yi, eff, sqrt(vi:+tausq))) return(ll - crit) } real scalar ML_profile_eff(real scalar eff, real colvector yi, real colvector vi, real scalar crit, real scalar tausq_init, real rowvector iteropts, string scalar hmethod, string scalar technique) { real scalar maxtausq, itol, maxiter maxtausq = iteropts[1] itol = iteropts[2] maxiter = iteropts[3] transmorphic S S = optimize_init() optimize_init_evaluator(S, &ML_est()) optimize_init_evaluatortype(S, "d2") optimize_init_params(S, tausq_init) optimize_init_argument(S, 1, yi) optimize_init_argument(S, 2, vi) optimize_init_argument(S, 3, eff) optimize_init_technique(S, technique) optimize_init_singularHmethod(S, hmethod) optimize_init_tracelevel(S, "none") real scalar tausq_ll, rc, ll tausq_ll = optimize(S) rc = optimize_result_returncode(S) if(rc) exit(error(rc)) ll = optimize_result_value(S) if(tausq_ll < 0) ll = sum(lnnormalden(yi, eff, sqrt(vi))) return(ll - crit) } real scalar ML_skov(real scalar b, real colvector yi, real colvector vi, real colvector wi, real rowvector params, real scalar crit, real rowvector iteropts, string scalar hmethod, string scalar technique) { // unpack iteropts and params real scalar maxtausq, itol, maxiter maxtausq = iteropts[1] itol = iteropts[2] maxiter = iteropts[3] // unpack params (ML values of eff, tausq, ll) real scalar eff, tausq, ll eff = params[1] tausq = params[2] ll = params[3] // maximize LL for fixed b transmorphic S S = optimize_init() optimize_init_evaluator(S, &ML_est()) optimize_init_evaluatortype(S, "d2") optimize_init_params(S, tausq) optimize_init_argument(S, 1, yi) optimize_init_argument(S, 2, vi) optimize_init_argument(S, 3, b) optimize_init_technique(S, technique) optimize_init_singularHmethod(S, hmethod) optimize_init_tracelevel(S, "none") real scalar tausq_b, rc, ll_b tausq_b = optimize(S) rc = optimize_result_returncode(S) if(rc) exit(error(rc)) ll_b = optimize_result_value(S) if(tausq_b < 0) { tausq_b = 0 ll_b = sum(lnnormalden(yi, b, sqrt(vi))) } real colvector wi_b wi_b = 1:/(vi:+tausq_b) // (unsigned, positive) likelihood statistic at b real scalar sll if (abs(ll_b - ll) <= itol) sll = 0 // in case ll, ll_b are very close (within itol) and rounding results in a negative value else { sll = sqrt(2*(ll - ll_b)) // calculate u for Skovgaard correction (always positive) real scalar u u = U(yi, wi, wi_b, eff, b) // Improved (Skovgaard-corrected) signed likelihood statistic if (callersversion()>=16.0) { // make use of ln1p, for possible improvement in numerical stability when u very close to sll sll = sll + (1/sll)*ln1p( (sll-u) / (-sll)) } else sll = sll + (1/sll)*ln(u/sll) } sll = sign(eff - b)*sll return(sll - crit) } real scalar U(real colvector yi, real colvector wi, real colvector wi_b, real scalar eff, real scalar b) { // Expected (I) & observed (J) information, evaluated at ML estimate real matrix Imat, Jmat Imat = Jmat = J(2, 2, 0) Imat[1,1] = Jmat[1,1] = sum(wi) Imat[2,2] = .5*quadcross(wi, wi) Jmat[1,2] = Jmat[2,1] = quadcross(wi, yi:-eff, wi) Jmat[2,2] = -Imat[2,2] + quadcrossdev(yi, eff, wi:^3, yi, eff) // Observed (J) information under constraint eff = b, corresponding to tausq real scalar Jtsq Jtsq = -.5*quadcross(wi_b, wi_b) + quadcrossdev(yi, b, wi_b:^3, yi, b) // S and q real matrix S, Sinvq real colvector q S = (sum(wi_b), (eff-b)*quadcross(wi_b, wi_b) \ 0, .5*quadcross(wi_b, wi_b)) q = ((eff-b)*sum(wi_b) \ -.5*sum(wi - wi_b)) Sinvq = luinv(S)*q real scalar u u = abs(Sinvq[1,1]) * sqrt(abs(det(Jmat))) * abs(det(S)) / (sqrt(abs(Jtsq)) * abs(det(Imat))) return(u) } /* REML */ void REML(string scalar varlist, string scalar touse, real scalar hlevel, real rowvector iteropts, string scalar hmethod, string scalar technique) { // setup real colvector yi, se, vi, wi varlist = tokens(varlist) st_view(yi=., ., varlist[1], touse) if(length(yi)==0) exit(error(2000)) st_view(se=., ., varlist[2], touse) vi = se:^2 wi = 1:/vi // Initialize real scalar eff0, tausq eff0 = mean(yi, wi) // effect size with zero tausq tausq = max((0, quadvariance(yi) - mean(vi))) // Initialize tausq using Hedges estimator // Iterative tau-squared using REML transmorphic S S = optimize_init() optimize_init_evaluator(S, &REML_est()) optimize_init_evaluatortype(S, "d0") optimize_init_params(S, tausq) optimize_init_argument(S, 1, yi) optimize_init_argument(S, 2, vi) optimize_init_technique(S, technique) optimize_init_singularHmethod(S, hmethod) optimize_init_tracelevel(S, "none") real scalar rc, ll tausq = optimize(S) rc = optimize_result_returncode(S) if(rc) exit(error(rc)) ll = optimize_result_value(S) if(tausq < 0) { tausq = 0 st_numscalar("r(ll_negtsq)", ll) ll = sum(lnnormalden(yi, eff0, sqrt(vi))) - 0.5*ln(sum(wi)) } st_numscalar("r(tausq)", tausq) st_numscalar("r(converged)", optimize_result_converged(S)) st_numscalar("r(ll)", ll) // Variance of tausq (using inverse Fisher information) real scalar tsq_var tsq_var = optimize_result_V(S) st_numscalar("r(tsq_var)", tsq_var) // Confidence interval for tausq using likelihood profiling real scalar maxtausq, itol, maxiter, crit maxtausq = iteropts[1] itol = iteropts[2] maxiter = iteropts[3] crit = ll - (invchi2(1, hlevel/100)/2) real scalar tsq_lci, rc_tsq_lci, tsq_uci, rc_tsq_uci rc_tsq_lci = mm_root(tsq_lci=., &REML_profile_tausq(), 0, tausq - itol, itol, maxiter, yi, vi, crit) st_numscalar("r(tsq_lci)", tsq_lci) st_numscalar("r(rc_tsq_lci)", rc_tsq_lci) rc_tsq_uci = mm_root(tsq_uci=., &REML_profile_tausq(), tausq + itol, 10*maxtausq, itol, maxiter, yi, vi, crit) st_numscalar("r(tsq_uci)", tsq_uci) st_numscalar("r(rc_tsq_uci)", rc_tsq_uci) } void REML_est(todo, tausq, yi, vi, lnf, S, H) { real colvector wi real scalar eff tausq = max((0, tausq)) wi = 1:/(vi:+tausq) eff = mean(yi, wi) lnf = sum(lnnormalden(yi, eff, sqrt(vi:+tausq))) - 0.5*ln(sum(wi)) // Note: using d2debug reveals discrepancy in the Hessian using my formulae below // not sure if I've done it correctly, so using d0 instead (as do -metaan- and -metareg- to be fair) // if(todo>=1) { // S = -0.5*sum(wi) + 0.5*mean(wi, wi) + 0.5*quadcrossdev(yi, eff, wi:^2, yi, eff) // if(todo>=2) { // H = 0.5*quadcross(wi, wi) - mean(wi:^2, wi) + 0.5*(mean(wi, wi)^2) - quadcrossdev(yi, eff, wi:^3, yi, eff) // } // } } real scalar REML_profile_tausq(real scalar tausq, real colvector yi, real colvector vi, real scalar crit) { real colvector wi real scalar eff, ll tausq = max((0, tausq)) wi = 1:/(vi:+tausq) eff = mean(yi, wi) ll = sum(lnnormalden(yi, eff, sqrt(vi:+tausq))) - 0.5*ln(sum(wi)) return(ll - crit) } /* Confidence interval for tausq estimated using approximate Gamma distribution for Q */ /* based on paper by Biggerstaff and Tweedie (Stat Med 1997; 16: 753-768) */ // Point estimate of tausq is simply the D+L estimate void BTGamma(string scalar varlist, string scalar touse, string scalar wtvec, real scalar hlevel, real rowvector iteropts) { // Setup real colvector yi, se, vi, wi varlist = tokens(varlist) st_view(yi=., ., varlist[1], touse) if(length(yi)==0) exit(error(2000)) st_view(se=., ., varlist[2], touse) vi = se:^2 wi = 1:/vi real scalar maxtausq, itol, maxiter, quadpts maxtausq = iteropts[1] itol = iteropts[2] maxiter = iteropts[3] quadpts = iteropts[4] // Estimate variance of tausq real scalar k, eff, Q, c, d, tausq_m, tausq, Q_var, tsq_var k = length(yi) eff = mean(yi, wi) // I-V common-effect estimate Q = crossdev(yi, eff, wi, yi, eff) // standard Q heterogeneity statistic c = sum(wi) - mean(wi,wi) // c = S1 - (S2/S1) d = cross(wi,wi) - 2*mean(wi:^2,wi) + (mean(wi,wi)^2) tausq_m = (Q - (k-1))/c // untruncated D+L tausq // Variance of Q and tausq (based on untruncated tausq) Q_var = 2*(k-1) + 4*c*tausq_m + 2*d*(tausq_m^2) tsq_var = Q_var/(c^2) st_numscalar("r(tsq_var)", tsq_var) // Find confidence limits for tausq real scalar tsq_lci, rc_tsq_lci, tsq_uci, rc_tsq_uci rc_tsq_lci = mm_root(tsq_lci=., &Gamma_crit(), 0, maxtausq, itol, maxiter, tausq_m, k, c, d, .5 + hlevel/200) st_numscalar("r(tsq_lci)", tsq_lci) st_numscalar("r(rc_tsq_lci)", rc_tsq_lci) rc_tsq_uci = mm_root(tsq_uci=., &Gamma_crit(), tsq_lci + itol, 10*maxtausq, itol, maxiter, tausq_m, k, c, d, .5 - hlevel/200) st_numscalar("r(tsq_uci)", tsq_uci) st_numscalar("r(rc_tsq_uci)", rc_tsq_uci) // Find and return new weights real scalar EQ, VQ, lambda, r, se_eff EQ = (k-1) + c*tausq_m VQ = 2*(k-1) + 4*c*tausq_m + 2*d*(tausq_m^2) lambda = EQ/VQ r = lambda*EQ real colvector wsi real rowvector params real scalar i wsi = wi for(i=1; i<=k; i++) { params = (vi[i], lambda, r, c, k) wsi[i] = integrate(&BTIntgrnd(), 0, ., quadpts, params) } wi = wi*gammap(r, lambda*(k-1)) :+ wsi // update weights st_store(st_viewobs(yi), wtvec, wi) // write new weights to Stata } real scalar Gamma_crit(real scalar tausq, real scalar tausq_m, real scalar k, real scalar c, real scalar d, real scalar crit) { real scalar lambda, r, limit, ans lambda = ((k-1) + c*tausq)/(2*(k-1) + 4*c*tausq + 2*d*(tausq^2)) r = ((k-1) + c*tausq)*lambda limit = lambda*(c*tausq_m + (k-1)) ans = gammap(r, limit) - crit return(ans) } real rowvector BTIntgrnd(real rowvector t, real rowvector params) { real scalar s, lambda, r, c, k, ans s = params[1,1] // vi[i] > 0 lambda = params[1,2] // lambda = E(Q)/Var(Q) [N.B. the inverse of this is used in Henmi & Copas] r = params[1,3] // r = [E(Q)^2]/Var(Q) c = params[1,4] // c = f(weights) k = params[1,5] // k = no. studies > 1 ans = (c:/(s:+t)) :* gammaden(r, 1/lambda, 1-k, c*t) return(ans) } /* Henmi and Copas method */ // Point estimate of tausq is simply the D+L estimate void HC(string scalar varlist, string scalar touse, real scalar level, real rowvector iteropts) { // Setup real colvector yi, se, vi, wi varlist = tokens(varlist) st_view(yi=., ., varlist[1], touse) if(length(yi)==0) exit(error(2000)) st_view(se=., ., varlist[2], touse) vi = se:^2 real scalar itol, maxiter, quadpts itol = iteropts[1] maxiter = iteropts[2] quadpts = iteropts[3] real scalar k, eff, Q, W1, W2, W3, W4, tausq, VR, SDR k = length(yi) wi = 1:/vi eff = mean(yi, wi) // I-V common-effect estimate Q = crossdev(yi, eff, wi, yi, eff) // standard Q heterogeneity statistic W1 = sum(wi) W2 = mean(wi, wi) W3 = mean(wi:^2, wi) W4 = mean(wi:^3, wi) tausq = max((0, (Q - (k-1))/(W1 - W2))) // truncated D+L VR = 1 + tausq*W2 SDR = sqrt(VR) // Coefficients of 1 and (x^2) for the following functions: // EQ(x) = conditional mean of Q given R=x // VQ(x) = conditional variance of Q given R=x // finv(x) = inverse function of f(Q). // All three functions are linear combinations of 1 and (x^2), // so all can be represented by a single function, f. real scalar aEQ, bEQ aEQ = (k - 1) + tausq*(W1 - W2) - (tausq^2)*(W3 - W2^2)/VR bEQ = (W3 - W2^2)*(tausq/VR)^2 real scalar aVQ, bVQ aVQ = 2*(k - 1) + 4*tausq*(W1 - W2) + 2*(tausq^2)*(W1*W2 - 2*W3 + W2^2) aVQ = aVQ - 4*(tausq^2)*(W3 - W2^2)/VR aVQ = aVQ - 4*(tausq^3)*(W4 - 2*W2*W3 + W2^3)/VR aVQ = aVQ + 2*(tausq^4)*(1/VR^2)*(W3 - W2^2)^2 bVQ = 4*(tausq^2)*((1/VR^2))*(W3 - W2^2) bVQ = bVQ + 4*(tausq^3)*(1/VR^2)*(W4 - 2*W2*W3 + W2^3) bVQ = bVQ - 2*(tausq^4)*2*(1/VR^3)*(W3 - W2^2)^2 real scalar afinv, bfinv afinv = (k-1) - (W1/W2 - 1) bfinv = (W1/W2 - 1) real rowvector params params = (aEQ, bEQ, aVQ, bVQ, afinv, bfinv, SDR) // Find quantile of approximate distribution // (u_alpha/2 in Henmi & Copas) real scalar t, rc_t rc_t = mm_root(t=., &Eqn(), 0, 2, itol, maxiter, quadpts, level, params) if (rc_t > 0) exit(error(498)) st_numscalar("r(crit)", SDR*t) // Find test statistic (u) and p-value real scalar u, p u = eff/sqrt((tausq*W2 + 1)/W1) p = 2*integrate(&HCIntgrnd(), abs(u)/SDR, 40, quadpts, (abs(u)/SDR, params)) st_numscalar("r(p)", p) st_numscalar("r(u)", u) } // N.B. Integration should be from x to infinity, // but we only integrate up to 40 since the integrand's value is indistinguishable from zero at this point. // To see this, note that the integrand is the product of a cumulative Gamma function ==> between 0 and 1 // and a standard normal density which is indistinguishable from zero at ~40. // (thanks to Ben Jann for pointing this out) real scalar Eqn(real scalar x, real scalar quadpts, real scalar level, real rowvector params) { real scalar ans ans = integrate(&HCIntgrnd(), x, 40, quadpts, (x, params)) return(ans - (.5 - level/200)) } real rowvector HCIntgrnd(real rowvector r, real rowvector params) { real scalar t, aEQ, bEQ, aVQ, bVQ, afinv, bfinv, SDR t = params[1] aEQ = params[2] bEQ = params[3] aVQ = params[4] bVQ = params[5] afinv = params[6] bfinv = params[7] SDR = params[8] real rowvector ans ans = gammap((f(r*SDR, aEQ, bEQ):^2):/f(r*SDR, aVQ, bVQ), f(r/t, afinv, bfinv):/(f(r*SDR, aVQ, bVQ):/f(r*SDR, aEQ, bEQ))) :* normalden(r) if(t==0) ans = normalden(r) return(ans) } real rowvector f(real rowvector x, real scalar a, real scalar b) { return(a :+ b*(x:^2)) } end