*! version 4.0.9 15mar2022 MJC

/*
History
15mar2022 version 4.0.9 - bug fix, when survsim_msm called with all exponential, NI fails - now fixed by doing integration
                          analytically for that special case of all CR hazards constant
12jan2022 version 4.0.8 - build error; now fixed
02dec2021 version 4.0.7 - bug fix with msm simulation, would occur if directory of installed files had a space in it
                        - missing .mata files in pkg file
13oct2021 version 4.0.6 - bug fix with varnames in chazard() or logchazard() not parsed correctly with survsim user; now fixed
                        - error check for negative hazard() added
                        - Now tidies up after itself if it errors out
02sep2021 version 4.0.5 - bug fix; mixture + tde() would error out, now fixed
24aug2020 version 4.0.4 - re-syncing with predictms for merlin model simulation
- bug fix in missing value warning text; now fixed
25may2020 version 4.0.3 - efficiency of _msm improved
                        - bug fix: Stata variables in user-defined hazard functions not parsed correctly; now fixed
30mar2020 version 4.0.2 - bug fix; missing value issue when t = 0 in chq function, now fixed
                        - bug fix; indexing issue occurred in varied circumstances, now fixed
29mar2020 version 4.0.1 - bug fix; error in distribution(gompertz) when specified in a hazard#() -> now fixed
                        - help file edits
26mar2020 version 4.0.0 - computation time reduced by about 70% for multi-state model simulation
                        - reset added to hazard#() for semi-Markov models
                        - syntax changed for user-defined functions, #t now {t} for more robust parsing
                        - support for general multiple timescales, including:
                                - {t0} can be used in user() and tdefunction() for time of entry into initial state for associated transition
                        - cr subclass removed, now part of general msm syntax
                        - bug fix; variables used directly in user() with a multi-state model were not passed to Mata -> now fixed
                        - bug fix; if the same variable was used in multiple hazard#()s, parsing would only pick up the first -> now fixed
24mar2020 version 3.1.0 - general Markov multi-state simulation now added
                        - transmatrix() option added
                        - tolerance changed to 0 from 1e-08 in mm_root()
21mar2020 version 3.0.0 - complete re-write of core code
                        - can now simulate from a fitted merlin survival model
                        - can now simulate from cause-specific hazards competing risks models, with either parametric or user-defined hazard functions, with a new syntax
                        - distribution() now required for parameteric models
                        - ltruncated() added, can be number or varname, synced with all except merlin model simulation
                        - maxtime() can now be a number or a varname
                        - help files re-written, hugely improved
                        - added error report on missing values when general algorithm used
                        - parsing made more robust
18dec2013 version 2.0.3 - tde() with logcumh() or cumh() caused an error. Now fixed.
10oct2013 version 2.0.2 - bug fix -> exact option added to confirm vars
08jul2013 version 2.0.1 - minor bug fix in mixture models
04jan2013 version 2.0.0 - maxtime() added to specify maximum generated survival time and event indicator
                        - n() removed so must set obs
                        - loghazard() and hazard() added for user-defined hazard functions, simulated using quadrature and root finding
                        - default centol() changed to 1E-08
                        - varnames can now appear in loghazard()/hazard() allowing time-dependent covariates
                        - mixture models now use Brent method -> much more reliable than NR, and allows tdes
                        - cumhazard() and logcumhazard() now added which just use root finding
15Nov2011 version 1.1.2 - Fixed bug when generating covariate tempvars with competing risks.
20sep2011 version 1.1.1 - Exponential distribution added.
10sep2011 version 1.1.0 - Added Gompertz distribution. Time-dependent effects available for all models except mixture. showerror option added.
09sep2011 version 1.0.1 - Time dependent effects now allowed for standard Weibull.
*/

