/****************************************************************************** *! version 4.0.2 # Ian White # 21apr2022 skip pi option if no Sigma fix pbest problem in Stata12: variable names were lost version 4.0.1 # Ian White # 07apr2022 change invt() to invtttail() to allow running under Stata12 version 4.0 # Ian White # 07apr2022 new release to UCL, github and SSC version 3.6 # Ian White # 8feb2022 warnings fixed warning() was parsed too late in Estimate warning() is now passed to Replay, mvmeta_estimate, mvmeta_bubble, bubble references dropped to disused argument warnings mvmeta_estimate now has arguments warningest and warningnoest noposdef check respects warningest removed unnecessary qui/noi on mvmeta_estimate version 3.5.2 # Ian White # 20jan2022 improved handling of checks for positive semi-definite use _getcovcorr instead of varcheck respect psdcrit and psdcheck run mvmeta_estimate, noestimate to get checks output nicely noposdef suppresses fe estimation version 3.5.1 # Ian White # 21dec2021 bubble updates: bug fix - bubble was ignored with noestimate because Replay wasn't called add warning in replay mode that bubble uses last mvmeta run omit statement of variables chosen when there are only 2 version 3.5.0 # Ian White # 15nov2021 fixed small bug to mvmeta_lmata.ado this affected results in ~5th sig figure in models without covariates and without the longparm option new bscovariance(exchangeable) version 3.4.1 # Ian White # 15nov2021 correct penalty() option so that global is set in mvmeta_estimate and hence cleared at end version 3.4.0 # Ian White # 11nov2021 added warning in Replay if Riley correlations are larger than 0.95 new penalty() option modifies ll and rll to penalise near-singular between-studies correlation matrix need an error if other methods are used version 3.3.0 # Ian White # 14oct2021 better error message if sencode is not installed added error if wrong ciscale() added pi option Replay-only options not explicitly parsed in Estimate Replay tables in common format - title, hline, table, hline qscalar and pi added to help file move variance check after listing of variables found wt(nostudies) suppresses study-specific weights all warnings controlled by new warning() option via three routines warn* covariance variable may be missing for all studies if no study has both point estimates version 3.2.5 # Ian White # 13oct2021 improved and documented the bubble option version 3.2.4 # Ian White # 8mar2021 RELEASED TO UCL 13OCT2021 enabled coeflegend option on replay pbest: added note when probabilities are loaded into memory default is to report all ranks (new option bestonly request just the best rank) neater output version 3.2.3 # Ian White # 22jan2021 fixed bug that ignored constraints but constraints may still not be correct REML wrong with constraints on fixed parameters? needs checking in test file version 3.2.2 # Ian White # 9apr2020 fixed bug in bubble: only top half of confidence regions appeared [line required c(l), as c(L) seems to be default] version 3.2.1 # Ian White # 11mar2020 fixed bug in bubble: confidence regions didn't appear unless variables were in alphabetical order version 3.2.0 # Ian White # 6apr2018 links corrected to UCL website RELEASED TO UCL AND SSC version 3.1.4 Ian White 11feb2016 improvments in i2 option: clearer note where terms are dropped in computing CI first column is widened if necessary to fit varnames version 3.1.3 Ian White 22jul2015 observation number causing error is stored in $MVMETA_obserror (for Stephen K) check that matrix argument in bscov() is symmetric removed mvmeta_forest version 3.1.2 Ian White 20jul2015 bug fix: Estimate calls mm1 not mm2 (hence i2 now works again) mvmeta_estimate labels e(ydata) using value labels, if id is numeric and labelled i2 with mm2 estimation reported mm1 results (hence discrepant tau^2 values) - now uses mm2 results - new ad hoc method for CI (see help file) identifying errors: mvmeta_mufromsigma reports errors using correct observation number or ID mvmeta_estimate quits after finding first error improved use of exit codes: 198: syntax error 459: something that should be true of your data is not 497: my programming error 498: other statistical errors version 3.1.1 Ian White 13jul2015 bug fixes: mvmeta_wt failed with missing data & wscorr(riley) - augmentation procedure fixed deleted last call to mvmeta_wt from Estimate new tutorial on getting the data in, minor changes to main tutorial version 3.1 Ian White 1jul2015 bug fixes: e(wscorr) wasn't returned hence Riley output wasn't printed differently randfix failed when e(V_fixed) not returned - now just ignored pbest now reports in data order (previously sorted string id) forest() now works when study names include spaces riley model now switches off: testsigma qscalar i2 randfix matrix mm (=mm2) is now default; mm1 calls componentwise mm NB bubble and forest are undocumented - for a future SJ paper version 3.0.4 Ian White 16jun2015 new fe suboption for suppress() fixed mvmeta_mufromsigma to respect riley method - was previously giving wrong answers for likelihood and fixed parameters - wscorr(riley) fixed now fails; previously it behaved like wscorr(0) fixed version 3.0.3 Ian White 12jun2015 bug fix: iterate() wasn't working version 3.0.2 Ian White 9jun2015 on external website problems when univariate or unstructured MM failed now captures problems sensibly better (shorter) warnings given no left-over matrices version 3.0.1 Ian White 8jun2015 on external website 25may2015 fixed bug in pbest, saving() 14may2015 bug fix in mvmeta_mufromsigma: bscov(equals) now returns correct [r]loglik new meanrank suboption in pbest (and mcse corrected?) v3.0.1 11may2015 updated mvmeta.pkg bug fixes: e(Q), e(Qa), e(Qb) always returned from mm method i2 option: now works for mm2 (same as mm) testsigma now works again (e(ll0) restored) ml options including iter now work options dropped: keepmat() corr() problems arising from uv estimation: new options: suppress([uv] [mm]) uviterate(default 20) commonparm => suppress(uv) REML fails with #obs = #parms, so (a) checks and issues warning, (b) captures uv failure mvmeta_bos renamed mvmeta_wt mvmeta_{estimate bubble wt forest} and bubble brought into this prog this version seems correct & stable v3.0.0 25apr2015 programming changes: ynames no longer passed as a global - use ylist new e() matrices: ydata, Sdata, Xdata, b_fixed, V_fixed, b_uv, V_uv, wt 24apr2015: ydata row names are not :estimate bos, bubble, randfix, forest - all updated functionality changes: no id() option on Replay, only on pbest - NB pbest is prediction type (can be out of sample) BUGS -mvmeta, i2- fails when estimated variance is 0 after -cov(propto M)- [e.g. H:\meta\incoherence\thrombolytics\mvmeta4.do] TO DO e(ll) etc. is inconsistent: reml & ml return e(ll) = maximum [restricted] loglik and e(ll0) = same for fixed-effect model mm returns nothing fixed returns e(ll0) = maximum loglik and e(rl0) = maximum restricted loglik Forest - still needs a lot more work weights on by default (Dan) default should either add univariate weights or suppress uv summary (Dan) how is prediction interval calculated? favours() - needs to be subgraph specific (for Berkey) nouv seems to leave a gap don't display those flagged as augmented Bubble Could plot against an EB estimate of the missing data??? PROGRAMMING NOTES This program uses method d0, not lf, because the REML log-likelihood does not have the required form as a sum of individual contributions. So I have to handle covariates through the program. -- could use lf for ML only The estimate and variance matrices must be numbered 1..`n' obs loops over 1/`N' (all obs); i loops over 1/`n' (obs in use); r, s loop over 1/`p' (outcome variables). Programs are defined in this order: mvmeta Estimate Replay pbest drawbeta mvmeta_wt mvmeta_bubble bubble [mvmeta_forest] mvmeta_estimate varcheck mvmeta_getdata mata selectpart() WARNINGS mm method doesn't cope with augmented data sets (i.e. with large variances) CERTIFICATION SCRIPTS H:\meta\mv\ado\mvmeta\mvmeta_cscript.do H:\meta\mv\regpaper\cscript.do - analyses as in SJ paper OPTIONS REMOVED report FUTURE PLANS lincom option? Kenward-Roger df for tests and \cint s? Compute empirical Bayes estimates? i2(white|cochran|jackson)?? ******************************************************************************/ *=============================== MASTER PROGRAM =============================== program define mvmeta version 9 // originally written under Stata 9, but updates are only tested under Stata 12+ if replay() { if "`e(cmd)'" != "mvmeta" error 301 Replay `0' } else Estimate `0' end *============================== ESTIMATE PROGRAM ============================== program define Estimate, eclass syntax [anything] [if] [in], /// [reml ml Fixed MM2 mm1 notrunc Vars(varlist) start(passthru) WSCORr(string) LONGParm EQuations(string) noCONStant BSCOVariance(string) COMMONparm /// Model options SUPPress(string) /// additional analysis options augment augquiet missest(real 0) missvar(real 1E4) /// Augmentation options - NOW UNDOCUMENTED maxse(real 5) psdcrit(passthru) noposdef mata no2pi TAULog /// Tuning options SHOWStart id(varname) WARNing(string) /// Output options for Estimate cholnames nolog quick debug noESTimates noJOINTCHeck mmfix constraints(passthru) COLLinear network(string) wscorrforce noWARNings UVITERate(int 20) ITERate(string) penalty(real 0) /// Undocumented Estimate options * /// ml and Replay options ] global MVMETA_obserror // new 21jul2015 ****************** PARSE ***************** mlopts mlopts replayopts, `options' if "`warning'"=="" local warning error if !inlist("`warning'", "error", "text", "off") { di as error "syntax: warning(error|text|off)" exit 198 } local warn warn`warning' * OUTCOMES AND VARIANCES tokenize "`anything'" local ystub `1' local Sstub `2' if "`Sstub'"=="" { di as error "Syntax: mvmeta estimate-stub variance-stub ..." exit 198 } * COVARIATES (AND CONSTANT) mac shift 2 if "`collinear'"=="" { _rmcoll `*', nocons local xvars `r(varlist)' } else local xvars `*' local isregression = "`xvars'"!="" * EQUATIONS tokenize "`equations'", parse(",") local neqs 0 while "`1'"!="" { gettoken eqy rest : 1, parse(":") gettoken colon eqx : rest, parse(":") if "`eqx'"!="" local isregression 1 foreach thisy of varlist `eqy' { local ++neqs local eqy`neqs' `thisy' if "`eqx'"!="" unab eqx`neqs' : `eqx' } mac shift 2 } * ESTIMATION METHOD if wordcount("`reml' `ml' `fixed' `mm1' `mm2'") > 1 { di as error "Please specify only one of reml, ml, fixed, mm1, mm2" exit 198 } local bsest `reml'`ml'`fixed'`mm1'`mm2' if "`bsest'"== "" local bsest reml * CHECK CORRELATIONS if "`wscorr'"=="riley" { if inlist("`bsest'","mm1","mm2") { di as error "Sorry, we don't yet have a method of moments for Riley's overall correlation model." exit 498 } local suppress `suppress' fe mm } else if "`wscorr'"!="" { cap assert abs(`wscorr')<=1 if _rc { di as error "Correlation (`wscorr') outside range [-1,1]." exit 459 } } else if !mi("`wscorrforce'") { `warn' "wscorrforce ignored - wscorr(#) not specified" local wscorrforce } * SHORT, LONG OR COMMON PARAMETERISATION if "`commonparm'"=="commonparm" local parmtype common else if "`longparm'"=="longparm" | "`xvars'"!="" | "`equations'"!="" | "`pbest'"!="" local parmtype long else local parmtype short if "`parmtype'"=="common" { local collinear collinear local suppress `suppress' uv } * notrunc not allowed without method mm -> move if "`trunc'"=="notrunc" & !inlist("`bsest'","mm1","mm2","fixed") { `warn' "notrunc option ignored with method `bsest'" local trunc } if "`debug'"=="" local ifdebug * local ifdebugnoi = cond("`debug'"=="","qui","noi") di as text "Note: using method " as result "`bsest'" if !mi("`id'") local idopt id(`id') foreach word of local suppress { if "`word'"=="uv" local suppressuv suppressuv if "`word'"=="mm" local suppressmm suppressmm if "`word'"=="fe" local suppressfe suppressfe } if !mi("`posdef'") local suppressfe suppressfe // 20jan2022 if `penalty'<0 { di as error "penalty() must be non-negative" exit 198 } if `penalty' di as text "Note: log likelihood enhanced by penalty term " as res "`penalty' * log(det(corr(Sigma)))" ****************** END OF PARSING ***************** ****************** IDENTIFY Y AND S VARIABLES TO USE, AND p ***************** if "`vars'"~="" local ylist0 `vars' /* change by SK to accomodate long variable list */ else local ylist0 `ystub'* local p 0 cap desc `ylist0' if _rc { di as error "No variables found starting with `ystub'" exit 459 } foreach yvar of varlist `ylist0' { if index("`yvar'","`ystub'")!=1 { di as error "Variable `yvar' in vars() does not start with `ystub'" exit 198 } local yend = substr("`yvar'",1+length("`ystub'"),.) if "`yend'"=="" { `warn' "variable `yvar' not used (looking for variable names `ystub'+suffix)" continue } summ `yvar' `if' `in', meanonly if r(N)==0 { `warn' "variable `yvar' not used (no non-missing values)" continue } * Found a valid yvar local ++p local yend_`p' `yend' local yvar_`p' `ystub'`yend' local yends `yends' `yend' local ylist `ylist' `ystub'`yend' cap confirm var `Sstub'`yend'`yend' if _rc { di as error "Variance `Sstub'`yend'`yend' not found for estimate `ystub'`yend'" exit 459 } * code not needed because zero variances are checked for later (and returned with id) * cap assert `Sstub'`yend'`yend'>0 * if _rc { * di as error "Variance `Sstub'`yend'`yend' must be >0" * exit 459 * } * Now see if we have any covariates for it local xvars_`p' `xvars' forvalues eq=1/`neqs' { if "`yvar_`p''"=="`eqy`eq''" local xvars_`p' `xvars_`p'' `eqx`eq'' } if `isregression' { di as text "Note: regressing " as result "`yvar_`p''" as text " on " _c if !mi("`xvars_`p''") di as result "`xvars_`p''" else di as result "(nothing)" } if "`collinear'"=="" & !mi("`xvars_`p''") { qui count if !mi(`yvar_`p'') if r(N)<=1 { // bug in _rmcoll - it fails with only 1 obs if mi("`constant'") local newxvars else local newxvars : word 1 of "`xvars_`p'" } else { _rmcoll `xvars_`p'' if !mi(`yvar_`p''), `constant' local newxvars = r(varlist) if "`newxvars'" == "." local newxvars } if "`xvars_`p''" != "`newxvars'" { local newxvarsdisp = cond(!mi("`newxvars'"), "`newxvars'", "(nothing)") `warn' "Collinearity detected: now regressing `yvar_`p'' on `newxvarsdisp'" local xvars_`p' `newxvars' } } if "`xvars_`p''"=="" & "`constant'"=="noconstant" { di as error "No covariates and no constant for outcome `p'" exit 498 } * find yvar label local yvarlab : var label `yvar_`p'' if mi("`yvarlab'") local yvarlab `yvar_`p'' if mi(`"`ylabels'"') local ylabels `"`"`yvarlab'"'"' else local ylabels `"`ylabels' `"`yvarlab'"'"' } if `p'==0 { di as error "No variables found starting with `ystub'" exit 111 } if !`isregression' { if `p'>1 local plural s di as text "Note: using variable`plural' " as result "`ylist'" } forvalues r=1/`p' { cap confirm var `Sstub'`yend_`r''`yend_`r'' if _rc { di as error "Variance `Sstub'`yend_`r''`yend_`r'' not found for estimate `ystub'`yend_`r''" exit 459 } } **************** SORT OUT COVARIANCE STRUCTURES ******************* * (moved from end of parsing because it needs `p') if "`bsest'"=="`fixed'" { if !mi("`bscovariance'") { di as error "fixed option used: bscovariance() ignored" local bscovariance } local jointcheck none } else { tempname sigma0 if "`bscovariance'"=="" local bscovariance unstructured local sigma0exp = word("`bscovariance'",2) if substr("`bscovariance'",1,3)=="uns" { local bscovariance unstructured } else if substr("`bscovariance'",1,4)=="prop" { if "`sigma0exp'"=="" { di as error "Syntax: bscovariance(proportional matrix)" exit 198 } mat `sigma0' = `sigma0exp' local bscovariance proportional } else if substr("`bscovariance'",1,2)=="eq" { if "`sigma0exp'"=="" { di as error "Syntax: bscovariance(equals matrix)" exit 198 } mat `sigma0' = `sigma0exp' local bscovariance equals } else if substr("`bscovariance'",1,4)=="corr" { if "`sigma0exp'"=="" { di as error "Syntax: bscovariance(correlation matrix)" exit 198 } mat `sigma0' = `sigma0exp' local bscovariance correlation } else if substr("`bscovariance'",1,4)=="exch" { local rho = word("`bscovariance'",2) if mi("`rho'") { * correlation is to be estimated, so this really is bscov(exch) local bscovariance exchangeable } else { * correlation is given, so this is shorthand for bscov(prop ..) cap confirm number `rho' if _rc { di as error "Syntax: bscovariance(exch #)" exit 198 } if `rho'<-1 | `rho'>1 { di as error "Syntax: bscovariance(exch #) with -1<=#<=1" exit 198 } local bscovariance proportional mat `sigma0' = (1-`rho')*I(`p')+`rho'*J(`p',`p',1) local sigma0exp `=1-`rho''*I(`p')+`=`rho''*J(`p',`p',1) } } else { di as error "bscovariance(`bscovariance') not available" exit 198 } if "`bscovariance'" != "unstructured" local jointcheck none } * check sigma0 matrices if !mi("`sigma0exp'") { if rowsof(`sigma0')!=`p' | colsof(`sigma0')!=`p' { di as error "bscovariance() option: `sigma0exp' is not `p'x`p'" exit 198 } if !issymmetric(`sigma0') { di as error "bscovariance() option: `sigma0exp' is not symmetrical" exit 198 } } ****************** ESTIMATION SAMPLE ***************** marksample touse * 9dec2014 - only mark out if x=. and y!=. tempvar hasyvalues hasxnoy gen `hasyvalues'=0 // indicates any non-missing y-value gen `hasxnoy' = 0 // indicates any missing x-value with non-missing y-value forvalues r=1/`p' { foreach xvar in `xvars_`r'' { qui replace `hasxnoy' = 1 if mi(`xvar') & !mi(`yvar_`r'') } qui replace `hasyvalues' = 1 if !mi(`yvar_`r'') } qui count if `hasyvalues'==0 & `touse' if r(N) { local observations = cond(r(N)>1,"observations","observation") local have = cond(r(N)>1,"have","has") `warn' r(N) " `observations' `have' all outcomes missing, and `have' been dropped from analysis" qui replace `touse' = 0 if `hasyvalues'==0 } qui count if `hasxnoy'==1 & `touse' if r(N) { local observations = cond(r(N)>1,"observations","observation") local have = cond(r(N)>1,"have","has") `warn' r(N) " `observations' `have' observed outcome and missing covariate, and `have' been dropped from analysis" qui replace `touse' = 0 if `hasxnoy'==1 } qui count if `touse' local n = r(N) local N = _N di as text "Note: " as result `n' as text " observations on " as result `p' as text " variables" ****************** CHECK VARIANCES ***************** if "`bscovariance'"=="unstructured" { forvalues r=1/`p' { forvalues s = `=`r'+1' / `p' { qui count if !missing(`yvar_`r'') & !missing(`yvar_`s'') & `touse' if r(N)<=1 { if r(N)==0 local problem have no jointly observed values if r(N)==1 local problem have only 1 jointly observed value if "`jointcheck'"=="" { di as error "`yvar_`r'' and `yvar_`s'' `problem' - consider alternatives to bscov(unstructured)" exit 459 } else `warn' "`yvar_`r'' and `yvar_`s'' `problem' - consider alternatives to bscov(unstructured)" } } } } **************** RUN SECONDARY ESTIMATIONS *************** ereturn clear // Needed to drop any old e-results (since we don't always do -ereturn post-) local optsforall `constant' parmtype(`parmtype') `augment' missest(`missest') missvar(`missvar') `augquiet' `taulog' `debug' `trunc' `mmfix' `idopt' `mlopts' `psdcrit' `posdef' warningest(`warning') // added psdcrit posdef 14jan2022; added warningest 8deb2022 forvalues r=1/`p' { if `r'>1 local xvarslist `xvarslist', if mi("`xvars_`r''") local xvarslist `xvarslist' . // missing is awkward else local xvarslist `xvarslist' `xvars_`r'' local nxvars_`r' = wordcount(`"`xvars_`r''"') + ("`constant'"!="noconstant") } * first run it with noestimates to check var-cov matrices and return `ifdebug' di as input _new "Checking data" mvmeta_estimate `ystub' `Sstub' if `touse', yends(`yends') xvarslist(`xvarslist') wscorr(`wscorr') `wscorrforce' `optsforall' noestimates warningnoest(`warning') if "`estimates'"!="noestimates" { ******** RUN FIXED-EFFECT ESTIMATION ******* if "`suppressfe'"!="suppressfe" { `ifdebug' di as input _new "Running fixed-effect estimation" `ifdebugnoi' mvmeta_estimate `ystub' `Sstub' if `touse', yends(`yends') xvarslist(`xvarslist') wscorr(`wscorr') `wscorrforce' bsest(fixed) `optsforall' warningnoest(off) if r(success) { tempname b_fixed V_fixed mat `b_fixed' = r(b) mat `V_fixed' = r(V) local Qscalar_chi2 = r(Qscalar) local Qscalar_df = r(Qscalardf) `ifdebug' di as input "Results of fixed-effect estimation:" _c `ifdebug' mat l `b_fixed', title(b_fixed) `ifdebug' mat l `V_fixed', title(V_fixed) if "`bsest'"=="reml" local ll0 = r(rl0) // [re]ll for FE model - needed for testsigma option else if "`bsest'"=="ml" local ll0 = r(ll0) } else `warn' "Fixed-effect estimation failed" } else `ifdebug' di as input _new "Estimation of fixed-effect model suppressed" ******** RUN UNIVARIATE ESTIMATIONS ******* if "`suppressuv'"!="suppressuv" { tempname b_uv V_uv thisb thisV zeroV forvalues r=1/`p' { `ifdebug' di as input _new "Running univariate estimation for outcome `yvar_`r''" qui count if `touse' & !mi(`yvar_`r'') if r(N)==`nxvars_`r'' & "`bsest'"=="reml" & !mi("`debug'") `warn' "univariate REML estimation may fail for outcome `yvar_`r'', since #obs= #parms" if "`bscovariance'"=="equals" { local sigma0uv = `sigma0'[`r',`r'] local bscovopt bscovariance(equals) sigma0(`sigma0uv') } else local bscovopt bscovariance(unstructured) `ifdebugnoi' mvmeta_estimate `ystub' `Sstub' if `touse' & !mi(`yvar_`r''), yends(`yend_`r'') xvarslist(`xvars_`r'') bsest(`bsest') `bscovopt' `cholnames' `showstart' `optsforall' iterate(`uviterate') warningnoest(off) if r(success) { mat `thisb' = r(b) mat `thisb' = `thisb'[1,1..`nxvars_`r''] mat `thisV' = r(V) mat `thisV' = `thisV'[1..`nxvars_`r'',1..`nxvars_`r''] } else { local failedvars `failedvars' `yvar_`r'' mat `thisb' = J(1,`nxvars_`r'',.) mat `thisV' = J(`nxvars_`r'',`nxvars_`r'',.) * name rows and cols local colnames foreach thing in `xvars_`r'' { local colnames `"`colnames' "`thing'""' } if mi("`constant'") local colnames `"`colnames' "_cons""' mat colnames `thisb' = `colnames' mat colnames `thisV' = `colnames' mat rownames `thisV' = `colnames' mat coleq `thisb' = "`yvar_`r''" mat coleq `thisV' = "`yvar_`r''" mat roweq `thisV' = "`yvar_`r''" } mat `b_uv' = nullmat(`b_uv'),`thisb' if `r'==1 mat `V_uv' = `thisV' else { mat `zeroV' = J(rowsof(`V_uv'),colsof(`thisV'),0) mat rownames `zeroV' = `: rowfullnames `V_uv'' mat colnames `zeroV' = `: colfullnames `thisV'' mat `V_uv' = (`V_uv', `zeroV' \ `zeroV'', `thisV') } } if !mi("`failedvars'") `warn' "univariate estimation failed for outcomes " as result "`failedvars'" `ifdebug' di as input "Results of univariate estimations:" _c `ifdebug' mat l `b_uv', title(b_uv) `ifdebug' mat l `V_uv', title(V_uv) } else `ifdebug' di as input _new "Estimation of univariate models suppressed" ******** RUN UNSTRUCTURED-SIGMA MM FOR ALL METHODS, BECAUSE I^2 NEEDS IT ******* if "`suppressmm'"!="suppressmm" { `ifdebug' di as input _new "Running unstructured-sigma mm" `ifdebugnoi' mvmeta_estimate `ystub' `Sstub' if `touse', yends(`yends') xvarslist(`xvarslist') wscorr(`wscorr') `wscorrforce' bsest(mm1) bscovariance(unstructured) `start' `cholnames' `showstart' `optsforall' warningnoest(off) if r(success) { tempname Q Qa Qb mat `Q' = r(Q) mat `Qa' = r(Qa) mat `Qb' = r(Qb) if inlist("`bsest'","ml","reml") & "`bscovariance'"=="unstructured" & mi("`start'") { tempname mmSigma mat `mmSigma' = r(Sigma) local start start(`mmSigma') } } else { if inlist("`bsest'","reml","ml") { `warn' "unstructured-variance method of moments failed - I2 statistic is not available" } if "`bscovariance'"=="unstructured" & inlist("`bsest'","mm2") { di as error "Error: method of moments failed" exit 459 } } } else `ifdebug' di as input _new "Estimation of unstructured-sigma method of moments suppressed" } else { di as text "Note: " as result "noestimates" as text " option - model not fitted" local bsest } **************** RUN MAIN ESTIMATION *************** // - must be last, as mvmeta_estimate implies -ereturn clear- // - this is run even with noestimates option, since it also writes data matrices `ifdebug' di as input _new "Running main estimation" mvmeta_estimate `ystub' `Sstub' if `touse', yends(`yends') xvarslist(`xvarslist') wscorr(`wscorr') `wscorrforce' bsest(`bsest') bscovariance(`bscovariance') sigma0(`sigma0') `start' `cholnames' `showstart' `optsforall' `estimates' sigma0exp(`sigma0exp') iterate(`iterate') `constraints' warningnoest(off) penalty(`penalty') // ereturn the main answers `ifdebug' di as input "Ereturning the main answers" foreach scalar in `r(scalarlist)' { ereturn scalar `scalar' = r(`scalar') } tempname temp foreach matrix in `r(matrixlist)' { if inlist("`matrix'","b","V") continue // b and V have already been ereturned by ml mat `temp' = r(`matrix') ereturn matrix `matrix' = `temp' } foreach local in `r(locallist)' { ereturn local `local' = r(`local') } ************ ERETURN DESCRIPTIONS AND TIDY UP *************** foreach thing in parmtype bsest constant cholnames wscorr wscorrforce network id ystub Sstub ylabels { ereturn local `thing' `"``thing''"' } forvalues r=1/`p' { ereturn local xvars_`r' `xvars_`r'' } ereturn local yvars "`ylist'" if "`bsest'"=="fixed" ereturn local bscovariance "(none)" else if "`sigma0exp'"=="" ereturn local bscovariance "`bscovariance'" else ereturn local bscovariance "`bscovariance' `sigma0exp'" ereturn local cmdline mvmeta `0' foreach thing in ll0 Qscalar_chi2 Qscalar_df { if !mi("``thing''") ereturn scalar `thing' = ``thing'' } ereturn scalar dims = `p' if "`estimates'"=="noestimates" { // some things have been missed ereturn scalar N = `N' } foreach thing in Q Qa Qb b_fixed V_fixed b_uv V_uv _wt { if !mi("``thing''") ereturn matrix `thing' = ``thing'' } // tidy up ereturn local cmd "mvmeta" // call replay if estimation has happened, or if bubble is called // 21dec2021 local 0 ,`replayopts' syntax, [bubble bubble2(string) *] if !mi("`bubble'`bubble2'") | "`estimates'"!="noestimates" Replay, `estimates' `debug' `replayopts' hasestimated warning(`warning') end *========================== END OF ESTIMATE PROGRAM =========================== *========================== START OF REPLAY PROGRAM =========================== program define Replay // PARSE syntax, [SHOWChol SHOWAll eform EFORM2(string) noUNCertainv print(string) Level(cilevel) dof(string) i2 i2fmt(string) ciscale(string) NCchi2 noESTimates TESTsigma pbest(string) RANDFix RANDFix2(varlist) Qscalar debug cformat(passthru) pformat(passthru) sformat(passthru) WT WT2(string) FORest FORest2(string asis) BUbble BUbble2(string asis) COEFLegend id(varname) PI PI2(string) hasestimated warning(string)] if !mi("`debug'") di "Starting Replay program" if !mi("`id'") { di as error "id() option not allowed in replay mode" exit 198 } if "`warning'"=="" local warning error if !inlist("`warning'", "error", "text", "off") { di as error "syntax: warning(error|text|off)" exit 198 } local warn warn`warning' foreach bit in `print' { if "`bit'"=="bscov" local printbscov on else if "`bit'"=="bscorr" local printbscorr on else di as error "option `bit' ignored in print(`print')" } if !mi("`eform2'") local eform eform(`eform2') else if !mi("`eform'") { if "`showall'"=="" local eform eform(exp(Coef)) else local eform eform([exp](Coef)) } // eform2 can now be ignored if index("`dof'","n") { local dof = subinstr("`dof'","n","e(N)",.) local dof = `dof' } if "`dof'"!="" local dofopt dof(`dof') local p = e(dims) local bsest = "`e(bsest)'" local bscov = word("`e(bscovariance)'",1) local neqs_mean = e(neqs_mean) local neq neq(`neqs_mean') /*2013-03-06*/ if "`uncertainv'" == "nouncertainv" & inlist("`bsest'","fixed","mm1","mm2") { di as text "(option nouncertainv ignored with estimation method `bsest')" local uncertainv } if "`showall'"=="showall" { if "`uncertainv'" == "nouncertainv" { di as text "(option showall ignored with option nouncertainv)" local showall } else if e(neqs_aux)==0 { di as text "(option showall ignored - no auxiliary parameters)" local showall } else { local neqs_tot = e(neqs_mean) + e(neqs_aux) local neq neq(`neqs_tot') } } if "`e(parmtype)'"=="long" & "`eform'"!="" & "`showall'"=="showall" { di as text "(option eform ignored with option showall and parmtype long)" local eform } forvalues r=1/`p' { local yvar_`r' : word `r' of `e(yvars)' } local levelopt level(`level') tempname b V V2 hold if "`ciscale'"=="" local ciscale sd local ciscale = lower("`ciscale'") if "`ciscale'"=="sd" local ciscalename SD else if "`ciscale'"=="logsd" local ciscalename log(SD) else if "`ciscale'"=="logh" local ciscalename log(H) else { di as error "Syntax: ciscale(sd|logsd|logh)" exit 198 } if "`i2fmt'"=="" local i2fmt %6.0f if "`debug'"=="" local ifdebug * if "`showchol'"=="showchol" { `warn' "showchol option has been renamed showall" local showall showall } if mi(e(bsest)) { // if nothing to display, don't display anything local estimates noestimates local wt } if "`e(wscorr)'"=="riley" { foreach opt in i2 qscalar testsigma randfix { if !mi("``opt''") di as error "Option `opt' is not available with wscorr(riley)" local `opt' } } // END OF PARSING // PRINT TABLE OF RESULTS if "`estimates'"!="noestimates" { // PRINT OUTPUT HEADER di as text _newline "Multivariate meta-analysis" di as text "Variance-covariance matrix = " as result e(bscovariance) di as text "Method = " as result "`bsest'" _c di _col(48) as text "Number of dimensions " _col(72) "=" as result %6.0f e(dims) if ("`bsest'"=="ml" ) di as text "Log likelihood = " as result e(ll) _c if ("`bsest'"=="reml" ) di as text "Restricted log likelihood = " as result e(ll) _c di _col(48) as text "Number of observations " _col(72) "=" %6.0f as result e(N) if "`dof'"!="" di _col(48) as text "Degrees of freedom " _col(72) "=" _col(72) %6.0f as result `dof' // PRINT ESTIMATES mat `b' = e(b) mat `V' = e(V) if ("`bsest'"=="ml" | "`bsest'"=="reml") & "`uncertainv'" == "nouncertainv" { mat `b' = `b'[1,1..`e(nparms_mean)'] mat `V2' = invsym(`V') mat `V2' = `V2'[1..`e(nparms_mean)',1..`e(nparms_mean)'] mat `V2' = invsym(`V2') _estimates hold `hold' ereturn post `b' `V2', `dofopt' *ereturn scalar k_eform = `neqs_mean' ereturn display, `eform' `levelopt' `cformat' `pformat' `sformat' // _coef_table fails here, not sure why di as text "Note: these standard errors ignore uncertainty in Sigma." _estimates unhold `hold' /* brings back the main results */ } else { _estimates hold `hold', copy ereturn post `b' `V', `dofopt' *ereturn scalar k_eform = `neqs_mean' if "`eform'"!="" & "`showall'"=="showall" { // here, _coef_table is better than ereturn display because it only exponentiates the first e(k_eform) equations - but older versions don't work with e(cmd) unset cap _coef_table, `eform' `levelopt' `cformat' `pformat' `sformat' `neq' if _rc { di as error "_coef_table failed - you may be using a version older than 3.0.0" di as error "All parameters will be exponentiated" ereturn display, `eform' `levelopt' `cformat' `pformat' `sformat' `neq' `coeflegend' } else { _coef_table, `eform' `levelopt' `cformat' `pformat' `sformat' `neq' `coeflegend' } di as text "Variance parameters are not exponentiated" } else ereturn display, `eform' `levelopt' `cformat' `pformat' `sformat' `neq' `coeflegend' _estimates unhold `hold' /* brings back the main results */ } if "`e(parmtype)'"=="common" & `p'>1 { di as text "The above coefficients also apply to the following equations:" if "`e(constant)'"!="noconstant" local constantvar _cons forvalues r=2/`p' { di as result _col(5) "`yvar_`r''" as text ": `e(xvars_`r')' `constantvar'" } } } // PRINT SIGMA if ("`estimates'"!="noestimates" | "`printbscov'`printbscorr'"!="") & "`bsest'"!="fixed" { if "`printbscov'`printbscorr'"=="" local printbscorr on if !mi("`printbscorr'") { tempname Corr SD all cap mat `Corr'=corr(e(Sigma)) local bsprob = (_rc>0) cap mat `SD'=vecdiag(cholesky(diag(vecdiag(e(Sigma)))))' if _rc local bsprob 1 else mat colnames `SD'="SD" if `bsprob' { `warn' "can't convert between-studies covariance matrix to correlation matrix" _c * di " (probably because one or more variances are zero or negative)" _c local printbscov on local printbscorr } } if "`printbscov'"=="on" { di as text _newline "Estimated between-studies covariance matrix Sigma" local linelength 49 di as txt "{hline `linelength'}" _c mat l e(Sigma), noheader } if "`printbscorr'"=="on" { mat `all'=(`SD', `Corr') * set above-diag elements to missing forvalues r=1/`p' { local s0=`r'+2 local s1=`p'+1 forvalues s=`s0'/`s1' { mat `all'[`r',`s']=. } } if "`e(wscorr)'"=="riley" { di as text _newline "Estimated between-studies SDs and OVERALL correlation matrix" local linelength 60 } else { if word(e(bscovariance),1)=="equals" di as text _newline "Specified" _c // 20jan2022 else di as text _newline "Estimated" _c di as text " between-studies SDs and correlation matrix" local linelength 52 } di as txt "{hline `linelength'}" _c mat l `all', noheader } if "`e(wscorr)'"=="riley" { // warnings re Riley tempname corr mat `corr' = corr(e(Sigma)) forvalues r=1/`p' { forvalues s=`r'/`p' { if `s'<=`r' continue if `corr'[`r',`s']>0.95 di as error "Warning: overall Riley correlation is near 1: estimates may be unstable" if `corr'[`r',`s']<-0.95 di as error "Warning: overall Riley correlation is near -1: estimates may be unstable" } } } di as txt "{hline `linelength'}" } // PRINT WEIGHTS AND BOS if !mi("`wt'`wt2'") { if "`e(parmtype)'"=="common" di as error "Option wt() ignored - incompatible with commonparm" else mvmeta_wt, `wt2' `debug' } // PRINT I-SQUARED if "`i2'"=="i2" { local error 0 foreach mat in Q Qa Qb { cap confirm matrix e(`mat') local error = cond(_rc,1,`error') } if `error' di as error "Method of moments failed when mvmeta ran, so I^2 cannot be computed" else { // SET UP local z = invnorm((100+`level')/200) foreach mat in Q Qa Qb Sigma { tempname `mat' mat ``mat'' = e(`mat') } // COMPUTE TAU AND I^2, WITH CIs IF POSSIBLE forvalues r=1/`p' { local s2`r' = `Qa'[`r',`r']/`Qb'[`r',`r'] if `s2`r''==. local s2problem yes if inlist("`bsest'", "fixed", "mm1", "mm2") { local thisQa = `Qa'[`r',`r'] if "`bsest'"=="mm2" { // ad hoc method 16jul2015 local thisQ = `Qa'[`r',`r'] + `Qb'[`r',`r'] * `Sigma'[`r',`r'] local reconstructed `"reconstructed "' } else local thisQ = `Q'[`r',`r'] cap heterogi `thisQ' `thisQa', `levelopt' `ncchi2' if _rc { local i2est`r' = max(0,100*(1 - `thisQa'/`thisQ')) local i2low`r' . local i2upp`r' . local heterogifailed yes } else { local i2est`r' = 100*r(I2) local meth = cond("`ncchi2'"=="ncchi2",2,1) local i2low`r' = 100*r(lb_I2_M`meth') local i2upp`r' = 100*r(ub_I2_M`meth') } // COMPUTE TAU foreach thing in est`r' low`r' upp`r' { local bssd`thing' = sqrt( `i2`thing'' * `s2`r'' / (100 - `i2`thing'') ) } local i2footnote "Note: I^2 computed from `reconstructed'Q matrix and its degrees of freedom" if "`heterogifailed'"=="yes" { cap heterogi if _rc==199 local heterogineeded yes } } else { // REML/ML: COMPUTE CIs FOR BETWEEN-STUDY SDs if "`bscov'"=="proportional" { // -nlcom- only works if you express sqrt(a*tau^2) as sqrt(a)*tau tokenize "`e(Sigma`r'`r')'", parse("*") local number `1' confirm number `1' assert "`3'" == "([tau]_b[_cons])^2" local tau sqrt(`number')*[tau]_b[_cons] local tausq `e(Sigma`r'`r')' } else if "`bscov'"=="unstructured" { local tausq `e(Sigma`r'`r')' // DROP ANY ZERO TERMS cap nlcom sqrt(`tausq') if _rc { tokenize "`tausq'", parse("+") local tausqnew local tausqdrop while "`1'"!="" { cap nlcom `1' if !_rc { if "`tausqnew'"!="" local tausqnew `tausqnew' + local tausqnew `tausqnew' `1' } else { if "`tausqdrop'"!="" local tausqdrop `tausqdrop' + local tausqdrop `tausqdrop' `1' } mac shift 2 } `ifdebug' di as text "Debug: term " as result "`tausqdrop'" as text " dropped in calculating i2 for `yvar_`r''" if abs((`tausq')/(`tausqnew')-1) < 1E-5 { * it's ok local tausq `tausqnew' local droppedterms yes local droppedterm`r' " *" } else di as error "nlcom failed and couldn't be fixed" } local tau sqrt(`tausq') } else local tau sqrt(`e(Sigma`r'`r')') local tau`r' `tau' // for correlations later // HANDLE SCALES (PART 1) if "`ciscale'"=="logsd" local transform log(`tau') else if "`ciscale'"=="logh" local transform 0.5*log(1+(`tausq')/`s2`r'') else if "`ciscale'"=="sd" local transform `tau' cap nlcom `transform' if _rc | ("`bscov'"=="equals") { local bssdest`r' = `tau' local bssdlow`r' . local bssdupp`r' . local nlcomfails yes if ("`bscov'"=="equals") local nlcomfails equals } else { mat `b' = r(b) mat `V' = r(V) local bssdest`r' = `b'[1,1] local bssdlow`r' = `b'[1,1] - `z' * sqrt(`V'[1,1]) local bssdupp`r' = `b'[1,1] + `z' * sqrt(`V'[1,1]) foreach thing in bssdest`r' bssdlow`r' bssdupp`r' { // HANDLE SCALES (PART 2) if "`ciscale'"=="logsd" local `thing' = exp(``thing'') else if "`ciscale'"=="logh" { if ``thing''<0 local `thing' 0 else local `thing' = sqrt(`s2`r''*(exp(2*``thing'')-1)) } else if "`ciscale'"=="sd" { if ``thing''<0 local `thing' 0 else local `thing' = ``thing'' } } } // COMPUTE I^2 foreach thing in est`r' low`r' upp`r' { local i2`thing' = 100*(`bssd`thing'')^2 / ((`bssd`thing'')^2 + `s2`r'') } local cifootnote "Note: CIs computed on `ciscalename' scale" local i2footnote "Note: I^2 computed from estimated between-studies variance" /// _new _skip(6) "and typical within-studies variances" } } * FOR OUTPUT local ylength 0 forvalues r=1/`p' { local ylength = max(`ylength',length("`yvar_`r''")) } local ylength1 = max(`ylength',15) local col2 _col(`=`ylength1'+2') local col3 _col(`=`ylength1'+13') local col4 _col(`=`ylength1'+24') local col5 _col(`=`ylength1'+35') local col6 _col(`=`ylength1'+44') local col7 _col(`=`ylength1'+55') local line4 as text _dup(`=`ylength1'+34') "{c -}" local line5 as text _dup(`=`ylength1'+43') "{c -}" local line7 as text _dup(`=`ylength1'+63') "{c -}" // OUTPUT BETWEEN-STUDY SDs AND I^2, WITH CIs if "`e(cholnames)'"=="cholnames" di as error "Sorry, i2 option is not available with cholnames" di as text _newline "Approximate confidence intervals for between-studies SDs and I^2:" di `line7' di as text "Variable" `col2' " SD" `col3' "[`level'% Conf. Interval]" `col5' " I^2" `col6' "[`level'% Conf. Interval]" _new `line7' forvalues r=1/`p' { di as text "`yvar_`r''" as result `col2' `bssdest`r'' `col3' `bssdlow`r'' `col4' `bssdupp`r'' `col5' `i2fmt' `i2est`r'' `col6' `i2fmt' `i2low`r'' `col7' `i2fmt' `i2upp`r'' "`droppedterm`r''" } di `line7' if "`nlcomfails'"=="yes" di as text "-nlcom- failed to produce some estimates: this usually occurs when one or more basic variance parameters is zero" if "`nlcomfails'"=="equals" di as text "Note: CIs are not appropriate with bscovariance(equals)" di as text "`i2footnote'" if "`cifootnote'"!="" di as text "`cifootnote'" if "`heterogineeded'"=="yes" di as error "To get CIs, please install heterogi using {stata ssc install heterogi}" if "`s2problem'"=="yes" di as text "Note: one or more values of I^2 weren't computed because Qa and/or Qb was missing" if "`droppedterms'"=="yes" di as text "* CI for I^2 ignores zero variance components" // COMPUTE AND OUTPUT BETWEEN-STUDY CORRELATIONS, WITH CIs local ylength1 = max(2*`ylength'+4,15) local col2 _col(`=`ylength1'+2') local col3 _col(`=`ylength1'+13') local col4 _col(`=`ylength1'+24') local col5 _col(`=`ylength1'+35') local col6 _col(`=`ylength1'+44') local col7 _col(`=`ylength1'+55') local line4 as text _dup(`=`ylength1'+34') "{c -}" local line5 as text _dup(`=`ylength1'+43') "{c -}" local line7 as text _dup(`=`ylength1'+63') "{c -}" if inlist("`bsest'","reml","ml") & `p'>1 & "`bscov'"=="unstructured" { di as text _newline "Between-study correlations:" _newline `line4' _newline "Variables" `col2' "Correl." `col3' "[`level'% Conf. Interval]" _new `line4' forvalues r=1/`p' { local rplus1 = `r'+1 local rminus1 = `r'-1 forvalues s=1/`rminus1' { local sdsd (`tau`r''*`tau`s'') local covar (`e(Sigma`s'`r')') if abs(`covar'/`sdsd')<1-1E-5 { qui nlcom log( (`sdsd'+`covar') / (`sdsd'-`covar') ) mat `b' = r(b) mat `V' = r(V) local est = `b'[1,1] local lower = `b'[1,1] - `z' * sqrt(`V'[1,1]) local upper = `b'[1,1] + `z' * sqrt(`V'[1,1]) foreach thing in est lower upper { local `thing' = (exp(``thing'')-1) / (exp(``thing'')+1) } } else { local est = (`covar'/`sdsd') local lower . local upper . } di as text _col(1) "`yvar_`s'' & `yvar_`r''" as result `col2' `est' `col3' `lower' `col4' `upper' } } di `line4' di as text "Note: CI computed on log((1+corr)/(1-corr)) scale" } } } if "`testsigma'"=="testsigma" { if inlist("`bsest'","reml","ml") { local lrt = 2*(e(ll)-e(ll0)) if abs(`lrt')0 local pvalue = `pvalue'/2 di if "`bsest'"=="reml" di as text "Restricted likelihood " _c else di as text "Likelihood " _c di as text "ratio test for heterogeneity: LRT = " as result %6.2f `lrt' as text " (d.f. = " as result `df' as text ") P = " as result %5.3f `pvalue' if `lrt'>0 di as text "P-value halved: see {help j_chibar:help j_chibar}" } else { di _newline as error "testsigma option is only available with reml or ml estimation methods" } } if "`qscalar'"=="qscalar" { di _new as text "Scalar Q test for heterogeneity: Q = " as result %6.2f e(Qscalar_chi2) _c di as text " (d.f. = " as result e(Qscalar_df) _c di as text ") P = " as result %5.3f chi2tail(e(Qscalar_df),e(Qscalar_chi2)) } if !mi("`randfix'`randfix2'") { cap confirm matrix e(V_fixed) local noVfixed = _rc>0 if "`e(bsest)'"=="" { di as error "option randfix ignored - no estimation done" } else if "`e(bsest)'"=="fixed" { di as error "option randfix ignored - not appropriate with fixed-effect analysis" } else if `noVfixed' { di as error "option randfix ignored - mvmeta didn't estimate the fixed-effect model" } else { di _new as text "Multivariate R statistic" di "Measures ratio of std errors in current vs. fixed-effect analysis" local yvars = e(yvars) if mi("`randfix2'") local randfix `yvars' else local randfix `randfix2' * check specified variables are part of e(yvars) local diff : list randfix - yvars if "`diff'"!="" { di as error "Variable(s) `diff' in randfix() are not outcome variables" exit 198 } local nparms_mean = `e(nparms_mean)' tempname hold varfixed varrand select mat `varfixed' = e(V_fixed) mat `varfixed' = `varfixed'[1..`nparms_mean',1..`nparms_mean'] mat `varrand' = e(V) mat `varrand' = `varrand'[1..`nparms_mean',1..`nparms_mean'] // find appropriate subvariance matrix if "`e(parmtype)'"=="short" local rownames: rownames `varfixed' else if "`e(parmtype)'"=="long" local rownames: roweq `varfixed' else exit 497 foreach varname of varlist `randfix' { local rownames : subinstr local rownames "`varname'" "1", all } foreach varname of varlist `yvars' { local rownames : subinstr local rownames "`varname'" "0", all } local rownames : subinstr local rownames " " ", ", all mat `select' = (`rownames') // calculate and output mata: selectpart("`select'","`varfixed'","`varfixed'") mata: selectpart("`select'","`varrand'","`varrand'") `ifdebug' mat l `varfixed', title(Variance-covariance under fixed-effect analysis) `ifdebug' mat l `varrand', title(Variance-covariance under current analysis) local dim = rowsof(`varfixed') local dimrand = rowsof(`varrand') if `dim' != `dimrand' { di as error "Dimensions of varfixed and varrand differ:" mat l `varfixed', title(Variance-covariance under fixed-effect analysis) mat l `varrand', title(Variance-covariance under current analysis) exit 459 } local detfixed = sqrt(det(`varfixed')) local detrand = sqrt(det(`varrand')) local rstat = (`detrand'/`detfixed')^(1/`dim') local col 49 local line as text _dup(78) "{c -}" di `line' di as text "Outcomes considered" _col(`col') as result "`randfix'" di as text "Number of parameters" _col(`col') as result `dim' di as text "Sqrt determinant of variance (`e(bsest)') " _col(`col') as result `detrand' di as text _col(30) "(fixed) " _col(`col') as result `detfixed' di as text "R statistic = ratio" _c if `dim'>1 di as text " ^ (1/" `dim' ")"_c di as text _col(`col') as result `rstat' di `line' } } *** PBEST if "`pbest'"!="" { if "`e(parmtype)'"!="long" di as error /// "Sorry, pbest() is only available after running mvmeta with the longparm option" else pbest `pbest' } *** BUBBLE PLOT if !missing(`"`bubble'`bubble2'"') { mvmeta_bubble, `bubble2' `debug' `hasestimated' warning(`warning') } *** FOREST PLOT if !missing(`"`forest'`forest2'"') { * mvmeta_forest, `forest2' `debug' di as error "Sorry, forest plot is not yet available" } *** PREDICTION INTERVAL if !mi("`pi'`pi2'") { cap confirm matrix e(Sigma) if "`e(bsest)'"=="fixed" di as error "Prediction intervals not reported - not meaningful after a fixed-effect model" else if _rc di as error "Prediction intervals not reported - Sigma was not estimated" else if !mi("`pi2'") mvmeta_pi, `pi2' else if !mi("`pi'") mvmeta_pi } end *=========================== END OF REPLAY PROGRAM ============================ *========================== START OF PBEST PROGRAM ============================ program define pbest // Parse syntax anything [if] [in], [ /// REPs(int 1000) zero gen(string) seed(int -1) format(string) /// documented id(varname) PREDict BESTonly saving(string) replace clear bar line /// documented CUMulative TABDISPoptions(string) mcse MEANrank /// documented zeroname(string) STRIPprefix(string) rename(string) all /// undocumented title(passthru) note(passthru) *] local minmax `anything' if !inlist("`minmax'","min","max") { di as error "Syntax: pbest(min|max, [options])" exit 198 } if `seed'!=-1 set seed `seed' if mi("`zeroname'") local zeroname zero marksample touse if !mi("`bestonly'") { if !mi("`all'") { di as error "Please don't specify both bestonly and all" exit 498 } if !mi("`line'") { di as text "Option bestonly ignored - cannot be specified with line" local bestonly } } if mi("`bestonly'") local all all * parse id variable if mi("`id'") & !mi(e(id)) { cap confirm var `e(id)' if !_rc local id `e(id)' } if mi("`id'") { tempvar id gen `id' = _n local idname _id } else { local idname `id' } cap isid `id' if _rc { di as error "Error in pbest: variable `id' does not uniquely identify the observations exit 459 } * also need a numeric id variable: changed v3.0.5 1jul2015 cap confirm numeric var `id' if _rc { tempvar idnum cap noi sencode `id', gen(`idnum') if _rc==199 di "Please install sencode using {stata ssc install sencode}" if _rc exit _rc } else local idnum `id' if !mi("`meanrank'") local all all forvalues r=1/`e(dims)' { foreach var in `e(xvars_`r')' { cap confirm var `var' if _rc { di as error "pbest: `var' not found" exit _rc } if _rc exit _rc cap assert !mi(`var') if `touse' if _rc==9 di as error "pbest: `var' must be non-missing `if' `in'" if _rc exit _rc } } * renaming local i 0 while !mi("`rename'") { local ++i gettoken thisrename rename : rename, parse(",") gettoken renamein`i' renameout`i' : thisrename, parse("=") local renameout`i' : subinstr local renameout`i' "=" "" local renamein`i' = trim(subinstr("`renamein`i''", ",","",1)) local renameout`i' = trim(subinstr("`renameout`i''", "=", "",1)) } local renamen `i' // retrieve stored results local p = e(dims) forvalues r=1/`p' { local yvar_`r' : word `r' of `e(yvars)' } // initialise counters if mi("`clear'") tempvar pbest rank treat else { local pbest _Pbest local rank _Rank local treat _Treat } tempvar best rbest pred mvn local rmin = cond("`zero'"=="zero", 0, 1) local smax = cond("`all'"=="all", `p' + ("`zero'"=="zero"), 1) forvalues s=1/`smax' { forvalues r=`rmin'/`p' { qui gen `pbest'`r'_`s'=0 if `touse' if `r'>0 { // 8apr2022: subinstr fails in Stata12 if "from" is missing if mi("`stripprefix'") local thischarold `yvar_`r'' else local thischarold : subinstr local yvar_`r' "`stripprefix'" "" } else local thischarold `zeroname' local thischarnew forvalues i=1/`renamen' { // rename if "`thischarold'" == "`renamein`i''" local thischarnew `renameout`i'' } if mi("`thischarnew'") local thischarnew `thischarold' if `s'==1 local trtlabel `trtlabel' `r' "`thischarnew'" if "`all'"=="all" local thischarnew `s':`thischarnew' char `pbest'`r'_`s'[varname] `thischarnew' } } tempname hold bstar b2 V2 // initialise save-file if "`saving'" != "" { if substr("`saving'",length("`saving'")-3,.)!=".dta" local saving `saving'.dta if mi("`replace'") { noi confirm new file `saving' } qui levelsof `idnum' if `touse', local(idlevels) local rmin = cond("`zero'"=="zero",0,1) forvalues r=`rmin'/`p' { local postfilevars `postfilevars' pred`r' } tempname pbestpost postfile `pbestpost' `idname' rep `postfilevars' using `saving', `replace' } // do parametric bootstraps -> creates variables `pbest'* forvalues rep=1/`reps' { drawbeta _estimates hold `hold', copy mat `bstar'=r(bstar) mat `b2' = `bstar'[1,1..`e(nparms_mean)'] mat `V2'= invsym(e(V)) mat `V2' = `V2'[1..`e(nparms_mean)',1..`e(nparms_mean)'] mat `V2' = invsym(`V2') ereturn post `b2' `V2' if "`zero'"=="zero" gen `pred'0 = 0 forvalues r=1/`p' { qui predict `pred'`r' if `touse', eq(`yvar_`r'') label var `pred'`r' "Prediction for outcome `r'" } _estimates unhold `hold' if "`predict'"=="predict" { * add random error with variance e(Sigma) drawnorm `mvn'_1 - `mvn'_`p', cov(e(Sigma)) forvalues r=1/`p' { qui replace `pred'`r' = `pred'`r' + `mvn'_`r'[1] if `touse' } drop `mvn'_1-`mvn'_`p' } if "`saving'" != "" { foreach i in `idlevels' { local postexps (`i') (`rep') forvalues r = `rmin'/`p' { local postexp = `pred'`r'[`i'] local postexps `postexps' (`postexp') } post `pbestpost' `postexps' } } forvalues s=1/`smax' { qui gen `best' = . // sth best value qui gen `rbest' = . // which value is sth best forvalues r=`rmin'/`p' { qui replace `rbest' = `r' if `pred'`r'==`minmax'(`pred'`r', `best') & `touse' & !mi(`pred'`r') qui replace `best' = `pred'`r' if `pred'`r'==`minmax'(`pred'`r', `best') & `touse' & !mi(`pred'`r') } forvalues r=`rmin'/`p' { // `pbest'`r'_`s' = whether rth trt is sth best qui replace `pbest'`r'_`s' = `pbest'`r'_`s'+100/`reps' if `rbest'==`r' & `touse' qui replace `pred'`r' = . if `rbest' ==`r' } drop `best' `rbest' } drop `pred'* } * Finish save file if "`saving'" != "" postclose `pbestpost' // OUTPUT preserve * header di as text _new "Estimated probabilities (%) of each treatment " _c if `smax'>1 di "having each rank" else di "being the best" di as text "- assuming the " as result "`minmax'imum" as text " parameter is the best" di as text "- using " as result `reps' as text " draws" _c if "`saving'" != "" di as text " (written to file " as result "`saving'" as text ")" _c di di "- allowing for parameter uncertainty" _c if "`predict'"=="predict" di " and between-studies variation" _c * reshape into a long data set qui keep if `touse' keep `idnum' `pbest'* *if "`id'"=="" { * tempvar id * gen `id'=_n * label var `id' "id" *} qui count local multid = r(N)>1 forvalues r=`rmin'/`p' { local pbestlist `pbestlist' `pbest'`r' } qui reshape long `pbestlist', i(`idnum') j(`rank') string qui replace `rank' = substr(`rank',2,.) // removes initial "_" qui destring `rank', replace qui reshape long `pbest', i(`idnum' `rank') j(`treat') * MC error if !mi("`mcse'") { tempvar mcse gen `mcse' = sqrt(`pbest'*(100-`pbest')/`reps') local dispvars `pbest' `mcse' di as text _n "- figures are estimated probability (upper), Monte Carlo error (lower)" _c } else local dispvars `pbest' * label variables and values char `rank'[varname] Rank char `treat'[varname] Treat char `idnum'[varname] `idname' char `pbest'[varname] Pbest label var `rank' "Rank" label var `treat' "Treatment" label var `idnum' "`idname'" label var `pbest' "Pbest" label def `treat' `trtlabel' label val `treat' `treat' if "`format'"=="" local format %6.1f // best for table format `pbest' `mcse' `format' forvalues s=1/`smax' { if `s'==1 local text "Best" else if `s'==`smax' local text "Worst" else if `s'==2 local text `s'nd else if `s'==3 local text `s'rd else local text `s'th label def `rank' `s' "`text'", add } label val `rank' `rank' * output mean rank & SUCRA local recordtype _Recordtype gen `recordtype' = 0 if !mi("`meanrank'") { sort `idnum' `treat' `rank' tempvar newrecord meanrank qui expand 2 if `rank'==`smax', gen(`newrecord') qui replace `rank'=`smax'+10 if `newrecord' // for meanrank qui replace `recordtype'=1 if `newrecord' drop `newrecord' qui expand 2 if `rank'==`smax', gen(`newrecord') qui replace `rank'=`smax'+11 if `newrecord' // for sucra qui replace `recordtype' = 2 if `newrecord' drop `newrecord' label def `recordtype' 0 "Probability" 1 "Mean rank" 2 "SUCRA" label val `recordtype' `recordtype' qui replace `pbest'=. if `recordtype' sort `idnum' `treat' `rank' qui by `idnum' `treat': gen `meanrank' = sum(`rank'*`pbest')/sum(`pbest') qui replace `pbest' = `meanrank' if `recordtype'==1 qui replace `pbest' = (`smax' - `meanrank') / (`smax' - 1) if `recordtype'==2 label def `rank' `=`smax'+10' "MEAN RANK" `=`smax'+11' "SUCRA", add if !mi("`mcse'") { tempvar meanrank2 qui by `idnum' `treat': gen `meanrank2' = sum(`rank'^2*`pbest')/sum(`pbest') replace `mcse' = sqrt(`meanrank2' - `meanrank'^2)/sqrt(`reps') if `recordtype' replace `mcse' = `mcse' / (`smax' - 1) if `recordtype'==2 drop `meanrank2' } // NB `smax' = #compared because meanrank option => all option } * tabulate /* old local byid by(`idnum') // could make this -if `multid'-? if `smax'==1 local tabcmd tabdisp `idnum' `treat', c(`dispvars') `tabdispoptions' else local tabcmd tabdisp `rank' `treat', `byid' c(`dispvars') `tabdispoptions' */ // could make this -if `multid'-? if `multid' local byid by(`idnum') local tabcmd tabdisp `rank' `treat', `byid' c(`dispvars') `tabdispoptions' `tabcmd' // GRAPH format `pbest' %6.0f // best for axis label if !mi("`bar'") { if !mi("`cumulative'") local stack stack if `smax'>1 { local overrank over(`rank') local legendtitle title("Rank") } else local legendtitle title("Treatment") if `multid' local byid by(`idnum', `title' `note') else local byid `title' `note' local graphcmd graph bar `pbest', `overrank' over(`treat') `byid' asy ytitle("Probability (%)") `stack' legend(`legendtitle') `options' } else if !mi("`line'") { if `multid' di as error "Graphs for multiple records will be overlaid" if !mi("`cumulative'") { sort `idnum' `treat' `rank' local pbestcum `pbest'cum qui by `idnum' `treat': gen `pbestcum' = sum(`pbest') if !`recordtype' local cumulative Cumulative local pbestvar `pbestcum' } else local pbestvar `pbest' if mi("`note'") local note note("") if !mi("`mcse'") { // I THINK THE MCSEs ARE WRONG WITH CUMULATIVE OPTION?? if !mi("`cumulative'") qui replace `mcse' = `pbestcum' * (100-`pbestcum') / `reps' if !`recordtype' tempvar pbestlow pbestupp qui gen `pbestlow'=`pbestvar'-1.96*`mcse' if !`recordtype' qui gen `pbestupp'=`pbestvar'+1.96*`mcse' if !`recordtype' local rspike (rspike `pbestlow' `pbestupp' `rank' if !`recordtype', pstyle(p1)) } local graphcmd twoway (line `pbestvar' `rank' if !`recordtype', c(L) pstyle(p1)) `rspike', by(`treat', `title' `note' imargin(medium) legend(off)) ytitle("`cumulative' Probability (%)") xlabel(minmax,val) `options' // imargin(medium) avoids "worst" of one column overlapping and "best" of next column } `graphcmd' format `pbest' `format' // FINISH OFF if !mi("`clear'") { // keep variables in memory di as text "Rank probabilities have been loaded into memory" global F7 `tabcmd' di as text "tabulate command is stored in F7" global F8 `graphcmd' di as text "graph command is stored in F8" restore, not exit } else if "`gen'"!="" { // optionally generate variables restore forvalues s=1/`smax' { forvalues r=`rmin'/`p' { local rname = cond(`r'>0,"`yvar_`r''","zero") local sname = cond(`s'>1,"`s'","`minmax'") rename `pbest'`r'_`s' `gen'`sname'_`rname' label var `gen'`sname'_`rname' "`rname' has rank `sname'" } } } else drop `pbest'* // not needed? end *============================= END OF PBEST PROGRAM ============================ *========================== START OF DRAWBETA PROGRAM (for pbest) ========================== program define drawbeta, rclass syntax, [tweak(real 1E-12)] tempname b V chol matrix `b' = e(b) matrix `V' = e(V) local p = colsof(`b') return scalar colsofb = colsof(`b') capture matrix `chol' = cholesky(`V') local rc = c(rc) if `rc' == 506 { matrix `chol' = cholesky((`tweak' * trace(`V')/`p') * I(`p') + `V') } else if `rc' > 0 { di as err "Cholesky decomposition failed with error " `rc' exit `rc' } tempname e bstar matrix `e' = J(1, `p', 0) forvalues i = 1 / `p' { matrix `e'[1, `i'] = invnormal(uniform()) } matrix `bstar' = `b' + `e' * `chol'' return matrix bstar = `bstar' end *=========================== END OF DRAWBETA PROGRAM =========================== *=========================== START OF MVMETA_WT PROGRAM =========================== program define mvmeta_wt /* v1.2 22apr2015 uses e(V_uv) id() no longer allowed (taken from e(ydata)) v1.1.1 11mar2015 bug fixes not sortpreserve (interferes with clear option) copes if 1st equation is empty mvmeta call to create matrices now has option -fixed- to avoid jointcheck returning error v1.1 9feb2015 no run option new syntax results as variables or matrices v1.0 3feb2015 try to make it handle covariates drop abs matrices aren't left lying around made output sensible and agreeing as far as poss across methods To do: Note: previous -borrow- option reported (Vuv/Vmv-1) cf here it correctly reports (1-Vmv/Vuv) Requires: sencode [from SSC] Note for mvmeta: this will be a weight option for mvmeta, because both wt and bos can be understood as sharing the weight between studies and between sources of information within studies. Syntax will be [nowt wt(suboptions)]. Default will be wt(sd) unless nowt is specified. */ // FIND ERETURNED RESULTS local n = e(N) local p = e(dims) * find covariates local hasx 0 local qsum0 0 forvalues r=1/`p' { local yvar_`r' : word `r' of `e(yvars)' local yvarsquoted `"`yvarsquoted' "`yvar_`r''""' local xvars_`r'=e(xvars_`r') if !mi(e(xvars_`r')) { local q`r' = wordcount(e(xvars_`r')) local hasx 1 if !mi("`eqs'") local eqs `eqs', local eqs `eqs' `yvar_`r'':`xvars_`r'' } else local q`r' 0 if mi(e(constant)) local q`r' = `q`r'' + 1 local qsum`r' = `qsum`=`r'-1'' + `q`r'' * Outcome `r' has `q`r'' covariates - will be useful later } if !mi(e(constant)) local constant = e(constant) // PARSE syntax, [ /// SD RV DPC /// the 3 methods, which can't be specified together DETails /// to output a table of the full SD (SD method only) /// or of the separate SEs (RV method only) Format(passthru) /// format for all methods CLear /// loads the data for the table into memory (SD method only) Keepmat(name) /// to save the matrices (all methods - 3 matrices for SD, 1 for RV & DPC) UNscaled /// also unscaled details (SD method only) Wide /// wide format (SD + details only) mat(name) /// weights matrix to be saved: undocumented noSTUDies /// no study weights (SD method only) debug /// undocumented ] if !mi("`debug'") di as input "Running: mvmeta_wt `0'" if wordcount("`sd' `rv' `dpc'")>1 { di as error "Please specify only one method out of: `sd' `rv' `dpc'" exit 198 } if mi("`sd'`rv'`dpc'") local sd sd // SD is default if !mi("`clear'") & (mi("`sd'") | mi("`details'")) di as error "mvmeta_wt: option clear ignored (only relevant with sd and details options)" if !mi("`wide'") & (mi("`sd'") | mi("`details'")) di as error "mvmeta_wt: option wide ignored (only relevant with sd and details options)" if !mi("`unscaled'") & mi("`sd'") di as error "mvmeta_wt: option unscaled ignored (only relevant with sd option)" if mi("`format'") & !mi("`sd'`rv'") local format format(%6.1f) if mi("`e(V_uv)'") & !mi("`rv'") { di as error "mvmeta_wt: option rv is not allowed since univariate results were not found" exit 459 } tempname V XViX Z W w wV tot xpart total direct borrowed // CREATE DATA MATRICES tempname hold ymat Smat Xmat mat `ymat' = e(ydata) mat `Smat' = e(Sdata) mat `Xmat' = e(Xdata) forvalues i = 1/`n' { mat `ymat'`i'=`ymat'[`i',1...] local first = `p'*(`i'-1)+1 local last = `p'*`i' mat `Smat'`i'=`Smat'[`first'..`last',1...] mat `Xmat'`i'=`Xmat'[`first'..`last',1...]' local studyname : rowname `ymat'`i' local studynames `"`studynames' `"`studyname'"'"' } // GENERAL CALCULATIONS forvalues i=1/`n' { mat `V'`i' = `Smat'`i' * "augment" missing values local Vmax 0 forvalues r=1/`p'{ local Vmax = max(`Vmax',`V'`i'[`r',`r']) } forvalues r=1/`p'{ forvalues s=1/`p'{ if mi(`V'`i'[`r',`s']) mat `V'`i'[`r',`s'] = cond(`r'==`s',`Vmax'*1E6,0) } } if "`e(bsest)'"!="fixed" { mat `V'`i' = `V'`i' + e(Sigma) if "`e(wscorr)'"=="riley" { mat `V'`i' = cholesky(diag(vecdiag(`V'`i'))) mat `V'`i' = `V'`i'*corr(e(Sigma))*`V'`i' * now need to re-set augmented off-diag elements to 0 for Riley method forvalues r=1/`p'{ forvalues s=1/`p'{ if `s'==`r' continue if mi(`Smat'`i'[`r',`r']) | mi(`Smat'`i'[`s',`s']) mat `V'`i'[`r',`s'] = 0 } } } } if !mi("`debug'") mat l `V'`i', title(Total variance matrix for study `i') mat `XViX'`i' = `Xmat'`i' * invsym(`V'`i') * `Xmat'`i'' if `i'==1 mat `XViX'sum = `XViX'`i' else mat `XViX'sum = `XViX'sum + `XViX'`i' } mat `Z' = invsym(`XViX'sum) // SCORE DECOMPOSITION METHOD: WEIGHTS, BoS, OPTIONAL DETAILS if !mi("`sd'") { forvalues i=1/`n' { forvalues r=1/`p' { local rprev = `r'-1 local start = `qsum`rprev''+1 local end = `qsum`r'' local rbit `start'..`end' mat `Z'part = `Z'[`rbit',1...] mat `Z'part2 = `Z'[`rbit',`rbit'] mat `xpart' = `Xmat'`i'[`rbit',`r'] * total contribution mat `total'_ir = vecdiag(`Z'part * `XViX'`i' * `Z'part') if `r'==1 mat `total'_i = `total'_ir else mat `total'_i = `total'_i , `total'_ir * direct contribution mat `direct'_ir = vecdiag(`Z'part2 * `xpart' * `xpart'' * `Z'part2 * (1 / `V'`i'[`r',`r'])) if `r'==1 mat `direct'_i = `direct'_ir else mat `direct'_i = `direct'_i , `direct'_ir } if `i'==1 mat `total' = `total'_i else mat `total' = `total' \ `total'_i if `i'==1 mat `direct' = `direct'_i else mat `direct' = `direct' \ `direct'_i } mat rownames `total' = `studynames' mat rownames `direct' = `studynames' mat `borrowed' = `total' - `direct' * now have: `direct', `borrowed' and `total' which are studies x parameters foreach source in total direct borrowed { if e(parmtype)=="short" { mat colnames ``source'' = `yvarsquoted' mat coleq ``source'' = "Overall" } mat ``source''_sum = J(1,`n',1) * ``source'' // row sum mat rownames ``source''_sum = "`=upper("`source'")'" mat ``source''_scaled = 100 * ``source'' * invsym(diag(`total'_sum)) mat ``source''_scaled_sum = J(1,`n',1) * ``source''_scaled // row sum mat rownames ``source''_scaled_sum = "`=upper("`source'")'" if !mi("`keepmat'") { if mi("`unscaled'") mat `keepmat'`source' = ``source''_scaled else mat `keepmat'`source' = ``source'' } } * now also have `direct'_sum, `borrowed'_sum and `total'_sum which are 1 x parameters * and `direct'_scaled, `borrowed'_scaled and `total'_scaled which are studies x parameters * Main output: study weights and BoS mat `total'_scaled_withtotal = `total'_scaled_sum \ `borrowed'_scaled_sum \ `direct'_scaled_sum if "`studies'"!="nostudies" mat `total'_scaled_withtotal = `total'_scaled \ `total'_scaled_withtotal di as text _n "Weights using score decomposition method" if "`studies'"!="nostudies" di as text "Study weights and borrowing of strength (% of total for each parameter)" else di as text "Borrowing of strength (% of total for each parameter)" di "{hline 71}" _c mat l `total'_scaled_withtotal, noheader `format' if !mi("`details'") { // now output as a data set if mi("`clear'") preserve qui { drop _all set obs `n' * gen id = "" gen id = _n /*1jul2015*/ forvalues i=1/`n' { local studyname : word `i' of `studynames' * replace id = "`studyname'" in `i' label def id `i' "`studyname'", add /*1jul2015*/ } label val id id /*1jul2015*/ tempname onecol forvalues k=1/`qsum`p'' { foreach source in direct borrowed total { mat `onecol' = ``source''[1...,`k'] svmat `onecol', names(`source'`k') rename `source'`k'1 `source'`k' } local outcome`k' : coleq `onecol' local covariate`k' : colnames `onecol' local parameter`k' : colfullnames `onecol' } reshape long direct borrowed total, i(id) j(parm) local parmids outcome covariate parameter foreach parmid in `parmids' { gen `parmid' = "" forvalues k=1/`qsum`p'' { replace `parmid' = "``parmid'`k''" if parm==`k' } cap noi sencode `parmid', replace if _rc==199 noi di "Please install sencode using {stata ssc install sencode}" if _rc exit _rc } drop parm egen totsum = sum(total), by(parameter) foreach source in direct borrowed total { gen scaled`source'=100*`source'/totsum rename `source' unscaled`source' } drop totsum reshape long scaled unscaled, i(id `parmids') j(source) string * next 4 lines choose order "total, borrowed, direct" gen sourcenum = 1*(source=="total") + 2*(source=="borrowed") + 3*(source=="direct") replace source=upper(source) sencode source, gsort(sourcenum) replace drop sourcenum } local var = cond(mi("`unscaled'"),"scaled","unscaled") di as text _n "Details of `var' weights:" _c if !mi("`wide'") { local parmid qui tab outcome if r(r)>1 local parmid outcome qui tab covariate if r(r)>1 local parmid `parmid' covariate if "`parmid'" == "outcome covariate" local parmid parameter local cmd table id `parmid' source, c(sum `var') `format' row } if mi("`wide'") { local cmd table id covariate outcome, c(sum `var') by(source) `format' row } `cmd' if !mi("`clear'") { di as text "Tabulated data are now in memory. Press F9 to recall table command." global F9 `cmd' exit } restore } if !mi("`mat'") { mat `mat' = `total'_scaled } di as text "{hline 71}" } // RELATIVE VARIANCES METHOD: BoS ONLY if !mi("`rv'") { tempname Vmv Vuv Vuvq rvmat all mat `Vmv' = e(V) mat `Vuv' = e(V_uv) mat `Vuv'=`Vuv'[1..`qsum`p'',1..`qsum`p''] local names : rowfullnames `Vuv', quoted foreach source in direct borrowed total { mat ``source''_sum = J(1,`qsum`p'',.) mat colnames ``source''_sum = `names' mat ``source''_scaled_sum = J(1,`qsum`p'',.) mat colnames ``source''_scaled_sum = `names' mat rownames ``source''_scaled_sum = "`=upper("`source'")'" } mat `rvmat' = J(3,`qsum`p'',.) mat `rvmat'_details = J(2,`qsum`p'',.) mat rownames `rvmat'_details = "Univariate" "Multivariate" mat rownames `rvmat' = "TOTAL" "BORROWED" "DIRECT" forvalues i=1/`qsum`p'' { local vmv = `Vmv'[`i',`i'] local vuv = `Vuv'[`i',`i'] mat `rvmat'_details[1,`i'] = sqrt(`vuv') mat `rvmat'_details[2,`i'] = sqrt(`vmv') mat `rvmat'[1,`i'] = 1 // total mat `rvmat'[2,`i'] = 1 - float(`vmv'/`vuv') // borrowed; float avoids tiny negative values mat `rvmat'[3,`i'] = `vmv'/`vuv' // direct } mat `rvmat' = 100 * `rvmat' mat colnames `rvmat' = `names' mat colnames `rvmat'_details = `names' di as text _n "Weights using relative variances method" di "Borrowing of strength (% of total for each parameter)" di as text "{hline 53}" _c mat l `rvmat', noheader `format' if !mi("`details'") { di _n as text "Standard errors used in the relative variance method:" _c mat l `rvmat'_details, noheader } if !mi("`keepmat'") { if !mi("`details'") mat `keepmat'se = `rvmat'_details mat `keepmat' = `rvmat' } di as text "{hline 53}" } // COMPUTE DATA-POINT COEFFICIENTS if !mi("`dpc'") { forvalues i=1/`n' { mat `W'`i' = `Z' * `Xmat'`i' * invsym(`V'`i') local thisid : word `i' of `studynames' mat coleq `W'`i' = `"`thisid'"' if `i'==1 mat `W'=`W'`i'' else mat `W' = `W' \ `W'`i'' } if e(parmtype)=="short" { mat colnames `W' = `yvarsquoted' mat coleq `W' = "Overall" } di as text _n "Weights using data-point coefficients method" di as text "Coefficients for parameters by data points" di as text "{hline 70}" _c mat l `W', noheader `format' if !mi("`keepmat'") mat `keepmat' = `W' di as text "{hline 70}" } end *=========================== END OF MVMETA_WT PROGRAM =========================== *=========================== START OF MVMETA_BUBBLE PROGRAM =========================== /* version 1.6 Ian White 22apr2015 uses e(ydata) etc. not current data no id() option 24apr2015 - data construction outsourced to mvmeta_getdata version 1.5 Ian White 11feb2015 */ program define mvmeta_bubble // PARSE syntax, /// [Variables(string) noMv noUv STUDylabel(passthru) UVLABel(passthru) MVLABel(passthru) clear MLABel /// hasestimated debug warning(string) * /// to be undocumented ] if mi("`hasestimated'") di as text "Bubble option: using data from last mvmeta run" // in replay mode if "`warning'"=="" local warning error if !inlist("`warning'", "error", "text", "off") { di as error "syntax: warning(error|text|off)" exit 198 } local warn warn`warning' * allows all bubble options if !mi("`debug'") di as input "Call: mvmeta_bubble `0'" if mi(e(bsest)) { local uv nouv local mv nomv local wt } if mi("`variables'") { local variables = word(e(yvars),1)+" "+word(e(yvars),2) if wordcount(e(yvars)) > 2 di as text "Bubble option: no variables specified, displaying " as result word(e(yvars),1) as text " and " as result word(e(yvars),2) } cap unab variables : `variables' tokenize "`variables'" if mi("`2'") | !mi("`3'") { di as error "Bubble option: variables() suboption must contain exactly two variables" exit 198 } local y `1' local x `2' local yend = subinstr("`y'",e(ystub),"",1) local xend = subinstr("`x'",e(ystub),"",1) // check if there are covariates local hasx 0 forvalues i = 1/`e(dims)' { if !mi("`e(xvars_`i')'") local hasx 1 } if (`hasx' | e(parmtype)=="common") { if "`uv'"!="nouv" | "`mv'"!="nomv" { if `hasx' local text "had covariates" else if e(parmtype)=="common" local text "used commonparm option" else local text `warn' "bubble plot can't show summaries - last mvmeta run `text'" local mv nomv local uv nouv } } // LOAD DATA if mi("`clear'") preserve if mi("`clear'") tempvar type else local type _type mvmeta_getdata `y' `x', type(`type') corr `uv' `mv' `studylabel' `uvlabel' `mvlabel' if !mi("`mlabel'") local mopts mopts(mlab(_id)) cap confirm variable corr_`y'_`x' local corrvar = cond(_rc,"corr_`x'_`y'","corr_`y'_`x'") bubble `y' `x' se_`y' se_`x' `corrvar', `mopts' group(`type') `debug' `clear' `options' warning(`warning') if !mi("`debug'") di as input "Ending mvmeta_bubble" end *=========================== END OF MVMETA_BUBBLE PROGRAM =========================== *=========================== START OF BUBBLE PROGRAM =========================== /* version 1.4.3 Ian White 5mar2015 colby renamed group lpattern and lwidth can be specified by group version 1.4.2 23dec2014 uses inbuilt color scheme via pstyle() colourbystudy renamed colby version 1.4.1 26jun2014 - can detect mis-matched mean and sd version 1.4 26mar2012 - COLOURBystudy option gives different colors! version 1.3 26mar2012 - variables correctly saved with clear option version 1.2 17jan2011 - treats missing values as known zeroes (or other value) - new options missval() and stagger() version 1.1 14Sep2010 - plots 1st vs 2nd not vice versa Simple usage: bubble ymean xmean ysd xsd corr */ program define bubble // PARSE syntax varlist(min=2 max=5) [if] [in], [pct(numlist) n(int 36) clear /// main options MISSval(string) STAGger(string) /// missing value options GRoup(varname) COLors(string) LPatterns(string) LWidths(string) /// by-group options MLABel(passthru) MSymbol(passthru) MSIZe(passthru) MSTYle(passthru) /// marker options LSTYle(passthru) /// line options warning(string) /// output options lopts(string) mopts(string) cov debug eform *] // graph options if !mi("`debug'") di as input "Call: bubble `0'" tokenize "`varlist'" local ymean `1' local xmean `2' if "`warning'"=="" local warning error if !inlist("`warning'", "error", "text", "off") { di as error "syntax: warning(error|text|off)" exit 198 } local warn warn`warning' tempvar ysd xsd corr qui { if mi("`3'") gen `ysd' = 1 if !mi(`ymean') else gen `ysd' = `3' if !mi("`cov'") replace `ysd' = sqrt(`ysd') if mi("`4'") gen `xsd' = 1 if !mi(`xmean') else gen `xsd' = `4' if !mi("`cov'") replace `xsd' = sqrt(`xsd') if mi("`5'") gen `corr' = 0 if !mi(`ymean',`xmean') else gen `corr' = `5' if !mi("`cov'") replace `corr' = `corr'/ (`xsd'*`ysd') } if !mi("`debug'") { di as text "Data to be plotted (ymean, xmean, ysd, xsd, corr):" l `ymean' `xmean' `ysd' `xsd' `corr' } // MISSING VALUES tokenize "`missval'" local xmissval = cond("`1'"!="","`1'","0") local ymissval = cond("`2'"!="","`2'","`xmissval'") tokenize "`stagger'" local xstagger `1' local ystagger = cond("`2'"!="","`2'","`1'") if "`pct'"=="" local pct 50 if !mi("`group'") & mi("`msymbol'") local msymbol msymbol(S) * bubble-specific options if !mi("`colors'") forvalues i=1/`=wordcount("`colors'")' { if word("`colors'",`i')!="=" local col = word("`colors'",`i') if !inlist("`col'","",".") { local lcol`i' lcol(`col') local mcol`i' mcol(`col') } } if !mi("`lpatterns'") forvalues i=1/`=wordcount("`lpatterns'")' { if word("`lpatterns'",`i')!="=" local lpatt = word("`lpatterns'",`i') if !inlist("`lpatt'","",".") local lpatt`i' lpatt(`lpatt') } if !mi("`lwidths'") forvalues i=1/`=wordcount("`lwidths'")' { if word("`lwidths'",`i')!="=" local lwid = word("`lwidths'",`i') if !inlist("`lwid'","",".") local lwid`i' lwid(`lwid') } * common options local lopts `lopts' `lstyle' local mopts `mopts' `mstyle' `mlabel' `msymbol' `msize' // PRESERVE if "`clear'"=="" { preserve tempvar x y } else { local x `xmean' local y `ymean' } // MISSING VALUES qui drop if mi(`xmean') & mi(`ymean') tempvar ismiss ismisssum /* these warnings are not needed when bubble is called by mvmeta_bubble foreach z in x y { // warnings (picks up mis-matched mean and sd) qui count if mi(``z'mean') > mi(``z'sd') if r(N) `warn' "`=r(N)' observations with missing ``z'mean' and non-missing ``z'sd' qui count if mi(``z'mean') < mi(``z'sd') if r(N) `warn' "`=r(N)' observations with non-missing ``z'mean' and missing ``z'sd' } */ foreach z in x y { if "``z'stagger'"=="" { qui summ ``z'mean' local `z'stagger = (r(max)-r(min))/100 } gen `ismiss' = mi(``z'mean') gen `ismisssum' = sum(`ismiss') qui count if `ismiss' if r(N)>0 { di as result r(N) as text " missing values found for ``z'mean' and displayed at " as result "``z'mean'=" ``z'missval' qui replace ``z'mean' = ``z'missval' + ``z'stagger'*(`ismisssum'-(r(N)+1)/2) if `ismiss' qui replace ``z'sd' = 0 if `ismiss' qui replace `corr' = 0 if `ismiss' local `z'ismiss yes if mi("`eform'") local options `options' `z'label(``z'missval' "Missing", add) else local options `options' `z'label(`=exp(``z'missval')' "Missing", add) } drop `ismiss' `ismisssum' } // EXPAND marksample touse, novarlist qui keep if `touse' tempvar id z gen `id'=_n local nstudies=_N local nplus1=`n'+2 qui expand `nplus1' sort `id' qui by `id': gen _theta = (_n-1)*2*_pi/`n' if _n<_N // ADD CONTOUR VARIABLES local ly : var label `ymean' local lx : var label `xmean' if "`ly'"=="" local ly `ymean' if "`lx'"=="" local lx `xmean' local i 0 if !mi("`group'") { cap confirm numeric var `group' if _rc { cap noi sencode `group', replace // convert to numeric if _rc==199 di "Please install sencode using {stata ssc install sencode}" if _rc exit _rc } qui levelsof `group', local(grouplevels) } local l 0 foreach p of numlist `pct' { local ++i local a = sqrt(-2*log(1-`p'/100)) qui gen `y'`i' = `a'*sin(_theta) qui gen `x'`i' = `a'*cos(_theta) qui replace `x'`i' = (`corr')*`y'`i' - sqrt(1-(`corr')^2)*`x'`i' qui replace `x'`i'=`xmean'+`xsd'*`x'`i' qui replace `y'`i'=`ymean'+`ysd'*`y'`i' local s 0 if !mi("`group'") { foreach level in `grouplevels' { local ++s local cond `group'==`level' local graphlist `graphlist' /// (line `y'`i' `x'`i' if `cond', pstyle(p`s') c(l) cmissing(n) `lcol`s'' `lpatt`s'' `lwid`s'' `lopts') /// (scatter `ymean' `xmean' if `cond' & _theta==0, pstyle(p`s') `mcol`s'' `mopts') local ++l if `i'==1 local legendorder `legendorder' `l' if `i'==1 local legendlabel `legendlabel' label(`l' "`:label (`group') `level''") local ++l } } else { local graphlist `graphlist' (line `y'`i' `x'`i', cmissing(n) `lopts' `lcol`i'' `lpatt`i'' `lwid`s'' `mcol`i'') local legend `legend' label(`i' "`p'%") local legendorder `legendorder' `i' } label var `y'`i' "`ly'" label var `x'`i' "`lx'" } if !mi("`group'") { local legendopt legend(order(`legendorder') `legendlabel') foreach p of numlist `pct' { local note `note' `p'% } local note note("Showing `note' confidence region(s)") } else { local i1=`i'+1 local graphlist `graphlist' (scatter `ymean' `xmean', pstyle(p`i1') `mcol`i1'' `mopts') local legendopt legend(`legend' order(`legendorder') title("Probability")) } if !mi("`ylab'") local options ytitle(`"`ylab'"') `options' if !mi("`xlab'") local options xtitle(`"`xlab'"') `options' if !mi("`eform'") { foreach var of varlist `ymean' `xmean' `y'* `x'* { qui replace `var'=exp(`var') } local options `options' xscale(log) yscale(log) } // GRAPH local command twoway `graphlist', `note' `legendopt' `options' if "`clear'"=="clear" { global F9 `command' di as text "Bubble graph data loaded into memory" di as text "(command saved in F9)" } if !mi("`debug'") di as input `"`command'"' `command' if !mi("`debug'") di as input "Ending bubble" end *=========================== END OF BUBBLE PROGRAM =========================== program define dicmd noi di as input `"`0'"' `0' end *=========================== START OF MVMETA_ESTIMATE PROGRAM =========================== program define mvmeta_estimate, rclass // taken from mvmeta.ado, 2apr2015 *** PARSE syntax namelist(min=2 max=2) if, yends(string) parmtype(string) missest(string) missvar(string) [bsest(string) xvarslist(string) bscovariance(string) wscorr(string) sigma0(string) start(string) noconstant wscorrforce warningest(string) warningnoest(string) augment augquiet notrunc cholnames debug showstart mmfix taulog id(varname) noestimates ITERate(passthru) sigma0exp(string) penalty(string) psdcrit(real 0) noposdef *] mlopts mlopts, `options' // iterate is parsed separately as it is changed in 2nd -ml- run if mi("`debug'") local ifdebug * `ifdebug' di as input "mvmeta_estimate `0'" tokenize "`namelist'" local ystub `1' local Sstub `2' local p = wordcount("`yends'") marksample touse if `p'>1 & mi("`xvarslist'") exit 497 tokenize `"`xvarslist'"', parse(",") forvalues r=1/`p' { if "`1'"!="." local xvars_`r' `1' // . means missing if `r'<`p' & "`2'"!="," { di as error "mvmeta_estimate: error in xvarslist option" exit 497 } mac shift 2 } qui count if `touse' local n = r(N) local N = _N if "`constant'" != "noconstant" { tempvar xcons gen `xcons' = 1 local xconsname _cons } if "`estimates'"!="noestimates" & "`bsest'"=="" { di as error "mvmeta_estimate: bsest() required" exit 497 } foreach type in est noest { // warnings for estimate and no-estimate parts if "`warning`type''"=="" local warning`type' error if !inlist("`warning`type''", "error", "text", "off") { di as error "syntax: warning`type'(error|text|off)" exit 198 } } local warn warn`warningnoest' *** START OF CODE TAKEN FROM MVMETA.ADO * COUNT FIXED PARAMETERS local fixedparms 0 forvalues r=1/`p' { local yend_`r' = word("`yends'",`r') local yvar_`r' `ystub'`yend_`r'' local ylist `ylist' `yvar_`r'' local fixedparms`r' = wordcount("`xvars_`r'' `xconsname'") if "`parmtype'"=="common" & `fixedparms`r''!=`fixedparms1' { di as error "commonparm option: equations 1 and `r' may not have different numbers of variables" exit 198 } local fixedparmsum = `fixedparmsum' + `fixedparms`r'' } local fixedparms = cond("`parmtype'"=="common", `fixedparms1', `fixedparmsum') ****************** IDENTIFY COVARIANCE EXPRESSIONS ***************** tempvar covvar // stub for any covariance variables we have to create local wscorrunused 0 local wscorrused 0 forvalues r=1/`p'{ forvalues s = `r' / `p' { cap confirm var `Sstub'`yend_`r''`yend_`s'', exact local okrs = _rc==0 cap confirm var `Sstub'`yend_`s''`yend_`r'', exact local oksr = _rc==0 if `okrs' & (`s'==`r' | mi("`wscorrforce'")) local Svar_`r'_`s' `Sstub'`yend_`r''`yend_`s'' else if `oksr' & (`s'==`r' | mi("`wscorrforce'")) local Svar_`r'_`s' `Sstub'`yend_`s''`yend_`r'' else if "`wscorr'"=="riley" { qui gen `covvar'`r'`s' = 0 // just to satisfy the pos def check local Svar_`r'_`s' `covvar'`r'`s' } else if "`wscorr'"!="" { qui gen `covvar'`r'`s' = `wscorr'*sqrt(`Sstub'`yend_`r''`yend_`r''*`Sstub'`yend_`s''`yend_`s'') if `touse' local Svar_`r'_`s' `covvar'`r'`s' } else { local message "Neither `Sstub'`yend_`r''`yend_`s'' nor `Sstub'`yend_`s''`yend_`r'' found, and wscorr() not specified" qui count if !mi(`yvar_`r'',`yvar_`s'') if r(N)==0 { `warn' "`message' - not a problem as `yvar_`r'', `yvar_`s'' are never jointly observed" qui gen `covvar'`r'`s' = 0 // just to satisfy the pos def check local Svar_`r'_`s' `covvar'`r'`s' } else { di as error "Error: `message'" exit 459 } } if "`wscorr'"!="" & `s'!=`r' { if (`oksr'|`okrs') & mi("`wscorrforce'") local wscorrunused 1 else local wscorrused 1 } } } if "`warningnoest'"!="off" { if "`wscorr'"=="riley" { di as text "Note: using Riley's overall correlation model" _c if `wscorrunused' di " (ignoring covariances)" _c di } else if "`wscorr'"!="" { if `wscorrused' & `wscorrunused' `warn' "wscorr(`wscorr')" as text " used for only some covariances" if !`wscorrused' & `wscorrunused' `warn' "wscorr(" as result "`wscorr'" as text ") not used" if `wscorrused' & !`wscorrunused' di as text "Note: wscorr(" as result "`wscorr'" as text ") used for all covariances" } } ************* CONVERT VARIABLES TO MATRICES ************* tempname xi ymat Smat Xmat local i 0 tempname bottom * set up id as char tempvar idstring if mi("`id'") gen `idstring' = strofreal(_n) else { cap confirm string var `id' if !_rc gen `idstring' = `id' // string else { if !mi("`: value label `id''") decode `id', gen(`idstring') else gen `idstring' = strofreal(`id') } } forvalues obs=1/`N' { if `touse'[`obs'] { local ++i mat `ymat'`i'=J(1,`p',0) local thisid = `idstring'[`obs'] // 20jul2015 local identifier_text = cond(mi("`id'"),"observation ","") local identifier_result = cond(mi("`id'"),"`obs'", "`id'=`thisid'") mat rownames `ymat'`i' = "`thisid'" mat colnames `ymat'`i' = `ylist' mat `Smat'`i'=J(`p',`p',0) mat rownames `Smat'`i' = `ylist' mat roweq `Smat'`i' = "Study `i'" mat colnames `Smat'`i' = `ylist' forvalues r=1/`p' { mat `ymat'`i'[1,`r'] = `yvar_`r''[`obs'] forvalues s=`r'/`p' { * FILL MATRIX OF (CO)VARIANCES * detect errors local misscase = 1*missing(`yvar_`r''[`obs']) + 2*missing(`yvar_`s''[`obs']) local missdesc1 as result "`yvar_`r''" as text " is" local missdesc2 as result "`yvar_`s''" as text " is" local missdesc3 as result "`yvar_`r''" as text " and " as result "`yvar_`s''" as text " are" if !`misscase' & missing(`Svar_`r'_`s''[`obs']) { if `r'==`s' di as error "Error at `identifier_text'`identifier_result': var(`yvar_`r'') is missing but `yvar_`r'' is non-missing" else di as error "Error at `identifier_text'`identifier_result': cov(`yvar_`r'',`yvar_`s'') is missing but `yvar_`r'' and `yvar_`s'' are non-missing" global MVMETA_obserror `obs' exit 459 } if `misscase' & !missing(`Svar_`r'_`s''[`obs']) { if `s'>`r' & "`Svar_`r'_`s''"!="`covvar'`r'`s'" & "`wscorr'"!="riley" /// `warn' "at `identifier_text'" as result "`identifier_result'" /// as text ", " `missdesc`misscase'' /// as text " missing, so I'm ignoring non-missing " /// as result "`Svar_`r'_`s''" if `s'==`r' /// `warn' "at `identifier_text'" as result "`identifier_result'" /// as text ", " as result "`yvar_`r''" /// as text " is missing, so I'm ignoring non-missing " /// as result "`Svar_`r'_`r''" } local Svalue = `Svar_`r'_`s''[`obs'] if `r'==`s' & `Svalue'==0 { di as error "Error at `identifier_text'`identifier_result': zero variance in `Sstub'`yend_`r''`yend_`s''" global MVMETA_obserror `obs' exit 459 } mat `Smat'`i'[`r',`s']=`Svalue' mat `Smat'`i'[`s',`r']=`Svalue' } } * create X matrices for reml local Xrownames forvalues j=1/`p' { mkmat `xvars_`j'' `xcons' in `obs', matrix(`xi') mat colnames `xi' = `xvars_`j'' `xconsname' mat coleq `xi' = "`yvar_`j''" if `j'==1 { mat `Xmat'`i'=`xi'' } else if "`parmtype'"!="common" { local Xmatrows = rowsof(`Xmat'`i') local Xmatcols = colsof(`Xmat'`i') mat `Xmat'`i' = (`Xmat'`i', J(`Xmatrows',1,0)) mat `bottom' = (J(rowsof(`xi''),`Xmatcols',0), `xi'') mat roweq `bottom' = "`yvar_`j''" mat rownames `bottom' = `xvars_`j'' `xconsname' mat `Xmat'`i' = (`Xmat'`i' \ `bottom') } else mat `Xmat'`i' = (`Xmat'`i', `xi'') mat coleq `Xmat'`i' = "Study `i'" } mat colnames `Xmat'`i' = `ylist' } } ************* AUGMENT AND CHECK VARIANCES ************* local i 0 forvalues obs=1/`N' { if `touse'[`obs'] { local ++i local thisid = `idstring'[`obs'] // 20jul2015 local identifier_text = cond(mi("`id'"),"observation ","") local identifier_result = cond(mi("`id'"),"`obs'", "`id'=`thisid'") // OPTIONALLY AUGMENT MISSING VALUES IN MATRICES local augmented 0 if "`augment'"=="augment" { forvalues r=1/`p' { forvalues s=`r'/`p' { if (missing(`ymat'`i'[1,`r']) | missing(`ymat'`i'[1,`s'])) { if `s'==`r' { qui summ `Svar_`r'_`s'' if `touse' local Svalue = `missvar'*r(max) mat `Smat'`i'[`r',`s'] = `Svalue' if "`augquiet'"=="" di as text "Note at `identifier_text'" as result "`identifier_result'" as text ": setting var(" as result "`yvar_`r''" as text ") to " as result `Svalue' as text " because " as result "`yvar_`r''" as text " is missing" } else { if "`augquiet'"=="" di as text "Note at `identifier_text'" as result "`identifier_result'" as text ": setting cov(" as result "`yvar_`r''" as text "," as result "`yvar_`s''" as text ") to " as result "0" as text " because " as result "`yvar_`r''" as text " or " as result "`yvar_`s''" as text " is missing" mat `Smat'`i'[`r',`s']=0 mat `Smat'`i'[`s',`r']=0 } local augmented 1 } } if missing(`ymat'`i'[1,`r']) { if "`augquiet'"=="" di as text ""Note at `identifier_text'" as result "`identifier_result'" as text": setting missing " as result "`yvar_`r''" as text " to " as result "`missest'" mat `ymat'`i'[1,`r'] = `missest' local augmented 1 } } } // CHECK VAR-COV MATRICES ARE POSITIVE DEFINITE if `p'>1 & !matmissing(`Smat'`i') { * cap varcheck `Smat'`i', check(psd) `psdcrit' // removed 20jan2022 if `psdcrit'>0 local tol tol(`psdcrit') // added 20jan2022 cap _getcovcorr `Smat'`i', check(psd) `tol' // added 20jan2022 if _rc==506 { if `augmented' local augtext "augmented " else local augtext if "`posdef'"=="" { // error di as error "Error: `augtext'variance-covariance matrix not positive semi-definite at `identifier_text'`identifier_result':" _c mat l `Smat'`i' global MVMETA_obserror `obs' exit 459 } else { // just a problem `warn' "`augtext'variance-covariance matrix not positive semi-definite at `identifier_text'`identifier_result':" _c if "`warningnoest'" != "off" mat l `Smat'`i' `warn' "this may cause estimation to fail" _new } } else if _rc { di as error "Error in varcheck code" exit _rc } } } } if "`estimates'"!="noestimates" { local warn warn`warningest' ************* SET UP GLOBALS ************* * needed in all subroutines including mvmeta_bscov_*.ado, mvmeta_mufromsigma.ado, mvmeta_lmata.ado local things things ystub Sstub ymat Smat Xmat p n N xcons bsest ylist yends parmtype wscorr 2pi bscovariance sigma0 sigma0exp touse fixedparms quick taulog id penalty forvalues r=1/`p' { local things `things' yvar_`r' xvars_`r' forvalues s=`r'/`p' { local things `things' Svar_`r'_`s' } } foreach thing in `things' { global MVMETA_`thing' ``thing'' } ******* START ESTIMATION ********* tempname Sigma local nparms_aux 0 local neqs_aux 0 ******* FIXED-EFFECT METHOD ******* if "`bsest'"=="fixed" { tempname yhat mat `Sigma' = J(`p',`p',0) mvmeta_mufromsigma, sigma(`Sigma') yhat(`yhat') makeq `debug' // cap noi removed, 21jul2015 if !_rc { local rl0 = r(rl) // ll0 is ll for FE model - needed for testsigma option local ll0 = r(ll) local Qscalar = r(Qscalar) local Qscalardf = r(sum_p)-`fixedparms' local return_scalar rl0 ll0 Qscalar Qscalardf } local success = !missing(`Qscalar') } ******** KNOWN-SIGMA METHOD - NEEDED??? ******* else if "`bscovariance'"=="equals" { mvmeta_mufromsigma, sigma(`sigma0') if "`bsest'"=="reml" local ll = r(rl) if "`bsest'"=="ml" local ll = r(ll) matrix `Sigma' = `sigma0' local return_scalar ll local success = !matmissing(`Sigma') } ******* METHOD OF MOMENTS ******* else if inlist("`bsest'", "mm1", "mm2") { mvmeta_bscov_`bscovariance' if `touse', `bsest' `cholnames' `trunc' `mmfix' `debug' local success = !matmissing(r(Sigma)) if `success' { if r(negevals)>0 { local plural = cond(r(negevals)>1,"s","") local thesehave = cond(r(negevals)>1,"these have","this has") di as text _new "Note: Sigma has " as result r(negevals) as text " negative eigenvalue`plural'" _c if "`trunc'"=="" di as text " - `thesehave' been set to 0" else di as text _newline "Note: negative eigenvalue`plural' may lead to problems - if so, drop the notrunc option" } mat `Sigma' = r(Sigma) mvmeta_mufromsigma, sigma(`Sigma') `ifdebug' mat l `Sigma', title(MM: Sigma) * next 4 lines were previously only for mm2 tempname Qa Qb Q matrix `Qa' = r(Qa) matrix `Qb' = r(Qb) matrix `Q' = r(Q) } } ***************** REML / ML ***************** else if inlist("`bsest'","reml","ml") { if "`bscovariance'"!="equals" { ***************** HANDLE SPECIFED COVARIANCE STRUCTURE ***************** tempname init mvmeta_bscov_`bscovariance' if `touse', setup start(`start') `cholnames' matrix `init' = r(init) local eqlist2 `r(eqlist)' local nvarparms `r(nvarparms)' local neqs_aux = `r(neqs_aux)' ***************** OPTIONALLY DISPLAY STARTING VALUES FOR Sigma ***************** if "`showstart'"=="showstart" { di as text _newline "Starting value for beta's:" _c mat l r(binit), noheader di as text _newline "Starting value for Sigma (between-studies variance):" _c mat l r(Sigma), noheader if !mi("`debug'") { di as text _newline "All starting values:" _c mat l `init', noheader } } ****************** CREATE MEAN EQUATIONS ****************** if "`parmtype'"=="long" { /* one equation per mean: "long parameterisation" */ forvalues r=1/`p' { local eqlist1 `eqlist1' (`yvar_`r'': `xvars_`r'', `constant') } } else if "`parmtype'"=="short" { /* one equation comprising all means */ local eqlist1 (Overall_mean: forvalues r=1/`p' { local eqlist1 `eqlist1' `yvar_`r'' } local eqlist1 `eqlist1', nocons) } else if "`parmtype'"=="common" { /* common equation for each means */ local eqlist1 (`yvar_1': `xvars_1', `constant') } * MODIFICATION FOR QUICK METHOD if "`quick'"=="quick" matrix `init' = `init'[1,(`fixedparms'+1)...] else local eqlist1a `eqlist1' ****************** MAXIMISE ***************** `ifdebug' di as text "*** starting to maximise ..." /* collinear option makes the short parameterisation work even when some b's are collinear */ /* neq(`p') fails here, not sure why */ local mlcommand ml model d0 mvmeta_lmata `eqlist1a' `eqlist2' if `touse', /// obs(`n') collinear init(`init', copy) maximize nooutput `mlopts' /// `iterate' nopreserve missing `log' `constraints' `ifdebug' di as input `"`mlcommand'"' `mlcommand' `ifdebug' di as text "*** maximisation is finished ..." **************** FINISH OFF FOR QUICK METHOD *************** if "`quick'"=="quick" { mvmeta_bscov_`bscovariance' if `touse', postfit `cholnames' mat `Sigma'=r(Sigma) global MVMETA_quick // not needed? tempname varparms mat `init'=e(b) mvmeta_mufromsigma, sigma(`Sigma') `ifdebug' di as text "Results not allowing for uncertain Sigma..." `ifdebug' ereturn display, `cformat' `pformat' `sformat' mat `init'=(e(b),`init') `ifdebug' mat l `init', title(Starting values for 2nd -ml- run) `ifdebug' di as text "Estimating variance allowing for uncertain Sigma..." qui ml model d0 mvmeta_lmata `eqlist1' `eqlist2' if `touse', obs(`n') collinear init(`init') maximize nooutput `mlopts' nopreserve missing `log' iterate(0) `constraints' } local success = e(converged) } ****************** COMPUTE Sigma ***************** mvmeta_bscov_`bscovariance' if `touse', postfit `cholnames' mat `Sigma'=r(Sigma) local nparms_aux = r(nparms_aux) forvalues r=1/`p' { forvalues s=1/`r' { local Sigma`s'`r' `r(Sigma`s'`r')' } } **************** ERETURN FOR BOTH METHODS *************** forvalues r=1/`p' { forvalues s=1/`r' { local return_local `return_local' Sigma`s'`r' } } } ************ ERETURN RESULTS AND TIDY UP *************** foreach thing in `things' { global MVMETA_`thing' } * nparms_* refers to numbers of parameters * neqs_* refers to numbers of equations tempname b V mat `b'= e(b) mat `V'= e(V) local nparms_mean = colsof(`b') - `nparms_aux' local neqs_mean = cond("`parmtype'"=="long",`p',1) local k_eform = `neqs_mean' // system name picked up by _coef_table in Replay } // noestimates resumes here local warn warn`warningnoest' * make combined data matrices tempname ydata Sdata Xdata forvalues i=1/`n' { mat `ydata' = nullmat(`ydata') \ `ymat'`i' mat `Sdata' = nullmat(`Sdata') \ `Smat'`i' mat `Xdata' = nullmat(`Xdata') \ `Xmat'`i'' // NB transpose } local return_scalar `return_scalar' neqs_mean neqs_aux nparms_mean nparms_aux k_eform fixedparms success foreach thing of local return_scalar { if "``thing''"!="" { return scalar `thing' = ``thing'' local scalarlist `scalarlist' `thing' } } return local scalarlist `scalarlist' local return_matrix Sigma Q Qa Qb b V ydata Sdata Xdata foreach thing of local return_matrix { cap confirm matrix ``thing'' if !_rc { return matrix `thing' = ``thing'' local matrixlist `matrixlist' `thing' } } return local matrixlist `matrixlist' foreach thing of local return_local { if "``thing''"!="" { return local `thing' ``thing'' local locallist `locallist' `thing' } } return local locallist `locallist' `ifdebug' di as input "End of mvmeta_estimate" end *=========================== END OF MVMETA_ESTIMATE PROGRAM =========================== *=========================== START OF MVMETA_GETDATA PROGRAM =========================== /* Reads in data from e(ydata), e(Sdata) Estimation results are appended and are flagged by variable `type' Options nouv, nomv suppress this */ program define mvmeta_getdata syntax namelist, type(name) [corr nouv nomv STUDylabel(string) UVLABel(string) MVLABel(string) wt] if mi("`studylabel'") local studylabel "Studies" if mi("`uvlabel'") local uvlabel "Pooled uv" if mi("`mvlabel'") local mvlabel "Pooled mv" drop _all tempname ydata y Sdata S value mat `ydata' = e(ydata) mat `Sdata' = e(Sdata) if !mi("`wt'") { tempname wtdata mat `wtdata' = e(wt) mat coleq `wtdata' = "" } local newobs = e(N) + 2 qui set obs `newobs' qui gen _id = "" qui gen `type' = "" foreach var in `namelist' { qui gen `var' = . qui gen se_`var' = . if !mi("`wt'") qui gen wt_`var' = . if !mi("`corr'") { foreach var2 in `namelist' { if "`var2'" <= "`var'" continue qui gen corr_`var'_`var2' = . } } } forvalues i=1/`newobs' { if `i'<=e(N) { local first = e(dims)*(`i'-1)+1 local last = e(dims)*`i' mat `y' = `ydata'[`i',1...] local studyname : rowname `y' qui replace _id = "`studyname'" in `i' mat `S' = `Sdata'[`first'..`last',1...] mat roweq `S' = "" qui replace `type' = "`studylabel'" in `i' } if `i'==e(N)+1 { // uv estimate if !mi("`uv'") continue mat `y' = e(b_uv) mat `S' = e(V_uv) qui replace `type' = "`uvlabel'" in `i' } if `i'==e(N)+2 { // mv estimate if !mi("`mv'") continue mat `y' = e(b) mat `S' = e(V) qui replace `type' = "`mvlabel'" in `i' } if `i'>e(N) { mat coleq `y' = "" mat coleq `S' = "" mat roweq `S' = "" } foreach var in `namelist' { * mat `value' = `y'[1,"`var'"] * qui replace `var' = `value'[1,1] in `i' * mat `value' = `S'["`var'","`var'"] * qui replace se_`var' = sqrt(`value'[1,1]) in `i' * if `i'<=e(N) & !mi("`wt'") { * mat `value' = `wtdata'[`i',"`var'"] * qui replace wt_`var' = `value'[1,1] in `i' * } local r = colnumb(`y',"`var'") qui replace `var' = `y'[1,`r'] in `i' qui replace se_`var' = sqrt(`S'[`r',`r']) in `i' if `i'<=e(N) & !mi("`wt'") { qui replace wt_`var' = `wtdata'[`i',`r'] in `i' } } if !mi("`corr'") { foreach var in `namelist' { local r = colnumb(`y',"`var'") foreach var2 in `namelist' { if "`var2'">="`var'" continue * mat `value' = `S'["`var'","`var2'"] * qui replace corr_`var2'_`var' = `value'[1,1] / (se_`var'*se_`var2') in `i' local s = colnumb(`y',"`var2'") qui replace corr_`var2'_`var' = `S'[`r',`s'] / (se_`var'*se_`var2') in `i' } } } } if !mi("`mv'") qui drop if _n==e(N)+2 if !mi("`uv'") qui drop if _n==e(N)+1 end *=========================== END OF MVMETA_GETDATA PROGRAM =========================== *======================= START OF MVMETA_PI PROGRAM ==================== /* IW 14oct2021 now works also with long parameterisation added warning after metaregression nicer output IW 8apr2020 Form prediction intervals after mvmeta Assumes df = total sample size minus 2: this matches univariate advice */ program define mvmeta_pi syntax [, format(string) level(cilevel) xvar(string)] if mi("`format'") local format %9.0g if mi("`xvar'") local xvar _cons if e(cmd)!="mvmeta" { di as error "mvmeta was not the last command run" exit 301 } * detect meta-regression local metareg 0 forvalues i=1/`e(dims)' { if !mi(e(xvars_`i')) local metareg 1 } local plevel = (100+`level')/200 local col1 _col(1) forvalues i=2/6 { local pos = 1+12*(`i'-1) local col`i' `format' _col(`pos') } if `metareg' local metaregnote for the meta-regression intercepts di as txt _n "Table of prediction intervals `metaregnote'" _n "{hline 70}" di `col1' as txt "Outcome" `col2' "Estimate" `col3' "`level'% Confidence Int." `col5' "`level'% Prediction Int." tempname Sigma mat `Sigma'=e(Sigma) local df = e(N)-2 foreach yvar in `e(yvars)' { if e(parmtype)=="short" { local est = _b[`yvar'] local se = _se[`yvar'] } else { local est = [`yvar']_b[`xvar'] local se = [`yvar']_se[`xvar'] } local pme = invnorm(`plevel') * `se' * local pmp = invt(`df',`plevel') * sqrt(`se'^2 + `Sigma'["`yvar'","`yvar'"]) // old way for stata17 local r = colnumb(`Sigma',"`yvar'") local pmp = invttail(`df',1-`plevel') * sqrt(`se'^2 + `Sigma'[`r',`r']) // 7apr2022 to allow Stata v12 di `col1' as txt "`yvar'" `col2' as res `format' `est' `col3' `format' `est'-`pme' `col4' `format' `est'+`pme' `col5' `format' `est'-`pmp' `col6' `format' `est'+`pmp' } di as txt "{hline 70}" di as txt "Note: using N-2=" as res `df' as txt " degrees of freedom" end *=========================== END OF MVMETA_PI PROGRAM =========================== *======================= START OF MATA ROUTINE FOR RANDFIX ==================== mata: void selectpart(string scalar cname, string scalar Vinname, string scalar Voutname ) { c=st_matrix(cname) V=st_matrix(Vinname) Vpart = select(V,c') Vpart = select(Vpart,c) st_matrix(Voutname,Vpart) } end *======================= END OF MATA ROUTINE FOR RANDFIX ==================== *======================= ROUTINES FOR WARNINGS ==================== program define warnerror foreach thing in txt text result res { local 0 : subinstr local 0 "as `thing' " "as error ", all } di as error "Warning: " `0' end program define warntext di as text "Warning: " `0' end program define warnoff local y end *======================= END OF ROUTINES FOR WARNINGS ==================== /****************************************************************************** HISTORY OF VERSIONS 1 & 2 version 2.14 Ian White 1apr2015 data, uv estimates and fixed estimates all returned as e() version 2.13.1 Ian White 27mar2015 removed noestimates from forest & bubble calls version 2.13 Ian White 14mar2015 randfix - improved re-call to mvmeta; without arguments defaults to all y's pbest - new mcse suboption; checks for non-missing covariates; add isid check; bug fix - rename() got confused if new names intersected with old names help file: specifies bscov(exch 0.5) as an option in the covariances list. now advises against -augment- option with wscorr() randfix, nowt and wt documented fixed bugs: with if/in, gave warning about dropping the observations not satisfying if/in eform option worked on replay but not on initial call improved error message "Error: yb and yc have no jointly observed values" in sparse NMA with bscov(uns) 12mar2015: copes with only one observation and meta-regression 13mar2015: bug fix in mvmeta_makecovs: if one name contains another, covariance is wrongly found 13mar2015: equations(yvarlist:...) now allowed; printing of equations improved 14mar2015: bscov(prop|eq|corr) checks correct dimension of matrix version 2.12 10feb2015 weights and bos automatically reported unless nowt specified [how to get formats in?] new id option: passed to wt() pbest() and forest() new option wt replaces old bos with suboptions: sd - score decomposition method (default) ) use rv - relative variance method ) one of dpc - data-point coefficients ) these details - more info for SD & RV clear - saves details (with sd details) format - for tables gen() - generates variables for weights (with sd) keepmat() - creates matrices id() - id for output smaller: recoded eform[()] using eform eform2() noestimates (replaces dryrun) allowed in Estimate: it proceeds to end but doesn't set e(bsest); Replay picks this up and suppresses estimates & weights. version 2.11.3 2feb2015 add row equation names to X matrices version 2.11.2 22dec2014 bos option is now a revised subroutine borrow option removed: now use bos(rv) keepmat for X: improved rowname for constant forest option also available as a subroutine version 2.11.1 19dec2014 wscorrforce option makes wscorr() override existing corrs (useful for mvmeta_runuv and mvmeta_forest) version 2.11 9dec2014 bug fixed - observations with missing covariates were deleted, even if outcome was also missing now also gives warning when observations are deleted due to missing covariates version 2.10.1 8sep2014 bscov(exch #) gives correct description from line 504 new suboption: pbest(, tabdispoptions()) version 2.10.0 29jul2014 rename suboption for pbest version 2.9 10jun2014 drops variables `y'* if no non-missing data new bscovariance(exch #) - designed for network meta where variables may be dropped version 2.8.2 6jun2014 bug fix in pbest version 2.8.1 30may2014 Commonparm note "The above coefficients also apply to the following equations" suppressed in only 1 equation version 2.8 15may2014 Changes to pbest: moved to separate program and restructured; new suboption clear implemented Aurelio Tobias' graphs: new suboptions bar, line, cum, graph_options version 2.7.3 6feb2014 bug fix - bos and borrow work together (NB new options need to be inserted carefully) bos works for incomplete data (by simple augmentation) version 2.7.2 5feb2014 bug fix - bos works after wscorr() - fixed by new e(wscorr) (replaces e(wscorrtype)) bos corrected to work after wscorr(riley) version 2.7.1 3feb2014 (on BSU website) fixed bug which limited length of e(xvars_*) if run in v12 version 2.7 27jan2014 Changes to incorporate new bos option: bos has optional suboptions new dryrun option version 2.6 3jan2014 (on BSU website) Changes for consistency with network.ado: pbest fails with parmtype=common: better error message new undocumented network option identifies calls from -network- new stripprefix and zeroname suboptions for pbest pbest(,reps(1000)) now default version 2.5.6 20jul2013 added standard options: cformat pformat sformat version 2.5.5 21mar2013 - ON BSU WEBSITE Anna Chaimani's fix to bug that made pbest fail with >10 treatments version 2.5.4 20mar2013 - ON BSU WEBSITE New saving() suboption on pbest version 2.5.3 6mar2013 Corrected showall again version 2.5.2 28aug2012 tries to continue even if fixed-effect method fails version 2.5.1 2jul2012 corrected small bug in pbest(,gen()) version 2.5 7jun2012 Renamed ereturned parameters more systematically as nparms_* and neqs_* Corrected showall eform option disabled with showall and parmtype long - previously only the 1st eq was exponentiated. The problem here is that e(k_eform) must be set to get correct eform behaviour, but it is lost when Replay posts new estimation results, which it must do to use dof() option; and Replay isn't eclass so can't re-set e(k_eform) version 2.4.2 1feb2012 - ON BSU WEBSITE "Warning: method of moments failed - I2 statistic not available" not now in red: it worried users unnecessarily Notes on variables not jointly observed now appear only for bscov(uns) Formatting of output notes improved Clearer error message if randfix() specifies variables that aren't outcome variables fixed bug making borrow fail if randfix() was also specified version 2.4.1 13jan2012 abbreviation common for commonparm version 2.4 20dec2011 new commonparm option which forces same coeffs in each equation mvmeta syntax unchanged but global storage changed: MVMETA_longparm is now conveyed by MVMETA_parmtype = long/short/common version 2.3.4 20dec2011 new constraints option for constrained maximisation new all option for pbest reports all ranks, not just the best new collinear option suppresses check for collinearity - use with care! version 2.3.3 15nov2011 new mm2 option (undocumented) version 2.3.2 15jun2011 new mmfix option sets unestimated off-diagonal Sigma terms to 0 version 2.3.1 15jun2011 tidied up output: displays either yvarlist or regressions, not both fixed now implies nojointcheck jointcheck flags error/warning with 0/1 jointly observed value (was 0 only) version 2.3.0 2jun2011 same as version 2.3 but with borrow, randfix and qscalari removed PUBLISHED WITH SJ 2011;11:255--270 version 2.3 31may2011 - ON BSU WEBSITE corrected the constant terms in log L and log RL, as suggested by Antonio Gasparrini: affects L/RL (only) in models with missing outcomes and models with covariates i2: when a variance expression is shortened due to zero parameters in the cholesky decomposition, the shorter expression is now carried forward to compute the correlations (avoids failure) version 2.2 6may2011 pbest has new suboption predict, giving prediction (for new study) not estimation (of overall mean) clearer error when pbest is called after fitting with parmtype=short i2 refuses when e(Q) etc. weren't set no checks of jointly observed data with bscov(eq|prop) iterate(#) no longer fails with quick option new bscov(CORRelation ) version 2.1.3 24mar2011 -borrow- output clarifies it compares squared std errors version 2.1.2 10mar2011 fixed bug in i2 without -heterogi-: negative i2 wasn't truncated. output tables widen to accommodate long variable names fixed bug in randfix: failed after if/in version 2.1.1 8mar2011 corrected borrow option with bscov(equals ..) testsigma now divides P-values by 2 if LRT>0. Should Sigma0 be stored?? version 2.1 26feb2011 new options: Borrow Randfix(varlist) Qscalar eqnames in e(b) and e(V) after fixed/mm now match those after ml/reml dropped e(yends), e(xvars) new e(cmdline), e(ystub), e(Sstub), etc. version 2.0.5 8mar2011 same as version 2.1.1 but with borrow, randfix and qscalar removed - resubmitted to SJ version 2.0.4 1feb2011 added drawbeta.ado into this file mvmeta_mufromsigma gives useful error message with unidentified model collinearity check is now performed in the subset with observed outcome version 2.0.3 19jan2011 corrected minor bug in mvmeta_mufromsigma that truncated equation names if name lengths totalled >256 characeters version 2.0.2 6dec2010 enabled ncchi2 option to pass to -heterogi- version 2.0.1 i2 output changed from "yB/yC" to "yB & yC" testsigma option added named version 2.0 14oct2010 version 1.10 14oct2010 best renamed pbest version 1.9 8sep2010 -corr(riley) mm- now fails; previously ran as -corr(0) mm- MM moved to covariance files covariance files don't require p() corr() renamed as wscorr() (but undocumented corr() still works) bscov|bscorr changed to print(bscov|bscorr) covariance() renamed to bscovariance() e(Sigma11) etc. now includes zero terms, so se and i2 can give funny results new undocumented noestimates option (suppresses coeff table) se moved within i2 option, and table reformatted i2 option: tau^2 and I^2 have their CIs computed on the same scale new ciscale(SD|logSD|logH) option determines which scale this is new i2fmt option formats I^2 new undocumented nojointcheck option suppresses check for jointly observed values version 1.8.2 18aug2010 fixed bug in Replay call that ignored nouncertainv on estimation run version 1.8.1 2aug2010 fixed bug in local ylist0 arising with long variable list (thanks to Stephen Kaptoge) version 1.8 29jun2010 CI for I2: reml/ml: improve by computing it on scale of log(H)=-0.5*log(1-I^2) fixed/mm: implemented by a call to -heterogi- CI for between-studies SD (on log scale) and correlation (on Fisher's transformation scale) undocumented quick option (re/ml): mu estimated from sigma using closed form expression within loglik routine, not as separate parameters tidied up se output version 1.7.1 3jun2010 drops any observations with all missing outcomes version 1.7 27may2010 Allow structured variance matrices: new option covariance() to echo -xtmixed- syntax. Can have: covariance(unstructured) covariance(proportional matrixname) New loglik and mufromsigma routines handle missing data - matrices are augmented only with new augment option - positive definite check now restricted to fully non-missing variance matrices Only uses mata routine for loglik Some improved error messages bscorr - correlations changed to . above diagonal Internal changes: new external files mvmeta1_mufromsigma.ado (now using mata) mvmeta1_covariance_unstructured.ado mvmeta1_covariance_proportional.ado version 1.6.1 26mar2010 _coef_table version 1.2.5 fails, 3.0.0 works. Changed to use _coef_table only when needed (i.e. withoptions eform and showchol) and to revert to -ereturn display- if necessary version 1.6 4jan2010 corr(riley) option NB this is the last version named mvmeta1 version 1.5.2 18nov2009 missing option added to ml command - previously gave wrong results with missing estimates version 1.5.1 18nov2009 fixed se option again (!) - instead of ignoring terms < a fixed value, it ignores terms that make nlcom crash new undocumented selist option shows what se is calculating - may be useful for doing more extended nlcom version 1.5 9oct2009 extensively tidied up code nopreserve option on ml dof(string) option se option fixed new i2 option new mata option uses new mvmeta1_lmata.ado - 2-5 times faster? version 1.4 29sep2009 re-written as Estimate and Replay non-interactive -ml- allows all standard maximization options (via -mlopts-) version 1.3 25-29sep2009 (Stephen K's suggestions) deleted -preserve- command at start of setup - e(sample) is now correct obs with incomplete covariates dropped corrected reml expression with unequal equations (eq() option) fixed bug that caused crash with long variable names (cholesky equation names could be over 32 characters): cholesky parameters are now named 11,12 etc unless new cholnames is used version 1.2 11aug2009 maximizeoptions() replaces trace difficult iterate() version 1.1 10aug2009 eform correctly implemented Ignores eform option if long parameterisation is used without covariates redisplays results if typed without arguments! output now occurs in same place in program for all methods (incl. redisplay) new noposdef option allows non-positive-semidefinite var/cov matrices criterion for pos semidef changed from min eval<0 to min/max<-1E-8 (can change with new psdcrit() option) allows separate regressions via equations() option the only maximize options allowed are difficult iterate() version 1.0 12jun2009 allows covariates in reml/ml new noconstant option start defaults to mm program no longer returns e(Mu) /******************************************************************************