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