program define survsim
	version 14.2

	CheckObs

	local opts1											/// -fitted merlin model-
                        MODel(passthru)				///
														//

	local opts2											///	-user-defined-
                        LOGHazard(passthru)		///
                        Hazard(passthru)		///
                        LOGCHazard(passthru)		///
                        CHazard(passthru)		///
                        NODES(passthru)			///	-# nodes-
                        TDEFUNCtion(passthru)		///	-function of time to interact with time-dependent effects-
                        MIXture				///	-two-component mixture-
                        PMix(passthru)			///	-mixture parameter-
                        CR				///	-simulate competing risks-
                        NCR(string)			///	-number of competing risks-
                        Distribution(passthru)		///
                        Lambdas(passthru)		///
                        Gammas(passthru)		///
                        COVariates(passthru)		///	-baseline covariates, e.g, (sex 0.5 race -0.4)-
                        TDE(passthru)			///	-time dependent effects-
                        MARGinal			///
                                                        //

	local commonopts	 			///
                        MAXTime(string)			///	-right censoring time-
                        LTruncated(string)		///	-left truncation/delayed entry-
							//

	//==============================================================================================================//
	//parse opts that are common to all three settings

		syntax newvarname(min=1 max=3) , [`commonopts' *]
		local newvars `varlist'
		local stime : word 1 of `newvars'
		local died  : word 2 of `newvars'
		local hasmaxtime = "`maxtime'"!=""

		if `hasmaxtime' {

			if "`died'"=="" {
				di as error "2 new variable names required"
				exit 198
			}

			cap confirm number `maxtime'
			if _rc {
				cap confirm numeric variable `maxtime'
				if _rc {
					di as error "Invalid maxtime()"
					exit 198
				}
			}
			else {
				if `maxtime'<0 {
					di as error "maxtime() must be >0"
					exit 198
				}
			}

			tempvar maxtvar
			gen `maxtvar' = `maxtime'
			local maxtopt maxtime(`maxtvar')
		}

		if "`ltruncated'"!="" {

			cap confirm number `ltruncated'
			if _rc {
				cap confirm numeric variable `ltruncated'
				if _rc {
					di as error "Invalid ltruncated()"
					exit 198
				}
			}
			else {
				if `ltruncated'<0 {
					di as error "ltruncated() must be >0"
					exit 198
				}
			}

			tempvar ltvar
			gen `ltvar' = `ltruncated'
			local ltopt ltruncated(`ltvar')
		}

		local 0 , `options'

	//==============================================================================================================//
	//survsim_model

		syntax , [`opts1' *]

		if "`model'"!="" {
			capture survsim_model `newvars', 	`model' 		///
												`maxtopt'		//
			local rc = c(rc)
			if `rc' {
				cap drop _survsim_rc
				cap drop `stime'
				cap drop `died'
                                exit `rc'
			}
			else {
				RCREPORT 	_survsim_rc
				MISSREPORT 	`stime'
				GENEVENT 	`stime' `died' `maxtime'
                                exit
			}
		}

		local 0 , `options'

	//==============================================================================================================//
	//survsim_user

		syntax , [`opts2' *]

		local useropt "`loghazard'`hazard'`logchazard'`chazard'`mixture'"

		if "`useropt'"!="" {
			capture survsim_user `newvars', 	`useropt' 		///
												`nodes'			///
												`maxtopt' 		///
												`ltopt'			///
												`covariates'	///
												`tde'			///
												`tdefunction'	///
												`distribution'	///	-for mixture-
												`lambdas'		///
												`gammas'		///
												`pmix'			///
												`marginal'		//
            local rc = c(rc)
			if `rc' {
				cap drop _survsim_rc
				cap drop `stime'
				cap drop `died'
                                exit `rc'
			}
			else {
				RCREPORT 	_survsim_rc
				MISSREPORT 	`stime'
				GENEVENT 	`stime' `died' `maxtime'
                                exit
			}
		}

		local 0 , `options'

	//==============================================================================================================//
	//survsim_msm

		syntax , [HAZARD1(passthru) HAZARD2(passthru) TRANSMATrix(passthru) STARTSTATE(passthru) *]

		if "`hazard1'"!="" {

			capture program drop survsim_msm			// refreshes Mata functions

			survsim_msm `newvars', 		`transmatrix'	///
										`hazard1'		///
										`hazard2'		///
										`maxtopt' 		///
										`ltopt'			///
										`startstate'	///
										`nodes'			///
										`options' 		//
			exit

		}

		local 0 , `options' `distribution' `lambdas' `gammas' `covariates' `tde'

	//==============================================================================================================//
	//parametric distribution

		syntax , 	Distribution(string) 		///
					Lambdas(string) 			///
				[								///
					Gammas(string) 				///
					COVariates(string)			///
					TDE(string)					///
					`opts3'						///
				]								//

		local ld = length("`distribution'")
		if 		substr("exponential",1,max(1,`ld'))=="`distribution'" {
			local dist "exp"
		}
		else if substr("gompertz",1,max(3,`ld'))=="`distribution'" {
			local dist "gompertz"
		}
		else if substr("weibull",1,max(1,`ld'))=="`distribution'" {
			local dist "weibull"
		}
		else {
			di as error "Unknown distribution"
			exit 198
		}

		foreach l of numlist `lambdas' `gammas' {
			if `l'<0 {
				di as error "lambdas()/gammas() must be > 0"
				exit 198
			}
		}

		if "`dist'"=="exp" & "`gammas'"!="" {
			di as error "gammas cannot be specified with distribution(exponential)"
			exit 198
		}

	//==============================================================================================================//
	//baseline covariates and time-dependent effects

		if "`covariates'"!="" {

			tokenize `covariates'
			local ncovlist : word count `covariates'
			local ncovvars = `ncovlist'/2
			cap confirm integer number `ncovvars'
			if _rc>0 {
				di as error "Variable/number missing in covariates"
				exit 198
			}
			local ind = 1
			local error = 0
			forvalues i=1/`ncovvars' {
				cap confirm var ``ind'', exact
				if _rc {
					local errortxt "invalid covariates(... ``ind'' ``=`ind'+1'' ...)"
					local error = 1
				}
				cap confirm num ``=`ind'+1''
				if _rc {
					local errortxt "invalid covariates(... ``ind'' ``=`ind'+1'' ...)"
					local error = 1
				}
				tempvar vareffect`i'
				gen double `vareffect`i'' = ``ind''*``=`ind'+1''

				local ind = `ind' + 2
			}
			if `error' {
				di as error "`errortxt'"
				exit 198
			}
			local cov_linpred "`vareffect1'"
			if `ncovvars'>1 {
				forvalues k=2/`ncovvars' {
					local cov_linpred "`cov_linpred' + `vareffect`k''"
				}
			}
			local cov_linpred "* exp(`cov_linpred')"

		}

		if "`tde'"!="" {

			tokenize `tde'
			local ntde : word count `tde'
			local ntdevars = `ntde'/2
			cap confirm integer number `ntdevars'
			if _rc>0 {
				di as error "Variable/number missing in tde"
				exit 198
			}

			local ind = 1
			local error = 0
			forvalues i=1/`ntdevars' {
				cap confirm var ``ind'', exact
				if _rc {
					local errortxt "invalid tde(... ``ind'' ``=`ind'+1'' ...)"
					local error = 1
				}
				cap confirm num ``=`ind'+1''
				if _rc {
					local errortxt "invalid tde(... ``ind'' ``=`ind'+1'' ...)"
					local error = 1
				}
				tempvar tdeeffect`i'
				gen double `tdeeffect`i'' = ``ind''*``=`ind'+1''

				local ind = `ind' + 2
			}
			if `error' {
				di as error "`errortxt'"
				exit 198
			}
			local tde_linpred "`tdeeffect1'"
			if `ntdevars'>1 {
				forvalues k=2/`ntdevars' {
					local tde_linpred "`tde_linpred' + `tdeeffect`k''"
				}
			}
			local tde_linpred "+ `tde_linpred'"

		}

	//==============================================================================================================//
	//stime and died

		tempvar u
		qui gen double `u' = runiform()

		if 		"`dist'"=="exp" {
			if "`ltruncated'"!="" {
				local ltcontrib "* exp(-`lambdas' `cov_linpred' * `ltvar' ^ (1 `tde_linpred') / (1 `tde_linpred'))"
			}
			qui gen double `stime' 	= (-ln(`u' `ltcontrib')*(1 `tde_linpred')/(`lambdas' `cov_linpred'))^(1/(1 `tde_linpred'))
		}
		else if "`dist'"=="weibull" {
			if "`ltruncated'"!="" {
				local ltcontrib "* exp(-`lambdas' `cov_linpred' * `gammas' * `ltvar' ^ (`gammas' `tde_linpred') / (`gammas' `tde_linpred'))"
			}
			qui gen double `stime' 	= (-ln(`u' `ltcontrib')*(`gammas' `tde_linpred')/(`lambdas'*`gammas' `cov_linpred'))^(1/(`gammas' `tde_linpred'))
		}
		else{
			if "`ltruncated'"!="" {
				local ltcontrib "* exp(- `lambdas' `cov_linpred' / (`gammas' `tde_linpred') * (exp(`ltvar' * (`gammas' `tde_linpred') ) - 1))"
			}
			qui gen double `stime' 	= (1/(`gammas' `tde_linpred'))*log(1-(((`gammas' `tde_linpred')*log(`u' `ltcontrib'))/(`lambdas' `cov_linpred')))
		}

		MISSREPORT `stime'
		if "`maxtime'"!="" {
			GENEVENT `stime' `died' `maxtvar'
			qui replace `stime' = `maxtvar' if `died'==0
		}

