*! version 1.1.2 12Jun2018

// checked RMST DM and Mest with bootstrapping

// work on other scales
// make certain transformations the default
// check rel surv models work for M-estimation
// error check for factor variables??

// error checks for centiles
// test ifs....
// weights for extenal standardisation.

// add weights

program define stpm2_standsurv, rclass
	version 12.2
	syntax [if] [in],			 							///
		[													///
		ATVars(string)										/// list or stub for at variables
		CONTRASTVars(string)								/// list or stub for contrasts
		LINCOMVar(string)									///	name of new variables for linear combination
		USERFUNCTIONVar(string)								///	name of new variables for userfunction
		ATReference(integer 1)								/// reference at() - default 1
		TImevar(varname) 									/// timepoints for predictions
		CONtrast(string)	 								/// type of contrast
		TRansform(string)									/// Transformation for variance calculation
		CENTILE(numlist)									///	centiles of standardized curve
		HAZard												/// standardized hazard function
		RMST												/// restricted mean survival time
		SURVival											/// standardized survival function (default)
		CI 													/// request CI to be calculated
		PER(integer 1)										/// per option (multiply to give pys etc)
		FAILure												/// calculate failure function
		CENTVar(string)										/// name of new centile variable (default _centvar)
		CENTILEUpper(real -99)								/// starting value for upper bound of centile search
		NOdes(integer 30)									/// number of nodes for numerical integration
		MESTimation											/// use M-estimation
		LEVel(real `c(level)')								/// level for CIs
		SE													/// calculate standard error
		LINCOM(string)										///
		INDWeights(string)									/// multiply observations by weights
		VERBOSE												/// show what is happening (speed tests)
		RUNMATA												///
		USERFunction(string)								/// user defined function
		*													/// atn() options
		]

	tempvar touse_time touse_model xb expxb dxb  expxb0
	marksample touse, novarlist
// default to standardized survival function
	if wordcount("`centile' `hazard' `survival' `rmst' `failure'") == 0 local survival "survival"
	if wordcount("`hazard' `survival' `rmst' `failure'") + ("`centile'" != "") > 1 {
		di as error "You can only specify one of the survival, hazard, centile and rmst options"
		exits 198
	}
	if "`e(cmd)'" != "stpm2" {
		di as error "You need to fit an stpm2 model to use stpm2_standsurv"
		exit 198
	}
// check moremata installed
	capture findfile lmoremata.mlib
	if _rc {
		display in yellow "You need to install moremata to use stpm2_stand"
		display in yellow "Type {stata ssc install moremata}, or just click on the link"
        exit  198
	}
	summ _t if _d==1 & e(sample), meanonly
	local maxt `r(max)'

// Check no factor variables
	fvexpand `e(varlist)'
	if "`r(fvops)'" == "true" {
		di as error "stpm2_standsurv does not allow factor variables"
		di as error "Refit the model using dummy variables etc"
		exit 198
	}


// M-estimation for centiles
	if "`centile'" != "" local mestimation mestimation

// Extract at() options
	local optnum 1
	local end_of_ats 0
	local 0 ,`options'
	while `end_of_ats' == 0 {
		capture syntax [,] AT`optnum'(string) [*]
		if _rc {
			local N_at_options = `optnum' - 1
			local end_of_ats 1
			continue, break
		}
		else local 0 ,`options'
		local optnum = `optnum' + 1
	}
	local N_at_options = `optnum' - 1
	if "`0'" != "," {
		di as error "Illegal option: `0'"
		exit 198
	}
	local hasatoptions = `N_at_options' > 0
	if !`hasatoptions' local N_at_options 1

// Parse at() options
// Probably does not work with factor variables
	if `hasatoptions' > 0 {
		forvalues i = 1/`N_at_options' {
	// parse "if" suboption
			tokenize "`at`i''", parse(",")
			local at`i'opt  `1'
			local atoptif `3'
			local 0 ,`3'
			syntax ,[if2(string)]
			if `"`if'"' != "" & `"`atoptif'"' != "" {
				di as error "You can either use an if statement or the if suboptions" _newline ///
							"of the at() options"
				exit 198
			}
			tempvar touse_at`i'
			if `"`atoptif'"' == "" {
				gen byte `touse_at`i'' = `touse'
			}
			else {
				gen byte `touse_at`i'' = (`if2')
			}

			tokenize `at`i'opt'
			while "`1'"!="" {
				if "`1'" == "." continue, break
				fvunab tmpfv: `1'
				local 1 `tmpfv'
				cap confirm var `1'
				if _rc {
					di "`1' is not in the data set"
				}
				local at`i'vars `at`i'vars' `1'
				if "`2'" != "=" {
					cap confirm num `2'
					if _rc {
						di as err "invalid at(... `1' `2' ...)"
						exit 198
					}
					local at`i'_`1'_value `2'
					mac shift 2

				}
				else {
					cap confirm var `3'
					if _rc {
						di as err "`var' is not in the data set"
						exit 198
					}
					local at`i'_`1'_value .
					local at`i'_`1'_variable `3'
					mac shift 3
				}
			}
		}
	}
	else {
		tempvar touse_at1
		gen byte `touse_at1' = `touse'
	}

// Number of observations for each at() option
	local varsinmodel = subinstr(strtrim(stritrim("`e(varlist)' `e(tvc)'"))," ",",",.)
	forvalues i = 1/`N_at_options' {
		quietly count if `touse_at`i'' == 1 & !missing(`varsinmodel')
		local Nobs_predict_at`i' `r(N)'
		local touse_at_list `touse_at_list' `touse_at`i''
	}
// names of new variables
	if "`atvars'" == "" {
		forvalues i = 1/`N_at_options' {
			local at_varnames `at_varnames' _at`i'
		}
	}
	else {
		capture _stubstar2names double `atvars', nvars(`N_at_options')
		local at_varnames `s(varlist)'
		if _rc>0 {
			di as error "atvars() option should either give `N_at_options' new variable names " ///
				"or use the {it:stub*} option. The specified variable(s) probably exists."
			exit 198
		}
	}
	if "`contrastvars'" == "" {
		forvalues i = 1/`N_at_options' {
			if `i' == `atreference' continue
			local contrast_varnames `contrast_varnames' _contrast`i'_`atreference'
		}
	}
	else {
		capture _stubstar2names double `contrastvars', nvars(`=`N_at_options'-1')
		local contrast_varnames  `s(varlist)'
		if _rc>0 {
			di as error "contrastvars() option should either give `=`N_at_options'-1' new variable names " ///
				"or use the {it:stub*} option. The specified variable(s) probably exists."
			exit 198
		}
	}
	if "`lincomvar'" == "" {
		local lincom_varname _lincom
	}
	else local lincom_varname `lincomvar'

	if "`userfunctionvar'" == "" {
		local userfunction_varname _userfunc
	}
	else local userfunction_varname `userfunctionvar'

	if `atreference' != 1 {
		if !inrange(`atreference',1,`N_at_options') {
			di as error "atreference option out of range"
			exit 198
		}
	}
// Transform option
	if "`transform'" == "" local transform log
	if !inlist("`transform'","loglog","logit","log","none") {
		di as error "Transform options are none, log, loglog or logit"
		exit 198
	}
// Number of observations used in the model
	quietly gen `touse_model' = e(sample)
	quietly count if `touse_model' == 1
	local Nobs_model `r(N)'

// time variable
	if "`timevar'" == "" local timevar _t
	gen byte `touse_time' = `timevar' != .

// currently only work with scale(hazard)
if "`e(scale)'" != "hazard" {
	di as error "stpm2_standsurv only currently works with scale(hazard) models."
	exit 198
}

// Check contrast option
	if "`contrast'" != "" {
		if !inlist("`contrast'","difference","ratio","pchange") {
			di as error "contrast option should either be difference or ratio"
			exit 198
		}
	}

// **** CHECK THESE ARE NEEDED WITH DELTA METHOD *****
// Use meansurv for point estimates
	if "`indweights'" != "" {
		local weight_option meansurvwt(`indweights')
	}

// exception if using meansurv with no contrast (use predict, meansurv as faster)
	if "`survival'" != "" & "`mestimation'" == "" & "`contrast'" == "" & "`lincom'" == "" & "`userfunction'" == "" & "`runmata'" == "" local meansurv_nomata 1

	if "`verbose'" != "" di in yellow "Predicting point estimates using meansurv"
	if "`meansurv_nomata'" == "" {
		if "`survival'" != "" | "`hazard'" != "" | "`failure'" != "" {
			forvalues i = 1/`N_at_options' {
				tempvar S`i'
				if `hasatoptions' {
					if ("`at`i'opt'" == ".") local tempatopt
					else local tempatopt at(`at`i'opt')
				}
				quietly predict `S`i'' if `touse_at`i'', meansurv `tempatopt' timevar(`timevar') `weight_option'
				local Smean_list `Smean_list' `S`i''
			}
		}
		if "`hazard'" != "" {
			forvalues i = 1/`N_at_options' {
				tempvar h`i' f`i'
				if `hasatoptions' {
					if ("`at`i'opt'" == ".") local tempatopt
					else local tempatopt at(`at`i'opt')
				}
				quietly predict `f`i'' if `touse_at`i'', meanft `tempatopt' timevar(`timevar') `weight_option'
				local fmean_list `fmean_list' `f`i''
			}
		}

		if "`rmst'" != "" {
			forvalues i = 1/`N_at_options' {
				tempvar rmst`i'
				if `hasatoptions' {
					if ("`at`i'opt'" == ".") local tempatopt
					else local tempatopt at(`at`i'opt')
				}
				local rmst_newvarname `rmst`i''
				mata: rmst_stand()
				local rmstmean_list `rmstmean_list' `rmst`i''
			}
		}
		// Need rootfinder to get centiles of standardized curves
		if "`centile'" != "" {
			if `centileupper' == -99 local centileupper = `maxt'*2
			if "`centvar'" == "" local centvar _centvals
			confirm new var `centvar'
			qui gen `centvar' = .
			local i = 1
			foreach c of numlist `centile' {
				qui replace `centvar' = `c' in `i'
				local ++i
			}
			tempvar touse_centiles
			gen byte `touse_centiles' = `centvar' != .

			forvalues i = 1/`N_at_options' {
				tempvar c`i'
				if `hasatoptions' {
					if ("`at`i'opt'" == ".") local tempatopt
					else local tempatopt at(`at`i'opt')
				}
				local newcentvar `c`i''
				mata: stand_centile()
				local centile_list `centile_list' `c`i''
			}
		}

	// predict xb and dxb (need for CIs)
		if ("`ci'" != "" | "`se'" != "") & "`mestimation'" != "" {
			if "`verbose'" != "" di in yellow "Predicting xb and dxb"
			quietly predict `xb'  if e(sample), xb
			quietly predict `dxb' if e(sample) , dxb
			if `e(del_entry)' & "`mestimation'" != "" {
				qui _predict `expxb0' if _t0>0 & e(sample), eq(xb0)
				qui replace `expxb0' = exp(`expxb0') if _t0>0 & e(sample)
				qui replace `expxb0' = 0 if _t0==0 & e(sample)
			}
		}
		if "`verbose'" != "" di in yellow "Calling main mata program"
		mata: standsurv()
	}
	else {
		forvalues i = 1/`N_at_options' {
			if `hasatoptions' {
				if ("`at`i'opt'" == ".") local tempatopt
				else local tempatopt at(`at`i'opt')
			}
			local newv = word("`at_varnames'",`i')
			quietly predict `newv' if `touse_at`i'', meansurv `tempatopt' timevar(`timevar') `weight_option' `ci'
		}
	}
