/********************************************************************* *! version 4.0.2 # Ian White # 21apr2022 last regression doesn't contaminate ereturn-ed values version 4.0 # Ian White # 07apr2022 version number changed to match mvmeta new release to UCL, github and SSC version 0.23.1 Ian White 7apr2022 reordered options to match help file version 0.23 Ian White 5apr2022 add option countby: written by S Kaptoge (Apr2009, Aug2016) to allow extraction of study-specific counts of events, participants, and person-time of follow up by a variable specified in added option countby(varname) version 0.22.2 Ian White 5apr2022 extended usecoefs to work without eqname (also corrected help) version 0.22.1 Ian White 4apr2022 added dicmd subroutine version 0.22 Ian White 28mar2022 major edit to enable prefix syntax small changes to output debug option enhanced 20jan2022: added but commented out possible change to pos def checking version 0.21.1 Ian White 2nov2021 quietly suppresses all output version 0.21 Ian White 22jun2015 by(varname) changed to by(varlist) multiple equations rewritten bug fix - some covariances were omitted new parsing from e(b) new useeqs() mlogit no longer requires baseoutcome() new clear option updated to v11 (makes -mixed- work better) NB to set up post file I have chosen the expensive option of fitting the full model in the whole data set: - using covariates because (a) stcox returns no e(b) without covariates (b) mixed needs to be correctly parsed - using full data set because covariates might have pp in 1st study but user can specify a subset using learnif() improved augmentation - uses all outcome levels in all subsets (not in the specific subset) new options: infix, longnames version 0.20 30may2014 esave(junk) returns missing _e_junk instead of failing (Phil Perry's request from 2011) preserve before augment to avoid changing data set if crash while augmented robust and cluster(varname) now disable augmentation - I tried changing iweight to pweight, but it gives bad (non-large) variances multiple-equation commands supported - command is run without covariates in first study before loop in order to identify equation names and set up post file - mlogit requires baseoutcome() - warning: I've found an example where mlogit does strange things after augmentation version 0.19 20apr2011 - Stephen Kaptoge's edits to esave() version 0.18 18mar2011 ON WEBSITE removes unwanted output with hard option version 0.17 22dec2009 - new counts() option version 0.16 22dec2009 corrected bug that gave wrong _e_N_sub version 0.15 21dec2009 corrected bug that caused augmentation to fail for clogit when matched sets weren't numbered 1, 2, 3, ... version 0.14 5oct2009 fixed bug causing crash with too many stratifying levels augmentation for clogit puts augmented values into new groups - this works version 0.13 25sep2009 augmentation for clogit is stratified by group variable - but still fails fixed bugs stopping filenames with spaces fixed bug in _augment that caused "Obs. nos. out of range" error with very long variable list now an eclass command returning e(sample) from each separate command version 0.12 3feb2009 Augmentation respects a weight that has been stset Option saving(filename, replace) now works (previously gave obscure "invalid `clear'" error -stset- features are preserved even if the stcox model crashes version 0.11 26sep2008 changes from Stephen K: made eclass and enabled esave(N_sub N_fail) with logistic regressions version 0.10 22aug2008 AS PUBLISHED IN SJ 9(1):40-56 varcheck rewritten to cope with Stephen K's very strange example in Stata 9 coef, or, nohr identified as separate options to be included in replay for hard new pause option decided not to capture augmented analysis as I've never had problems version 0.9 8aug2008 gives error when usevars() contains variables not in xvarlist bug fixed when usevars() is just one variable (PP was missed) new esave - adds selected e() to dataset (similarly undocumented rsave, for Stephen K) nodetails also suppresses augmentation details noauglist suppresses listing of augmenting observations ppcmd allows alternative command (e.g. plogit) as alternative to augmentation ppfix(all) uses augmentation/ppcmd unconditionally hard captures results of first model fit & augments if error works with blogit (but augmentation doesn't) remove eqnames for any command, not just clogit augment() renamed augwt() nocons included as explicit option thus getting #parms right in _augment and disabling useconstant new output variables _ppfixed and _ppremains version 0.8 17mar2008 _getcovcorr replaced with varcheck (_getcovcorr's check option isn't in Stata 9) version 0.7.1 18dec2007 much improved augmentation checks for missing variables and zero variances early event for any st command zerovar, missvar dropped and replaced with error and missing value version 0.7 11dec2007 bugs fixed: single covariate -> perfect prediction always detected uselist default wrong much improved augmentation checks for missing variables and zero variances early event for any st command zerovar, missvar dropped and replaced with error and missing value version 0.6 19oct2007 uses augmented procedure (option augment() instead of hessadd()) version 0.5 18oct2007 attempts to detect when Stata has used a generalised inverse for V: if detected, it augments the Hessian by `hessadd'*var(X) before inverting (version 10 only) missing beta's are replaced by zero with variance `zerovar', covariance 0 NOTE: OK if there are no events in any category*study, but fails if there are no individuals in the reference category (because a new reference category is chosen) version 0.4 24sep2007 enable weights version 0.3 5jul2007 rename use, cons as usevars, useconstant; rename outputs as bvar1, Vvar1var2 etc. not b1, V11 etc.; rename singular option as zerovar and improve warning message version 0.2 28jun2007 handle clogit row/colnames: e(b), e(V) labelled y:x1, y:x2 etc. instead of x1, x2 etc.; append option; if/in bug removed Ideas: make it a prefix command? would still need to use weights, strata from main command PROBLEMS: if y is string, -mvmeta_make logit y x, ...- doesn't fail! (just creates empty data) *********************************************************************/ prog def mvmeta_make, eclass version 11 local mvoptions /// SAVing(string) replace append clear /// save-file options USEVars(varlist) USEConstant USEEqs(string) LEARNif(string) /// what results are stored USECOEfs(string) esave(string) counts(string) COLLapse(string) /// what results are stored NAMEs(namelist min=2 max=2) infix(string) LONGnames keepmat /// how results are stored noDETails pause /// output PPFix(string) AUGwt(real 0.01) noAUGList PPCmd(string) /// perfect prediction behaviour hard /// estimation debug dryrun rsave(string) countby(varname) // undocumented options local cmdoptions STrata(passthru) noCONstant coef or nohr /// command options requiring special treatment robust cluster(passthru) BAseoutcome(passthru) /// command options requiring special treatment * // other command options * WHICH SYNTAX? * if first character after mvmeta_make is "," then we have prefix syntax, otherwise we have classic syntax if substr("`0'",1,1)=="," { local syntype prefix di as text "Prefix syntax detected" } else if substr("`0'",1,1)==" " { di as text "Program error" exit 498 } else { local syntype classic di as text "Classic (non-prefix) syntax detected" } * PARSE CLASSIC (NON-PREFIX) SYNTAX * detect multiple commas gettoken left right : 0, parse(",") bind local right: subinstr local right "," "" gettoken left right : right, parse(",") bind if !mi("`right'") di as error "Probable syntax error: multiple commas found" if "`syntype'"=="classic" { syntax anything(equalok) [if] [in] [fweight aweight pweight iweight], /// by(varlist) [`mvoptions' `cmdoptions'] } * PARSE PREFIX SYNTAX else { gettoken prefixpart anything : 0, parse(":") local anything = trim(subinstr("`anything'",":","",1)) * parse the mvmeta_make options local 0 `prefixpart' syntax, by(varlist) [`mvoptions'] * parse the regcmd: might be prefix (main command after colon) or mixed (main command before ||) local prefixcmds bootstrap jackknife permute "mi estimate" while inlist(word("`anything'",1), "bootstrap", "jackknife", "permute", "mi") { // while a prefix command is detected gettoken one anything : anything, parse(":") local prefix `prefix' `one': local anything: subinstr local anything ":" "" } if word("`anything'",1) == "mixed" { local strpos = strpos("`anything'","||") local postfix = substr("`anything'",`strpos',.) local anything = substr("`anything'",1,`strpos'-1) } local 0 `anything' syntax anything(equalok) [if] [in] [fweight aweight pweight iweight], /// [`cmdoptions'] } * if postfix doesn't contain a comma, append one _parse comma lhs rhs : postfix if mi("`rhs'") local postfix `postfix', if !mi("`debug'") { di as input "Initial results of parsing:" foreach thing in prefix anything by if in weight postfix options { di as text _col(5) "`thing'" as result _col(16) "``thing''" } } local options `options' `robust' `cluster' `baseoutcome' *********************** PARSING *********************** * Find yvar and xvarlist gettoken regcmd regbody : anything, parse(" ") unabcmd `regcmd' local anything `prefix' `anything' if substr("`regcmd'",1,2)=="st" local xvarlist `regbody' else if "`regcmd'"=="mvreg" { gettoken yvar xvarlist : regbody, parse("=") local xvarlist : subinstr local xvarlist "=" "" } else gettoken yvar xvarlist : regbody if "`regcmd'"=="blogit" { gettoken nvar xvarlist : xvarlist } if "`regcmd'"=="mlogit" & mi("`baseoutcome'") { di as error `"mvmeta_make mlogit: baseoutcome() required"' exit 198 // otherwise -mlogit- might choose different baseoutcome in different studies } if inlist("`regcmd'","mixed","xtmixed") { gettoken xvarlist junk : xvarlist, parse("|") * local regbody `yvar' `xvarlist' } if inlist("`regcmd'","mlogit","mvreg") & !mi("`useeqs'") { di "useeqs() ignored" local useeqs } if !mi("`yvar'") unab yvar : `yvar' unab xvarlist : `xvarlist' // gives error message if factor variables used if "`debug'"=="debug" { di as text "y variables: " _col(21) as result "`yvar'" di as text "x variables: " _col(21) as result "`xvarlist'" } if "`weight'"!="" local wtexp [`weight'`exp'] if "`names'"!="" { tokenize "`names'" local y `1' local S `2' } else { local y y local S S } if "`keepmat'"=="keepmat" { local ymat `y' local Smat `S' } else { tempname ymat Smat } if "`details'"=="nodetails" local qui qui if !mi("`debug'") local noidicmd dicmd tempname onevalue post if "`hard'"=="hard" local hardcap capture * augmentation options if !mi("`robust'`cluster'") | "`weight'"=="pweight" { // Robust SE if "`ppfix'"=="" { local ppfix none di as error `"Robust standard error requested - ppfix set to none"' } else if inlist("`ppfix'", "check", "all") { di as error `"Robust standard error requested - ppfix(`ppfix') not allowed"' exit 198 } } if inlist("`regcmd'","regress","mvreg","mixed","xtmixed") { if !inlist("`ppfix'","","none") di as error "ppfix(`ppfix') not available with `regcmd'" local ppfix none } if "`ppfix'"=="" local ppfix check if !inlist("`ppfix'","none","check","all") { di as error `"Please specify ppfix(none|check|all)"' exit 198 } if "`ppcmd'"!="" { gettoken ppcmd1 ppcmd2 : ppcmd, parse(",") local ppcmd2=substr("`ppcmd2'",2,.) } if "`ppcmd'"=="" & `augwt'<=0 { local ppfix none } if "`append'"=="append" & "`replace'"=="replace" { di as error `"Can't have append and replace"' exit 198 } local eformopts `coef' `or' `hr' * check if user has saving(file, replace) tokenize `"`saving'"', parse(",") if "`3'"=="replace" { local replace replace local saving `1' } * sort out clear local savingorig `saving' if mi("`saving'") { if mi("`clear'") { di as error "Please specify saving() or clear" exit 198 } if !mi("`append'") { di as error "append ignored without saving()" local append } if !mi("`replace'") { di as error "replace ignored without saving()" local replace } if !mi("`clear'") { tempfile saving } } if "`useconstant'" == "useconstant" & "`constant'"=="noconstant" { di as error "Can't have both useconstant and noconstant options" exit 198 } if "`useconstant'" == "useconstant" & "`regcmd'"=="stcox" { di as error "useconstant option not allowed - stcox doens't have a constant" exit 198 } if !mi("`usecoefs'") { if !mi("`usevars'") di as error "usecoefs() specified: usevars() ignored" if !mi("`useconstant'") di as error "usecoefs() specified: useconstant ignored" if !mi("`useeqs'") di as error "usecoefs() specified: useeqs() ignored" if !mi("`learnif'") di as error "usecoefs() specified: learnif() ignored" } else if mi("`learnif'") local learnif 1 // by default, use all data to learn coefficients *********************** END OF PARSING *********************** if !mi("`debug'") { di as input "Final results of parsing:" foreach thing in anything if in wtexp postfix strata constant eformopts options { di as text _col(5) "`thing' " as result _col(16) "``thing''" } } *********************** START ANALYSIS: IF/IN AND BYVARS *********************** marksample touse, novarlist markout `touse' `yvar' `xvarlist' * Make summary by-variable local byvarlist `by' if wordcount("`by'")>1 { tempvar byvarname qui egen `byvarname' = group(`byvarlist') if `touse', label } else local byvarname `by' qui levelsof `byvarname' if `touse', local(byvarnamelevels) local nby : word count `bylevels' cap confirm numeric var `byvarname' local ischarbyvarname = _rc>0 * Create postfile expression for by-variables foreach byvar of local byvarlist { cap confirm string var `byvar' if _rc==0 { // string variables: find length tempvar length gen `length'=length(`byvar') summ `length', meanonly local bypost `bypost' str`=r(max)' } else { // numeric variables: any label? local label`byvar' : value label `byvar' local labels2save `labels2save' `label`byvar'' } local bypost `bypost' `byvar' } if !mi("`labels2save'") { tempfile labelsfile label save `labels2save' using `labelsfile' } *********************** SORT OUT USEVARS *********************** if "`usevars'"=="" { local usevars `xvarlist' } local usevars : list uniq usevars local p : word count `usevars' * check `usevars' are within `xvarlist' local errorvars : list usevars - xvarlist if !mi("`errorvars'") { di as error `"`errorvars' found in usevars() but not in xvarlist"' exit 498 } * assign `var`i'' locals forvalues r=1/`p' { local var`p' : word `r' of `usevars' } if "`useconstant'" == "useconstant" { local ++p local var`p' _cons local usevars `usevars' _cons } *** SORT OUT COUNTS *** foreach name in records subjects failures { foreach count in `counts' { if "`count'"==substr("`name'",1,length("`count'")) local counts2 `counts2' `name' } } local counts `counts2' local counts2 *********************** IDENTIFY COEFFICIENTS *********************** if !mi("`usecoefs'") { di as text "Using coefficients: " _c local ncoefs 0 foreach usecoef of local usecoefs { local ++ncoefs local pos = strpos("`usecoef'","]") if `pos' local eq = substr("`usecoef'",2,`pos'-2) local colon = cond(`pos',":","") local var = substr("`usecoef'",`pos'+1,.) local coef`ncoefs' `eq'`colon'`var' local coefname`ncoefs' `infix'`eq'`infix'`var' if `pos' di as result "[`eq']`var' " _c else di as result "`var' " _c } } else { * initial run to identify equations without pp problems * stcox and mixed are problems * by default this fits FULL command to ALL data - learnif() simplifies di as text "Learning equation names..." if "`regcmd'"=="stcox" local estimate estimate `noidicmd' qui `anything' if `touse' & `learnif' `wtexp' `postfix' /// `strata' `constant' `eformopts' `options' `estimate' tempname b mat `b' = e(b) local beqs : coleq `b' local beqsuniq : list uniq beqs local bnames : colname `b' local beq1 : word 1 of `beqs' local ncoefs 0 di as text "Using coefficients: " _c * special commands if inlist("`regcmd'","mlogit","mvreg") { local useeqs1 `beqsuniq' if inlist("`regcmd'","mlogit") { // drop equation for base outcome local dropeq : label (`yvar') `e(baseout)' local useeqs1 : list useeqs1 - dropeq } local useeqs2 } else { local useeqs1 `beq1' if "`useeqs'"=="*" local useeqs2 : list beqsuniq - beq1 else local useeqs2 `useeqs' } local neqs : word count `useeqs1' `useeqs2' * coefficients in first type of equation (use all x's) foreach useeq of local useeqs1 { foreach usevar of local usevars { local ++ncoefs local coef`ncoefs' `useeq':`usevar' // for matrix subscripting if `neqs'==1 & mi("`longnames'") { // short name when possible local coefname`ncoefs' `infix'`usevar' // variable name in postfile di as result "`usevar' " _c } else { local coefname`ncoefs' `infix'`useeq'`infix'`usevar' // variable name in postfile di as result "[`useeq']`usevar' " _c } } } * coefficients in second type of equation (use everything in last fit) foreach useeq of local useeqs2 { forvalues j=1/`=colsof(`b')' { local thiseq : word `j' of `beqs' if "`useeq'" == "`thiseq'" { // jth parm is from equation `useeq' local thisname : word `j' of `bnames' if substr("`thisname'",1,2)=="o." continue local ++ncoefs local coef`ncoefs' `thiseq':`thisname' local coefname`ncoefs' `infix'`thiseq'`infix'`thisname' di as result "[`thiseq']`thisname' " _c } } } } di if !mi("`dryrun'") { di as text "Dryrun requested - no analysis done" exit } *********************** SET UP POSTFILE *********************** if "`debug'"=="debug" di as input "Setting up postfile..." * identify coefficients and variances for post file forvalues i = 1 / `ncoefs' { local bvars `bvars' `y'`coefname`i'' local bpostvars `bpostvars' (`coef`i'') forvalues j = `i' / `ncoefs' { local Svars `Svars' `S'`coefname`i''`coefname`j'' } } * ereturned and returned quantities foreach e in `esave' { local evars `evars' _e_`e' } foreach r in `rsave' { local rvars `rvars' _r_`r' } if "`countby'"~="" { // SK code capture confirm string variable `countby' local strcount = _rc==0 qui levelsof `countby', local(countbylevls) foreach levl in `countbylevls' { local addcount ns_`countby'_`levl' nf_`countby'_`levl' pt_`countby'_`levl' local countsby `countsby' `addcount' } capture confirm variable _rec _nrec, exact local uniquecond = cond(_rc == 0, "& _rec == _nrec", "") /* distinct observations in case-cohort data declared by stsetcco */ } /* SK generate time variable for computing person-time below, taking account of delayed entry */ tempvar timevar qui gen double `timevar' = . qui if "`:char _dta[_dta]'" == "st" { replace `timevar' = _t - _t0 local idvar: char _dta[st_id] } if "`append'"=="append" { cap confirm file "`saving'" if _rc confirm file "`saving'.dta" tempfile postfile } else local postfile `saving' postfile `post' `bypost' `bvars' `Svars' `evars' `rvars' `counts' `countsby' /// _ppfixed _ppremains using "`postfile'", `replace' if !mi("`debug'") { foreach thing in bvars Svars evars rvars counts { di as text "`thing':" as result _col(21) "``thing''" } } *********************** MAIN ANALYSIS *********************** if !mi("`debug'") di as input `"Basic command is: `anything' `if' `in' `wtexp' `postfix' `strata' `constant' `eformopts' `options'"' tempvar esample individual gen `esample' = 0 gen `individual' = _n foreach level of local byvarnamelevels { `qui' di _newline if !mi("`debug'") di as input `"level `level'"' * identify this level if `ischarbyvarname' local levelvalue `""`level'""' else local levelvalue `level' summ `individual' if `byvarname' == `levelvalue', meanonly local first = r(min) // an obs in this by-group if mi("`first'") { di as error "Program error in mvmeta_make has led to wrongly seeking `byvarname'==`levelvalue'" exit 498 } local bylevels local bydesc foreach byvar of local byvarlist { local bylevel = `byvar'[`first'] cap confirm string variable `byvar' if !_rc { // string local bylevel `"`"`bylevel'"'"' local bylevel2 `bylevel' } else { // numeric local bylevel2 : label (`byvar') `bylevel' } local bylevels `bylevels' (`bylevel') local bydesc `bydesc' `byvar'==`bylevel2' } * di as input `"-> `bydesc'"' // not silenced by quietly di `"{input}-> `bydesc'"' // silenced by quietly local pptofix 0 local fix fix forvalues pass=0/1 { local ppfound 0 if "`ppfix'"=="all" { * always adjust for perfect prediction: bypass initial model fit local fix avoid local pptofix 1 } * DO THE MAIN REGRESSION local regrc 0 if `pptofix'==0 { // WITHOUT FIXING PP `noidicmd' `hardcap' `qui' `anything' if `byvarname'==`levelvalue' & `touse' `wtexp' `postfix' `strata' `constant' `eformopts' `options' if "`hard'"=="hard" { local regrc = _rc if `regrc' { di as error `"`bydesc': `regcmd' has returned error code `regrc': this may indicate perfect prediction"' local ppfound 1 } else if "`details'"!="nodetails" { cap `regcmd', `eformopts' `options' if !_rc `regcmd', `eformopts' `options' else { cap `regcmd', `eformopts' if !_rc `regcmd', `eformopts' else { cap `regcmd' if !_rc `regcmd' else di "Can't redisplay `regcmd' results" } } } } } else if `pptofix'==1 { // WITH FIXING PP if "`ppcmd'"!="" { di as error `"`bydesc': attempting to `fix' perfect prediction by using `ppcmd'"' `qui' `ppcmd1' `regbody' if `byvarname'==`levelvalue' & `touse' `wtexp', `strata' `constant' `eformopts' `options' `ppcmd2' } else { if "`regcmd'"=="cox" | substr("`regcmd'",1,2)=="st" local st st di as error `"`bydesc': attempting to `fix' perfect prediction by augmenting with `augwt' observations per parameter"' if "`st'"=="st" di as error `" (counting baseline hazard as 1 parameter)"' preserve tempvar wtvar augvar local list = cond("`auglist'"!="noauglist", "list", "") if "`regcmd'"=="clogit" local groupvar groupvar(`e(group)') `qui' _augment if `touse' `wtexp', `list' `st' xvarlist(`xvarlist') /// yvar(`yvar') subgroup(`byvarname'==`levelvalue') timesweight(`augwt') /// wtvar(`wtvar') augvar(`augvar') `strata' `groupvar' `constant' local wtexp2 [iweight=`wtvar'] if "`st'"=="st" { qui stset _t _d if _st `wtexp2', time0(_t0) local wtexp2 } `qui' di as text _newline "Analysis after augmentation:" `noidicmd' `qui' `anything' if (`byvarname'==`levelvalue' & `touse')|`augvar' `wtexp2' `postfix' `strata' `constant' `eformopts' `options' restore } } if `regrc'==0 { mat `ymat'`level'=e(b) mat `Smat'`level'=e(V) if !mi("`debug'") { mat l `ymat'`level', title(e(b) for study `level') mat l `Smat'`level', title(e(V) for study `level') } * CHECK FOR PERFECT PREDICTION: /* *** change 20jan2022: reduce to the required matrix before checking for PP * fails with logit - need to add the equation names to `row', `col' tempname SS mat `SS'=J(`p',`p',.) forvalues r=1/`p' { local row : word `r' of `usevars' forvalues c=1/`p' { local col: word `c' of `usevars' mat `SS'[`r',`c'] = `Smat'`level'["`row'", "`col'"] } } if !mi("`debug'") mat l `SS', title("Variance matrix being checked") cap varcheck `SS', check(pd) */ cap varcheck `Smat'`level', check(pd) if _rc { `qui' di di as error `"`bydesc': perfect prediction or inestimable parameter"' di as error `" (V matrix is not positive definite)"' `qui' mat l `Smat'`level', title(V) local ppfound 1 } else { local missing forvalues j = 1 / `ncoefs' { cap mat `onevalue' = `ymat'`level'[1,"`coef`j''"] if _rc local missing `missing' `y'`coefname`j'' } if !mi(`"`missing'"') { `qui' di di as error `"`bydesc': perfect prediction or inestimable parameter"' di as error `" (parameters `missing' have been dropped)"' local ppfound 1 } } } * OPTIONALLY PREPARE FOR 2ND PASS if "`regcmd'"=="blogit" & `ppfound' & "`ppcmd'"=="" { di as error `"Sorry, augmentation is not available with blogit as it doesn't allow weights"' continue, break } if "`ppfix'"=="all" | `pass'==1 { if `ppfound' di as error `"`bydesc': perfect prediction has not been `fix'ed"' else di as text `"`bydesc': perfect prediction has been `fix'ed"' } if "`ppfix'"=="none" & `ppfound'==1 di as error `"`bydesc': perfect prediction fixing disabled"' if "`ppfix'"=="all" | "`ppfix'"=="none" | `ppfound'==0 continue, break if `pass'==0 local pptofix `ppfound' } if !mi("`debug'") di as input "Estimation completed for study `levelvalue'" * POST THE RESULTS local blist local Slist forvalues i = 1 / `ncoefs' { cap mat `onevalue' = `ymat'`level'[1,"`coef`i''"] if _rc { di as error `"`bydesc': missing `coefname`i''"' local bvalue . } else { local bvalue = `onevalue'[1,1] } local blist `blist' (`bvalue') forvalues j = `i' / `ncoefs' { cap mat `onevalue' = `Smat'`level'["`coef`i''", "`coef`j''"] if !_rc { local Svalue = `onevalue'[1,1] } else { di as error `"`bydesc': missing `S'`coefname`i'`coefname`j''"' local Svalue . } if `i'==`j' & `Svalue' == 0 { di as error `"`bydesc': missing `S'`coefname`i'`coefname`j''"' local Svalue . } local Slist `Slist' (`Svalue') } } local rlist foreach r in `rsave' { local rlist `rlist' (r(`r')) } local elist foreach e in `esave' { if inlist("`e(cmd)'", "logistic", "logit", "clogit") & "`e'"=="N_sub" { qui count if e(sample) /*21dec09*/ local elist `elist' (`=r(N)') /*21dec09*/ } else if inlist("`e(cmd)'", "logistic", "logit", "clogit") & "`e'"=="N_fail" { qui count if `e(depvar)'==1 & e(sample) local elist `elist' (`=r(N)') } else local elist `elist' (`=e(`e')') /* 30may2014: unavailable `e' evaluates as . not empty, hence doesn't crash */ } local clist foreach count in `counts' { if "`count'"==substr("records",1,length("`count'")) { qui count if e(sample) local clist `clist' (`=r(N)') } else if "`count'"==substr("subjects",1,length("`count'")) { tempvar isin gen `isin' = e(sample) if "`regcmd'"=="stcox" { // identify unique obs per id local id : char _dta[st_id] if "`id'"!="" { local sort : sortedby sort `id' `isin' qui by `id': replace `isin'=0 if _n<_N if "`sortedby'"!="" sort `sortedby' } } qui count if `isin' local clist `clist' (`=r(N)') } else if "`count'"==substr("failures",1,length("`count'")) { if "`regcmd'"=="stcox" local event _d else local event `e(depvar)' qui count if e(sample) & `event' local clist `clist' (`=r(N)') } else { di as error `"counts(`count') not allowed"' exit 198 } } if "`countby'"~="" { local countsby local depvar = cond("`e(cmd)'"=="cox", "_d", "`e(depvar)'") if !missing("`idvar'") { tempvar idsubgroup egen `idsubgroup' = tag(`idvar' `countby') if e(sample) `uniquecond' local idsubgcond = "& `idsubgroup' == 1" } qui foreach levl in `countbylevls' { local subgroup = cond(`strcount'==1, `"`countby'=="`levl'""', `"`countby'==`levl'"') count if e(sample) & `subgroup' `uniquecond' `idsubgcond' local ns = r(N) count if e(sample) & `depvar'==1 & `subgroup' `uniquecond' local nf = r(N) local pt = . if !missing("`timevar'") { summ `timevar' if e(sample) & `subgroup' local pt = r(sum) } local addcount (`ns') (`nf') (`pt') local countsby `countsby' `addcount' } capture drop `idsubgroup' } qui replace `esample' = e(sample) if `byvarname'==`levelvalue' & `touse' local postcommand post `post' `bylevels' `blist' `Slist' `elist' `rlist' `clist' `countsby' (`pptofix') (`ppfound') if "`debug'"=="debug" di as input `"Post command: `postcommand'"' `postcommand' `pause' } postclose `post' preserve if !mi("`collapse'") { // combine postfile with summary results cap collapse `collapse' if `touse', by(`byvarlist') if _rc { di as error "mvmeta_make: collapse option failed" } else { tempfile collapsedata sort `byvarlist' qui save `collapsedata' use `postfile', clear sort `byvarlist' qui merge 1:1 `byvarlist' using `collapsedata' cap assert _merge==3 if _rc di as error "Warning: merge was imperfect" drop _merge qui save `postfile', replace } } // Tidy up saved results (and append to existing results if requested) use "`saving'", clear if "`append'"=="append" append using `postfile' label var _ppfixed "Perfect prediction fixed?" label var _ppremains "Perfect prediction remains in this result?" if !mi("`labelsfile'") qui do `labelsfile' foreach byvar of local byvarlist { if !mi("`label`byvar''") label val `byvar' `label`byvar'' } `qui' di if !mi("`savingorig'") save "`saving'", replace if !mi("`clear'") { if !mi("`savingorig'") di as text "mvmeta_make results are also loaded into memory" else di as text "mvmeta_make results are now loaded into memory" restore, not ereturn post // removes ereturned results from last regression } else { restore ereturn post, esample(`esample') } // e-class ereturn local cmdline mvmeta_make `0' ereturn local cmd mvmeta_make end ****************************************************************************************** program define _augment * varlist changed to new options 22jun2015: xvarlist, yvar, subgroup version 11 syntax [if] [in] [fweight aweight pweight iweight], /// xvarlist(varlist) subgroup(string) /// [yvar(varlist) TOTALweight(real 0) TIMESweight(real 0) noPREServe list /// wtvar(string) augvar(string) noequal st strata(varlist) /// noconstant groupvar(string)] if mi("`yvar'`st'") { di as error "mvmeta_make: _augment call must have yvar() or st option" exit 498 } if "`wtvar'"=="" local wtvar _weight if "`augvar'"=="" local augvar _augment confirm new var `wtvar' `augvar' tempvar augment wt if "`weight'"!="" gen `wt' `exp' else gen `wt' = 1 marksample touseall tempvar touse gen `touse' = `touseall' & (`subgroup') // Sort out some locals if "`st'"=="" { qui levelsof `yvar' if `touseall', local(ylevels) // NB and not if `subgroup' qui tab `yvar' if `touseall' [iw=`wt'] local nylevels = r(r) local totw = r(N) } else { st_is 2 analysis qui replace `touse'=0 if _st==0 local yvar _t local nylevels 2 local st_wv : char _dta[st_wv] if "`st_wv'"!="" replace `wt'=`wt'*`st_wv' summ `wt' if _st & `touse', meanonly local totw = r(sum) local stvars _t _d summ _t if _st & `touse', meanonly local tmin = r(min)/2 local tmax = r(max)+r(min)/2 local ylevels `tmin' `tmax' } local Nold = _N // Count x's local nx : word count `xvarlist' // Set total weight of added observations if `totalweight'>0 & `timesweight'>0 { di as error `"_augment: Can't specify totalweight() and timesweight()"' exit 498 } if `totalweight' == 0 local totalweight = `nx' + ("`constant'"!="noconstant") if `timesweight'>0 local totalweight = `totalweight' * `timesweight' // Number of added observations local Nadd = `nx'*2*`nylevels' local Nnew = `Nold'+`Nadd' di as text "Adding " as result `Nadd' as text " pseudo-observations "_c if "`strata'"!="" di as text "per stratum " _c di as text "with total weight " as result `totalweight' qui { set obs `Nnew' gen `augment' = _n>`Nold' local thisobs `Nold' // Handle groupvar - intended for clogit if "`groupvar'"!="" { qui summ `groupvar' if `touse' local augxnum = r(max) // Counts x-groups } else local augxnum 0 foreach x of var `xvarlist' { local augxnum=`augxnum'+2 sum `x' [iw=`wt'] if `touse' & !`augment' if r(sd)==0 { * No variation in this cohort: use SD across all cohorts sum `x' [iw=`wt'] if !`augment' } local mean = r(mean) local sd = r(sd)*sqrt((r(sum_w)-1)/r(sum_w)) replace `x' = `mean' if `augment' local auginum 0 // Counts obs within x-groups foreach yval in `ylevels' { if "`equal'"=="noequal" { summ `wt' if `yvar'==`yval' & `touse' & !`augment', meanonly local py = r(sum)/`totw' } else local py = 1/`nylevels' foreach xnewstd in -1 1 { local ++thisobs local ++auginum replace `x' = `mean'+(`xnewstd')*`sd' in `thisobs' replace `yvar' = `yval' in `thisobs' replace `wt'= `totalweight' * `py' / (2*`nx') in `thisobs' if "`groupvar'"!="" replace `groupvar' = `augxnum' + (`auginum'==2|`auginum'==3) in `thisobs' } } } if "`st'"=="st" { replace _d=1 if `yvar'==`tmin' & `augment' replace _d=0 if `yvar'==`tmax' & `augment' replace _t0=0 if `augment' replace _st=1 if `augment' } } // Handle strata if "`strata'"!="" { qui { foreach stratvar of varlist `strata' { cap confirm string variable `stratvar' local isstring = _rc==0 levelsof `stratvar' if `touse' & !`augment', local(stratlevels) local nlevels : word count `stratlevels' forvalues i=2/`nlevels' { local Noldplus1=_N+1 expand 2 if `augment' & missing(`stratvar') local thislevel : word `i' of `stratlevels' if `isstring' replace `stratvar'="`thislevel'" in `Noldplus1'/l else replace `stratvar'=`thislevel' in `Noldplus1'/l } local thislevel : word 1 of `stratlevels' if `isstring' replace `stratvar'="`thislevel'" if `augment' & missing(`stratvar') else replace `stratvar'=`thislevel' if `augment' & missing(`stratvar') replace `wt'=`wt'/`nlevels' if `augment' } } } rename `wt' `wtvar' rename `augment' `augvar' if "`list'"=="list" { di as txt _newline "Listing of added records:" char `wtvar'[varname] weight label var `wtvar' "Weight" char `augvar'[varname] augment label var `augvar' "Augmented?" list `stvars' `strata' `groupvar' `xvarlist' `yvar' `wtvar' if `augvar', sep(4) subvarname } end prog def varcheck syntax anything, [check(string)] if colsof(`anything') != rowsof(`anything') { di as error `"varcheck: `anything' is not a square matrix"' exit 498 } * Remove any coeffs whose name starts "o." matdropo `anything' `anything' local dim = rowsof(`anything') * Check for missing values, adapted to cope with the following special case * In Stata 9, I have come across a matrix M of values 1.#QNAN that returns matmissing(M)=0 but matmissing(M+M)=1 tempname anything1 mat `anything1'=`anything'+`anything' if matmissing(`anything1') { di as error `"varcheck: `anything' has missing values"' exit 498 } if `dim' == 1 { * 1x1 matrix local mineigen = `anything'[1,1] } else { * larger matrix forvalues r=1/`dim' { if `anything'[`r',`r']<=0 { di as error: `"varcheck: `anything' has variance<=0"' exit 498 } forvalues s=1/`dim' { if (`anything'[`r',`s'])^2>`anything'[`r',`r']*`anything'[`s',`s'] { di as error: `"varcheck: `anything' has correlation outside [-1,1]"' exit 498 } } } if "`check'"!="" { tempname evecs evals mat symeigen `evecs' `evals' = `anything' local mineigen = `evals'[1,`dim'] } } if "`check'"=="psd" { if `mineigen'<0 { di as error `"varcheck: `anything' is not positive semi-definite"' mat l `evals' exit 506 } } else if "`check'"=="pd" { if `mineigen'<=0 { di as error `"varcheck: `anything' is not positive definite"' mat l `evals' exit 506 } } end prog def matpart args matin matselect matout * selects submatrix of matin indicated by matselect and stores it in matout mata { sel=st_matrix("`matselect'") matin=st_matrix("`matin'") matout=select(select(matin,sel),sel') st_matrix("`matout'",matout) } end prog def matdropo args matin matout * selects submatrix of matin whose rownames don't start "o." and stores it in matout local rowfullnames : rowfullnames `matin' tempname select foreach rowfullname of local rowfullnames { if strpos("`rowfullname'",":") { gettoken roweq rowname : rowfullname, parse(":") local rowname : subinstr local rowname ":" "" } local ok = strpos("`rowname'","o.")!=1 mat `select' = nullmat(`select'),(`ok') if `ok' local newrowfullnames `newrowfullnames' `rowfullname' } matpart `matin' `select' `matout' mat rownames `matout' = `newrowfullnames' mat colnames `matout' = `newrowfullnames' end prog def dicmd noi di as input `"`0'"' `0' end