*! version 1.0.1 20jan2022 MJC /* History MJC 20jan2022: version 1.0.1 - bug fix; equation names for scalar parameters not updated for v15 - now fixed - version statements added MJC 24nov2021: version 1.0.0 - rebasing, moved to github - minor edits/tidying MJC 31jul2017: version 0.1.2 - af prediction fixed for tvcs - matrix printed out when noorthog was used, now fixed MJC 26jun2017: version 0.1.1 - rcsgen2.ado was missing, xb bug fix in staft_pred.ado MJC 23jun2017: version 0.1.0 Development -> main help file done -> se's for predictions hard coded in Mata using delta method -> note TVC's are included in stpm2 starting values, but zeros are passed to initmat --> only need starting values for baseline splines, but using rcs vars that stpm2 calculates in ml equarions -> was a mistake in starting values when tvc's present -> greatly improved now -> df(1) needs log bknots extracted explicitly, as they're not in e(ln_bknots) from stpm2 when df=1 -> if no covariates and no tvc's are specified it fits stpm2 with a message saying this -> tvc added to a markout call, incase tvc vars not in varlist */ /* To do -> sync with other constraints -> initmat - override initial values but still do stpm2 iter(0) to get splines etc. -> label final spline vars -> check set obs missing problem */ program staft, eclass sortpreserve properties(st) version 15.1 if replay() { syntax [, DF(string) KNOTS(string) *] if "`df'`knots'" != "" { Estimate `0' ereturn local cmdline `"staft `0'"' } else { if "`e(cmd)'" != "staft" { error 301 } if _by() { error 190 } Replay `0' } exit } Estimate `0' ereturn local cmdline `"staft `0'"' end program Estimate, eclass version 15.1 st_is 2 analysis syntax [varlist(default=empty)] [fw pw iw aw] [if] [in] /// /// [, /// DF(passthru) /// -degrees of freedom for the baseline- KNOTS(passthru) /// -knots for the baseline- BKnots(passthru) /// -boundary knots for the baseline- /// BKNOTSTVC(passthru) /// -boundary knots for time-dependent effects- TVC(varlist) /// -varlist of time-dependent effects- DFTVC(passthru) /// -degrees of freedom for time-dependent effects- KNOTSTVC(passthru) /// -knots for time-dependent effects- NOORTHog /// -do not orthogonalise splines- /// /* Results display options */ /// EFORM /// -Exponentiate first ml equation- SHOWINITial /// -Show output from fitting initial value models- Level(cilevel) /// -Statistical significance level- /// /* Maximisation options */ /// * /// -ML options- /// ADAPT /// -undocumented- LININIT /// MLMETHOD(string) /// ] //===================================================================================================================================================// // Error checks and preliminaries // // Check stpm2 and rcsgen are installed capture which stpm2 if _rc >0 { display in yellow "You need to install the command stpm2. This can be installed using," display in yellow ". {stata ssc install stpm2}" exit 198 } capture which rcsgen if _rc >0 { display in yellow "You need to install the command rcsgen. This can be installed using," display in yellow ". {stata ssc install rcsgen}" exit 198 } // Weights //!!not done if "`weight'" != "" { display as err "weights must be stset" exit 101 } local wt: char _dta[st_w] local wtvar: char _dta[st_wv] if "`wt'" != "" { local fw fw(`wtvar') } // Marksample and mlopts marksample touse markout `touse' `tvc' qui replace `touse' = 0 if _st==0 | `touse'==. qui count if `touse' local N `r(N)' if `r(N)' == 0 { display in red "No observations" exit 2000 } qui count if `touse' & _d if `r(N)' == 0 { display in red "No failures" exit 198 } if "`mlmethod'"=="" { local mlmethod lf2 } if "`showinitial'"!="" { local noisily noisily } // Check time origin for delayed entry models local delentry = 0 qui summ _t0 if `touse', meanonly if r(max)>0 { display in green "Note: delayed entry models are being fitted" di "" local delentry = 1 tempvar t0ind index qui egen `index' = seq() if `touse'==1 qui gen byte `t0ind' = (_t0>0 & `touse'==1) qui replace `t0ind'=0 if `t0ind'==. } // temp timevars -> used in rcsgen calls from lnl evaluator tempvar lntxb qui gen double `lntxb' = log(_t) if `touse' if `delentry' { tempvar lnt0xb qui gen double `lnt0xb' = log(_t0) if `t0ind' } //if no covariates or tvc's then exit after stpm2 model local stpm2 = 0 if "`varlist'"=="" & "`tvc'"=="" { di in green "No covariates have been specified -> fitting stpm2 model" local noisily noisily //makes sure stpm2 is shown local stpm2 = 1 } //=======================================================================================================================================================// // Starting values if !`stpm2' { di as txt "Obtaining initial values:" local ncovs : word count `varlist' //starting values for covariate effects tempname initmat1 initmat2 qui `noisily' streg `varlist', dist(weib) time level(`level') mat `initmat1' = e(b) mat `initmat1' = `initmat1'[1,1..`ncovs'] //starting values for splines, and creates splines qui `noisily' stpm2, scale(h) `df' `knots' `noorthog' level(`level') `bknots' `bknotstvc' tvc(`tvc') `dftvc' `knotstvc' `knscale' failconvlininit local Nsplines : word count `e(rcsterms_base)' local ln_bknots `e(ln_bhknots)' //all log baseline knots including boundary knots if "`ln_bknots'"=="" { //this is empty when df(1) local ln_bknots `=log(`: word 1 of `e(boundary_knots)'')' local ln_bknots `ln_bknots' `=log(`: word 2 of `e(boundary_knots)'')' } if "`noorthog'"=="" { tempname rmat mat `rmat' = e(R_bh) local rmatopt rmatrix(`rmat') } local tvc `e(tvc)' if "`tvc'"!="" { foreach tvcvar in `tvc' { local rcsterms_`tvcvar' `e(rcsterms_`tvcvar')' local tvcvars `tvcvars' `rcsterms_`tvcvar'' local drcsterms_`tvcvar' `e(drcsterms_`tvcvar')' local dtvcvars `dtvcvars' `drcsterms_`tvcvar'' local ln_tvcknots_`tvcvar' `e(ln_tvcknots_`tvcvar')' local boundary_knots_`tvcvar' `e(boundary_knots_`tvcvar')' if "`nooorthog'"=="" { tempname R_`tvcvar' mat `R_`tvcvar'' = e(R_`tvcvar') } local ind = 1 foreach tvar in `e(rcsterms_`tvcvar')' { local dtvar : word `ind' of `e(drcsterms_`tvcvar')' constraint free constraint `r(free)' [xb][`tvar'] = [dxb][`dtvar'] local conslist `conslist' `r(free)' local `++ind' } } local Ntvcvars : word count `tvcvars' mat `initmat1' = `initmat1',J(1,`=2*`Ntvcvars'',0) if `delentry' { foreach tvcvar in `tvc' { local s0rcsterms_`tvcvar' : subinstr local rcsterms_`tvcvar' "_rcs" "_s0_rcs", all local s0tvcvars `s0tvcvars' `s0rcsterms_`tvcvar'' local ind = 1 foreach tvar in `e(rcsterms_`tvcvar')' { local s0tvar : word `ind' of `s0rcsterms_`tvcvar'' constraint free constraint `r(free)' [xb][`tvar'] = [s0xb][`s0tvar'] local conslist `conslist' `r(free)' local `++ind' } } foreach nontvcvar in `varlist' { constraint free constraint `r(free)' [xb][`nontvcvar'] = [s0xb][`nontvcvar'] local conslist `conslist' `r(free)' } //initmat needs non tvc covariate starting values and tvc zeros mat `initmat1' = `initmat1',`initmat1'[1,1..`ncovs'],J(1,`Ntvcvars',0) } local constraints constraints(`conslist') } //starting values for splines (skipping any tvc's to get intercept) tempname initmat2 mat `initmat2' = e(b) mat `initmat2' = `initmat2'[1,"xb:"] mat `initmat1' = `initmat1',`initmat2'[1,1..`Nsplines'],`initmat2'[1,colsof(`initmat2')] } else { //stpm2 and exit qui `noisily' stpm2, scale(h) `df' `knots' `noorthog' level(`level') `bknots' `bknotstvc' tvc(`tvc') `dftvc' `knotstvc' `knscale' exit } di "" //=======================================================================================================================================================// // ML equations local mleqns (xb: `varlist' `tvcvars', nocons) if "`tvc'"!="" { local mleqns `mleqns' (dxb: `dtvcvars',nocons) if `delentry' { local mleqns `mleqns' (s0xb: `varlist' `s0tvcvars', nocons) } } forvalues i=1/`Nsplines' { local mleqns `mleqns' /rcs`i' } local mleqns `mleqns' /cons //=======================================================================================================================================================// // Maximisation if "`adapt'"!="" { local dprolog derivprolog(staft_prolog()) } local initopt init(`initmat1',copy) if "`lininit'"!="" { quietly staft `varlist', df(1) `adapt' mat `initmat1' = e(b) local Nvars : word count `varlist' mat `initmat1' = `initmat1'[1,1..`=`Nvars'+1'],J(1,`=`Nsplines'-1',0),`initmat1'[1,colsof(`initmat1')] local initopt init(`initmat1',copy) } mata: staft_setup() di as txt "Fitting full model:" ml model `mlmethod' staft_lf2() `mleqns' /// if `touse' /// `wt', /// init(`initmat1',copy) /// `options' /// waldtest(0) /// search(off) /// userinfo(`staft_struct') /// `searchopt' /// `dprolog' /// `collinopt' /// `constraints' /// `nolog' /// maximize constraint drop `conslist' //Tidy up and create final _rcs splines tempvar lntxbfinal qui predictnl double `lntxbfinal' = log(_t * exp(-xb(xb))) if `touse' forvalues i=1/`Nsplines' { cap drop _rcs`i' cap drop _d_rcs`i' } qui rcsgen `lntxbfinal' if `touse', gen(_rcs) knots(`ln_bknots') dgen(_d_rcs) `rmatopt' ereturn local rcsterms `r(rcslist)' ereturn local drcsterms `r(drcslist)' if `delentry' & "`tvc'"!="" { tempvar lnt0xbfinal qui predictnl double `lnt0xbfinal' = log(_t0 * exp(-xb(s0xb))) if `touse' forvalues i=1/`Nsplines' { cap drop _s0_rcs`i' cap drop _d_s0_rcs`i' } qui rcsgen `lnt0xbfinal' if `touse', gen(_s0_rcs) knots(`ln_bknots') dgen(_d_s0_rcs) `rmatopt' ereturn local s0rcsterms `r(rcslist)' ereturn local ds0rcsterms `r(drcslist)' } ereturn local predict staft_pred ereturn local title "Restricted cubic spline accelerated failure time model" ereturn local cmd staft ereturn local varlist `varlist' ereturn local tvc `tvc' ereturn local ln_bknots `ln_bknots' ereturn local noorthog `noorthog' if "`noorthog'"=="" { ereturn matrix R = `rmat' } ereturn local delentry = `delentry' if "`tvc'"!="" { foreach tvcvar in `tvc' { ereturn local rcsterms_`tvcvar' `rcsterms_`tvcvar'' ereturn local drcsterms_`tvcvar' `drcsterms_`tvcvar'' ereturn local ln_tvcknots_`tvcvar' `ln_tvcknots_`tvcvar'' ereturn local boundary_knots_`tvcvar' `boundary_knots_`tvcvar'' if "`noorthog'"=="" { ereturn matrix R_`tvcvar' = `R_`tvcvar'' } } } ereturn scalar dev = -2*e(ll) ereturn scalar AIC = -2*e(ll) + 2 * e(rank) qui count if `touse' == 1 & _d == 1 ereturn scalar BIC = -2*e(ll) + ln(r(N)) * e(rank) Replay, level(`level') `showcons' `variance' `eform' end program Replay version 15.1 syntax [, Level(cilevel) EFORM] ml display, level(`level') `eform' end