// Warnings
	if "`centile'" != "" {
		foreach var in `at_varnames' {
			quietly count if (`var'>`maxt') & (`touse_centiles')
			if `r(N)'>0 {
				di as result "Warning: centile point estimate for `var' > maximum event time"
			}
			if "`ci'" != "" {
				quietly count if (`var'_uci>`maxt') & (`touse_centiles')
				if `r(N)'>0 {
					di as result "Warning: CI for centile for `var' > maximum event time"
				}
			}
		}
	}

// Return stuff
	return local varmethod=cond("`mestimation'" == "","delta-method","M-estimation")
end

// note need to run build_ado to put mata programs at the end of this file
mata





//======================================================================================================================================//
// Main structure for main analysis
struct stpm2_standardisation {
		real scalar		N_at_options, 		// number of at options
						hasatoptions,		// has at options
						Nobs_model,			// N in model
						calchazard,			// calculate hazard
						calcsurvival,		// calculate survival
						calcfailure,		// calculate failure function
						calcrmst,			// calculate rmst
						calccentile,		// calculate centiles
						deltamethod,		// use delta method
						mestimation,		// use M estimation
						verbose,			// display stage of estimation (speed tests)
						loopmax,			// loop for calculating centiles
						at_reference,		// reference at level
						hascontrast,		// contrast option
						haslincom,			// has linear combination
						hasuserfunction,	// has userfunction combination
						hasweights,			// as indweight option
						se,					// calculate standard error
						ci,					// indicator to calculate CIs
						per,				// multiply estimates by per
						level,				// level for CI (default 95%)
						dots,				// print dots as loop over t
						hasfailure,			// has the failure option,
						hascons,			// has a constant term in the model
						hastvc,				// has time-dependent effects
						hasdel_entry,		// has delayed entry
						hasbhazard,			// has bhazard option
						orthog,				// orthogonalisation of splines
						rcsbaseoff,			// rcsbaseoff option specified
						Ntvc,				// Number of tvc variables
						Nt,					// Number of time points to predict at
						Nvarlist,			// Number of main effects (excluding constant)
						hasvarlist,			// Has a varlist
						Nparameters,		// Number of parameters
						df,					// df for baseline
						Ncentiles,
						Nnodes,				// Number of nodes for integration for rmst
						j


		string scalar	touse_time,			// indicator for timevar
						touse_model,		// touse e(sample)
						touse_centiles,
						contrast,			// type of contrast
						transform,			// type of transformation for CI
						scale				// scale option for model

		real matrix 	Nobs_predict_at,	// Number of observations for each at() option
						Smean,				// Standardised survival
						hmean,				// Standardised hazard
						fmean,				// mean of f(t)
						rmstmean,			// rmst of standardized survival
						centmean,
						xb,					// linear predictor
						expxb,				// exp of linear predictor
						dxb,				// derivative of linear predictor
						expxb0,				// exp of linear predictor for delayed entry
						d,					// event indicator
						obst,				// observed time (_t)
						t,					// times to predict at
						bhazard,			// expected hazard rate
						touse_at,			// touse option for each at option
						X,					// design matrix
						Xdrcs,				// derivative for design matrix
						nodeweights,		// weights for numerical integration
						nodes,				// nodes for numerical integration
						beta,				// beta coefficients
						betarcs,			// spline coefficients
						dftvc,				// df for time-dependent effects
						V,					// Variance matrix
						tj,					// current value of time
						Si,					// Individual contribution to mean survival
						hi,					// Individual contribution to mean hazard
						fi,					// Contribution to density
						rmsti,				// contribution to RMST
						ResVcov,			// Variance
						Vest,				//
						lincom,				// linear combination of at options
						weights				// weight (for external standardization)

		string matrix	at_vars,			// Name of variables for each at() option.
						at_varnames,		// Names of new variables
						tvcnames,			// Names of tvc variables
						contrast_varnames,	// Names of contrast variables
						lincom_varname,		// Names of linear combination
						userfunction_varname // Name of user_functions

	transmorphic matrix	knots,				// array for knot locations
						Rmatrix,			// Array for Rmatrix
						X_at,				// Change in X matrix for at() options
						X_at_index,			// Index of covariate that needs to change
						X_at_tvc,			// Changed  tvc covariates for at() options
						X_at_t,				// change for at within loop
						Xrcstvc,			// tvc splines while looping over t
						Xdrcstvc,			// spline derivatives while looping over t
						Xdrcs_at_t,			// spline derivatives while looping over t
						Sti_nodes,			// survival function at nodes
						dxbi_nodes,
						Xrcsbase,			// Used when looping over t
						Xdrcsbase,			// Used when looping over t
						rcs_nodes

	pointer				GenA11,				// Generate A11 depending on standardisation option
						GenA22,				// Generate A22 depending on standardisation option
						GenA12,				// Generate A12 depending on standardisation option
						GenU,
						Stgen,				// generate survival depending on scale
						htgen,				// generate hazard depending on scale
						RMSTgen,			// generate RMST depending on scale
						Gen_Sderiv,			// derivative of survival function
						GenUbeta,
						Gen_fderiv,
						Gen_RMSTderiv
}


//======================================================================================================================================//
// Structure for Results
struct stpm2_standardisation_results {
	real matrix CI_at_lci,				// array for at CI
				CI_at_uci,
				at_est,
				se,
				CI_contrast_lci,
				CI_contrast_uci,
				contrast_est,
				CI_lincom_lci,
				CI_lincom_uci,
				lincom_est,
				userfunction_est,
				CI_userfunction_lci,
				CI_userfunction_uci

	pointer scalar userfunction
}