end

program CheckObs
	if _N==0 {
		di as error "No observations"
		exit 198
	}
end

program GENEVENT
	args stime died maxtime
	qui gen byte `died' = `stime'<`maxtime' if `stime'!=.
	qui replace `stime' = `maxtime' if `stime'>`maxtime' & `stime'!=.
end

program RCREPORT
	args rcvar
	qui su `rcvar' if `rcvar'==1, meanonly
	if r(N)>0 {
		di in yellow "Warning: `r(N)' survival times did not converge"
		di in yellow "         They have been set to the final iteration value"
		di in yellow "         You can identify them by _survsim_rc = 1"
	}
	qui su `rcvar' if `rcvar'==2, meanonly
	if r(N)>0 {
		di in yellow "Warning: `r(N)' survival times were below the lower limit"
		di in yellow "         You can identify them by _survsim_rc = 2"
	}
	qui su `rcvar' if `rcvar'==3, meanonly
	if r(N)>0 {
		di in yellow "Warning: `r(N)' survival times were above the upper limit of maxtime()"
		di in yellow "         They have been set to maxtime()"
		di in yellow "         You can identify them by _survsim_rc = 3"
	}
end

program MISSREPORT
	args time
	qui count if `time'==.
	if `r(N)'>0 {
		di in yellow "`r(N)' missing values generated in simulated survival times"
	}
end