*! version 4.4.1 17jan2023 MJC
/*
History - multistate/predictms
17jan2023: version 4.4.1 - chintpoints() option added to allow the user to specify the number of
quadrature points used in any numerical integrations
22dec2021: version 4.4.0 - Cox model is now supported with predictms (ci's not supported yet)
- added moremata check
02sep2021: version 4.3.1 - bug fix, graphms ignored extra options; now fixed
- bug fix, when rmst was called with simulate it would error out; now fixed
14mar2021: version 4.3.0 - accuracy and efficiency of numerical integration based predictions from an illness-death and extended
illness death model are improved
- support for family(hazard) now documented
- bug fix for numerical integration based predictions of illness-death and extended illness-death, failing
to pass the entry time appropriately would return a vector of 0s instead of appropriate probabilities;
now fixed
- bug fix; if reset or tsreset() was used with an illness-death model, then numerical integration was used
to calculate predictions instead of simulation, giving (obvious) incorrect results - this has been fixed
- bug fix when using a relative survival model within predictms, would attempt to call bhazard() when it
should only simulate from the relative survival function; now fixed
- warning added when passing a relative survival model to -predictms-, to ensure the expected rate has been
modelled.
09feb2021: version 4.2.1 - bug fix; check on transmatrix() skipped final row, so could erroneously exit telling the user the number
of transitions and models didn't match - now fixed
26jan2012: version 4.2.0 - interactive & freqat() options added to msboxes
16dec2020: version 4.1.0 - family(pwexponential) synced & doc'd
28oct2020: version 4.0.0 - delta method is now the default for ci calculation for methods that use numerical integration or the hybrid AJ estimator for predictions
-> this is substantially faster than the bootstrap
-> includes contrasts
-> bootstrap is still the only method available for large-sample simulation based predictions
- bootstrap option added to override the default delta method
- latent option now documented (added in 3.0.0 to invoke latent times simulation, overiding the default total hazard + multinomial draws method)
- log normal, log logistic and general log cumulative hazard scale models added
- memory leak fixed, resulting in speed improvements
11sep2020: version 3.0.3 - msaj option enter() renamed to ltruncated() for consistency with survsim & predictms
08sep2020: version 3.0.2 - rmst missed from standardise calculation; now fixed
- ci + visit or rmst error'ed; now fixed
- storage improved, so bootstrap cis slightly faster
01sep2020: version 3.0.1 - left over print statements when simulating; now removed
24aug2020: version 3.0.0 - predictms rewrite:
- Stata version 15.1 now required
- all models must be fitted using stmerlin or merlin
- transprob renamed to probability
- survival renamed to singleevent for transmat() shorthand
- ltruncated() introduced, which defines the starting time for all predictions
- rmst added to predict restricted mean survival time, i.e. total time spent in non-absorbing states
- hazard and survival options added to predict transition specific functions
- standard survival, competing risks, illness-death and extended illness-death predictions are now
calculated using numerical integration as the default method, rather than large-sample simulation
- the default method of simulation to obtain predictions uses Beyersmann et al., with a new -latent- option
to use the previous default of latent times
- each simulated dataset can now be saved
- ci now displays a bootstrap output so it can be monitored
- dist(exp) failed; now fixed
- ci's now calculated on log/atanh/logit scale and back transformed to ensure appropriate bounds
- graphms added, with help file
- help files substantially improved, new main multistate help file added
*/
/*
History - galahad
MJC 21mar2020: version 1.3.0 - some behind the scenes syncing with survsim
MJC 17mar2020: version 1.2.1 - bug fix; singleevent option failed -> now fixed
- help file example fixes for new transprob option
- bug fix; missing error check, Cox model is not currently supported
MJC 07mar2020: version 1.2.0 - hazard option added to predictms to predict each transition-specific hazard function
- survival changed to singleevent
- survival now predicts transition-specific survival functions
- transprob option now required to request transition probabilities
- bug fix: visit with cis incorrectly returned los predictions instead -> now fixed.
- bug fix: contrasts with visit, new variable names were incorrect -> now fixed
- bug fix: userlink, transformation not applied to point estimate when calculating CI -> now fixed
- bug fix: contrasts and ci, point estimates weren't passed -> now fixed
MJC 30jan2020: version 1.1.0 - bug fix displaying dots for multiple ats()
- aj wasn't synced properly - now working
- synced multiple timescales with all survival models
- added timevar to ret list for use with ghgraph
- ghgraph added, including help file
- bug fix; family(logchazard) wasn't synced - now fixed
MJC 28jan2020: version 1.0.0
*/
/*
History - multistate/predictms
MJC 13aug2019: version 2.3.2 - bug fix: missed an edit in 2.3.1 fix preventing aj to work in predictms (merlin function names); now fixed
- error check added for streg, d(ggamma) combined with aj - not supported
MJC 25jul2019: version 2.3.1 - bug fix: edits in merlin function names caused predictms to fail; now fixed
- error check added for bhtime with strcs models - not currently supported
MJC 03jun2019: version 2.3.0 - bug fix: stpm2 and strcs models failed, introduced in 2.2.0; now fixed
- streg, dist(lognormal) now allowed with clock-forward models
- bug fix: streg, dist(lognormal) with reset incorrectly assumed loglogistic; now fixed
MJC 28may2019: version 2.2.0 - bug fix: with novcv() when ci not specified introduced in 2.1.0; now fixed
- bug in error check for log normal requiring reset; now fixed
- Generalised gamma streg model now supported when used in the models() syntax
- bug fix in msaj when number of transitions was not equal to number dimension of transmatrix()
MJC 15apr2019: version 2.1.0 - novcv() option added to use mean vector instead of draws in CI calculations
MJC 19mar2019: version 2.0.1 - help file improved for msaj
MJC 23feb2018: version 2.0.0 - bug fix: los was incorrectly scaled to 1, introduced in 1.2.0. Now fixed.
- normal approximation for CIs is now the default. percentile option replaces normal.
- bug fix: level() was ignored in CI calculations, always assumed 95%. Now fixed.
- gen() removed
- bug fix: in some cases mm_root_vec produced an error when solving the root. Now fixed.
- all predictions now available in one call
- infinite at#()s added (limited to 50 for error check reasons)
- at() changed to at1()
- difference option added
- added atref(#) - default 1
- standardise option added to calculate standardised (population-averaged) predictions
- bug fix: with Stata 15 and streg (apart from dist(exp)), with at2(), the ancilary parameter constant was ignored. Now fixed.
- _time = 0 now allowed with normal ci's
- userfunction() added for user-defined prediction function, subroutines for probs and los
- userlink() added for normal cis of userfunction() - default identity, can also be log or logit
- cr option added to avoid specifying a transition matrix. For use with models() only.
- bug fix: reset now synced with root-finding updates for stpm2 and strcs models
- outsample added
- stms added
- enter() removed, now min(timevar or _t)
- aj added for Aalen-Johansen estimator with Markov models
- reversible transitions now allowed
MJC 17nov2017: version 1.2.1 - bug fix with simulation from strcs models with delayed entry
MJC 13nov2017: version 1.2.0 - tscale2() and time2() added for multiple timescales
- probabilities/los scaled to 1/t
release moved to website
MJC 21nov2016: verison 1.1.1 - bug fix
MJC 16nov2016: version 1.1.0 - some re-writes of source code for massive speed improvements when root-finding required
- trans#() syntax removed, now parsing is done automatically and using covariates expanded using msset
- model#() changed to models()
- improvements to help files
- covariates() added to msset to create transition-specific dummies, which must be used in model fitting (if required)
- default n() now 100,000, or 10,000 with ci
- survival added to calculate all predictions on a standard single event model
MJC 09aug2016: version 1.0.1 - gen() option to allow a stub for created variables, default is pred
MJC 12jul2016: version 1.0.0
Development
stms not allowed yet
MJC 19may2015: version 1.1.0 - error check on transmatrix() added
- separate models allowed through model#(name, [at()])
MJC 14may2015: verison 1.0.8 - added error check not allowing AFT weibull or exp
MJC 11may2015: verison 1.0.7 - timevar() added
- basic graph option added which creates stacked plots
MJC 11may2015: verison 1.0.6 - fixed bug with weibull and forward approach
- Error checks improved
MJC 09may2015: version 1.0.5 - fixed bug in from() which occurred when anything but from(1) was used
- fixed bug in forward calculations when enter>0
- fixed bug that only calculated predictions to states you could go to from first state
MJC 09may2015: version 1.0.4 - clock-forward approach now the default (simulations incorporate delayed entry), reset option added to use clock-reset
- only reset approach allowed with streg, dist(lnormal)
MJC 06may2015: version 1.0.3 - now synced with streg, dist(exp|weib|gompertz|llogistic|lnormal)
- normal approximation added for CI calculations
MJC 15apr2015: version 1.0.2 - when no ci's calculated it was using first draw from MVN, this has been fixed to be e(b)
MJC 01apr2015: version 1.0.1 - stpm2 simulation improved by creating and passing struct
- odds and normal scales added for stpm2 models
MJC 31mar2015: version 1.0.0
*/
program define predictms, sortpreserve properties(st) rclass
version 15.1
syntax , ///
///
[ ///
TRANSMatrix(string) /// -transition matrix-
SINGLEEVENT /// -a single event survival model was fitted-
CR /// -competing risks model was fitted-
///
PROBability /// -transition probabilities-
LOS /// -calculate length of stay in each state-
RMST /// -restricted mean survival time-
VISIT /// -prob. of ever visiting each state-
HAZard /// -transition-specific hazard functions-
SURVival /// -transition-specific survival functions-
USERFunction(string) /// -name of Mata function-
USERLink(string) /// -link function to apply to Mata function-
///
MODels(string) /// -est store objects-
///
FROM(numlist >0 asc int) /// -starting state for predictions-
TIMEvar(varname numeric) /// -prediction times-
LTruncated(numlist >=0 max=1) /// -time at which observations become at risk-
OBS(string) /// -Number of time points-
MINT(string) /// -minimum time to calculate predictions at-
MAXT(string) /// -maximum time to calculate predictions at-
///
SEED(string) /// -pass to set seed-
N(string) /// -sample size-
SIMulate /// -use simulation-
LATENT /// -use latent times simulation-
AJ /// -use AJ estimator-
///
M(numlist >=20 int max=1) /// -number of parameter samples from MVN-
CI /// -calculate confidence intevals-
BOOTstrap /// -use parametric bootstrap for cis-
PERCentile /// -use percentile based cis-
NOVCV(string) /// -transitions to skip VCV when calculating CIs-
Level(cilevel) /// -level for CIs-
///
DIFFerence ///
RATIO ///
ATREFerence(string) /// -reference at for contrasts-
///
STANDardise /// -standardise predictions-
STANDIF(string) ///
///
RESET /// -use clock reset approach in simulations-
TSCALE2(string) /// -transition models on a second timescale-
TIME2(string) /// -time to add to main timescale-
TSRESET(string) /// -transition-specific resets-
///
CHINTPoints(numlist max=1 >0 int) ///
///
OUTsample ///
SAVE(string) ///
///
///
SURVSIM(string) /// -not documented-
SURVSIMTOUSE(string) /// -not documented-
///
///
DEVCODE1(string) /// -not documented-
DEVCODE2(string) /// -not documented-
DEVCODE3(string) /// -not documented-
DEVCODE9(string) /// -not doc-
///
INTERACTIVE JSONPATH(string) ///
///
* /// -infinite ats-
]
//====================================================================//
// Preliminaries
capture which merlin
if _rc>0 {
display in yellow "You need to install the command merlin. This can be installed using,"
display in yellow ". {stata ssc install merlin}"
exit 198
}
capture which lmoremata.mlib
if _rc>0 {
display in yellow "You need to install the moremata package. This can be installed using,"
display in yellow ". {stata ssc install moremata}"
exit 198
}
local ats `options'
if "`probability'`los'`visit'`hazard'`survival'`userfunction'`rmst'"=="" {
di as error "At least one prediction type must be specified"
exit 198
}
if "`aj'"!="" {
if "`probability'"=="" {
local probability probability
}
if "`tscale2'"!="" {
di as error "aj not allowed with tscale2()"
exit 198
}
if "`reset'"!="" {
di as error "aj not allowed with reset"
exit 198
}
if "`n'"!="" {
di as error "n() not needed with aj"
exit 198
}
if "`visit'"!="" {
di as error "visit not allowed with aj"
exit 198
}
}
local K = 1
if "`standardise'"!="" {
if "`hazard'"!="" {
di as error "hazard and standardise not supported"
exit 198
}
if "`outsample'"!="" {
di as error "Not valid"
exit 1986
}
if "`userfunction'"!="" {
di as error "userfunction() not allowed with standardise"
exit 198
}
tempvar stdtouse
if "`standif'"!="" {
qui gen byte `stdtouse' = (`standif')
}
else {
qui gen byte `stdtouse' = 1
}
qui count if `stdtouse'
local K = r(N)
}
if "`outsample'"!="" {
local out = 1
}
else local out = 0
if `out' & "`timevar'"=="" {
di as error "timevar() must be specified with outsample"
exit 198
}
if "`visit'"!="" {
local simulate simulate
}
if ("`difference'"!="" | "`ratio'"!="") & "`ats'"=="" {
di as error "at least at2() needed"
exit 198
}
if "`seed'"!="" {
set seed `seed'
}
local survsimcall = "`survsim'"!=""
if "`singleevent'"!="" {
if "`transmatrix'"!="" {
di in yellow "transmatrix(`transmatrix') ignored"
}
tempname transmatrix
mat `transmatrix' = (.,1\.,.)
local latent latent //faster with single event
}
if "`cr'"!="" {
if "`models'"=="" {
di as error "models() required with cr"
exit 198
}
if "`transmatrix'"!="" {
di in yellow "transmatrix(`transmatrix') ignored"
}
local Ntrans : word count `models'
tempname transmatrix
mat `transmatrix' = J(`Ntrans'+1,`Ntrans'+1,.)
forvalues i=1/`Ntrans' {
mat `transmatrix'[1,`i'+1] = `i'
}
}
cap confirm matrix `transmatrix'
if _rc>0 & "`singleevent'"=="" & "`cr'"=="" {
di as error "transmatrix(`transmatrix') not found"
exit 198
}
mata: check_transmatrix()
if ("`tscale2'"!="" & "`time2'"=="" ) | ("`tscale2'"=="" & "`time2'"!="" ) {
di as error "tscale2() and time2() must both be specified"
exit 198
}
if "`tscale2'"!="" & "`reset'"!="" {
di as error "reset not allowed with tscale2()"
exit 198
}
if "`ci'"=="" & "`m'"!="" {
di as error "Cannot specify m() without ci"
exit 198
}
if "`save'"!="" {
if "`aj'"!="" {
di as error "save() not allowed with aj"
exit 198
}
local 0 `save'
syntax anything(name=savestub) , [REPlace]
local savereplace `replace'
}
if "`novcv'"!="" & "`models'"=="" {
di as error "novcv() requires models()"
exit 198
}
if "`models'"!="" {
//pase models()
local Nmodels : word count `models'
if `Ntrans'!=`Nmodels' {
di as error "Number of estimates objects in model() must be equal to number of transitions"
exit 198
}
forvalues i=1/`Ntrans' {
local model`i' : word `i' of `models'
}
forvalues i=1/`Ntrans' {
local 0 `model`i''
syntax name(id="model estimates required")
local modelests`i' `namelist'
}
}
if "`chintpoints'"=="" {
local chintpoints = 30
}
// Checks for interactive options
if "`jsonpath'" != "" & "`interactive'" == "" {
di as error "You have used the jsonpath option without using the interactive option."
exit 198
}
if "`interactive'" != "" {
if "`jsonpath'" != "" {
mata st_local("direxists",strofreal(direxists("`jsonpath'")))
if !`direxists' {
di as error "Folder `jsonpath' does not exist."
exit 198
}
mata st_local("jsonfile",pathjoin("`jsonpath'","msboxes_predictions.json"))
}
else {
local jsonfile msboxes_predictions.json
}
capture confirm file "`jsonfile'"
if !_rc {
capture erase "`jsonfile'"
if _rc {
display as error "`jsonfile' cannot be deleted'"
}
}
}
//====================================================================//
//core method choices
if "`visit'"!="" {
local simulate simulate
local bootstrap bootstrap
}
//====================================================================//
//parse ats
//check at#()
local atind = 1
while "`ats'"!="" {
if `atind'>50 {
di as error "at#() limit reached, or unrecognised option"
exit 198
}
local 0 , `ats'
syntax , [at`atind'(string) *]
local ats `options'
if !`out' {
local varcount : word count `at`atind''
local count = `varcount'/2 //!!add error check = integer
tokenize ``at`atind'''
while "`1'"!="" {
unab 1: `1'
cap confirm var `1'
if _rc {
di in red "invalid at`atind'(... `1' `2' ...)"
exit 198
}
forvalues i=1/`Ntrans' {
if "`1'"=="_trans`i'" {
di as error "Cannot specify _trans# variables in at`atind'()"
exit 198
}
}
cap confirm num `2'
if _rc {
di in red "invalid at`atind'(... `1' `2' ...)"
exit 198
}
mac shift 2
}
}
local atind = `atind'+1
}
if `atind'>1 {
local Nats = `atind'-1
}
else {
local Nats = 1
}
if `Nats'==1 & ("`difference'"!="" | "`ratio'"!="") {
di as error "difference/ratio require at least 2 at#() statements"
exit 198
}
if "`atreference'"=="" {
local atref = 1
}
else {
cap confirm integer number `atreference'
if _rc {
di as error "atreference() must be an integer"
exit 198
}
if `atreference'<1 {
di as error "atreference() must be >=1"
exit 198
}
local atref = `atreference'
}
if `atref'>`Nats' & `Nats'>1 {
di as error "atreference(#) must be an at#()"
exit 198
}
//==//
if "`from'"=="" local from 1
// prediction time variable
if "`timevar'"!="" & ("`mint'"!="" | "`maxt'"!="" | "`obs'"!="") {
di as error "timevar() cannot be specified with mint()/maxt()/obs()"
exit 198
}
if "`ltruncated'"=="" {
local ltruncated = 0
}
if "`timevar'"=="" {
if "`maxt'"=="" {
qui su _t, meanonly
local maxt = `r(max)'
}
if "`mint'"=="" {
local mint = 0
}
if "`obs'"=="" & "`aj'"=="" local obs = 20
else if "`obs'"=="" & "`aj'"!="" local obs = 500
local timevar _time
cap drop _time
cap range2 _time `mint' `maxt' `obs'
label var _time "Follow-up time"
//touse variable for predictions etc.
tempvar touse
qui gen byte `touse' = _n<= `obs'
}
else {
//touse variable for predictions etc.
tempvar touse
qui gen byte `touse' = `timevar'!=.
qui count if `touse'==1
local obs = `r(N)'
qui su `timevar', meanonly
if `r(max)'<`ltruncated' {
di as error "max(`timevar') must be > ltruncated()"
exit 198
}
if `r(min)'<`ltruncated' {
di as error "min(`timevar') must be >= ltruncated()"
exit 198
}
}
if "`aj'"!="" {
qui su `timevar' if `touse'
if `r(min)'!=`ltruncated' {
di as error "mininum of timevar() must be equal to ltruncated() when using aj"
exit 198
}
}
//get core stuff
//possible transitions from each state
local Nstates = colsof(`transmatrix')
if `Nstates'<2 {
di as error "Must be at least 2 possible states, including starting state"
exit 198
}
forvalues i=1/`Nstates' {
forvalues j=1/`Nstates' {
if (`transmatrix'[`i',`j']!=.) {
local row`i'trans `row`i'trans' `=`transmatrix'[`i',`j']'
local row`i'next `row`i'next' `j'
}
}
}
//check somewhere to go
foreach frm in `from' {
if "`row`frm'next'"=="" {
di as error "No possible next states from(`frm')"
exit 198
}
}
//====================================================================//
//CORE
// models framework
if "`models'"!="" {
forvalues i=1/`Ntrans' {
//error checks
cap estimates restore `modelests`i''
if _rc>0 {
di as error "model estimates `modelests`i'' not found"
exit 198
}
predictms_modelcheck
local familys `familys' `r(family)'
//get estimates and variances
tempname emat`i'
mat `emat`i'' = e(b)
if "`ci'"!="" {
tempname evmat`i'
mat `evmat`i'' = e(V)
}
local Nparams`i' = colsof(`emat`i'')
}
//to sync with out of sample predictions - all vars must be present
if `out' {
forvalues i=1/`Ntrans' {
cap estimates restore `modelests`i''
foreach corevar in `e(allvars)' {
cap gen `corevar' = 1 in 1
if !_rc {
local todrop `todrop' `corevar'
}
}
}
}
forval a=1/`Nats' {
forvalues i=1/`Ntrans' {
//handle outsample and extract standardise variables
predictms_model, trans(`i') ///
nparams(`Nparams`i'') ///
ntrans(`Ntrans') ///
at(`at`a'') ///
aind(`a') ///
`standardise'
if "`standardise'"!="" {
local at`a'stdvars`i' `r(stdvars)'
}
}
}
}
//stacked single model framework
else {
predictms_modelcheck
local familys `r(family)'
//get estimates and variances
tempname emat
mat `emat' = e(b)
if "`ci'"!="" {
tempname evmat
mat `evmat' = e(V)
}
local Nparams = colsof(`emat')
//to sync with out of sample predictions - all vars must be present
if `out' {
foreach corevar in `e(allvars)' {
cap gen `corevar' = 1 in 1
if !_rc {
local todrop `todrop' `corevar'
}
}
}
forval a=1/`Nats' {
forvalues i=1/`Ntrans' {
//handle outsample and extract standardise variables
predictms_model, trans(`i') ///
nparams(`Nparams`i'') ///
ntrans(`Ntrans') ///
at(`at`a'') ///
aind(`a') ///
`standardise'
local toupdate `toupdate' `r(toupdate)'
if "`standardise'"!="" {
local at`a'stdvars`i' `r(stdvars)'
}
}
}
}
if "`standardise'"!="" {
local stdcheck = 0
forval a=1/`Nats' {
forvalues i=1/`Ntrans' {
di as text "Transition `i', at`a'(): Standardising over -> `at`a'stdvars`i''"
local stdcheck = `stdcheck' + ("`at`a'stdvars`i''"=="")
}
}
if `stdcheck' {
di as error "No variables for standardising"
exit 1986
}
}
mata: predictms()
if "`interactive'" != "" {
mata: predictms_writejson2()
}
//=====================================================================================================================//
//finish
//tidy up
cap drop `todrop'
if `survsimcall' {
cap drop _time
exit
}
// return list
return matrix transmatrix = `transmatrix', copy
return local Nstates = `Nstates'
return local from `from'
return local timevar `timevar'
//Done
end
// double added to range
program define range2
version 3.1
if "`3'"=="" | "`5'"!="" { error 198 }
confirm new var `1'
if _N==0 {
if "`4'"=="" { error 198 }
set obs `4'
local o "`4'"
}
else {
if "`4'"!="" {
local o "`4'"
if `4' > _N { set obs `4' }
}
else local o=_N
}
gen double `1'=(_n-1)/(`o'-1)*((`3')-(`2'))+(`2') in 1/`o'
end
program predictms_modelcheck, rclass
if "`e(cmd)'"!="merlin" {
di as error "Only merlin models are supported"
exit 198
}
else {
if `e(Nlevels)'>1 {
di as error "Multilevel models are not supported"
exit 198
}
if "`e(failure1)'"=="" {
di as error "The merlin model must be a survival/time-to-event model"
exit 198
}
if "`e(family2)'"!="" {
di as error "Only univariate merlin models are currently supported"
exit 198
}
if "`e(bhazard1)'"!="" {
di as text "warning -> relative survival model detected; only valid when modelling the expected rate"
di as text " -> see " in smcl "{helpb exptorcs:help exptorcs}" as text " for more information"
}
}
return local family "`e(family1)'"
end
mata
void check_transmatrix()
{
tmat = st_matrix(st_local("transmatrix"))
tmat_ind = tmat:!=. //indicator matrix
//Error checks
if (max(diagonal(tmat_ind))>0) {
errprintf("All elements on the diagonal of transmatrix() must be coded missing = .\n")
exit(198)
}
row = 1
rtmat = rows(tmat)
trans = 1
while (row<=rtmat) {
for (i=1;i<=rtmat;i++) {
if (sum(tmat:==tmat[row,i])>1 & tmat[row,i]!=.) {
errprintf("Elements of transmatrix() are not unique\n")
exit(198)
}
if (tmat[row,i]!=. & tmat[row,i]!=trans){
errprintf("Elements of transmatrix() must be sequentially numbered from 1,...,K, where K = number of transitions\n")
exit(198)
}
if (tmat[row,i]!=.) trans++
}
row++
}
st_local("Ntrans",strofreal(trans-1))
}
end
program predictms_model, rclass
syntax [, TRANS(string) NPARAMS(string) Ntrans(real 1) AT(string) STANDardise aind(string)]
foreach corevar in `e(allvars)' {
if !strpos("`e(response1)'","`corevar'") {
local inat = 0
predictms_atparse, corevar(`corevar') ///
at(`at') ///
ntrans(`ntrans') ///
i(`trans')
local inat = r(inat)
local toupdate `toupdate' `r(toupdate)'
if !`inat' & "`standardise'"!="" {
predictms_stdparse `corevar' `trans' `ntrans'
if r(include) {
local stdvars `stdvars' `r(stdvar)'
local toupdate `toupdate' `r(toupdate)'
}
}
}
}
return local stdvars `stdvars'
return local toupdate `toupdate'
end
program predictms_atparse, rclass
syntax, [ corevar(string) ///
i(real 0) ///
Ntrans(real 1) ///
at(string)]
local inat = 0
tokenize `at'
while "`1'"!="" {
unab 1: `1'
if "`corevar'"=="`1'_trans`i'" {
local inat = 1
local toupdate `toupdate' `1'
}
else if "`corevar'"=="`1'" {
local inat = 1
}
mac shift 2
}
//this makes sure you can't use any of the _trans# vars in at1(), at2() or standardising
forvalues j=1/`ntrans' {
if "`corevar'"=="_trans`j'" {
local inat = 1
}
}
return scalar inat = `inat'
return local toupdate `toupdate'
end
program predictms_stdparse, rclass
args corevar i Ntrans
local include = 0
//see if _trans`i' is in it
local has_trans = strpos("`corevar'","_trans`i'")
if `has_trans' {
local strlength = strlen("_trans`i'")
local strlength2 = strlen("`corevar'")
local corevarstub = substr("`corevar'",1,`=`strlength2'-`strlength'')
//need to check that each var matches with *_trans`i', others not included in trans specific design matrix
local extracttrans = substr("`corevar'",`=`strlength2'-`strlength'+1',.)
if "`extracttrans'"=="_trans`i'" {
local include = 1
local toupdate `corevarstub'
}
}
else {
//could still be ?_transj
local hastransj = 0
forvalues j=1/`Ntrans' {
if `i'!=`j' {
if strpos("`corevar'","_trans`j'") {
local hastransj = 1
}
}
}
if !`hastransj' {
local include = 1
local corevarstub `corevar'
}
}
return local stdvar `corevarstub'
return local toupdate `toupdate'
return scalar include = `include'
end