//======================================================================================================================================//
// GetStuff() - get options and data and store in struct
function GetStuff() {
	struct stpm2_standardisation scalar S
	S.verbose = st_local("verbose") != ""
	if(S.verbose) display("Reading in things to set up structure")
	S.N_at_options = strtoreal(st_local("N_at_options"))
	S.hasatoptions = S.N_at_options >1
	S.touse_time = st_local("touse_time")
	S.touse_model = st_local("touse_model")
	S.Nobs_model = strtoreal(st_local("Nobs_model"))
	S.calchazard = st_local("hazard") != ""
	S.calcsurvival = st_local("survival") != ""
	S.calcfailure = st_local("failure") != ""
	S.calcrmst = st_local("rmst") != ""
	S.calccentile = st_local("centile") != ""
	S.deltamethod = st_local("mestimation") == ""
	S.mestimation = 1- S.deltamethod
	if(S.calccentile) S.loopmax = S.N_at_options
	else S.loopmax = 1
	S.Nobs_predict_at = J(1,S.N_at_options,.)
	S.at_vars = J(1,S.N_at_options,"")
	for(i=1;i<=S.N_at_options;i++) {
		S.Nobs_predict_at[1,i] = strtoreal(st_local("Nobs_predict_at"+strofreal(i)))
		S.at_vars[1,i] = st_local("at"+strofreal(i)+"vars")
	}
	S.at_reference = strtoreal(st_local("atreference"))
	S.at_varnames = tokens(st_local("at_varnames"))
	S.hascontrast = st_local("contrast") != ""
	S.haslincom = st_local("lincom") != ""
	S.hasuserfunction = st_local("userfunction") != ""
	S.hasweights = st_local("indweights") != ""
	S.contrast = st_local("contrast")
	S.ci = st_local("ci") != ""
	S.se = st_local("se") != ""
	S.level = strtoreal(st_local("level"))
	S.per = strtoreal(st_local("per"))
	S.transform = st_local("transform")
	S.hasfailure = st_local("failure") != ""
	S.hascons = st_global("e(noconstant)") == ""
	S.hastvc = st_global("e(tvc)") != ""
	S.hasdel_entry = st_numscalar("e(del_entry)")
	S.hasbhazard = st_global("e(bhazard)") != ""
	S.orthog = st_global("e(orthog)") != ""
    S.rcsbaseoff = st_global("rcsbaseoff") != ""
	S.df = st_numscalar("e(dfbase)")
	S.scale = st_global("e(scale)")
	if (S.hastvc) {
		S.tvcnames = tokens(st_global("e(tvc)"))
		S.Ntvc = cols(S.tvcnames)
	}
	else S.Ntvc = 0
	if(S.hascontrast) S.contrast_varnames = tokens(st_local("contrast_varnames"))
	if(S.haslincom)	S.lincom_varname = st_local("lincom_varname")
	if(S.hasuserfunction) S.userfunction_varname = st_local("userfunction_varname")
//======================================================================================================================================//
// Variables created in Stata
	if(S.calcsurvival | S.calchazard | S.calcfailure) S.Smean = st_data(.,(st_local("Smean_list")),S.touse_time)
	if(S.calchazard) {
		S.fmean = st_data(.,(st_local("fmean_list")),S.touse_time)
		S.hmean = S.fmean:/S.Smean
	}
	if(S.calcrmst) {
		S.rmstmean = st_data(.,(st_local("rmstmean_list")),S.touse_time)
	}
	if(S.calccentile) {
		S.touse_centiles = st_local("touse_centiles")
		S.centmean = st_data(.,(st_local("centile_list")),S.touse_centiles)
		S.Ncentiles = rows(S.centmean)
	}
	S.t = st_data(.,st_local("timevar"),S.touse_time)
	S.Nt = rows(S.t)
	if(S.ci | S.se) {
		if(S.mestimation) {
			S.xb = st_data(.,st_local("xb"),S.touse_model)
			S.expxb = exp(S.xb)
			S.dxb = st_data(.,st_local("dxb"),S.touse_model)
			S.d = st_data(.,"_d",S.touse_model)
			if(S.hasbhazard) {
				S.bhazard = st_data(.,st_global("e(bhazard)"),S.touse_model)
				S.obst = st_data(.,"_t",S.touse_model)
			}
			if(S.hasdel_entry) S.expxb0 = st_data(.,st_local("expxb0"),S.touse_model)

		}
		S.touse_at = st_data(.,st_local("touse_at_list"),S.touse_model)

		if(S.hasweights) S.weights = st_data(.,st_local("indweights"),S.touse_model)
		else S.weights = 1

		//======================================================================================================================================//
// get knot locations and df for each knot
		S.knots = asarray_create()
		S.dftvc = J(1,S.Ntvc,.)
		if(!S.rcsbaseoff) asarray(S.knots,"baseline",strtoreal(tokens(st_global("e(ln_bhknots)"))))
		if(S.hastvc) {
			for(i=1;i<=S.Ntvc;i++) {
				asarray(S.knots,S.tvcnames[i],strtoreal(tokens(st_global("e(ln_tvcknots_"+S.tvcnames[i]+")"))))
				S.dftvc[1,i] = st_numscalar("e(df_"+S.tvcnames[i]+")")
			}
		}

//======================================================================================================================================//
// get R matrices
		S.Rmatrix = asarray_create()
		if(S.orthog & !S.rcsbaseoff) asarray(S.Rmatrix,"baseline",st_matrix("e(R_bh)"))
		else asarray(S.Rmatrix,"baseline",J(0,0,.))
		if(S.hastvc) {
			for(i=1;i<=S.Ntvc;i++) {
				if(S.orthog) asarray(S.Rmatrix,S.tvcnames[i],st_matrix("e(R_"+S.tvcnames[i]+")"))
				else asarray(S.Rmatrix,S.tvcnames[i],J(0,0,.))
			}
		}

//======================================================================================================================================//
// Observed X matrix and derivative of spline functions
		S.Nvarlist = cols(tokens(st_global("e(varlist)")))
		S.hasvarlist = S.Nvarlist>0
		covariates = J(1,0,"")
		drcsvars = J(1,0,"")
		if(S.Nvarlist > 0) covariates = covariates, tokens(st_global("e(varlist)"))
			if(!S.rcsbaseoff) {
			covariates = covariates, tokens(st_global("e(rcsterms_base)"))
			drcsvars = drcsvars, tokens(st_global("e(drcsterms_base)"))
		}
		if(S.hastvc) {
			for(i=1;i<=S.Ntvc;i++) {
				covariates = covariates, tokens(st_global("e(rcsterms_"+S.tvcnames[i]+")"))
				drcsvars = drcsvars, tokens(st_global("e(drcsterms_"+S.tvcnames[i]+")"))
			}
		}
		S.X = st_data(.,covariates,S.touse_model) //
		if(S.hascons) S.X = S.X,J(S.Nobs_model,1,1)
		S.Xdrcs = st_data(.,drcsvars,S.touse_model)

		S.Nparameters = cols(S.X)
		S.Nparameters = cols(covariates) + S.hascons
		S.beta = st_matrix("e(b)")'[1..S.Nparameters,1]
		S.betarcs = S.beta[(S.Nvarlist+1)..(S.Nparameters - S.hascons)]
		S.V = st_matrix("e(V)")[1..S.Nparameters,1..S.Nparameters]

//======================================================================================================================================//
// changes to X matrix needed for at() options
		S.X_at = asarray_create()
		S.X_at_index = asarray_create()
		if(S.hastvc) S.X_at_tvc = asarray_create()
		for(i=1;i<=S.N_at_options;i++) {
			ati = "at" + strofreal(i)
			asarray(S.X_at_index,ati,J(1,0,.))
			asarray(S.X_at,ati,J(S.Nobs_model,0,.))
			// main effects
			for(j=1;j<=S.Nvarlist;j++) {
				if(subinword(S.at_vars[1,i],covariates[j],"") != S.at_vars[1,i]) {
					asarray(S.X_at_index,ati, (asarray(S.X_at_index,ati), j))
					atval_macro = st_local("at" + strofreal(i) + "_" + covariates[j]+"_value")
					if(atval_macro != ".") {
						asarray(S.X_at,ati,(asarray(S.X_at,ati), J(S.Nobs_model,1,strtoreal(atval_macro))))
					}
					else {
						asarray(S.X_at,ati,(asarray(S.X_at,ati), st_data(.,st_local("at" + strofreal(i) + "_" + covariates[j]+"_variable"),S.touse_model)))
					}
				}
			}
			// tvcs
			if(S.hastvc) {
				asarray(S.X_at_tvc,ati, J(S.Nobs_model,0,.))
				for(j=1;j<=S.Ntvc;j++) {
					if(subinword(S.at_vars[1,i],S.tvcnames[j],"") != S.at_vars[1,i]) {
						atval_macro = st_local("at" + strofreal(i) + "_" + S.tvcnames[j]+"_value")
						asarray(S.X_at_tvc, ati, (asarray(S.X_at_tvc, ati),J(S.Nobs_model,1,strtoreal(atval_macro))))
					}
					else {
						asarray(S.X_at_tvc,ati, (asarray(S.X_at_tvc,ati), st_data(.,S.tvcnames[j],S.touse_model)))
					}
				}
			}
		}
	}
// linear combination
	if(S.haslincom) {
		S.lincom = strtoreal(tokens(st_local("lincom")))
	}

// nodes for rmst
	if(S.calcrmst) {
		S.Nnodes = strtoreal(st_local("nodes"))
		gq(S.nodeweights,S.nodes)
	}
	if(S.verbose) display("Finished setting up structure")
	return(S)
}

//======================================================================================================================================//
// Get DeclarePointers - declare all pointers

function DeclarePointers(struct stpm2_standardisation scalar S)
{
	if(S.calcsurvival) {
		S.GenA11 = &GenA11_surv()
		S.GenA22 = &GenA22_surv()
		S.GenA12 = &GenA12_surv()
		S.GenU = &GenU_surv()
	}
	if(S.calcfailure) {
		S.GenA11 = &GenA11_surv()
		S.GenA22 = &GenA22_surv()
		S.GenA12 = &GenA12_surv()
		S.GenU = &GenU_failure()
	}
	else if(S.calchazard) {
		S.GenA11 = &GenA11_hazard()
		S.GenA22 = &GenA22_hazard()
		if(S.mestimation) S.GenA12 = &GenA12_hazard_mestimation()
		else S.GenA12 = &GenA12_hazard_deltamethod()
		S.GenU = &GenU_hazard()
	}
	else if(S.calcrmst) {
		S.GenA11 = &GenA11_surv()		// same as survival
		S.GenA22 = &GenA22_surv()		// same as survival
		S.GenA12 = &GenA12_rmst()
		S.GenU = &GenU_rmst()
	}
	else if(S.calccentile) {
		S.GenA11 = &GenA11_centile()
		S.GenA22 = &GenA22_surv()		// same as survival
		S.GenU = &GenU_centile()
		S.GenA12 = &GenA12_surv()		// same as survival
	}
// scale option
	if(S.scale == "hazard") {
		S.Stgen = &Stgen_hazard()
		S.htgen = &htgen_hazard()
		S.RMSTgen = &RMSTgen_hazard()
		S.Gen_Sderiv = &Gen_Sderiv_hazard()
		S.Gen_fderiv = &Gen_fderiv_hazard()
		S.Gen_RMSTderiv = &Gen_RMSTderiv_hazard()

		if(S.hasbhazard) {
			S.GenUbeta = &GenUbeta_hazard_rs()
		}
		else S.GenUbeta = &GenUbeta_hazard()
	}
	else if(S.scale == "odds") {
		S.Stgen = &Stgen_odds()
		S.htgen = &htgen_odds()
		S.RMSTgen = &RMSTgen_odds()
		S.Gen_Sderiv = &Gen_Sderiv_odds()
		if(S.hasbhazard) {
			S.GenUbeta = &GenUbeta_odds_rs()
		}
		else S.GenUbeta = &GenUbeta_odds()
	}
}



//======================================================================================================================================//
// main functions


//======================================================================================================================================//
// standsurv() - main function
void function standsurv()
{
	struct stpm2_standardisation scalar S
	struct stpm2_standardisation_results scalar R
	S = GetStuff()
	R = CreateResultsMat(S)
	DeclarePointers(S)
	if(S.ci) {
		if(S.deltamethod) standsurv_deltamethod(S,R)
		else if(S.mestimation) standsurv_Mestimation(S,R)
	}
	WriteResults(S,R)
}

