/*********************************************************************
*! 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