//======================================================================================================================================//
// standsurv() - Delta Method
void function standsurv_deltamethod(struct stpm2_standardisation scalar S, struct stpm2_standardisation_results scalar R)
{
	if(S.verbose) timer_clear()

	S.ResVcov = asarray_create("real")
	for(j=1;j<=S.Nt;j++) {
		if(S.verbose) display("Time point " + strofreal(j))
		S.tj = S.t[j]
		S.j = j
		lnt = ln(S.tj)
		if(S.verbose) timer_on(1)
		GenSplines(S,lnt)
		if(S.verbose) timer_off(1)
		if(S.verbose) timer_on(2)
		GenX_at(S)
		if(S.verbose) timer_off(2)
		if(S.calchazard) {
			GenDerivSplines(S,lnt)
			Gen_Xdrcs_at(S)
		}
		if(S.verbose) timer_on(3)
		GenUind(S)
		if(S.verbose) timer_off(3)
		if(S.verbose) timer_on(4)
		A12 = (*S.GenA12)(S)
		if(S.verbose) timer_off(4)
		asarray(S.ResVcov,j,(A12*S.V*A12'))
		if(S.verbose) timer_on(5)
		CalcCI_at(S,R,j)
		if(S.hascontrast) CalcContrasts(S,R,j)
		if(S.haslincom) CalcLincom(S,R,j)
		if(S.hasuserfunction) CalcUserfunction(S,R,j)
		if(S.verbose) timer_off(5)
	}

}

//======================================================================================================================================//
// standsurv() - M estimation
void function standsurv_Mestimation(struct stpm2_standardisation scalar S, struct stpm2_standardisation_results scalar R) {
// Get matrices needed for all measures	that do not depend on t
	Ubeta = (*S.GenUbeta)(S)
	(*S.GenA22)(S,A22,zeros)
	S.ResVcov = asarray_create("real")
	if(S.calcsurvival | S.calcrmst | S.calcfailure) (*S.GenA11)(S,A11)	// does not vary over t
	if(S.calccentile) S.Nt = S.Ncentiles
	for(j=1;j<=S.Nt;j++) {
		if(S.calccentile) S.tj = S.centmean[j,]
		else S.tj = S.t[j]
		lnt = ln(S.tj)
		GenSplines(S,lnt)
		if(S.calchazard | S.calccentile) {
			GenDerivSplines(S,lnt)
			Gen_Xdrcs_at(S)
		}
		GenX_at(S)
		GenUind(S)
		if(S.calchazard | S.calccentile) (*S.GenA11)(S,j,A11)
		// A12 and U depend on what is standardized.
		U = (*S.GenU)(S,j)
		A12 = (*S.GenA12)(S)
		// Now put everything together
		VarU = quadvariance((U,Ubeta))
		Ainv = luinv((A11, A12 \ zeros, A22))
		asarray(S.ResVcov,j,((Ainv*VarU*Ainv'):/S.Nobs_model)[1..S.N_at_options,1..S.N_at_options])
		// Calculate CI for at() options & contrasts
		CalcCI_at(S,R,j)
		if(S.hascontrast) CalcContrasts(S,R,j)
		if(S.haslincom) CalcLincom(S,R,j)
		if(S.hasuserfunction) CalcUserfunction(S,R,j)
	}
}

//======================================================================================================================================//
// Get U_beta (derivative of score function) as does not depend on t.
// Same for all options
// Need to add h_star....

function GenUbeta_hazard(struct stpm2_standardisation scalar S)
{
	Ubeta = (S.d :- S.expxb):*S.X
	dxb_select=(S.Nvarlist+1..S.Nparameters-S.hascons)
	Ubeta[,dxb_select] = Ubeta[,dxb_select] :+ (S.d:/S.dxb):*S.Xdrcs
	if(S.hasdel_entry) Ubeta = Ubeta + S.expxb0:*S.X
	return(Ubeta)
}

void function GenUbeta_odds(struct stpm2_standardisation scalar S)
{
	Ubeta = (S.d :- S.expxb):/(1:+S.expxb):*S.X
	dxb_select=(S.Nvarlist+1..S.Nparameters-S.hascons)
	Ubeta[,dxb_select] = Ubeta[,dxb_select] + (S.d:/S.dxb):*S.Xdrcs
	if(S.hasdel_entry) Ubeta = Ubeta + S.expxb0:/(1:+S.expxb0):*S.X
	return(Ubeta)
}

function GenUbeta_hazard_rs(struct stpm2_standardisation scalar S)
{
	Ubeta = (S.expxb:*(S.d:*S.dxb :/ (S.expxb:*S.dxb:+S.obst:*S.bhazard):-1)):*S.X
	dxb_select=(S.Nvarlist+1..S.Nparameters-S.hascons)
	Ubeta[,dxb_select] = Ubeta[,dxb_select] + (S.d:*S.expxb:/(S.expxb:*S.dxb:+S.obst:*S.bhazard)):*S.Xdrcs
	if(S.hasdel_entry) Ubeta = Ubeta + S.expxb0 :*S.X
	return(Ubeta)
}


function GenUbeta_odds_rs(struct stpm2_standardisation scalar S)
{
	Ubeta = -(S.expxb:*(-S.d:*S.dxb:+S.obst:*S.bhazard:+
				S.expxb:*(S.dxb:+S.t:*S.bhazard))):/
				((1:+S.expxb):*(S.t:*S.bhazard :+ S.expxb:*(S.dxb:+ S.obst:*S.bhazard))):*S.X
	dxb_select=(S.Nvarlist+1..S.Nparameters-S.hascons)
	Ubeta[,dxb_select] = Ubeta[,dxb_select] + S.d:*S.expxb:/(S.t:*S.bhazard :+ S.expxb:*(S.dxb :+ S.obst:*S.bhazard)):*S.Xdrcs
	if(S.hasdel_entry) Ubeta = Ubeta + S.expxb0:/(1:+S.expxb0):*S.X
	return(Ubeta)
}


//======================================================================================================================================//
// A11 and A22 do not depend on t for survival
real matrix function GenA11_surv(struct stpm2_standardisation scalar S,A11)
{
	A11 = -1*I(S.N_at_options)
}

real matrix function GenA22_surv(struct stpm2_standardisation scalar S,A22,zeros)
{
	A22 = -luinv(S.V):/ S.Nobs_model
	zeros = J(S.Nparameters,S.N_at_options,0)
}

//======================================================================================================================================//
// A11 varies by t for hazard, but not A22
real matrix function GenA11_hazard(struct stpm2_standardisation scalar S, real scalar j, A11)
{
//	A11 =  diag((S.Smean[j,],-mean(S.touse_at),-mean(S.touse_at)))
	A11 =  diag((S.Smean[j,],J(1,S.N_at_options,-1),J(1,S.N_at_options,-1)))
	index_start = S.N_at_options + 1
	index_stop = index_start + S.N_at_options - 1
	A11[1..S.N_at_options,index_start..index_stop] = (-I(S.N_at_options))
	index_start = 2*S.N_at_options + 1
	index_stop = index_start + S.N_at_options - 1
	A11[1..S.N_at_options,index_start..index_stop] = diag(S.hmean[j,])
}

real matrix function GenA22_hazard(struct stpm2_standardisation scalar S,A22,zeros)
{
	A22 = -luinv(S.V):/ S.Nobs_model
	zeros = J(S.Nparameters,S.N_at_options:*3,0)
}

//======================================================================================================================================//
// A11 varies by centile, but not A22
real matrix function GenA11_centile(struct stpm2_standardisation scalar S, real scalar j, A11)
{
	A11 = -diag(colsum(S.weights:*S.fi:*S.touse_at):/S.Nobs_predict_at)
}

//======================================================================================================================================//
// GenSplines function - generate baseline and tvc splines (at different t values for centile options).
// Within arrays for splines use "0" for survival, hazard and rmst "at"+i for centiles

void function GenSplines(struct stpm2_standardisation scalar S,lnt)
{
	if(!S.rcsbaseoff) S.Xrcsbase = asarray_create()
	if(S.hastvc) S.Xrcstvc = asarray_create("string",2)
	for(i=1;i<=S.loopmax;i++) {
		ati = "at" + strofreal(i)
		if(S.calccentile) asarray_index = ati
		else asarray_index = "0"
		if(!S.rcsbaseoff) {
			if(S.orthog) asarray(S.Xrcsbase,asarray_index,rcsgen_core(lnt[i],asarray(S.knots,"baseline"),0,asarray(S.Rmatrix,"baseline")))
			else asarray(S.Xrcsbase,asarray_index,rcsgen_core(lnt[i],asarray(S.knots,"baseline"),0))
		}
		if(S.hastvc) {
			for(k=1;k<=S.Ntvc;k++) {
				if(S.orthog) asarray(S.Xrcstvc,(asarray_index,S.tvcnames[k]),rcsgen_core(lnt[i],asarray(S.knots,S.tvcnames[k]),0,asarray(S.Rmatrix,S.tvcnames[k])))
				else asarray(S.Xrcstvc,(asarray_index,S.tvcnames[k]),rcsgen_core(lnt[i],asarray(S.knots,S.tvcnames[k]),0))
			}
		}
	}
}

//======================================================================================================================================//
// GenDerivSplines function - loops over at options and calculates derivative of spline functions
// Within arrays for splines use "0" for survival, hazard and rmst "at"+i for centiles
void function GenDerivSplines(struct stpm2_standardisation scalar S,real matrix lnt)
{
	S.Xdrcsbase = asarray_create()
	if(S.hastvc) S.Xdrcstvc = asarray_create("string",2)
	for(i=1;i<=S.loopmax;i++) {
		ati = "at" + strofreal(i)
		if(S.calccentile) asarray_index = ati
		else asarray_index = "0"
		if(S.orthog) asarray(S.Xdrcsbase,asarray_index,rcsgen_core(lnt[i],asarray(S.knots,"baseline"),1,asarray(S.Rmatrix,"baseline")))
		else asarray(S.Xdrcsbase,asarray_index,rcsgen_core(lnt[i],asarray(S.knots,"baseline"),1))
		if(S.hastvc) {
			for(k=1;k<=S.Ntvc;k++) {
				if(S.orthog) asarray(S.Xdrcstvc,(asarray_index,S.tvcnames[k]),rcsgen_core(lnt[i],asarray(S.knots,S.tvcnames[k]),1,asarray(S.Rmatrix,S.tvcnames[k])))
				else asarray(S.Xdrcstvc,(asarray_index,S.tvcnames[k]),rcsgen_core(lnt[i],asarray(S.knots,S.tvcnames[k]),1))
			}
		}

	}
}

//======================================================================================================================================//
// Different X matrices for each at() option
// returns S.X_at_t
void function GenX_at(struct stpm2_standardisation scalar S)
{
	S.X_at_t = asarray_create()
	for(i=1;i<=S.N_at_options;i++) {
		ati = "at" + strofreal(i)
		if(S.calccentile) asarray_index = ati
		else asarray_index = "0"

		X_tmp = J(S.Nobs_model,S.Nparameters,.)
		if(S.hasvarlist) X_tmp[,1..S.Nvarlist] = S.X[,1..S.Nvarlist]
		X_tmp[,asarray(S.X_at_index,ati)] = asarray(S.X_at,ati)

		Xstart = S.Nvarlist + 1
		Xstop = Xstart + S.df - 1

		if(!S.rcsbaseoff) X_tmp[,Xstart..Xstop] = J(S.Nobs_model,1,asarray(S.Xrcsbase,asarray_index))
		if(S.hastvc) {
			for(k=1;k<=S.Ntvc;k++) {
				Xstart = Xstop + 1
				Xstop = Xstart + S.dftvc[k] - 1
				X_tmp[,Xstart..Xstop] = asarray(S.Xrcstvc,(asarray_index,S.tvcnames[k])) # asarray(S.X_at_tvc,ati)[,k]
			}
		}
		if(S.hascons) X_tmp[,S.Nparameters] = J(S.Nobs_model,1,1)
		asarray(S.X_at_t,ati,X_tmp)
	}
}

//======================================================================================================================================//
// Different spline derivative matrix for each at() option
// returns S.Xdrcs_at_t
// could speed up (see function above)

void function Gen_Xdrcs_at(struct stpm2_standardisation scalar S) {
	S.Xdrcs_at_t = asarray_create()
	for(i=1;i<=S.N_at_options;i++) {
		ati = "at" + strofreal(i)
		if(S.calccentile) asarray_index = ati
		else asarray_index = "0"
		X_tmp = J(S.Nobs_model,0,.)
		if(!S.rcsbaseoff) X_tmp = X_tmp,J(S.Nobs_model,1,asarray(S.Xdrcsbase,asarray_index))
		if(S.hastvc) {
			for(k=1;k<=S.Ntvc;k++) {
				X_tmp = X_tmp,J(S.Nobs_model,1,asarray(S.Xdrcstvc,(asarray_index,S.tvcnames[k]))):*asarray(S.X_at_tvc,ati)[,k]
			}
		}
		asarray(S.Xdrcs_at_t,ati,X_tmp)
	}
}


//======================================================================================================================================//
// generates survival and hazard functions for each at option
// uses pointers to call appropraite survial/hazard functions
void function GenUind(struct stpm2_standardisation scalar S)
{
	if(!S.calcrmst) S.Si = J(S.Nobs_model,S.N_at_options,.)
	if(S.calchazard | S.calccentile) {
		S.hi = J(S.Nobs_model,S.N_at_options,.)
		S.fi = J(S.Nobs_model,S.N_at_options,.)
	}
	if(S.calcrmst) {
		S.rmsti = J(S.Nobs_model,S.N_at_options,.)
		S.Sti_nodes = asarray_create() // create as needed for A12
	}
	for(i=1;i<=S.N_at_options;i++) {
		ati = "at" + strofreal(i)
		if(S.calccentile) tt = S.tj[i]
		else tt = S.tj
		if(!S.calcrmst) S.Si[,i] = (*S.Stgen)(S,ati)
		if(S.calchazard | S.calccentile) {
			S.hi[,i] = (*S.htgen)(S,ati,tt)
			S.fi[,i] = S.Si[,i]:*S.hi[,i]
		}
		if(S.calcrmst) {
			S.rmsti[,i] = (*S.RMSTgen)(S,ati)
		}
	}
}

//======================================================================================================================================//
// survival and hazard function for each scale
// need to add other scales
function Stgen_hazard(struct stpm2_standardisation scalar S,string scalar ati)
{
	return(exp(-exp(asarray(S.X_at_t,ati)*S.beta)))
}
function Stgen_odds(struct stpm2_standardisation scalar S,string scalar ati)
{
	return((1 :+ exp(asarray(S.X_at_t,ati)*S.beta)):^(-1))
}

// Note that also creates S.Sti_nodes for use when calculating A12
real matrix function RMSTgen_hazard(struct stpm2_standardisation scalar S,string scalar ati)
{
	nodes_t = (0.5:*S.tj:*S.nodes :+ 0.5:*S.tj)
	xb_tmp = J(S.Nobs_model,S.Nnodes,.)
	S.rcs_nodes = J(S.Nparameters-S.Nvarlist-S.hascons,S.Nnodes,.)
	for(n=1;n<=S.Nnodes;n++) {
		tn = nodes_t[n]
		lntn = ln(tn)
		if(S.hasvarlist) X_tmp = S.X[,1..S.Nvarlist]
		else X_tmp = J(S.Nobs_model,0,.)
		if(S.hasvarlist) X_tmp[,asarray(S.X_at_index,ati)] = asarray(S.X_at,ati)
		if(!S.rcsbaseoff) {
			if(S.orthog) rcstmp = rcsgen_core(lntn,asarray(S.knots,"baseline"),0,asarray(S.Rmatrix,"baseline"))
			else rcstmp = rcsgen_core(lntn,asarray(S.knots,"baseline"),0)
			S.rcs_nodes[1..S.df,n] = rcstmp'
			X_tmp = X_tmp, J(S.Nobs_model,1,rcstmp)
		}
		p_index = S.df + 1
		if(S.hastvc) {
			for(k=1;k<=S.Ntvc;k++) {
				if(S.orthog) rcstmp = rcsgen_core(lntn,asarray(S.knots,S.tvcnames[k]),0,asarray(S.Rmatrix,S.tvcnames[k]))
				else rcstmp = rcsgen_core(lntn,asarray(S.knots,S.tvcnames[k]),0,)
				cols_select = (p_index..p_index + cols(rcstmp) - 1)
				S.rcs_nodes[cols_select,n] = rcstmp'
				X_tmp = X_tmp, J(S.Nobs_model,1,rcstmp):*asarray(S.X_at_tvc,ati)[,k]
				p_index = p_index + S.dftvc[1,k]
			}
		}
		if(S.hascons) X_tmp = X_tmp,J(S.Nobs_model,1,1)
		xb_tmp[,n] = (X_tmp*S.beta)
	}
	asarray(S.Sti_nodes,ati,exp(-exp(xb_tmp)))
	return(0.5:*S.tj:*quadrowsum(S.nodeweights':*asarray(S.Sti_nodes,ati),1))
}

function htgen_hazard(struct stpm2_standardisation scalar S,string scalar ati,tt) {
	return((asarray(S.Xdrcs_at_t,ati)*S.betarcs):/tt :* (exp(asarray(S.X_at_t,ati)*S.beta)))
}
function htgen_odds(struct stpm2_standardisation scalar S,string scalar ati,tt) {
	return((asarray(S.Xdrcs_at_t,ati)*S.betarcs):/tt :* (exp(asarray(S.X_at_t,ati)*S.beta)):*(1:+exp(asarray(S.X_at_t,ati)*S.beta)):^(-1))
}


//======================================================================================================================================//
// A12 matrix
// returns A12
// Survival

// generic functions for scales using pointers
real matrix function GenA12_surv(struct stpm2_standardisation scalar S)
{
	A12 = J(S.N_at_options,S.Nparameters,.)
	for(i=1;i<=S.N_at_options;i++) 	{
		A12[i,] = (*S.Gen_Sderiv)(S,i)
	}
	if(S.calcfailure) _negate(A12)
	return(A12)
}

real matrix function Gen_Sderiv_hazard(struct stpm2_standardisation scalar S,i) {
	ati = "at" + strofreal(i)
	return(colsum(S.weights:*S.Si[,i] :* log(S.Si[,i]):*asarray(S.X_at_t,ati):*S.touse_at[,i]):/S.Nobs_predict_at[,i])
}
real matrix  function Gen_Sderiv_odds(struct stpm2_standardisation scalar S,i) {
	ati = "at" + strofreal(i)
	return(mean(S.weights:*S.Si[,i] :* log(S.Si[,i]):*asarray(S.X_at_t,ati):*S.touse_at[,i]))
}

// Hazard
// odds need updating
// generic functions for scales using pointers

real matrix function GenA12_hazard_mestimation(struct stpm2_standardisation scalar S)
{
	A12 = J(S.N_at_options:*3,S.Nparameters,.)
	A12[1..S.N_at_options,1..S.Nparameters] = J(S.N_at_options,S.Nparameters,0)
	for(i=1;i<=S.N_at_options;i++) {
		ati = "at" + strofreal(i)
		findex = i :+ S.N_at_options
		Sindex = i :+ 2*S.N_at_options
		A12[findex,] = (*S.Gen_fderiv)(S,i)
		A12[Sindex,] = (*S.Gen_Sderiv)(S,i)
	}
	return(A12)
}



real matrix function GenA12_hazard_deltamethod(struct stpm2_standardisation scalar S)
{
	A12 = J(S.N_at_options,S.Nparameters,.)
	for(i=1;i<=S.N_at_options;i++) {
		ati = "at" + strofreal(i)
		A12[i,] =((*S.Gen_fderiv)(S,i):*S.Smean[S.j,i] - (*S.Gen_Sderiv)(S,i):*S.fmean[S.j,i]):/(S.Smean[S.j,i]:^2)
	}
	return(A12)
}

// derivative of f(t) for scale(hazard) models
real matrix function Gen_fderiv_hazard(struct stpm2_standardisation scalar S,i) {
	ati = "at" + strofreal(i)
	dxb_select=(S.Nvarlist+1..S.Nparameters-S.hascons)
	A12f = colsum(S.weights:*S.fi[,i]:*asarray(S.X_at_t,ati):*(log(S.Si[,i]) :+ 1):*S.touse_at[,i]):/S.Nobs_predict_at[,i]
	A12f[,dxb_select] = A12f[,dxb_select] :- colsum(S.weights:*(asarray(S.Xdrcs_at_t,ati):*S.Si[,i] :* log(S.Si[,i]):/S.tj):*S.touse_at[,i]):/S.Nobs_predict_at[,i]
	return(A12f)
}



// Just here for placement
// currently just a copy of hazard scale
function Gen_fderiv_odds(struct stpm2_standardisation scalar S,i) {
/*	ati = "at" + strofreal(i)
	dxb_select=(S.Nvarlist+1..S.Nparameters-S.hascons)
	A12f = mean(S.fi[,i]:*asarray(S.X_at_t,ati):*(1 - 2:*:/S.Si[,i])   :*(log(S.Si[,i]) :+ 1):*S.touse_at[,i])
	A12f[,dxb_select] = A12f[,dxb_select] :- mean((asarray(S.Xdrcs_at_t,ati):*S.Si[,i] :* log(S.Si[,i]):/S.tj):*S.touse_at[,i])
	return(A12f)
*/
}

// RMST
real matrix function GenA12_rmst(struct stpm2_standardisation scalar S)
{
	A12 = J(S.N_at_options,S.Nparameters,.)
	for(i=1;i<=S.N_at_options;i++) {
		A12[i,] = (*S.Gen_RMSTderiv)(S,i)
	}
	return(A12)
}

real matrix function Gen_RMSTderiv_hazard(struct stpm2_standardisation scalar S,real scalar i)
{
	ati = "at" + strofreal(i)
	A12i = J(1,S.Nparameters,.)
	rcs_index = tvc_index = 1
	for(k=1;k<=S.Nparameters;k++) {
		if(k<=S.Nvarlist | (k==S.Nparameters & S.hascons)) {
			A12i[1,k] = 0.5:*S.tj:*quadrowsum(S.nodeweights':*mean(S.weights:*asarray(S.Sti_nodes,ati):*log(asarray(S.Sti_nodes,ati)):*asarray(S.X_at_t,ati)[,k]:*S.touse_at[,i]),1)
		}
		else if (k>S.Nvarlist & k<=(S.Nvarlist +  S.df)){
			A12i[1,k] = 0.5:*S.tj:*quadrowsum(S.nodeweights':*mean(S.weights:*asarray(S.Sti_nodes,ati):*log(asarray(S.Sti_nodes,ati)):*S.rcs_nodes[rcs_index,]:*S.touse_at[,i]),1)
			rcs_index++
		}
		else {
			A12i[1,k] = 0.5:*S.tj:*quadrowsum(S.nodeweights':*mean(S.weights:*asarray(S.Sti_nodes,ati):*log(asarray(S.Sti_nodes,ati)):*S.rcs_nodes[rcs_index,]:*asarray(S.X_at_tvc,ati)[,tvc_index]:*S.touse_at[,i]),1)
			rcs_index++
			if (k>=(S.Nvarlist +  S.df + rowsum(S.dftvc[1..tvc_index]))) tvc_index++
		}
	}
	return(A12i)
}

//======================================================================================================================================//
// U matrix
// Survival
real matrix function GenU_surv(struct stpm2_standardisation scalar S,real scalar j)
{
	return((S.weights:*S.Si:*S.touse_at :- S.Smean[j,]))
}

real matrix function GenU_failure(struct stpm2_standardisation scalar S,real scalar j)
{
	return((S.Smean[j,] :- S.weights:*S.Si:*S.touse_at))
}

// Hazard
real matrix function GenU_hazard(struct stpm2_standardisation scalar S,real scalar j)
{
	return((J(S.Nobs_model,1,S.Smean[j,]:*S.hmean[j,] :- S.hmean[j,]),(S.weights:*S.fi :- S.fmean[j,]):*S.touse_at,(S.weights:*S.Si :- S.Smean[j,]):*S.touse_at))
}
// RMST
real matrix function GenU_rmst(struct stpm2_standardisation scalar S,real scalar j)
{
	return((S.weights:*S.rmsti:*S.touse_at :- S.rmstmean[j,]))
}

// centile
real matrix function GenU_centile(struct stpm2_standardisation scalar S,real scalar j)
{
	return((S.weights:*S.Si:*S.touse_at :- S.centmean[j,]))
}



// This file contins the following functions
//
// mm_root_vec()	- Brent's root finder
// epsilon_vec()	- check if root found to degree of accuracy
// gq()				- guassian quadrature nodes and weights
// rcsgen_core()	- generate restricted cubic splines

//======================================================================================================================================//
// mm_root_vec - adpated form Michael Crowther's adaption of Ben Jann mm_root.
real colvector mm_root_vec(transmorphic x,      // bj: will be replaced by solution
						 pointer(real scalar function) scalar f,
											  // Address of the function whose zero will be sought for
						 real scalar ax,      // Root will be sought for within a range [ax,bx]
						 real scalar bx,      //
						 | real scalar tol,   // Acceptable tolerance for the root value (default 0)
						   real scalar maxit, // bj: maximum # of iterations (default: 1000)
						   o1, o2, o3, o4, o5, o6, o7, o8, o9, o10)            // bj: additional args to pass on to f
{
    transmorphic  fs            // setup for f
    real colvector   a, b, c       // Abscissae, descr. see above
    //real scalar   fa, fb, fc    // f(a), f(b), f(c)
    real scalar   prev_step     // Distance from the last but one
    real scalar   tol_act       // Actual tolerance
    real scalar   p             // Interpolation step is calcu-
    real scalar   q             // lated in the form p/q; divi-
                                // sion operations is delayed
                                // until the last moment
    real scalar   new_step      // Step at this iteration
    real scalar   t1, cb, t2
    real scalar   itr

	real scalar nobs
	real colvector index
	real colvector   fa, fb, fc    // f(a), f(b), f(c)

    if (args()<5) tol = 0       // bj: set tolerance
    if (args()<6) maxit = 1000  // bj: maximum # of iterations

//  fs = mm_callf_setup(f, args()-6, `opts') // bj: prepare function call
    fs = mm_callf_setup(f, args()-6, o1, index) // bj: prepare function call
	nobs = rows(x)
	index = 1::nobs

	result = J(nobs,1,.)
	x = J(nobs,1,.)
    a = J(nobs,1,ax);  b = J(nobs,1,bx);  fa = mm_callf(fs, a);  fb = mm_callf(fs, b)
    c = a;  fc = fa

    //if ( fa==. ) return(0)      // bj: abort if fa missing
	tempindex = selectindex(fa:==.)
	nti = rows(tempindex)
	if (nti & cols(tempindex)) result[tempindex,] = J(nti,1,0)

	//remove tempindex as they are done
	index = select(index,fa:!=.)

	if (rows(index)) { //not done
		tempindex = select(index,((fa[index]:>0) :* (fb[index]:>0)))
		if (cols(tempindex)) {

			flag1 = abs(fa[tempindex]) :< abs(fb[tempindex])
			flag2 = 1:-flag1
			tempindex2 = select(tempindex,flag1)
			nti = rows(tempindex2)
			if (nti) {
				result[tempindex2] = J(nti,1,2)
				x[tempindex2] = a[tempindex2]
			}
			tempindex2 = select(tempindex,flag2)
			nti = rows(tempindex2)
			if (nti) {
				result[tempindex2] = J(nti,1,3)
				x[tempindex2] = b[tempindex2]
			}
			//update index
			index = select(index,x:==.)					//better way to do this?
			if (rows(index)==0) return(result)
		}
		tempindex = select(index,((fa[index]:<0) :* (fb[index]:<0)))
		if (cols(tempindex)) {

			flag1 = abs(fa[tempindex]) :< abs(fb[tempindex])
			flag2 = 1:-flag1
			tempindex2 = select(tempindex,flag1)
			nti = rows(tempindex2)
			if (nti) {
				result[tempindex2] = J(nti,1,2)
				x[tempindex2] = a[tempindex2]
			}
			tempindex2 = select(tempindex,flag2)
			nti = rows(tempindex2)
			if (nti) {
				result[tempindex2] = J(nti,1,3)
				x[tempindex2] = b[tempindex2]
			}
			//update index
			index = select(index,x:==.)					//better way to do this?
			if (rows(index)==0) return(result)
		}
	}
	else return(result)

	for (itr=1; itr<=maxit; itr++) {
		tempindex = index[selectindex(fb[index]:==.)]

		if (cols(tempindex)) result[tempindex] = J(rows(tempindex),1,0)

		//remove tempindex as they are done
		index = select(index,fb[index]:!=.)
		if (!cols(index)) return(result)

		tempindex = select(index,abs(fc[index]) :< abs(fb[index]))

		if (cols(tempindex)) {
			a[tempindex] = b[tempindex];  b[tempindex] = c[tempindex];  c[tempindex] = a[tempindex];         // best approximation
            fa[tempindex] = fb[tempindex];  fb[tempindex] = fc[tempindex];  fc[tempindex] = fa[tempindex]
		}

		tol_act = 2:*epsilon_vec(b[index]) :+ tol:/2
        new_step = (c[index]:-b[index]):/2

		flag1 = (abs(new_step):<=tol_act) :+ (fb[index]:==0)
		flag2 = (flag1:==0)
		tempindex = select(index,flag1)

		if (cols(tempindex)) {
			x[tempindex] = b[tempindex]
			result[tempindex] = J(rows(tempindex),1,0)
		}

		index = select(index,flag2)
		if (!cols(index) | !rows(index)) return(result)

		//update stuff
		tol_act = select(tol_act,flag2)
		new_step = select(new_step,flag2)

        // Decide if the interpolation can be tried
		prev_step = b[index]:-a[index]

		tempindex11 = (abs(prev_step) :>= tol_act) :* (abs(fa[index]) :> abs(fb[index]))
		tempindex = select(index,tempindex11)

		if (cols(tempindex)) {

			cb = c[tempindex] :- b[tempindex]

			p = q  = cb:*0					//fix

			flag1 = a[tempindex] :== c[tempindex]
			flag2 = 1:-flag1
			tempindex2 = select(tempindex,flag1)
			if (cols(tempindex2)) {
				t1 = fb[tempindex2]:/fa[tempindex2]
				p[selectindex(flag1)] = select(cb,flag1) :* t1
				q[selectindex(flag1)] = 1:- t1
			}
			tempindex2 = select(tempindex,flag2)
			if (cols(tempindex2)) {
				q[selectindex(flag2)] = fa[tempindex2]:/fc[tempindex2]; t1 = fb[tempindex2]:/fc[tempindex2]; t2 = fb[tempindex2]:/fa[tempindex2]
				p[selectindex(flag2)] = t2 :* ( select(cb,flag2) :* q[selectindex(flag2)] :* (q[selectindex(flag2)] :- t1) :- (b[tempindex2]:-a[tempindex2]):*(t1:-1) )
                q[selectindex(flag2)] = (q[selectindex(flag2)]:-1) :* (t1:-1) :* (t2:-1)
			}
			flag1 = p:>0
			flag2 = 1:-flag1
			tempindex = selectindex(flag1)
			if (cols(tempindex)) q[tempindex] = -q[tempindex]
			tempindex = selectindex(flag2)
			if (cols(tempindex)) p[tempindex] = -p[tempindex]

			tempindex = (p :< (0.75:*cb:*q:-abs(select(tol_act,tempindex11):*q):/2))  :* (p :< abs(select(prev_step,tempindex11):*q:/2))
			if (cols(tempindex)) {
				//update tempindex11
				tempindex22 = select(selectindex(tempindex11),tempindex)
				if (cols(tempindex22) & rows(tempindex22)) new_step[tempindex22] = p[selectindex(tempindex)]:/q[selectindex(tempindex)]
			}

		}

		tempindex = selectindex(abs(new_step) :< tol_act)

		if (rows(tempindex)) {
			flag1 = new_step[tempindex] :> 0
			flag2 = 1:-flag1
			tempindex2 = select(tempindex,flag1)
			if (rows(tempindex2)) new_step[tempindex2] = tol_act[tempindex2]
			tempindex2 = select(tempindex,flag2)
			if (rows(tempindex2)) new_step[tempindex2] = -tol_act[tempindex2]
        }

        a[index] = b[index];  fa[index] = fb[index]                   // Save the previous approx.
        b[index] = b[index] + new_step

		fb[index] = mm_callf(fs, b[index]) // Do step to a new approxim.

		tempindex1 = select(index,((fb[index]:>0) :* (fc[index]:>0)))
		tempindex2 = select(index,((fb[index]:<0) :* (fc[index]:<0)))

		if (cols(tempindex1)) {
			c[tempindex1] = a[tempindex1]
			fc[tempindex1] = fa[tempindex1]
		}
		if (cols(tempindex2)) {
			c[tempindex2] = a[tempindex2]
			fc[tempindex2] = fa[tempindex2]
		}
    }
	x[index] = b[index]
	result[index] = J(rows(index),1,0)
    return(result)                             // bj: convergence not reached
}


//======================================================================================================================================//
// epsilon_vec -
real colvector epsilon_vec(real colvector y)
{
	res = J(rows(y),1,.)
	for (i=1;i<=rows(y);i++) res[i] = epsilon(y[i])
	return(res)
}


//======================================================================================================================================//
// gq() - gauss/Legenrde nodes and weights

void gq(real matrix weights, real matrix nodes)
{
	n =  strtoreal(st_local("nodes"))

        i = range(1,n,1)'
        i1 = range(1,n-1,1)'
        muzero = 2
        a = J(1,n,0)
        b = i1:/sqrt(4 :* i1:^2 :- 1)
        A= diag(a)
        for(j=1;j<=n-1;j++){
                A[j,j+1] = b[j]
                A[j+1,j] = b[j]
        }
        symeigensystem(A,vec,nodes)
        weights = (vec[1,]:^2:*muzero)'
        weights = weights[order(nodes',1)]
        nodes = nodes'[order(nodes',1)']
}

//======================================================================================================================================//
// rcsgen_core function - calculate splines with provided knots

real matrix rcsgen_core(	real colvector variable,	///
							real rowvector knots, 		///
							real scalar deriv,|			///
							real matrix rmatrix			///
						)
{
	real scalar  Nobs, Nknots, kmin, kmax, interior, Nparams
	real matrix splines, knots2
	//======================================================================================================================================//
	// Extract knot locations

	Nobs 	= rows(variable)
	Nknots 	= cols(knots)
	if(Nknots==0) {
		if (deriv==0) splines = variable
		else splines = J(Nobs,1,0)
		Nparams = 1
	}
	else {
		kmin 	= knots[1,1]
		kmax 	= knots[1,Nknots]

		if (Nknots==2) interior = 0
		else interior = Nknots - 2
		Nparams = interior + 1

		splines = J(Nobs,Nparams,.)

		//======================================================================================================================================//
		// Calculate splines

		if (Nparams>1) {
			lambda = J(Nobs,1,(kmax:-knots[,2..Nparams]):/(kmax:-kmin))
			knots2 = J(Nobs,1,knots[,2..Nparams])
		}

		if (deriv==0) {
			splines[,1] = variable
			if (Nparams>1) {
				splines[,2..Nparams] = (variable:-knots2):^3 :* (variable:>knots2) :- lambda:*((variable:-kmin):^3):*(variable:>kmin) :- (1:-lambda):*((variable:-kmax):^3):*(variable:>kmax)
			}
		}
		else if (deriv==1) {
			splines[,1] = J(Nobs,1,1)
			if (Nparams>1) {
				splines[,2..Nparams] = 3:*(variable:-knots2):^2 :* (variable:>knots2) :- lambda:*(3:*(variable:-kmin):^2):*(variable:>kmin) :- (1:-lambda):*(3:*(variable:-kmax):^2):*(variable:>kmax)
			}
		}
		else if (deriv==2) {
			splines[,1] = J(Nobs,1,0)
			if (Nparams>1) {
				splines[,2..Nparams] = 6:*(variable:-knots2) :* (variable:>knots2) :- lambda:*(6:*(variable:-kmin)):*(variable:>kmin) :- (1:-lambda):*(6:*(variable:-kmax)):*(variable:>kmax)
			}
		}
		else if (deriv==3) {
			splines[,1] = J(Nobs,1,0)
			if (Nparams>1) {
				splines[,2..Nparams] = 6:*(variable:>knots2) :- lambda:*6:*(variable:>kmin) :- (1:-lambda):*6:*(variable:>kmax)
			}
		}
	}

	//orthog
	if (args()==4) {
		real matrix rmat
		rmat = luinv(rmatrix)
		if (deriv==0) splines = (splines,J(Nobs,1,1)) * rmat[,1..Nparams]
		else splines = splines * rmat[1..Nparams,1..Nparams]
	}
	return(splines)
}


//======================================================================================================================================//
// stand_centile - obtain standardized centiles

// Define structure
struct struct_centile_stand
{
	string scalar	atoption,
					wtoption,
					touse_centiles,
					centvar

	real scalar		t_start,		// lower starting time
					t_stop,			// upper starting time
					hasweights,
					Ncentiles		// Number of centiles listed

	real matrix		centiles
}

// Main program
function stand_centile() {
	struct struct_centile_stand scalar C
	C = CentileGetStruct()
	callBrent = mm_root_vec(CentResults=J(C.Ncentiles,1,.),&centile_calc(),C.t_start,C.t_stop,1E-8,1000,C,index)
	(void) st_addvar("double",st_local("newcentvar"))
	st_store(.,st_local("newcentvar"),C.touse_centiles,CentResults)
}

// CentileGetStruct()
// Information need for calculation of centiles
function CentileGetStruct() {
	struct struct_centile_stand scalar C
	C.centvar = st_local("centvar")
	C.centiles = 1 :- (strtoreal(tokens(st_local("centile"))):/100)'
	C.Ncentiles = rows(C.centiles)
	C.atoption = st_local("tempatopt")
	C.touse_centiles = st_local("touse_centiles")
	C.t_start = 0
	C.t_stop = strtoreal(st_local("centileupper"))
	C.hasweights = st_local("indweights") != ""
	if(C.hasweights) C.wtoption = "meansurvwt(" + st_local("indweights") + ")"
	else C.wtoption = ""
	return(C)
}

// function called by mm_root_vec()
function centile_calc(tmpc,struct struct_centile_stand scalar C, index) {

	(void) st_addvar("double","_tmp_t")
	(void) st_addvar("double","_tmp_centiles")

	readrows = (1,rows(tmpc))
	st_store(readrows,"_tmp_t",.,tmpc)
	stata("qui predict double ms_cent, meansurv timevar(_tmp_t) " + C.atoption + C.wtoption)
	st_store(readrows,"_tmp_centiles",.,C.centiles[index])

	stata("qui replace ms_cent = ms_cent - (_tmp_centiles)")
	update_cent = st_data(readrows,"ms_cent",.)
	stata("capture drop _tmp_t ms_cent _tmp_centiles")
	return(update_cent)
}


//======================================================================================================================================//
// rmst - calculate rmst at all timevar points.
// structure
struct struct_rmst_stand  {
	string scalar	atoption,			// current at option
					wtoption,			// external weights option
					touse_time,
					newvarname

	real scalar		Nt,
					Nnodes,
					hasweights


	real matrix		t,
					nodes,
					nodeweights,
					nodes_t,
					rmst
}

// main function
function rmst_stand() {
	struct struct_rmst_stand scalar R
	R = RMSTGetStruct()
	RMST_calcstand(R)
	(void) st_addvar("double",R.newvarname)
	st_store(.,R.newvarname,R.touse_time,R.rmst)
}

function RMSTGetStruct() {
	struct struct_rmst_stand scalar R
	R.newvarname = st_local("rmst_newvarname")
	R.touse_time = st_local("touse_time")
	R.t = st_data(.,st_local("timevar"),R.touse_time)
	R.Nt = rows(R.t)
	R.Nnodes = strtoreal(st_local("nodes"))
	gq(R.nodeweights,R.nodes)
	R.atoption = st_local("tempatopt")
	if(R.hasweights) R.wtoption = "meansurvwt(" + st_local("indweights") + ")"
	else R.wtoption = ""
	return(R)
}

function RMST_calcstand(struct struct_rmst_stand scalar R)
{
	R.rmst = J(R.Nt,1,.)
	for(j=1;j<=R.Nt;j++) {
		tj = R.t[j]
		nodes_t = (0.5:*tj:*R.nodes :+ 0.5:*tj)
		(void) st_addvar("double",tmpt_var = st_tempname())
		ms_var = st_tempname()
		st_store((1,R.Nnodes),tmpt_var,.,nodes_t)
		stata("predict double "+ms_var+", meansurv timevar("+tmpt_var+")" + R.atoption + R.wtoption)
		R.rmst[j] = 0.5:*tj:*quadcolsum(R.nodeweights:*st_data((1,R.Nnodes),ms_var,.))
		stata("capture drop" + invtokens(tmpt_var,ms_var))
	}
}

//======================================================================================================================================//
// Calculate CIs, transformations and write results
//
// CreateResultsMat() - results matrix
// CalcCI_at() - calculate CIS for each at() option
// CalcContrasts() calculate contrast and CIs
// CalcFailure() - 1 - S(t)
// WriteResults() - write results to new variables

//======================================================================================================================================//
// Create Empty Results Matices
function CreateResultsMat(struct stpm2_standardisation scalar S)
{
	struct stpm2_standardisation_results scalar R
	if(S.calccentile) N = S.Ncentiles
	else N=S.Nt

	if(S.ci) R.CI_at_lci = R.CI_at_uci = J(N,S.N_at_options,.)
	if(S.se) R.se = J(N,S.N_at_options,.)
	R.at_est = J(N,S.N_at_options,.)

	// set-up for contrasts
	if(S.hascontrast) {
		if(S.ci) R.CI_contrast_lci = R.CI_contrast_uci = J(N,S.N_at_options,.)
		R.contrast_est = J(N,S.N_at_options,.)
	}
	if(S.haslincom) {
		if(S.ci) R.CI_lincom_lci = R.CI_lincom_uci = J(N,1,.)
		R.lincom_est = J(N,1,.)
	}
	if(S.hasuserfunction) {
		if(S.ci) R.CI_userfunction_lci = R.CI_userfunction_uci = J(N,1,.)
		R.userfunction_est = J(N,1,.)
		stata("mata: pf = &"+st_local("userfunction")+"()")
		external pf
		R.userfunction = pf
	}
	return(R)
}





//======================================================================================================================================//
// Calculate CIs for at() option
real matrix function CalcCI_at(struct stpm2_standardisation scalar S,struct stpm2_standardisation_results scalar R,j)
{
	// transform to scale to calculate CIs
	if(S.calcsurvival) R.at_est[j,] = S.Smean[j,]:*S.per
	else if(S.calcfailure) R.at_est[j,] = (1 :- S.Smean[j,]):*S.per
	else if(S.calchazard) R.at_est[j,] = (S.fmean[j,]:/S.Smean[j,]):*S.per
	else if(S.calccentile) R.at_est[j,] = S.centmean[j,]:*S.per
	else if(S.calcrmst) R.at_est[j,] = S.rmstmean[j,]:*S.per

	S.Vest = S.per:^2:*(asarray(S.ResVcov,j)[1..S.N_at_options,1..S.N_at_options])
	if(S.transform == "none") {
		est_trans = R.at_est[j,]
		dtransform = I(S.N_at_options)
	}
	else if(S.transform == "log") {
		dtransform = diag(1:/R.at_est[j,])
		est_trans = log(R.at_est[j,])
	}
	else if(S.transform == "logit") {
		dtransform = diag(1:/(R.at_est[j,]:*(1:-R.at_est[j,])))
		est_trans = logit(R.at_est[j,])
	}
	else if(S.transform == "loglog") {
		dtransform = diag(1:/(R.at_est[j,]:*ln(R.at_est[j,]))) // Update this
		est_trans = log(-log(R.at_est[j,]))
	}


	Vest_trans = dtransform' * S.Vest * dtransform
	// return CIs of mean survival for each at() option
	theta = invnormal(1-(1-S.level/100)/2)*sqrt(diagonal(Vest_trans))
	if(S.ci) {
		at_lci = est_trans - theta'
		at_uci = est_trans + theta'

		if(S.transform == "none") {
			R.CI_at_lci[j,] = (at_lci)
			R.CI_at_uci[j,] = (at_uci)
		}
		else if(S.transform == "log") {
			R.CI_at_lci[j,] = exp(at_lci)
			R.CI_at_uci[j,] = exp(at_uci)
		}
		else if(S.transform == "logit") {
			R.CI_at_lci[j,] = invlogit(at_lci)
			R.CI_at_uci[j,] = invlogit(at_uci)
		}
		else if(S.transform == "loglog") {
			R.CI_at_lci[j,] = exp(-exp(at_lci))
			R.CI_at_uci[j,] = exp(-exp(at_uci))
		}
	}
	if(S.se) R.se[j,] = sqrt(diagonal(Vest_trans))'
}

//======================================================================================================================================//
// Calculate contrasts and CIs when using contrast() option.
function CalcContrasts(struct stpm2_standardisation scalar S,struct stpm2_standardisation_results scalar R,j)
{
	if(S.contrast == "difference") {
		dcontrast_dtransform = I(S.N_at_options)
		dcontrast_dtransform[,S.at_reference] = J(S.N_at_options,1,-1)
		dcontrast_dtransform[S.at_reference,S.at_reference] = 0
		R.contrast_est[j,] = R.at_est[j,] :- R.at_est[j,S.at_reference]
	}
	if(S.contrast == "ratio") {	// calculate on log scale
		dcontrast_dtransform = I(S.N_at_options):/R.at_est[j,]
		dcontrast_dtransform[,S.at_reference] = J(S.N_at_options,1,-1:/R.at_est[j,S.at_reference])
		dcontrast_dtransform[S.at_reference,S.at_reference] = 0
		R.contrast_est[j,] = ln(R.at_est[j,]) :- ln(R.at_est[j,S.at_reference])
	}
	if(S.contrast == "pchange") {	// calculate on absolute scale
		dcontrast_dtransform = 	-I(S.N_at_options) :/( R.at_est[j,S.at_reference])
		dcontrast_dtransform[,S.at_reference] = (R.at_est[j,])':/(R.at_est[j,S.at_reference]:^2)
		dcontrast_dtransform[S.at_reference,S.at_reference] = 0
		R.contrast_est[j,] = ( R.at_est[j,S.at_reference] :-R.at_est[j,]) :/ (R.at_est[j,S.at_reference])
	}
	Vcont = dcontrast_dtransform*S.Vest*dcontrast_dtransform'

	if(S.contrast == "difference" | S.contrast == "pchange") {
		R.CI_contrast_lci[j,] = R.contrast_est[j,] - invnormal(1-(1-S.level/100)/2)*sqrt(diagonal(Vcont)')
		R.CI_contrast_uci[j,] = R.contrast_est[j,] + invnormal(1-(1-S.level/100)/2)*sqrt(diagonal(Vcont)')
	}
	else if(S.contrast == "ratio") {
		R.CI_contrast_lci[j,] = exp((R.contrast_est[j,]) - invnormal(1-(1-S.level/100)/2)*sqrt(diagonal(Vcont)'))
		R.CI_contrast_uci[j,] = exp((R.contrast_est[j,]) + invnormal(1-(1-S.level/100)/2)*sqrt(diagonal(Vcont)'))
		R.contrast_est[j,] = exp(R.contrast_est[j,])
	}
}

//======================================================================================================================================//
// Calculate contrasts and CIs when using lincom() option.
function CalcLincom(struct stpm2_standardisation scalar S,struct stpm2_standardisation_results scalar R,j)
{
	R.lincom_est[j] = S.lincom*R.at_est[j,]'
	lincom_se = sqrt(S.lincom*S.Vest*S.lincom')
	R.CI_lincom_lci[j] = 	R.lincom_est[j] - invnormal(1-(1-S.level/100)/2)*lincom_se
	R.CI_lincom_uci[j] = 	R.lincom_est[j] + invnormal(1-(1-S.level/100)/2)*lincom_se
}

//======================================================================================================================================//
// Calculate contrasts and CIs when using userfunction option.
function CalcUserfunction(struct stpm2_standardisation scalar S,struct stpm2_standardisation_results scalar R,j)
{
	R.userfunction_est[j] = (*R.userfunction)(R.at_est[j,])

	if(S.ci) {
		Var_user = J(S.N_at_options,1,.)
		eps = 1.0x-1a
		hh = J(1,S.N_at_options,0)
		for(i=1;i<=S.N_at_options;i++) {
			grad = J(1,S.N_at_options,.)
			for(k=1;k<=S.N_at_options;k++) {
				hh[1,k] = eps
				grad[k] = ((*R.userfunction)(R.at_est[j,] + hh) - (*R.userfunction)(R.at_est[j,] - hh)) / (2*eps)
				hh[1,k] = 0
			}
			Var_user[i] = grad*S.Vest*grad'
		}
		R.CI_userfunction_lci[j] = 	R.userfunction_est[j] - invnormal(1-(1-S.level/100)/2)*sqrt(Var_user[1])
		R.CI_userfunction_uci[j] = 	R.userfunction_est[j] + invnormal(1-(1-S.level/100)/2)*sqrt(Var_user[1])
	}
}



//======================================================================================================================================//
// Store results in Stata
// Need to change Smean to make more general

void function WriteResults(struct stpm2_standardisation scalar S,struct stpm2_standardisation_results scalar R)
{
	if(!S.ci) {
		if(S.calcsurvival) R.at_est = S.Smean
		else if(S.calcfailure) R.at_est = 1 :- S.Smean
		else if(S.calchazard) R.at_est = S.fmean:/S.Smean
		else if(S.calccentile) R.at_est = S.centmean
		else if(S.calcrmst) R.at_est = S.rmstmean
	}

	if(S.calccentile) touse_tmp = S.touse_centiles
	else touse_tmp = S.touse_time

	for(i=1;i<=S.N_at_options;i++) {
		(void) st_addvar("double",S.at_varnames[1,i])
		st_store(.,S.at_varnames[1,i],touse_tmp,R.at_est[,i])

		if(S.ci) {
			(void) st_addvar(("double","double"),(S.at_varnames[1,i]+"_lci",S.at_varnames[1,i]+"_uci"))
			st_store(.,S.at_varnames[1,i]+"_lci",touse_tmp,R.CI_at_lci[,i])
			st_store(.,S.at_varnames[1,i]+"_uci",touse_tmp,R.CI_at_uci[,i])
		}
		if(S.se) {
			(void) st_addvar(("double"),S.at_varnames[1,i]+"_se")
			st_store(.,S.at_varnames[1,i]:+"_se",touse_tmp,R.se[,i])
		}
	}
	if(S.hascontrast) {
		cont_index = 1
		for(i=1;i<=S.N_at_options;i++) {
			if(i != S.at_reference) {  // do not write for reference
				(void) st_addvar("double",S.contrast_varnames[1,cont_index])
				st_store(.,S.contrast_varnames[1,cont_index],touse_tmp,R.contrast_est[,i])
				if(S.ci) {
					(void) st_addvar(("double","double"),(S.contrast_varnames[1,cont_index]+"_lci",S.contrast_varnames[1,cont_index]+"_uci"))
					st_store(.,S.contrast_varnames[1,cont_index]+"_lci",touse_tmp,R.CI_contrast_lci[,i])
					st_store(.,S.contrast_varnames[1,cont_index]+"_uci",touse_tmp,R.CI_contrast_uci[,i])
				}
				cont_index = cont_index + 1
			}
		}
	}
	if(S.haslincom) {
		(void) st_addvar("double",S.lincom_varname)
		st_store(.,S.lincom_varname,touse_tmp,R.lincom_est)
		if(S.ci) {
			(void) st_addvar(("double","double"),(S.lincom_varname+"_lci",S.lincom_varname+"_uci"))
			st_store(.,S.lincom_varname+"_lci",touse_tmp,R.CI_lincom_lci)
			st_store(.,S.lincom_varname+"_uci",touse_tmp,R.CI_lincom_uci)
		}
	}

	if(S.hasuserfunction) {
		(void) st_addvar("double",S.userfunction_varname)
		st_store(.,S.userfunction_varname,touse_tmp,R.userfunction_est)
		if(S.ci) {
			(void) st_addvar(("double","double"),(S.userfunction_varname+"_lci",S.userfunction_varname+"_uci"))
			st_store(.,S.userfunction_varname+"_lci",touse_tmp,R.CI_userfunction_lci)
			st_store(.,S.userfunction_varname+"_uci",touse_tmp,R.CI_userfunction_uci)
		}
	}
}

end