*! version 1.7.6 18Jan2023 /* History PL 12mar2019: meansurv with no timevar opton and td effects. PL 18jul2018: fixed bug with rcsbaseoff and meansurv PL 16jun2018: allow at(x1=x2) for example PL 10oct2017: Fix for using option for lifelost PL 04oct2017: Previous fix broke use of if statement PL 17sep2017: Fixed so meansurv now works with Stata version < 15 again. PL 10sep2017: added meanhazard option PL 11aug2017: corrected bugs in stpm2_centpred (did nt work with df(1) or with no covariates) PL 24may2017: fixed bug when standardising with cure models PL 05mar2017: tidied lifelost code PR 14jul2016: Fixed bug in rmst with tmin>0. Disallowed rsdst with tmin>0 (calc of rsdst not worked out). PL 20jun2016: Update centile option - now calls stpm2_centpred and used Brents algorithm. PR 21dec2015: Fixed bug in tmax() option, also deleted an erroneous blank line. PR 09OCT2015: option rmst, meansurv recoded to allow tmax() to contain a numlist. PR 02SEP2015: Corrected bug in rmst option - could fail if tmax() not specified. PL 08APR2015: Fixed covpat as use of expr_query not compatible with Stata 12 PR 12NOV2014: Corrected minor bug with rmst option PR 30JUL2014: Allowed meansurv option with rmst for covariate-adjusted RMST calc with meansurv. PL 29JUL2014: fixed bug in meansurv option that did not allow the use of factor variables PL 03JAN2014: Added meansurvwt() option to allow external standardization when using meansurv. PL 09DEC2013: Added if2() option to increase possibilities when using predictnl PL 21NOV2013: fixed another bug when using factor variables with interactions. PL 28AUG2013: fixed bug when using hrnum etc with factor variables with continuous interactions. PL 20JAN2013: added predictions for loss in expectation of life for relative survival models. PL 20JAN2013: added predictions for log link. PL 19JUN2012: fixed bug with centile option when df() and dftvc() had different values. PL 12Feb2012: added failure option (1-S(t)). PR 09sep2011: rmst now works with factor variables. PL 08sep2011: fixed bug in centiles for uncured in cure models / PR rmst update PL 15aug2011: if timevar variable = 0 then survival is predicted to be 1 rather than missing. PR 26may2011: undocumented option nliter() to set option iterate() in predictnl for rmst. If not set, predictnl's default iterate(100) is used. For simulations, it may be good to take nliter(0) or nliter(some small integer) to avoid long run times. PR 09apr2011: hidden option power() of _rmst subroutine not implemented as option for predict PR 25feb2011: add rmst and abc options for mean survival time and area between curves PR 25feb2011: add stdp option to various predict options PL 29NOV2010: fixed bug when using meansurv option with a null model. PL 07Nov2010: corrected bug that caused error if tempory folder had spaces. PL 01Nov2010: introduced factor variables. PL 02sep2010: added Therese Anderssons additions for prediction after fitting a cure model PL 15mar2010: added reverse option to spline generation PL 01feb2010: change to Stata 11. PL 12nov2009: added tvc option PL 21sep2009: ensure predictions work with new rcsbaselineoff option in stpm2. PL 08Sep2009: modification of msurvpop to work in Stata 11. PL 12Mar2009: changed to using Newton-Raphson method for estimating centiles. PL 11Mar2009: fixed problem with strings>244 characters for long varlists. PL 11dec2008: changed to using rcsgen for spline functions. PL 20aug2008: Added hdiff1 and hdiff2 options for difference in survival curves. PL 20aug2008: Added hdiff1 and hdiff2 options for difference in survival curves. PL 09aug2008: Fixed bug for some predictions involving time-dependent effects. */ program stpm2_pred, sortpreserve version 11.1 syntax newvarname [if] [in], [Survival Hazard XB XBNOBaseline DXB DZDY HRNumerator(string) HRDenominator(string) MEANSurv /// CENtile(string) CUMHazard CUMOdds NORmal MARTingale DEViance DENSity AT(string) ATVAR(string) /// ZEROs FAILure /// noOFFset SDIFF1(string) SDIFF2(string) HDIFF1(string) HDIFF2(string) tvc(varname) /// CI LEVel(real `c(level)') TIMEvar(varname) STDP PER(real 1) /// CENTOL(real 0.0001) Cure UNCured STARTUNC(real -1) /// RMst RSDst TMAx(string) TMIn(string) n(int 1000) POWer(real 1) ABc hr0(string) first(int 1) /// NLiter(string) CENTITER(int 100) /// MERGEBY(string) DIAGAge(name) DIAGYear(name) BY(string) MAXYear(int 2050) /// ATTAge(name) ATTYear(name) ID(name) MAXAge(int 99) NOBS(int 100) MAXT(real 10) /// SURVProb(name) LIFelost OBSSurv GRpd STub(string) EXPSurv(name) /// NODEs(int 50) TINF(int 50) TCOND(int 0) /// IF2(string) noCOLLAPSEMEANSURV *] local newvarname `varlist' if "`if'" != "" & "`if2'" != "" { di in green "WARNING: if statement has been replaced by expression in if2() option" } if "`if2'" != "" { local if `if2' } marksample touse, novarlist local 0 `", `options'"' syntax [,MEANSURVWT(passthru) MEANHAZard MEANFT USIng(string)] qui count if `touse' if r(N)==0 { error 2000 /* no observations */ } /* Check Options */ /* First check rmst options */ if "`tmax'" != "" { if "`rmst'`rsdst'`abc'" == "" { display as error "tmax() valid only with rmst, rsdst or abc" exit 198 } if "`rsdst'" != "" local rmst rmst } else local tmax . // ! PR: allow meansurv with rmst to permit rmst adjusted for covariates if "`meansurv'"!="" & "`rmst'"!="" { local meansurvrmst meansurv local meansurv } // ! end PR insert if "`tmin'" != "" { if "`rmst'`rsdst'`abc'" == "" { display as error "tmin() valid only with rmst, rsdst or abc" exit 198 } confirm number `tmin' if `tmin' < 0 { display as error "tmin() may not be negative" exit 198 } if `tmin' > 0 & "`rsdst'" != "" { display as error "rsdst not supported for tmin() > 0" exit 198 } } else local tmin 0 if "`hr0'" != "" { if "`rmst'`abc'" == "" { display as error "hr0() valid only with rmst or abc" exit 198 } cap confirm var `hr0' if c(rc) { cap confirm number `hr0' if c(rc) { display as error "hr0() must be a numeric constant or variable" exit 198 } if `hr0' <= 0 { display as error "hr0() must exceed 0" exit 198 } } } if "`hrdenominator'" != "" & "`hrnumerator'" == "" { display as error "You must specifiy the hrnumerator option if you specifiy the hrdenominator option" exit 198 } if "`sdiff2'" != "" & "`sdiff1'" == "" { display as error "You must specifiy the sdiff1 option if you specifiy the sdiff2 option" exit 198 } if "`hdiff2'" != "" & "`hdiff1'" == "" { display as error "You must specifiy the hdiff1 option if you specifiy the hdiff2 option" exit 198 } local hratiotmp = substr("`hrnumerator'",1,1) local sdifftmp = substr("`sdiff1'",1,1) local hdifftmp = substr("`hdiff1'",1,1) if wordcount(`"`survival' `hazard' `failure' `meansurv' `meanhazard' `meanft' `hratiotmp' `sdifftmp' `hdifftmp' `centile' `xb' `xbnobaseline' `dxb' `dzdy' `martingale' `deviance' `cumhazard' `cumodds' `normal' `density' `tvc' `cure' `rmst' `obssurv' `lifelost'"') > 1 { display as error "You have specified more than one option for predict" exit 198 } if wordcount(`"`survival' `hazard' `failure' `meansurv' `meanhazard' `meanft' `hrnumerator' `sdiff1' `hdifftmp' `centile' `xb' `xbnobaseline' `dxb' `dzdy' `martingale' `deviance' `cumhazard' `cumodds' `normal' `density' `tvc' `cure' `rmst' `abc' `obssurv' `lifelost'"') == 0 { display as error "You must specify one of the predict options" exit 198 } if `per' != 1 & "`hazard'" == "" & "`hdiff1'" == "" & "`meanhazard'" == ""{ display as error "You can only use the per() option in combination with the hazard, hdiff1()/hdiff2() or meanhazard options." exit 198 } if "`stdp'" != "" & "`ci'" != "" { display as error "You can not specify both the ci and stdp options." exit 19 } if "`stdp'" != "" & /// wordcount(`"`xb' `dxb' `xbnobaseline' `rmst' `abc' `hrnumerator' `hdiff1' `sdiff1'"') == 0 { display as error "The stdp option cannot be used with this predict option." exit 198 } if "`ci'" != "" & /// wordcount(`"`survival' `hazard' `failure' `hrnumerator' `sdiff1' `hdiff1' `centile' `xb' `dxb' `xbnobaseline' `tvc' `meansurv' `meanhazard' `cure' `rmst' `abc' `obssurv' `lifelost'"') == 0 { display as error "The ci option can not be used with this predict option." exit 198 } if "`zeros'" != "" & ("`meansurv'" != "" | "`meanhazard'" != "" | "`meanft'" != "") { display as error "You can not specify the zero option with the `meansurv'`meanhazard'`meanft' option." exit 198 } if "`zeros'" != "" & "`tvc'" != "" { display as error "You can not specify the zero option with the tvc option." exit 198 } if "`zeros'" != "" & ("`hrnumerator'" != "" | "`hdiff1'" != "" | "`sdiff1'" != "") { display as error "You can not specify the zero option with the hrnumerator, hdiff or sdiff options." exit 198 } if "`zeros'" != "" & ("`obssurv'" != "" | "`lifelost'" != "") { display as error "You can not specify the zero option with the obssurv or lifelost options." exit 198 } if ("`at'" != "" | "`atvar'" != "") & "`hrnumerator'" != "" { display as error "You can not use the at() or atvar() options with the hrnumerator option" exit 198 } if ("`at'" != "" | "`atvar'" != "") & "`sdiff1'" != "" { display as error "You can not use the at() option with the sdiff1 and sdiff2 options" exit 198 } if ("`at'" != "" | "`atvar'" != "") & "`hdiff1'" != "" { display as error "You can not use the at() option with the hdiff1 and hdiff2 options" exit 198 } if "`at'" != "" & ("`obssurv'" != "" | "`lifelost'" != "") { display as error "You can not specify the at option with the obssurv or lifelost options." exit 198 } if "`timevar'" != "" & ("`obssurv'" != "" | "`lifelost'" != "") { display as error "You can not specify the timevar option with the obssurv or lifelost options." exit 198 } if ("`uncured'" != "" | "`cure'" != "") & "`e(cure)'" == "" { display as error "You can only use the cure and uncured options after fitting a cure model" exit 198 } if "`uncured'" != "" & "`survival'" == "" & "`hazard'" == "" & "`centile'" == "" & "`meansurv'" == "" { display as error "You must specify the survival, hazard or centile option if you specify uncured" exit 198 } if ("`meansurv'" == "" & "`meanhazard'" == "" & "`meanft'" == "") & "`meansurvwt'" != "" { display as error "You must specify the meansurv() option when using meansurvwt()." exit 198 } /* call Stmeancurve if meansurv option specified */ if "`meansurv'" != "" | "`meanhazard'" !="" | "`meanft'" != "" { Stmeancurve `newvarname' if `touse', timevar(`timevar') at(`at') atvar(`atvar') `collapsemeansurv' `meansurvwt' `meanhazard' `meanft' `ci' `offset' `uncured' level(`level') per(`per') exit } /* options for obssurv and lifelost */ if "`obssurv'" != "" | "`lifelost'" != "" { if "`diagage'" == "" { local diagage agediag } capture confirm variable `diagage' if _rc>0 { display as error "Age at diagnosis variable, `diagage', does not exist. Specify using the diagage() option." exit 198 } if "`diagyear'" == "" { local diagyear yeardiag } capture confirm variable `diagyear' if _rc>0 { display as error "Year at diagnosis variable, `diagyear', does not exist. Specify using the diagyear() option." exit 198 } if "`mergeby'" == "" { display "You need to specify the mergeby() when using the lifelost/obssurv options to merge in the population mortality file." exit 198 } if "`using'" == "" { display "You need to specify the using() option to give the population mortality filename." exit 198 } capture qui desc using `"`using'"', varlist local using_varlist `r(varlist)' if _rc == 601 { display "File specified in using() option does not exist." exit 198 } foreach mvar in `mergeby' { if `"`: list posof `"`mvar'"' in using_varlist'"' == "0" { display as error "Variable `mvar' specified in mergeby() option does not exist in using dataset." exit 198 } } capture confirm variable _merge if _rc==0 { display as error "Variable _merge exists and needs to be dropped to use the lifelost/obssurv options." exit 198 } if "`attage'" == "" { local attage _age } novarabbrev capture confirm variable `attage' if _rc == 0 { display as error "Attained age variable `attage' already exists in dataset. Delete or rename. See the attage() option." exit 198 } if "`attyear'" == "" { local attyear _year } novarabbrev capture confirm variable `attyear' if _rc == 0 { display as error "Attained year variable `attyear' already exists in dataset. Delete or rename. See the attyear() option." exit 198 } } /* call Stobssurv if obssurv option specified */ if "`obssurv'" != "" { Stobssurv `newvarname' if `touse', timevar(`timevar') `offset' mergeby(`mergeby') diagage(`diagage') diagyear(`diagyear') /// by(`by') attage(`attage') attyear(`attyear') id(`id') maxage(`maxage') maxyear(`maxyear') `ci' level(`level') /// survprob(`survprob') using(`"`using'"') `grpd' nobs(`nobs') maxt(`maxt') expsurv(`expsurv') exit } /* call Stlifelost if lifelost option specified */ if "`lifelost'" != "" { Stlifelost `newvarname' if `touse', `offset' mergeby(`mergeby') diagage(`diagage') diagyear(`diagyear') /// by(`by') attage(`attage') attyear(`attyear') maxage(`maxage') maxyear(`maxyear') `ci' level(`level') /// survprob(`survprob') using(`"`using'"') `grpd' stub(`stub') tinf(`tinf') nodes(`nodes') tcond(`tcond') exit } /* call _rmst (formerly stpm2_mttf) for SE and CI if rmst option specified */ /* call _rmst (formerly stpm2_mttf) for SE and CI if rmst option specified */ if ("`rmst'" != "") & ("`meansurvrmst'" != "") { // rmst meansurv allows >1 values of tmax // !! ?? add numlist parse of tmax() - done anyway by _rmst_meansurv if "`at'" != "" local At at(`at') cap confirm var `newvarname' if c(rc) != 0 { _rmst_meansurv if `touse', tmin(`tmin') tmax(`tmax') generate(`newvarname') /// `At' `zeros' n(`n') `rsdst' if "`ci'`stdp'" == "" { exit } } cap drop `newvarname' if "`ci'" != "" { local predictnl_opts , ci(`newvarname'_lci `newvarname'_uci) } else if "`stdp'" != "" { local predictnl_opts , se(`newvarname'_se) } if "`nliter'"!="" local nli iterate(`nliter') else local nli qui predictnl `newvarname' = predict(rmst tmin(`tmin') tmax(`tmax') /// `At' `zeros' n(`n') `rsdst' meansurv) `nli' if `touse' /// `predictnl_opts' /* force */ exit } if "`rmst'" != "" { if !missing(`tmax') { confirm number `tmax' if `tmax' <= 0 { display as error "tmax() must be positive" exit 198 } } // else local tmax if "`at'" != "" local At at(`at') if ("`meansurvrmst'" != "") local meansurv meansurv cap confirm var `newvarname' if c(rc) != 0 { _rmst if `touse', tmin(`tmin') tmax(`tmax') generate(`newvarname') /// `At' `zeros' n(`n') power(`power') `rsdst' `meansurv' if "`ci'`stdp'" == "" { exit } } cap drop `newvarname' if "`ci'" != "" { local predictnl_opts , ci(`newvarname'_lci `newvarname'_uci) } else if "`stdp'" != "" { local predictnl_opts , se(`newvarname'_se) } if "`nliter'"!="" local nli iterate(`nliter') else local nli if !missing(`tmax') local Tmax tmax(`tmax') qui predictnl `newvarname' = predict(rmst tmin(`tmin') `Tmax' /// `At' `zeros' n(`n') power(`power') `rsdst' `meansurv') `nli' if `touse' /// `predictnl_opts' force exit } if "`abc'" != "" { // if "`at'" != "" local At at(`at') if "`hrdenominator'" != "" local hrd hrdenominator(`hrdenominator') if `first' local wt wt else local wt if "`tmin'" != "" local tmin tmin(`tmin') if "`hr0'" != "" local hr0 hr0(`hr0') cap confirm var `newvarname', exact if c(rc) != 0 { // `newvarname' does not exist qui _abc `newvarname' if `touse', `hr0' hrnumerator(`hrnumerator') `hrd' /// tmax(`tmax') `tmin' `At' `zeros' n(`n') `wt' cap drop _weight if "`ci'`stdp'" == "" { *cap erase _predict_weights.dta exit } } cap drop `newvarname' if "`ci'" != "" { local predictnl_opts ci(`newvarname'_lci `newvarname'_uci) } else if "`stdp'" != "" { local predictnl_opts se(`newvarname'_se) } cap predictnl double `newvarname' = predict(abc `hr0' hrnumerator(`hrnumerator') `hrd' /// tmax(`tmax') `tmin' `At' `zeros' n(`n') first(0)) if `touse', `predictnl_opts' level(`level') local rc = c(rc) // Tidy up weight variable and temporary file created by _abc cap drop _weight cap erase _predict_weights.dta if `rc' { noi di as err "could not estimate SE or CI - predictnl failed to converge" exit 498 } else exit } /* calculate midt for centile option */ summ _t, meanonly local midt = (r(max) - r(min))/2 /* calculate startt for centile option if uncured option is specified*/ /* if `startunc' == -1 { summ _t, meanonly local startt = (r(max) - r(min))/8 } else { local startt = `startunc' } */ /* store time-dependent covariates and main varlist */ local etvc `e(tvc)' local main_varlist `e(varlist)' /* dydx option of old version of stpm */ if "`dzdy'" != "" { local dxb dxb } /* generate ocons for use when orthogonalising splines */ tempvar ocons gen `ocons' = 1 /* Use _t if option timevar not specified */ tempvar t lnt if "`timevar'" == "" { qui gen double `t' = _t if `touse' qui gen double `lnt' = ln(_t) if `touse' } else { qui gen double `t' = `timevar' if `touse' qui gen double `lnt' = ln(`timevar') if `touse' } /* Check to see if nonconstant option used */ if "`e(noconstant)'" == "" { tempvar cons qui gen `cons' = 1 if `touse' } /* Preserve data for out of sample prediction */ tempfile newvars preserve /* Calculate new spline terms if timevar option specified */ if "`timevar'" != "" & "`e(rcsbaseoff)'" == "" { capture drop _rcs* _d_rcs* if "`e(orthog)'" != "" { tempname rmatrix matrix `rmatrix' = e(R_bh) local rmatrixopt rmatrix(`rmatrix') } qui rcsgen `lnt' if `touse', knots(`e(ln_bhknots)') gen(_rcs) dgen(_d_rcs) `e(reverse)' `rmatrixopt' `e(nosecondder)' `e(nofirstder)' } /* calculate new spline terms if timevar option or hrnumerator option is specified */ if "`timevar'" != "" | "`hrnumerator'" != "" | "`sdiff1'" != "" | "`hdiff1'" != "" { foreach tvcvar in `e(tvc)' { if (("`hrnumerator'" != "" | "`sdiff1'" != "" | "`hdiff1'" != "") & "`timevar'" == "") | "`e(rcsbaseoff)'" != "" { capture drop _rcs_`tvcvar'* _d_rcs_`tvcvar'* } if "`e(orthog)'" != "" { tempname rmatrix_`tvcvar' matrix `rmatrix_`tvcvar'' = e(R_`tvcvar') local rmatrixopt rmatrix(`rmatrix_`tvcvar'') } qui rcsgen `lnt' if `touse', gen(_rcs_`tvcvar') knots(`e(ln_tvcknots_`tvcvar')') dgen(_d_rcs_`tvcvar') `e(reverse)' `rmatrixopt' `e(nosecondder)' `e(nofirstder)' if "`hrnumerator'" == "" & "`sdiff1'" == "" & "`hdiff1'" == "" { forvalues i = 1/`e(df_`tvcvar')'{ qui replace _rcs_`tvcvar'`i' = _rcs_`tvcvar'`i'*`tvcvar' if `touse' qui replace _d_rcs_`tvcvar'`i' = _d_rcs_`tvcvar'`i'*`tvcvar' if `touse' } } } } /* zeros */ if "`zeros'" != "" { local tmptvc `e(tvc)' foreach var in `e(varlist)' { _ms_parse_parts `var' if `"`: list posof `"`r(name)'"' in at'"' == "0" { qui replace `r(name)' = 0 if `touse' if `"`: list posof `"`r(name)'"' in tmptvc'"' != "0" { forvalues i = 1/`e(df_`r(name)')' { qui replace _rcs_`r(name)'`i' = 0 if `touse' qui replace _d_rcs_`r(name)'`i' = 0 if `touse' } } } } } /* Out of sample predictions using at() */ if "`at'" != "" { tokenize `at' while "`1'"!="" { if "`1'" == "." continue, break fvunab tmpfv: `1' local 1 `tmpfv' _ms_parse_parts `1' if "`r(type)'"!="variable" { display as error "level indicators of factor" /* */ " variables may not be individually set" /* */ " with the at() option; set one value" /* */ " for the entire factor variable" exit 198 } cap confirm var `1' if _rc { di "`1' is not in the data set" } if "`2'" != "=" { if _rc { cap confirm num `2' if _rc { di as err "invalid at(... `1' `2' ...)" exit 198 } } qui replace `1' = `2' if `touse' local shift 2 } else { cap confirm var `3' if _rc { di as err "`var' is not in the data set" exit 198 } qui replace `1' = `3' if `touse' local shift 3 } if `"`: list posof `"`1'"' in etvc'"' != "0" { local tvcvar `1' if "`e(orthog)'" != "" { tempname rmatrix_`tvcvar' matrix `rmatrix_`tvcvar'' = e(R_`tvcvar') local rmatrixopt rmatrix(`rmatrix_`tvcvar'') } capture drop _rcs_`tvcvar'* _d_rcs_`tvcvar'* qui rcsgen `lnt' if `touse', knots(`e(ln_tvcknots_`tvcvar')') gen(_rcs_`tvcvar') dgen(_d_rcs_`tvcvar') `e(reverse)' `rmatrixopt' `e(nosecondder)' `e(nofirstder)' forvalues i = 1/`e(df_`tvcvar')'{ qui replace _rcs_`tvcvar'`i' = _rcs_`tvcvar'`i'*`tvcvar' if `touse' qui replace _d_rcs_`tvcvar'`i' = _d_rcs_`tvcvar'`i'*`tvcvar' if `touse' } } mac shift `shift' } } /* Add offset term if exists unless no offset option is specified */ if "`e(offset1)'" != "" & /* !! PR */ "`offset'" != "nooffset" { local addoff "+ `e(offset1)'" } /* check ci and stdp options */ if "`ci'" != "" & "`stdp'" != "" { display as error "Only one of the ci and se options can be specified" exit 198 } /* Deviance and Martingale Residuals */ if "`deviance'" != "" | "`martingale'" != "" { tempvar cH res qui predict `cH' if `touse', cumhazard timevar(`t') `offset' gen double `res' = _d - `cH' if `touse' if "`deviance'" != "" { gen double `newvarname' = sign(`res')*sqrt( -2*(`res' + _d*(ln(_d -`res')))) if `touse' } else rename `res' `newvarname' } /* Failure (1-S(t)) */ else if "`failure'" != "" { qui predict `newvarname' if `touse', s timevar(`t') `offset' ci qui replace `newvarname' = 1 - `newvarname' if `touse' tempvar tmpSt qui gen double `tmpSt' = 1 - `newvarname'_uci if `touse' qui replace `newvarname'_uci = 1 - `newvarname'_lci if `touse' qui replace `newvarname'_lci = `tmpSt' if `touse' } /* Cumulative Hazard */ else if "`cumhazard'" != "" { tempvar S predict `S' if `touse', s timevar(`t') `offset' gen double `newvarname' = -ln(`S') if `touse' } /* Cumulative Odds */ else if "`cumodds'" != "" { tempvar S predict `S' if `touse', s timevar(`t') /* !! pr */ `offset' gen double `newvarname' = (1 -`S')/`S' if `touse' } /* Standard Normal Deviate */ else if "`normal'" != "" { tempvar S predict `S' if `touse', s timevar(`t') /* !! pr */ `offset' gen double `newvarname' = -invnormal(`S') if `touse' } /* density */ else if "`density'" != "" { tempvar S h predict `S' if `touse', s timevar(`t') /* !! pr */ `offset' predict `h' if `touse', h timevar(`t') /* !! pr */ `offset' gen double `newvarname' = `S'*`h' if `touse' } /* linear predictor */ else if "`xb'" != "" { if "`ci'" != "" { local prednlopt ci(`newvarname'_lci `newvarname'_uci) } else if "`stdp'" != "" { local prednlopt se(`newvarname'_se) } qui predictnl double `newvarname' = xb(xb) `addoff' if `touse', `prednlopt' level(`level') } /* derivative of linear predictor */ else if "`dxb'" != "" { if "`ci'" != "" { local prednlopt ci(`newvarname'_lci `newvarname'_uci) } else if "`stdp'" != "" { local prednlopt se(`newvarname'_se) } qui predictnl double `newvarname' = xb(dxb) if `touse', `prednlopt' level(`level') } /* linear predictor exluding spline terms */ else if "`xbnobaseline'" != "" { if "`ci'" != "" { local prednlopt ci(`newvarname'_lci `newvarname'_uci) } else if "`stdp'" != "" { local prednlopt se(`newvarname'_se) } /* commented out for now - ignores the constant (gamma_0) - may be needed later if "`e(noconstant)'" == "" { local xbnobhpred [xb][_cons] } */ foreach var in `e(varlist)' { if "`xbnobhpred'" == "" { local xbnobhpred [xb][`var']*`var' } else { local xbnobhpred `xbnobhpred' + [xb][`var']*`var' } if `"`: list posof `"`var'"' in etvc'"' != "0" { forvalues i = 1/`e(df_`var')' { local xbnobhpred `xbnobhpred' + [xb][_rcs_`var'`i']*_rcs_`var'`i' } } } * if "`e(noconstant)'" != "" { * local xbnobhpred = subinstr("`xbnobhpred'","+","",1) * } predictnl double `newvarname' = `xbnobhpred' if `touse', `prednlopt' level(`level') } /* tvc option */ else if "`tvc'" != "" { if "`ci'" != "" { local prednlopt ci(`newvarname'_lci `newvarname'_uci) } else if "`stdp'" != "" { local prednlopt se(`newvarname'_se) } local tvcpred [xb][`tvc']*`tvc' if `"`: list posof `"`tvc'"' in etvc'"' != "0" { forvalues i = 1/`e(df_`tvc')' { local tvcpred `tvcpred' + [xb][_rcs_`tvc'`i']*_rcs_`tvc'`i' } } predictnl double `newvarname' = (`tvcpred') if `touse', `prednlopt' level(`level') } /* Survival Function */ else if "`survival'" != "" & "`uncured'" == "" { tempvar sxb if "`ci'" != "" { tempvar sxb_lci sxb_uci local prednlopt ci(`sxb_lci' `sxb_uci') } if "`e(scale)'" != "theta" { qui predictnl double `sxb' = xb(xb) `addoff' if `touse', `prednlopt' level(`level') } /* predict on ln(-ln S(t)) scale for theta */ else if "`e(scale)'" == "theta" { qui predictnl double `sxb' = ln(ln(exp(xb(ln_theta))*exp(xb(xb)`addoff')+1)/exp(xb(ln_theta))) if `touse', `prednlopt' level(`level') } /* Transform back to survival scale */ if "`e(scale)'" == "hazard" { qui gen double `newvarname' = exp(-exp(`sxb')) if `touse' if "`ci'" != "" { qui gen `newvarname'_lci = exp(-exp(`sxb_uci')) if `touse' qui gen `newvarname'_uci = exp(-exp(`sxb_lci')) if `touse' } } else if "`e(scale)'" == "odds" { qui gen double `newvarname' = (1 +exp(`sxb'))^(-1) if `touse' if "`ci'" != "" { qui gen `newvarname'_lci = (1 +exp(`sxb_uci'))^(-1) if `touse' qui gen `newvarname'_uci = (1 +exp(`sxb_lci'))^(-1) if `touse' } } else if "`e(scale)'" == "normal" { qui gen double `newvarname' = normal(-`sxb') if `touse' if "`ci'" != "" { qui gen `newvarname'_lci = normal(-`sxb_uci') if `touse' qui gen `newvarname'_uci = normal(-`sxb_lci') if `touse' } } else if "`e(scale)'" == "log" { qui gen double `newvarname' = 1 - exp(`sxb') if `touse' if "`ci'" != "" { qui gen `newvarname'_lci = 1 - exp(`sxb_uci') if `touse' qui gen `newvarname'_uci = 1 - exp(`sxb_lci') if `touse' } } else if "`e(scale)'" == "theta" { qui gen double `newvarname' = exp(-exp(`sxb')) if `touse' if "`ci'" != "" { qui gen `newvarname'_lci = exp(-exp(`sxb_lci')) if `touse' qui gen `newvarname'_uci = exp(-exp(`sxb_uci')) if `touse' } } qui replace `newvarname' = 1 if `t' == 0 & `touse' if "`ci'" != "" { qui replace `newvarname'_lci = 1 if `t' == 0 & `touse' qui replace `newvarname'_uci = 1 if `t' == 0 & `touse' } } /* Hazard Function */ else if "`hazard'" != "" & "`uncured'" == "" { tempvar lnh if "`ci'" != "" { tempvar lnh_lci lnh_uci local prednlopt ci(`lnh_lci' `lnh_uci') } if "`e(scale)'" == "hazard" { qui predictnl double `lnh' = -ln(`t') + ln(xb(dxb)) + xb(xb) `addoff' if `touse', `prednlopt' level(`level') } if "`e(scale)'" == "odds" { qui predictnl double `lnh' = -ln(`t') + ln(xb(dxb)) + (xb(xb)`addoff') -ln(1+exp(xb(xb)`addoff')) if `touse', `prednlopt' level(`level') } if "`e(scale)'" == "normal" { qui predictnl double `lnh' = -ln(`t') + ln(xb(dxb)) + ln(normalden(xb(xb)`addoff')) - ln(normal(-(xb(xb)`addoff'))) if `touse', `prednlopt' level(`level') } if "`e(scale)'" == "log" { qui predictnl double `lnh' = -ln(`t') + ln(xb(dxb)) + xb(xb) -ln(1-exp(xb(xb))) `addoff' if `touse', `prednlopt' level(`level') } if "`e(scale)'" == "theta" { qui predictnl double `lnh' = -ln(`t') + ln(xb(dxb)) + xb(`xb') - ln(exp(xb(ln_theta))*exp(xb(xb)`addoff') + 1) if `touse', `prednlopt' level(`level') } /* Transform back to hazard scale */ qui gen double `newvarname' = exp(`lnh')*`per' if `touse' if "`ci'" != "" { qui gen `newvarname'_lci = exp(`lnh_lci')*`per' if `touse' qui gen `newvarname'_uci = exp(`lnh_uci')*`per' if `touse' } } /* Predict Hazard Ratio */ else if "`hrnumerator'" != "" { tempvar lhr if `"`ci'"' != "" { tempvar lhr_lci lhr_uci local predictnl_opts ci(`lhr_lci' `lhr_uci') } else if "`stdp'" != "" { tempvar lhr_se local predictnl_opts se(`lhr_se') } forvalues i=1/`e(dfbase)' { local dxb1 `dxb1' [xb][_rcs`i']*_d_rcs`i' local dxb0 `dxb0' [xb][_rcs`i']*_d_rcs`i' local xb1_plus `xb1_plus' [xb][_rcs`i']*_rcs`i' local xb0_plus `xb0_plus' [xb][_rcs`i']*_rcs`i' if `i' != `e(dfbase)' { local dxb0 `dxb0' + local dxb1 `dxb1' + local xb1_plus `xb1_plus' + local xb0_plus `xb0_plus' + } } /* use Parse_list to select appropriate values of factor variables */ Parse_list, listname(hrnumerator) parselist(`hrnumerator') tokenize `r(retlist)' while "`1'"!="" { cap confirm var `2' if _rc { if "`2'" == "." { local 2 `1' } else { cap confirm num `2' /* if _rc { di as err "invalid hrnumerator(... `1' `2' ...)" exit 198 } */ } } if "`xb10'" != "" & "`2'" != "0" & `: list posof `"`1'"' in main_varlist' != 0{ local xb10 `xb10' + } if "`xb1_plus'" != "" & "`2'" != "0" & `: list posof `"`1'"' in main_varlist' != 0{ local xb1_plus `xb1_plus' + } if "`2'" != "0" & `: list posof `"`1'"' in main_varlist' != 0 { local xb10 `xb10' [xb][`1']*`2' local xb1_plus `xb1_plus' [xb][`1']*`2' } if `"`: list posof `"`1'"' in etvc'"' != "0" & "`2'" != "0" { if "`e(rcsbaseoff)'" == "" | (`: list posof `"`1'"' in etvc'>1) { local dxb1 `dxb1' + } local xb10 `xb10' + local xb1_plus `xb1_plus' + forvalues i=1/`e(df_`1')' { local dxb1 `dxb1' [xb][_rcs_`1'`i']*_d_rcs_`1'`i'*`2' local xb10 `xb10' [xb][_rcs_`1'`i']*_rcs_`1'`i'*`2' local xb1_plus `xb1_plus' [xb][_rcs_`1'`i']*_rcs_`1'`i'*`2' if `i' != `e(df_`1')' { local dxb1 `dxb1' + local xb10 `xb10' + local xb1_plus `xb1_plus' + } } } mac shift 2 } if "`hrdenominator'" != "" { /* use Parse_list to select appropriate values of factor variables */ Parse_list, listname(hrdenominator) parselist(`hrdenominator') tokenize `r(retlist)' while "`1'"!="" { cap confirm var `2' if _rc { if "`2'" == "." { local 2 `1' } else { /* cap confirm num `2' if _rc { di as err "invalid hrdenominator(... `1' `2' ...)" exit 198 } */ } } if "`2'" != "0" & `: list posof `"`1'"' in main_varlist' != 0 { local xb10 `xb10' - [xb][`1']*`2' if "`e(rcsbaseoff)'" == "" & `: list posof `"`1'"' in main_varlist' != 0 { local xb0_plus `xb0_plus' + [xb][`1']*`2' } else if `: list posof `"`1'"' in main_varlist' != 0 { local xb0_plus `xb0_plus' [xb][`1']*`2' } } if `"`: list posof `"`1'"' in etvc'"' != "0" & "`2'" != "0" { if "`e(rcsbaseoff)'" == "" | (`: list posof `"`1'"' in etvc'>1) { local dxb0 `dxb0' + } local xb0_plus `xb0_plus' + local xb10 `xb10' - forvalues i=1/`e(df_`1')' { local dxb0 `dxb0' [xb][_rcs_`1'`i']*_d_rcs_`1'`i'*`2' local xb10 `xb10' [xb][_rcs_`1'`i']*_rcs_`1'`i'*`2' local xb0_plus `xb0_plus' [xb][_rcs_`1'`i']*_rcs_`1'`i'*`2' if `i' != `e(df_`1')' { local dxb0 `dxb0' + local xb10 `xb10' - local xb0_plus `xb0_plus' + } } } mac shift 2 } } if "`e(noconstant)'" == "" { local xb0_plus `xb0_plus' + [xb][_cons] local xb1_plus `xb1_plus' + [xb][_cons] } if "`e(scale)'" =="hazard" { qui predictnl double `lhr' = ln(`dxb1') - ln(`dxb0') + `xb10' if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" =="odds" { qui predictnl double `lhr' = ln(`dxb1') - ln(`dxb0') + `xb10' - /// ln(1+exp(`xb1_plus')) + ln(1+exp(`xb0_plus')) /// if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" =="normal" { qui predictnl double `lhr' = ln(`dxb1') - ln(`dxb0') + /// ln(normalden(`xb1_plus')) - ln(normalden(`xb0_plus')) - /// ln(normal(-(`xb1_plus'))) + ln(normal(-(`xb0_plus'))) /// if `touse', `predictnl_opts' level(`level') } if "`e(scale)'" == "log" { qui predictnl double `lhr' = ln(`dxb1') - ln(`dxb0') + `xb10' /// -ln(1-exp(`xb1_plus')) + ln(1-exp(`xb0_plus')) /// if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" =="theta" { qui predictnl double `lhr' = ln(`dxb1') - ln(`dxb0') + `xb10' /// -ln(exp(xb(ln_theta))*exp(`xb1_plus') + 1) + ln(exp(xb(ln_theta))*exp(`xb0_plus') + 1) /// if `touse', `predictnl_opts' level(`level') } qui gen double `newvarname' = exp(`lhr') if `touse' if `"`ci'"' != "" { qui gen double `newvarname'_lci=exp(`lhr_lci') if `touse' qui gen double `newvarname'_uci=exp(`lhr_uci') if `touse' } else if "`stdp'" != "" { qui gen double `newvarname'_se = `lhr_se' * `newvarname' } } /* Predict Difference in Hazard Functions */ else if "`hdiff1'" != "" { if `"`ci'"' != "" { local predictnl_opts "ci(`newvarname'_lci `newvarname'_uci)" } else if "`stdp'" != "" { local predictnl_opts se(`newvarname'_se) } forvalues i=1/`e(dfbase)' { local dxb1 `dxb1' [xb][_rcs`i']*_d_rcs`i' local dxb0 `dxb0' [xb][_rcs`i']*_d_rcs`i' local xb1_plus `xb1_plus' [xb][_rcs`i']*_rcs`i' local xb0_plus `xb0_plus' [xb][_rcs`i']*_rcs`i' if `i' != `e(dfbase)' { local dxb0 `dxb0' + local dxb1 `dxb1' + local xb1_plus `xb1_plus' + local xb0_plus `xb0_plus' + } } /* use Parse_list to select appropriate values of factor variables */ Parse_list, listname(hdiff1) parselist(`hdiff1') tokenize `r(retlist)' while "`1'"!="" { cap confirm var `2' if _rc { if "`2'" == "." { local 2 `1' } else { /* cap confirm num `2' if _rc { di as err "invalid hdiff1(... `1' `2' ...)" exit 198 } */ } } if "`xb1_plus'" != "" & "`2'" != "0" & `: list posof `"`1'"' in main_varlist' != 0 { local xb1_plus `xb1_plus' + } if "`2'" != "0" & `: list posof `"`1'"' in main_varlist'!=0 { local xb1_plus `xb1_plus' [xb][`1']*`2' } if `"`: list posof `"`1'"' in etvc'"' != "0" & "`2'" != "0" { if "`e(rcsbaseoff)'" == "" { local dxb1 `dxb1' + } local xb1_plus `xb1_plus' + local xb10 `xb10' + forvalues i=1/`e(df_`1')' { local dxb1 `dxb1' [xb][_rcs_`1'`i']*_d_rcs_`1'`i'*`2' local xb1_plus `xb1_plus' [xb][_rcs_`1'`i']*_rcs_`1'`i'*`2' if `i' != `e(df_`1')' { local dxb1 `dxb1' + local xb1_plus `xb1_plus' + } } } mac shift 2 } if "`hdiff2'" != "" { /* use Parse_list to select appropriate values of factor variables */ Parse_list, listname(hdiff2) parselist(`hdiff2') tokenize `r(retlist)' while "`1'"!="" { cap confirm var `2' if _rc { if "`2'" == "." { local 2 `1' } else { /* cap confirm num `2' if _rc { di as err "invalid hdiff2(... `1' `2' ...)" exit 198 } */ } } if "`2'" != "0" { if "`e(rcsbaseoff)'" == "" & `: list posof `"`1'"' in main_varlist' != 0 { local xb0_plus `xb0_plus' + [xb][`1']*`2' } else if `: list posof `"`1'"' in main_varlist' != 0 { local xb0_plus `xb0_plus' [xb][`1']*`2' } } if `"`: list posof `"`1'"' in etvc'"' != "0" & "`2'" != "0" { if "`e(rcsbaseoff)'" == "" { local dxb0 `dxb0' + } local xb0_plus `xb0_plus' + forvalues i=1/`e(df_`1')' { local dxb0 `dxb0' [xb][_rcs_`1'`i']*_d_rcs_`1'`i'*`2' local xb0_plus `xb0_plus' [xb][_rcs_`1'`i']*_rcs_`1'`i'*`2' if `i' != `e(df_`1')' { local dxb0 `dxb0' + local xb0_plus `xb0_plus' + } } } mac shift 2 } } if "`e(noconstant)'" == "" { local xb0_plus `xb0_plus' + [xb][_cons] local xb1_plus `xb1_plus' + [xb][_cons] } if "`e(scale)'" =="hazard" { qui predictnl double `newvarname' = (1/`t' * (`dxb1')*exp(`xb1_plus') - 1/`t' * (`dxb0')*exp(`xb0_plus'))*`per' /// if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" =="odds" { qui predictnl double `newvarname' = (1/`t' *(`dxb1')*exp(`xb1_plus')/((1 + exp(`xb1_plus'))) - /// 1/`t' *(`dxb0')*exp(`xb0_plus')/((1 + exp(`xb0_plus'))))*`per' /// if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" =="normal" { qui predictnl double `newvarname' = (1/`t' *(`dxb1')*normalden(`xb1_plus')/normal(-(`xb1_plus')) - /// 1/`t' *(`dxb0')*normalden(`xb0_plus')/normal(-(`xb0_plus')))*`per' /// if `touse', `predictnl_opts' level(`level') } if "`e(scale)'" =="log" { qui predictnl double `newvarname' = (1/`t' * (`dxb1')*exp(`xb1_plus')/(1-exp(`xb1_plus')) - /// 1/`t' * (`dxb0')*exp(`xb0_plus')/(1-exp(`xb0_plus')))*`per' /// if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" =="theta" { qui predictnl double `newvarname' = (1/`t' *((`dxb1')*exp(`xb1_plus'))/((exp([ln_theta][_cons])*exp(`xb1_plus') + 1)) - /// 1/`t' *((`dxb0')*exp(`xb0_plus'))/((exp([ln_theta][_cons])*exp(`xb0_plus') + 1)))*`per' /// if `touse', `predictnl_opts' level(`level') } } /* Predict Difference in Survival Curves */ else if "`sdiff1'" != "" { if `"`ci'"' != "" { local predictnl_opts "ci(`newvarname'_lci `newvarname'_uci)" } else if "`stdp'" != "" { local predictnl_opts se(`newvarname'_se) } forvalues i=1/`e(dfbase)' { local xb1_plus `xb1_plus' [xb][_rcs`i']*_rcs`i' local xb0_plus `xb0_plus' [xb][_rcs`i']*_rcs`i' if `i' != `e(dfbase)' { local xb1_plus `xb1_plus' + local xb0_plus `xb0_plus' + } } /* use Parse_list to select appropriate values of factor variables */ Parse_list, listname(sdiff1) parselist(`sdiff1') tokenize `r(retlist)' while "`1'"!="" { cap confirm var `2' if _rc { if "`2'" == "." { local 2 `1' } else { /* cap confirm num `2' if _rc { di as err "invalid sdiff1(... `1' `2' ...)" exit 198 } */ } } if "`xb1_plus'" != "" & "`2'" != "0" & `: list posof `"`1'"' in main_varlist' != 0 { local xb1_plus `xb1_plus' + } if "`2'" != "0" & `: list posof `"`1'"' in main_varlist' != 0 { local xb1_plus `xb1_plus' [xb][`1']*`2' } if `"`: list posof `"`1'"' in etvc'"' != "0" & "`2'" != "0" { local xb1_plus `xb1_plus' + forvalues i=1/`e(df_`1')' { local xb1_plus `xb1_plus' [xb][_rcs_`1'`i']*_rcs_`1'`i'*`2' if `i' != `e(df_`1')' { local xb1_plus `xb1_plus' + } } } mac shift 2 } if "`sdiff2'" != "" { /* use Parse_list to select appropriate values of factor variables */ Parse_list, listname(sdiff2) parselist(`sdiff2') tokenize `r(retlist)' while "`1'"!="" { cap confirm var `2' if _rc { if "`2'" == "." { local 2 `1' } else { /* cap confirm num `2' if _rc { di as err "invalid sdiff2(... `1' `2' ...)" exit 198 } */ } } if "`2'" != "0" & `: list posof `"`1'"' in main_varlist' != 0 { local xb0_plus `xb0_plus' + [xb][`1']*`2' } if `"`: list posof `"`1'"' in etvc'"' != "0" & "`2'" != "0" { local xb0_plus `xb0_plus' + forvalues i=1/`e(df_`1')' { local xb0_plus `xb0_plus' [xb][_rcs_`1'`i']*_rcs_`1'`i'*`2' if `i' != `e(df_`1')' { local xb0_plus `xb0_plus' + } } } mac shift 2 } } if "`e(noconstant)'" == "" { local xb0_plus `xb0_plus' + [xb][_cons] local xb1_plus `xb1_plus' + [xb][_cons] } if "`e(scale)'" =="hazard" { qui predictnl double `newvarname' = exp(-exp(`xb1_plus')) - exp(-exp(`xb0_plus')) if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" =="odds" { qui predictnl double `newvarname' = 1/(exp(`xb1_plus')+1) - 1/(exp(`xb0_plus')+1) if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" =="normal" { qui predictnl double `newvarname' = normal(-(`xb1_plus')) - normal(-(`xb0_plus')) if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" == "log" { qui predictnl double `newvarname' = (1 - exp(`xb1_plus')) - (1-exp(`xb0_plus')) if `touse', `predictnl_opts' level(`level') } else if "`e(scale)'" =="theta" { qui predictnl double `newvarname' = (exp([ln_theta][_cons])*exp(`xb1_plus') + 1)^(-1/exp([ln_theta][_cons])) /// -(exp([ln_theta][_cons])*exp(`xb0_plus') + 1)^(-1/exp([ln_theta][_cons])) /// if `touse', `predictnl_opts' level(`level') } } /* estimate cure, survival of uncured or hazard of uncured */ else if "`cure'" != "" | "`uncured'" != "" { local xblist [xb][_cons] local rcslist local drcslist tempvar temp if "`ci'" != "" { local prednlopt ci(`temp'_lci `temp'_uci) } foreach var in `e(varlist)' { local xblist `xblist' + [xb][`var']*`var' } if "`cure'" != "" { /*if cure is specified this is what we want to estimate*/ qui predictnl double `temp' = `xblist' if `touse', `prednlopt' level(`level') qui gen double `newvarname' = exp(-exp(`temp')) if `touse' /*we model on log(-log) scale*/ if "`ci'" != "" { qui gen double `newvarname'_lci = exp(-exp(`temp'_uci)) if `touse' qui gen double `newvarname'_uci = exp(-exp(`temp'_lci)) if `touse' } } else { /*continue, estimate survival or hazard of uncured or Predicted survival time among uncured for a given centile*/ forvalues i = 1/`e(dfbase)' { if "`rcslist'" == "" local rcslist [xb][_rcs`i']*_rcs`i' /*create a list of the sum of all spline variables*/ else local rcslist `rcslist' + [xb][_rcs`i']*_rcs`i' } foreach var in `e(tvc)' { forvalues i = 1/`e(df_`var')' { local rcslist `rcslist' + [xb][_rcs_`var'`i']*_rcs_`var'`i' } } forvalues i = 1/`e(dfbase)' { if "`drcslist'" == "" local drcslist [xb][_rcs`i']*_d_rcs`i' /*need derivatives of rcslist to calculate hazard and for centile*/ else local drcslist `drcslist' + [xb][_rcs`i']*_d_rcs`i' } foreach var in `e(tvc)' { forvalues i = 1/`e(df_`var')' { local drcslist `drcslist' + [xb][_rcs_`var'`i']*_d_rcs_`var'`i' } } local pi exp(-exp(`xblist')) /*we need cure for estimation of survival and hazard*/ local exprcs exp(`rcslist') /*we need exp of the sum of all spline variables for estimation of survival and hazard*/ /*predicted survival of uncured*/ if "`survival'" != "" & "`uncured'" != "" { tokenize `e(boundary_knots)' local lastknot = `2' qui predictnl double `temp' = ln(-(ln(`pi'^(`exprcs') - `pi') - ln(1 - `pi'))) if `touse', `prednlopt' level(`level') qui gen double `newvarname' = exp(-exp(`temp')) if `touse' qui replace `newvarname' = 0 if `newvarname' == . & `t'>=`lastknot' & `touse' if "`ci'" != "" { qui gen double `newvarname'_lci = exp(-exp(`temp'_uci)) if `touse' qui gen double `newvarname'_uci = exp(-exp(`temp'_lci)) if `touse' } } /*predicted hazard of uncured*/ else if "`hazard'" != "" & "`uncured'" != "" { qui predictnl double `temp' = ln(-ln(`pi')*((`drcslist')/`t')*`exprcs'*`pi'^(`exprcs'))- ln(`pi'^(`exprcs') - `pi') if `touse', `prednlopt' level(`level') qui gen double `newvarname' = exp(`temp') if `touse' if "`ci'" != "" { qui gen double `newvarname'_lci = exp(`temp'_lci) if `touse' qui gen double `newvarname'_uci = exp(`temp'_uci) if `touse' } } /* Predicted survival time among uncured for a given centile */ else if "`centile'" != "" & "`uncured'" != ""{ tempvar centilevar gen `centilevar' = 1 - `centile'/100 if `touse' stpm2_centpred `newvarname' if `touse', centile(`centilevar') centol(`centol') offset(`offset') `uncured' `ci' level(`level') } } } /* Predicted survival time for a given centile */ /* Updated to useing Brents root finder (version 1.6.4) */ /* calls stpm2_centpred : a separate ado file */ else if "`centile'" != "" & "`uncured'" == "" { tempvar centilevar gen `centilevar' = 1 - `centile'/100 if `touse' stpm2_centpred `newvarname' if `touse', centile(`centilevar') centol(`centol') offset(`offset') `ci' level(`level') } /* restore original data and merge in new variables */ local keep `newvarname' if "`ci'" != "" { local keep `keep' `newvarname'_lci `newvarname'_uci } else if "`stdp'" != "" { local keep `keep' `newvarname'_se } keep `keep' qui save `"`newvars'"' restore merge 1:1 _n using `"`newvars'"', nogenerate noreport end /* meansurv added to stpm2_pred as sub program */ * 10March2009 - added averages for models on odds, probit or theta scales. program Stmeancurve, sortpreserve version 10.0 syntax newvarname [if] [in],[TIMEvar(varname) AT(string) ATVAR(string) noCOLLAPSEMEANSURV /// meansurvwt(varname) ci noOFFSET UNCURED MEANHAZard MEANFT LEVEL(real `c(level)') per(integer 1)] marksample touse, novarlist local newvarname `varlist' tempvar t lnt touse_time preserve if "`meanhazard'" == "" & "`meanft'" == "" local meansurv meansurv // use timevar option or _t // do not use if if "`timevar'" == "" { qui gen `t' = _t qui gen double `lnt' = ln(_t) } else { qui gen double `t' = `timevar' qui gen double `lnt' = ln(`timevar') } /* index which time units are selected */ gen `touse_time' = `t' != . if "`meanhazard'" != "" | "`meanft'" != "" { qui replace `touse_time' = 0 if `t' == 0 } /* generate ocons for use when orthogonalising splines */ tempvar ocons gen `ocons' = 1 /* Calculate new spline terms */ if "`e(rcsbaseoff)'" == "" { capture drop _rcs* _d_rcs* if "`e(orthog)'" != "" { tempname rmatrix matrix `rmatrix' = e(R_bh) local rmatrixopt rmatrix(`rmatrix') } else local rmatrixopt "" qui rcsgen `lnt' if `touse_time', knots(`e(ln_bhknots)') gen(_rcs) dgen(_d_rcs) `e(reverse)' `rmatrixopt' `e(nosecondder)' `e(nofirstder)' } if "`e(tvc)'" != "" capture drop _rcs_* _d_rcs_* foreach tvcvar in `e(tvc)' { if "`e(orthog)'" != "" { tempname rmatrix_`tvcvar' matrix `rmatrix_`tvcvar'' = e(R_`tvcvar') local rmatrixopt rmatrix(`rmatrix_`tvcvar'') } else local rmatrixopt "" qui rcsgen `lnt' if `touse_time', gen(_rcs_`tvcvar') knots(`e(ln_tvcknots_`tvcvar')') dgen(_d_rcs_`tvcvar') `e(reverse)' `rmatrixopt' `e(nosecondder)' `e(nofirstder)' } /* Out of sample predictions using at() */ if "`at'" != "" { tokenize `at' while "`1'"!="" { if "`1'" == "." continue, break unab 1: `1' cap confirm var `1' if _rc { di "`1' is not in the data set" } if "`2'" != "=" { if _rc { cap confirm num `2' if _rc { di as err "invalid at(... `1' `2' ...)" exit 198 } } qui replace `1' = `2' if `touse' local shift 2 } else { cap confirm var `3' if _rc { di as err "`var' is not in the data set" exit 198 } qui replace `1' = `3' if `touse' local shift 3 } mac shift `shift' } } if "`atvar'" != "" { // strip spaces before and after equals sign local stillspaces 1 while `stillspaces' !=0 { local newatvar: subinstr local atvar " =" "=", all local newatvar: subinstr local newatvar "= " "=", all if "`newatvar'" == "`atvar'" { local stillspaces 0 } else { local atvar `newatvar' } } tokenize `atvar' while "`1'" != "" { capture replace `1' if `touse' if _rc>0 { di "Error in atvar() option ... `atvar'" exit 198 } mac shift 2 } } foreach tvcvar in `e(tvc)' { local rcstvclist `rcstvclist' `e(rcsterms_`tvcvar')' local drcstvclist `drcstvclist' `e(drcsterms_`tvcvar')' } if "`e(scale)'" == "theta" local theta = exp([ln_theta][_cons]) qui count if `touse_time' local Nt `r(N)' local touse_msurvpop `touse' if "`collapsemeansurv'" == "" { tempvar duppattern dupfirst dupsum tmpid newtouse gen `tmpid' = _n covpat `e(varlist)' `e(tvc)' if `touse', gen(`duppattern') quietly bysort `duppattern' : gen `dupfirst' = _n == 1 if !missing(`duppattern') quietly bysort `duppattern' : gen `dupsum' = _N if !missing(`duppattern') quietly gen `newtouse' = (`dupfirst' == 1) local touse_msurvpop `newtouse' if "`meansurvwt'" != "" { tempvar ismiss capture bysort `duppattern' (`meansurvwt'): assert `meansurvwt'[1] == `meansurvwt'[_N] if _rc>0 { di as error "meansurvwt not constant within unique combinations of covariates" di as error "check correct or use nocollapsemeansurv option." } } sort `tmpid' } tempname tmpb tmpV matrix `tmpb' = e(b) matrix `tmpV' = e(V) _ms_findomitted `tmpb' `tmpV' mata: msurvpop() /* restore original data and merge in new variables */ qui replace `newvarname' = 1 if `t' == 0 & `touse_time' local keep `newvarname' if "`ci'" != "" { if "`meanhazard'" == "" { qui replace `newvarname'_lci = 1 if `t' == 0 & `touse_time' qui replace `newvarname'_uci = 1 if `t' == 0 & `touse_time' qui replace `newvarname'_lci = 0 if `newvarname'_lci <0 & `touse_time' qui replace `newvarname'_uci = 1 if `newvarname'_uci > 1 & `touse_time' } local keep `keep' `newvarname'_lci `newvarname'_uci } else if "`stdp'" != "" { local keep `keep' `newvarname'_se } keep `keep' tempfile newvars qui save `"`newvars'"' restore merge 1:1 _n using `"`newvars'"', nogenerate noreport end /* program Parse_list converts the options give in hrnum, hrdenom, hdiff1, hdiff2, sdiff1 and sdiff2 to factor notation _varlist contains the variables listed in the predict option */ program define Parse_list, rclass syntax, listname(string) parselist(string) tokenize `parselist' local etvc `e(tvc)' local main_varlist `e(varlist)' while "`1'"!="" { fvunab tmpfv: `1' local 1 `tmpfv' _ms_parse_parts `1' if "`r(type)'"!="variable" { display as error "level indicators of factor" /* */ " variables may not be individually set" /* */ " with the `listname'() option; set one value" /* */ " for the entire factor variable" exit 198 } cap confirm var `2' if _rc { cap confirm num `2' if _rc { if "`2'" != "." { di as err "invalid `listname'(... `1' `2' ...)" exit 198 } } } local _varlist `_varlist' `1' local `1'_value `2' mac shift 2 } _ms_extract_varlist `e(varlist)', noomitted local varlist_omitted `r(varlist)' /* check if any tvc variables not in varlist_omitted */ local tvconly foreach tvcvar in `e(tvc)' { mata: st_local("addtvconly",strofreal(subinword(st_local("varlist_omitted"),st_local("tvcvar"),"")==st_local("varlist_omitted"))) if `addtvconly' { local tvconly `tvconly' `tvcvar' } } /* loop over all variables in model */ foreach var in `varlist_omitted' `tvconly' { _ms_parse_parts `var' local vartype `r(type)' local intmult foreach parse_var in `_varlist' { /* check parse_var in model */ /* _ms_extract_varlist `parse_var' if "`r(varlist)'" == "" { display as error "`parse_var' is not included in the model" exit 198 } */ * NOW SEE IF MODEL VARIABLE IS LISTED IN PARSE_VAR _ms_parse_parts `var' local invar 0 if "`r(k_names)'" == "" { if "`r(name)'" == "`parse_var'" { local invar 1 } } else { forvalues i = 1/`r(k_names)' { if "`r(name`i')'" == "`parse_var'" { local invar 1 } } } if `invar' { if "`vartype'" == "variable" { local retlist `retlist' `var' ``parse_var'_value' } else if "`vartype'" == "factor" { if `r(level)' == ``parse_var'_value' { local retlist `retlist' `var' 1 } else { _ms_extract_varlist ``parse_var'_value'.`parse_var' } } else if "`vartype'" == "interaction" { if strpos("`var'","`parse_var'") >0 { _ms_parse_parts `var' forvalues i = 1/`r(k_names)' { if "`r(name`i')'" == "`parse_var'" { if "`r(op`i')'" == "``parse_var'_value'" { local intmult `intmult'*1 } else if "`r(op`i')'" == "c" { local intmult `intmult'*`parse_var' } else { local intmult `intmult'*0 } } } } else { local intmult `intmult'*0 } } else if "`vartype'" == "product" { display "products not currently available" } } } if "`vartype'" == "interaction" { // & `invar' { local intmult: subinstr local intmult "*" "" if `intmult' != 0 { local retlist `retlist' `var' `intmult' } } return local retlist `retlist' } end program define _rmst, sortpreserve /* Mean or restricted mean or restriced SD of time to failure for an stpm2 model, based on covariate patterns or an `at()' specification. */ version 11.1 syntax [if] [in], [ at(string) GENerate(string) N(int 1001) /// POWer(real 1) RSDst TMAx(real 0) TMIn(real 0) ZEROs ] local cmdline `e(cmdline)' if (`tmin' < 0) local tmin 0 if (`tmax' <= 0) local tmax . local survival survival // Use original variable names saved by stpm2 local varlist `e(varnames)' quietly { marksample touse replace `touse' = 0 if _st == 0 tempvar t s if missing(`tmax') { noi di as txt "[tmax() not specified, estimating unrestricted mean survival time]" // Find t corresponding to 99.999999 centile for range of t local cmax 99.999999 predict `t', centile(`cmax') sum `t', meanonly local tmax = r(max) drop `t' } local n_orig = _N if `n' > `n_orig' { tempvar mark_orig gen byte `mark_orig' = 1 } range `t' `tmin' `tmax' `n' tempvar integral /* `power' option is undocumented and not implemented for rsdst */ if `power' != 1 { tempvar lnt gen `lnt' = cond(`power' == 0, ln(`t') - ln(`tmax'), /// (`t'^`power' - 1) / `power' - (`tmax'^`power' - 1) / `power') } gen `integral' = . /* identify covariate patterns and compute mean survival for them */ if "`zeros'" == "" { if "`at'" != "" { // Identify covariates `s(omitted)' in `varlist' but not in `at'. CheckAt "`at'" "`varlist'" local covpatvarlist `s(omitted)' } else local covpatvarlist `varlist' } if "`covpatvarlist'" != "" { tempvar covpat covpat `covpatvarlist' if `touse', generate(`covpat') sum `covpat', meanonly local ncovpat = r(max) noi di as txt _n "Processing " as res `ncovpat' as txt " covariate patterns ..." } else { local ncovpat 1 local covpat 1 } forvalues i = 1 / `ncovpat' { local At tokenize `covpatvarlist' while "`1'" != "" { sum `1' if (`covpat' == `i') & (`touse' == 1), meanonly local mm = r(mean) local At `At' `1' `mm' mac shift } *noi di in red `"predict double `s' if `touse', timevar(`t') at(`At' `at') `survival' `meansurv' `zeros'"' predict double `s' /*if `touse'*/, timevar(`t') at(`At' `at') `survival' `zeros' replace `s' = 1 if `t'==0 if `power' != 1 { // work with f(t) - f(tmax), so that on power-transformed scale, tmax (effective) = 0. replace `s' = 1 - `s' // distribution function integ `s' `lnt' if !missing(`lnt') replace `integral' = cond(`power' == 0, ln(`tmax'), (`tmax'^`power'-1)/`power') - r(integral) if `covpat' == `i' } else { // !! PR bug fix 14jul2016: get correct integrated survival between tmin > 0 and tmax = t* if `tmin' > 0 { // direct integration of survival function integ `s' `t' local ex = r(integral) } else { // via integration of distribution function, F(t) = 1-S(t) replace `s' = cond(`t' == 0, 0, 1 - `s') integ `s' `t' local ex = `tmax' - r(integral) } if "`rsdst'" != "" { // supported only if tmin = 0 replace `s' = `s' * `t' integ `s' `t' local ex2 = `tmax'^2 - 2 * r(integral) replace `integral' = sqrt(`ex2' - `ex'^2) if `covpat' == `i' } else replace `integral' = `ex' if `covpat' == `i' } drop `s' if mod(`i', 10) == 0 noi di as txt `i', _cont } if !missing("`generate'") { cap confirm var `generate', exact if (c(rc) == 0) replace `generate' = `integral' else rename `integral' `generate' /* local txt = cond("`rsdst'" != "", "SD", "mean") if (missing("`cmax'")) lab var `generate' "restricted `txt' at time `tmax'" else lab var `generate' "`txt' time" */ } if `n' > `n_orig' { drop if missing(`mark_orig') drop `mark_orig' } } end * version 1.1 PCL 28JUL2014 (now allows for factor variables) * Based on PR 9-Jan-94. Based on covariate-pattern counter in lpredict.ado. program define covpat, sortpreserve version 11 syntax [varlist(default=none fv)] [if] [in], Generate(string) confirm new var `generate' tempvar keep fvexpand `varlist' local expvarlist `r(varlist)' // The following is changes as expr_query did not work in Stata 12. /* local expvars: subinstr local expvarlist " " "+", all expr_query `expvars' local uniquevars `r(varnames)' */ local uniquevars foreach var in `expvarlist' { _ms_parse_parts `var' if "`r(k_names)'" == "" { local uniquevars `uniquevars' `r(name)' } else { forvalues i = 1/`r(k_names)' { local uniquevars `uniquevars' `r(name`i')' } } } local uniquevars: list uniq uniquevars quietly { marksample keep sort `keep' `uniquevars' gen long `generate' = . by `keep' `uniquevars': replace `generate'=cond(_n==1 & `keep',1,.) replace `generate' = sum(`generate') replace `generate'=. if `generate'==0 } end program define _rmst_meansurv, sortpreserve version 11.1 /* Mean or restricted mean or restriced SD of time to failure for an stpm2 model, based on covariate patterns or an `at()' specification. Assumes stpm2 model has been fitted. */ syntax [if] [in], [ at(string) GENerate(string) N(int 1001) meansurv /// POWer(real 1) RSDst TMIn(real 0) TMAx(numlist >0 ascending) ZEROs ] local cmdline `e(cmdline)' local tmin 0 local ntstar : word count `tmax' local tbig : word `ntstar' of `tmax' // Use original variable names saved by stpm2 local varlist `e(varnames)' quietly { marksample touse replace `touse' = 0 if _st == 0 tempvar t s integral local n_orig = _N local n_plus = `n' + `ntstar' - 1 if `n_plus' > `n_orig' { tempvar mark_orig gen byte `mark_orig' = 1 set obs `n_plus' } range `t' `tmin' `tbig' `n' if `ntstar' > 1 { // Augment t with tstar values excluding highest local n1 = `n' + 1 local ntstar1 = `ntstar' - 1 tokenize `tmax' forvalues j = 1/`ntstar1' { replace `t' = ``j'' in `n1' local ++n1 } } gen `integral' = . /* identify covariate patterns and compute mean survival for them */ if "`zeros'" == "" { if "`at'" != "" { // Identify covariates `s(omitted)' in `varlist' but not in `at'. CheckAt "`at'" "`varlist'" local covpatvarlist `s(omitted)' } else local covpatvarlist `varlist' } predict double `s', timevar(`t') at(`At' `at') meansurv `zeros' replace `s' = cond(`t' == 0, 0, 1 - `s') // eval integral on each (0,tstar] interval tokenize `tmax' forvalues j = 1/`ntstar' { local ts ``j'' integ `s' `t' if `t'<=(`ts' * 1.000001) local ex = `ts' - r(integral) if "`rsdst'" != "" { replace `s' = `s' * `t' integ `s' `t' if `t'<=`ts' local ex2 = `ts'^2 - 2 * r(integral) replace `integral' = sqrt(`ex2' - `ex'^2) in `j' } else replace `integral' = `ex' in `j' } drop `s' if !missing("`generate'") { cap confirm var `generate', exact if (c(rc) == 0) replace `generate' = `integral' else rename `integral' `generate' } if `n_plus' > `n_orig' { drop if missing(`mark_orig') drop `mark_orig' } *stop } end program define CheckAt, sclass version 11 /* Checks whether the `at' list specifies all variables in varlist. If not, provides in s(omitted) a list of omitted variables */ args at varlist sret clear local omitted `varlist' tokenize `at' while "`1'"!="" { fvunab v: `1' local 1 `v' _ms_parse_parts `1' if "`r(type)'"!="variable" { display as error "level indicators of factor" /* */ " variables may not be individually set" /* */ " with the at() option; set one value" /* */ " for the entire factor variable" exit 198 } cap confirm var `2' if _rc { cap confirm num `2' if _rc { di as err "invalid at(... `1' `2' ...)" exit 198 } } ChkIn `1' "`varlist'" local omitted : list omitted - 1 macro shift 2 } sreturn local omitted `omitted' end program define ChkIn, sclass version 9.2 * Returns s(k) = index # of target variable v in varlist, or 0 if not found. args v varlist sret clear local k: list posof "`v'" in varlist sret local k `k' if `s(k)' == 0 { di as err "`v' is not a valid covariate" exit 198 } end // mata program to obtain mean survival mata: void msurvpop() { // Transfer data from Stata newvar = st_local("newvarname") touse = st_local("touse_msurvpop") touse_time = st_local("touse_time") tvcrcslist = st_local("rcstvclist") Nt = strtoreal(st_local("Nt")) hascons = st_global("e(noconstant)") == "" uncured = st_local("uncured") != "" meansurv_opt = st_local("meansurv") != "" meanhazard_opt = st_local("meanhazard") != "" per = strtoreal(st_local("per")) meanft_opt = st_local("meanft") != "" t = st_data( ., st_local("t"), touse_time) if (st_global("e(rcsbaseoff)") == "") { rcsbase = st_data( ., tokens(st_global("e(rcsterms_base)")), touse_time) if(meanhazard_opt | meanft_opt) drcsbase = st_data( ., tokens(st_global("e(drcsterms_base)")), touse_time) } else rcsbase = J(Nt,0,0) if (tvcrcslist != "") { rcstvc = st_data( ., tokens(tvcrcslist), touse_time) drcstvc = st_data( ., tokens(st_local("drcstvclist")), touse_time) } else rcstvc = I(0) if (st_global("e(varlist)") != "") { x = st_data(.,tokens(st_global("e(varlist)")),touse) } else x = J(1,0,.) if (st_global("e(tvc)") != "") { xtvc = st_data(.,tokens(st_global("e(tvc)")),touse) } else xtvc = J(1,0,.) tvcvar = tokens(st_global("e(tvc)")) ntvc = cols(tvcvar) Nvarlist = cols(tokens(st_global("e(varlist)"))) // Beta matrix tmpb = st_local("tmpb") beta = st_matrix(tmpb)'[1..cols(rcsbase)+cols(rcstvc) + Nvarlist + hascons,1] betarcs = beta[(Nvarlist+1)..(Nvarlist + cols(rcsbase)+cols(rcstvc)),1] scale = st_global("e(scale)") if (scale == "theta") theta = strtoreal(st_local("theta")) // check whether to include offset offset_name = st_global("e(offset1)") if (offset_name != "" & st_local("offset") != "nooffset") offset = st_data(.,offset_name,touse) else offset = 0 startstop = J(ntvc,2,.) tvcpos = J(ntvc,1,.) tmpstart = 1 // Loop over number of time observations for (i=1;i<=ntvc;i++){ startstop[i,1] = tmpstart tmpntvc = cols(tokens(st_global("e(rcsterms_"+tvcvar[1,i]+")"))) startstop[i,2] = tmpstart + tmpntvc - 1 tmpstart = startstop[i,2] + 1 } if(st_local("collapsemeansurv")!="") { Nx = rows(x) Nobs = Nx covfreq = J(Nx,1,1) } else { covfreq = st_data(.,tokens(st_local("dupsum")),touse) Nx = rows(covfreq) Nobs = sum(covfreq) } if(st_local("meansurvwt") == "") { wt = J(Nx,1,1) } else { wt = st_data(.,st_local("meansurvwt"),touse) } // Loop over all selected observations meansurv = J(Nt,1,0) if(meanhazard_opt | meanft_opt) meanft = J(Nt,1,0) if (hascons) addcons = J(Nt,1,1) else addcons = J(Nt,0,0) if(uncured) { pi_cure = exp(-exp((x[,])*(beta[1..Nvarlist]) :+ hascons :* beta[rows(beta)])) start_exprcs = Nvarlist :+ 1 stop_exprcs = rows(beta) :- 1 if(!hascons) stop_exprcs = stop_exprcs :+ 1 } for (i=1;i<=Nx;i++) { tmprcs = J(Nt,0,.) tmpdrcs = J(Nt,0,.) for (j=1;j<=ntvc;j++) { tmprcs = tmprcs, rcstvc[,startstop[j,1]..startstop[j,2]]:*xtvc[i,j] tmpdrcs = tmpdrcs, drcstvc[,startstop[j,1]..startstop[j,2]]:*xtvc[i,j] } if(!uncured) { if (scale == "hazard") { meansurv = meansurv :+ wt[i,]:*covfreq[i,]:*exp(-exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*beta :+ offset)):/Nobs if(meanhazard_opt | meanft_opt) { meanft = meanft + wt[i,]:*covfreq[i,]:*1:/t :*((drcsbase,tmpdrcs)*betarcs):*exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*beta - exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*beta :+ offset)):/Nobs } } else if (scale == "odds") meansurv = meansurv :+ wt[i,]:*covfreq[i,]:*((1 :+ exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*beta :+ offset)):^(-1)):/Nobs else if (scale == "normal") meansurv = meansurv :+ wt[i,]:*covfreq[i,]:* normal(-(J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*beta :+ offset):/Nobs else if (scale == "theta") meansurv = meansurv :+ wt[i,]:*covfreq[i,]:* ((theta:*exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*beta :+ offset) :+ 1):^(-1/theta)):/Nobs } else if(uncured) { exprcs = exp((rcsbase,tmprcs) * beta[start_exprcs..stop_exprcs] :+ offset) meansurv = meansurv :+ wt[i,]:*covfreq[i,]:* ((pi_cure[i]:^(exprcs) :- pi_cure[i]):/(1 :- pi_cure[i])):/Nobs } } (void) st_addvar("double",newvar) if(meansurv_opt) st_store(., newvar, touse_time, meansurv) else if(meanhazard_opt) { meanhaz_result = meanft:/meansurv st_store(., newvar, touse_time, per:*meanhaz_result) } else if(meanft_opt) st_store(., newvar, touse_time, meanft) // Calculate CI using delta method if(st_local("ci") != "") { Nparams = rows(beta) stata("_ms_omit_info " + tmpb) omitvars = st_matrix("r(omit)")[,1..Nparams] V = st_matrix("e(V)")[1..Nparams,1..Nparams] G = J(Nt,0,.) for(k=1;k<=Nparams;k++) { if(omitvars[k]==1) { G = (G,J(Nt,1,1)) continue } h = sqrt(epsilon(1))*beta[k] newbeta1 = beta newbeta1[k] = beta[k] + h/2 newbeta1rcs = newbeta1[(Nvarlist+1)..(Nvarlist + cols(rcsbase)+cols(rcstvc)),1] newbeta2 = beta newbeta2[k] = beta[k] - h/2 newbeta2rcs = newbeta2[(Nvarlist+1)..(Nvarlist + cols(rcsbase)+cols(rcstvc)),1] f1 = J(Nt,1,0) f2 = J(Nt,1,0) if(uncured) { pi_cure1 = exp(-exp((x[,])*(newbeta1[1..Nvarlist]) :+ hascons :* newbeta1[rows(beta)])) pi_cure2 = exp(-exp((x[,])*(newbeta2[1..Nvarlist]) :+ hascons :* newbeta2[rows(beta)])) } if(meanhazard_opt) { f1b = J(Nt,1,0) f2b = J(Nt,1,0) } // loop over unique x's for (i=1;i<=Nx;i++) { tmprcs = J(Nt,0,.) tmpdrcs = J(Nt,0,.) for (j=1;j<=ntvc;j++) { tmprcs = tmprcs, rcstvc[,startstop[j,1]..startstop[j,2]]:*xtvc[i,j] tmpdrcs = tmpdrcs, drcstvc[,startstop[j,1]..startstop[j,2]]:*xtvc[i,j] } if(!uncured) { if (scale == "hazard") { f1 = f1 :+ wt[i,]:*covfreq[i,]:*exp(-exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta1 :+ offset)):/Nobs f2 = f2 :+ wt[i,]:*covfreq[i,]:*exp(-exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta2 :+ offset)):/Nobs } else if(scale == "odds") { f1 = f1 :+ wt[i,]:*covfreq[i,]:*((1 :+ exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta1 :+ offset)):^(-1)):/Nobs f2 = f2 :+ wt[i,]:*covfreq[i,]:*((1 :+ exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta2 :+ offset)):^(-1)):/Nobs } else if(scale == "normal") { f1 = f1 :+ wt[i,]:*covfreq[i,]:* normal(-(J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta1 :+ offset):/Nobs f2 = f2 :+ wt[i,]:*covfreq[i,]:* normal(-(J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta2 :+ offset):/Nobs } else if(scale == "theta") { f1 = f1 :+ wt[i,]:*covfreq[i,]:* ((theta:*exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta1 :+ offset) :+ 1):^(-1/theta)):/Nobs f2 = f2 :+ wt[i,]:*covfreq[i,]:* ((theta:*exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta2 :+ offset) :+ 1):^(-1/theta)):/Nobs } } else if(uncured) { exprcs1 = exp((rcsbase,tmprcs) * newbeta1[start_exprcs..stop_exprcs] :+ offset) exprcs2 = exp((rcsbase,tmprcs) * newbeta2[start_exprcs..stop_exprcs] :+ offset) f1 = f1 :+ wt[i,]:*covfreq[i,]:*((pi_cure1[i]:^(exprcs1) :- pi_cure1[i]):/(1 :- pi_cure1[i])):/Nobs f2 = f2 :+ wt[i,]:*covfreq[i,]:*((pi_cure2[i]:^(exprcs2) :- pi_cure2[i]):/(1 :- pi_cure2[i])):/Nobs } if(meanhazard_opt) { if (scale == "hazard") { f1b = f1b :+ wt[i,]:*covfreq[i,]:*1:/t :*((drcsbase,tmpdrcs)*newbeta1rcs):*exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta1 - exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta1 :+ offset)):/Nobs f2b = f2b :+ wt[i,]:*covfreq[i,]:*1:/t :*((drcsbase,tmpdrcs)*newbeta2rcs):*exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta2 - exp((J(Nt,1,x[i,]),rcsbase, tmprcs,addcons)*newbeta2 :+ offset)):/Nobs } } } if(meansurv_opt) G = (G,(log(f1) :- log(f2)):/h) else if(meanhazard_opt)G = (G,(log(f1b:/f1) :- log(f2b:/f2)):/h) } //se = sqrt(diagonal(G*V*G')) se = sqrt(rowsum((G*V):*G)) // This is quicker // calculate confidence interval level = strtoreal(st_local("level")) z = abs(invnormal((1-(level/100))/2)) (void) st_addvar("double",newvar+"_lci") (void) st_addvar("double",newvar+"_uci") if(meansurv_opt) { st_store(., newvar+"_lci", touse_time, exp(log(meansurv) - z*se)) st_store(., newvar+"_uci", touse_time, exp(log(meansurv) + z*se)) } else if(meanhazard_opt){ per st_store(., newvar+"_lci", touse_time, exp(log(meanhaz_result) - z*se):*per) st_store(., newvar+"_uci", touse_time, exp(log(meanhaz_result) + z*se):*per) } } } end **************************************************** * Predict observed survival (relative survival times expected survival) from a relative survival model * Added by Therese Andersson *************************** program Stobssurv, sortpreserve version 11.1 syntax newvarname [if] [in],[TIMEvar(varname) noOFFSET MERGEBY(string) DIAGAge(name) DIAGYear(name) BY(string) GRpd ATTAge(name) /// ATTyear(name) ID(name) MAXAge(int 99) MAXYear(int 2050) Nobs(int 100) CI LEVel(real `c(level)') MAXT(real 10) survprob(name) USIng(string) EXPSURV(string)] marksample touse, novarlist local newvarname `varlist' if "`e(scale)'" != "hazard" | "`e(bhazard)'" == ""{ display as error "You can only use the obssurv option after fitting a relative survival model with the scale(hazard) option" exit 198 } if "`grpd'"=="" & "`by'" != "" { display as error "You can not specify the by option without the grpd option" exit 198 } *Set all variable names needed for strs, if they have not been specified in options ********************** local atstep = `maxt'/`nobs' if "`id'" == "" { local id id } if "`survprob'" == "" { local survprob prob } if "`by'" != "" { local byexp by(`by') } * Now preserve the data and recalculate the splines if timevar option is used *************************** preserve tempvar t lnt /* use timevar option or _t */ if "`timevar'" == "" { qui gen `t' = _t if `touse' qui gen double `lnt' = ln(_t) if `touse' } else { qui gen double `t' = `timevar' if `touse' qui gen double `lnt' = ln(`timevar') if `touse' } /* Calculate new spline terms if timevar option is specified */ if "`timevar'" != "" & "`e(rcsbaseoff)'" == "" { drop _rcs* _d_rcs* if "`e(orthog)'" != "" { tempname rmatrix matrix `rmatrix' = e(R_bh) local rmatrixopt rmatrix(`rmatrix') } qui rcsgen `lnt' if `touse', knots(`e(ln_bhknots)') gen(_rcs) dgen(_d_rcs) `e(reverse)' `rmatrixopt' `e(nosecondder)' `e(nofirstder)' } /* calculate new spline terms if timevar option is specified */ if "`timevar'" != "" { foreach tvcvar in `e(tvc)' { capture drop _rcs_`tvcvar'* _d_rcs_`tvcvar'* if "`e(orthog)'" != "" { tempname rmatrix_`tvcvar' matrix `rmatrix_`tvcvar'' = e(R_`tvcvar') local rmatrixopt rmatrix(`rmatrix_`tvcvar'') } qui rcsgen `lnt' if `touse', gen(_rcs_`tvcvar') knots(`e(ln_tvcknots_`tvcvar')') dgen(_d_rcs_`tvcvar') `e(reverse)' `rmatrixopt' `e(nosecondder)' `e(nofirstder)' forvalues i = 1/`e(df_`tvcvar')'{ qui replace _rcs_`tvcvar'`i' = _rcs_`tvcvar'`i'*`tvcvar' if `touse' qui replace _d_rcs_`tvcvar'`i' = _d_rcs_`tvcvar'`i'*`tvcvar' if `touse' } } } * We need the predictions xb(xb), to get the relative survival *************************** tempvar sxb if "`ci'" != "" { tempvar sxb_uci sxb_lci local prednlopt ci(`sxb_lci' `sxb_uci') } qui predictnl double `sxb' = xb(xb)`addoff' if `touse', `prednlopt' level(`level') * Calculate the expected survival using strs ************************************************** tempfile grpdat inddat tempvar St_star tempfail temptime gen `tempfail'=0 gen `temptime'=`maxt' qui stset `temptime', failure(`tempfail'=1) id(`id') qui strs using "`using'" if `touse', br(0(`atstep')`maxt') mergeby(`mergeby') diagage(`diagage') diagyear(`diagyear') /// savgroup(`grpdat', replace) savind(`inddat', replace) attage(`attage') attyear(`attyear') maxage(`maxage') /// notables survprob(`survprob') `byexp' ederer1 qui stset `t', failure(`tempfail'=1) id(`id') if "`grpd'" != "" { if "`by'" != "" { joinby `by' using `grpdat', nolabel unmatched(both) assert _merge==3 if `touse' assert _merge==3 | _merge==1 drop _merge } else { cross using `grpdat' } sort `id' end gen `St_star'=cp_e1 if `touse' } if "`grpd'" == "" { qui merge 1:m `id' using `inddat', keepusing(p_star end) assert(match master) drop _merge sort `id' end by `id': gen double `St_star' = exp(sum(ln(p_star))) if `touse' } qui replace _t = float(_t) if float(_t)==`maxt' qui keep if end>=_t qui by `id': keep if _n==1 * Use expected and relative survival to obtain observed survival ***************************** qui gen double `newvarname' = exp(-exp(`sxb'))*`St_star' if `touse' if "`ci'" != "" { qui gen `newvarname'_lci = exp(-exp(`sxb_uci'))*`St_star' if `touse' qui gen `newvarname'_uci = exp(-exp(`sxb_lci'))*`St_star' if `touse' } * restore original data and merge in new variables ******************************* local keep `id' `newvarname' if "`ci'" != "" { local keep `keep' `newvarname'_lci `newvarname'_uci } if "`expsurv'" != "" { gen `expsurv'= `St_star' local keep `keep' `expsurv' } keep `keep' tempfile newvars qui save `newvars' restore merge 1:1 `id' using `newvars', nogenerate noreport end *********************************************************************************************** * Program for calculating loss in expectation of life, ADDED Feb 2012 by Therese Andersson * Nov 2012, changed the way the numerical integration is done, to speed up calculation *********************************************************************************************** program Stlifelost, sortpreserve version 11.1 syntax newvarname [if] [in],[ noOFFSET MERGEBY(string) DIAGAge(name) DIAGYear(name) MAXAge(int 99) ATTAge(name) ATTyear(name) survprob(name) USIng(string) /// BY(string) MAXYear(int 2050) Nodes(int 50) Tinf(int 50) Tcond(int 0) GRpd Stub(string) CI LEVel(real `c(level)') ] marksample touse, novarlist local newvarname `varlist' /*variables names for mean observed and mean expected survival, if stub option used*/ local meanexp `stub'exp local meanobs `stub'obs * CHECKS ********** if "`e(scale)'" != "hazard" | "`e(bhazard)'" == ""{ display as error "You can only use the lifelost option after fitting a relative survival model with the scale(hazard) option" exit 198 } if "`grpd'"=="" & "`by'" != "" { display as error "You can not specify the by option without the grpd option" exit 198 } if `nodes' < 0 { display as err "You need non-negative number of nodes" exit 198 } if `tcond' < 0 { display as err "Tcond must be non-negative" exit 198 } if `tinf' < 0 { display as err "Tinf must be non-negative" exit 198 } *Set all variable names needed for merging on p_star, if they have not been specified in options ********************** if "`survprob'" == "" { local survprob prob } ******************************************************************** * Now preserve the data and loop over each year until tinf * * and merge on interval specific expected survival at every year * * calculate cumulative expected survival at each year * ******************************************************************** preserve tempvar obsn qui gen `obsn' = _n tempvar S_star_`tcond' qui gen `S_star_`tcond''=1 if `touse' local b=`tcond'+1 forvalues i=`b'/`tinf' { tempvar S_star_`i' `attage' `attyear' local j=`i'-1 qui gen `attage'=floor(min(`diagage'+ `j',`maxage')) if `touse' qui gen `attyear'=floor(min(`diagyear' + `j', `maxyear')) if `touse' sort `mergeby' qui merge m:1 `mergeby' using `"`using'"', nolabel keepusing(`survprob') /// keep(1 3) noreport sorted capture assert _merge==3 if `touse' if _rc>0 { display as error "Records fail to match with the population file (`using'.dta)." exit 459 } *qui gen `p_star_`i''=`survprob' if `touse' qui gen `S_star_`i''=`S_star_`j''*`survprob' if `touse' drop `attage' `attyear' `survprob' _merge } /*if using grouped expected survival, calculate mean according to Ederer1*/ if "`grpd'" != "" { if "`by'"!= "" { forvalues i=`b'/`tinf' { tempvar S_star_`i'_grp qui bysort `by': egen `S_star_`i'_grp' =mean(`S_star_`i'') qui replace `S_star_`i''= `S_star_`i'_grp' if `touse' qui drop `S_star_`i'_grp' } } else { forvalues i=`b'/`tinf' { tempvar S_star_`i'_grp qui egen `S_star_`i'_grp'=mean(`S_star_`i'') qui replace `S_star_`i''= `S_star_`i'_grp' if `touse' qui drop `S_star_`i'_grp' } } } ****************************************************************** * Find the nodes and the weights for the numerical integration * * using the guassquad command and then calculate the time points * ****************************************************************** tempname weightsmat nodesmat mata gq("`weightsmat'","`nodesmat'") forvalues i=1/`nodes' { // loop over all nodes to create the time points local t`i'= (`tinf'-`tcond')*0.5*el(`nodesmat',`i',1)+(`tinf'+`tcond')*0.5 tempvar tvar`i' qui gen `tvar`i''=`t`i'' if `touse' } **************************************************************************** * Calculate cumulative expected survival at every time point of interest, * * and multiply with the weights, and do the integration (summation) * **************************************************************************** forvalues i=1/`nodes' { tempvar S`i' S_W_`i' local floort=floor(`t`i'') local ceilt=ceil(`t`i'') local dist=`t`i''-`floort' qui gen `S`i''=`S_star_`floort''-(`S_star_`floort''-`S_star_`ceilt'')*`dist' if `touse' qui gen `S_W_`i''= `S`i''*el(`weightsmat',`i',1) if `touse' } /*Sum up to get the mean expected survival*/ local SW_list `S_W_1' forvalues i=2/`nodes' { local SW_list `SW_list' `S_W_`i'' } qui egen `meanexp' = rowtotal(`SW_list') if `touse' qui replace `meanexp' = `meanexp'*(`tinf'-`tcond')*0.5 if `touse' ****************************************************************************************** * Loop over each time point and save R(t)S*(t) (the integrand) for each time point * * estimate the mean observed survival and CI, then calculate loss in expectation of life * ****************************************************************************************** local predictstat (predict(survival timevar(`tvar1'))*`S_W_1') forvalues i=2/`nodes' { local predictstat `predictstat' + (predict(survival timevar(`tvar`i''))*`S_W_`i'') } tempvar sxb if "`ci'" != "" { local prednlopt ci(`meanobs'_lci `meanobs'_uci) local levelci level(`level') } if `tcond'==0 { qui predictnl double `meanobs' = 0.5*`tinf'*(`predictstat') if `touse', `prednlopt' `levelci' } else { display as text "Estimating conditional loss in expectation of life" tempvar t_cond qui gen `t_cond'=`tcond' qui predictnl double `meanobs' = 0.5*(`tinf'-`tcond')*(`predictstat')/(predict(survival timevar(`t_cond'))) if `touse', `prednlopt' `levelci' } qui gen `newvarname' = `meanexp'-`meanobs' if `touse' if "`ci'" != "" { gen `newvarname'_lci= `meanexp'-`meanobs'_uci if `touse' gen `newvarname'_uci= `meanexp'-`meanobs'_lci if `touse' } * Restore original data and merge in new variables ************************************************************ // local keep `newvarname' local keep `id' `newvarname' if "`ci'" != "" { local keep `keep' `newvarname'_lci `newvarname'_uci } if "`stub'" != "" { local keep `keep' `meanexp' `meanobs' if "`ci'" != "" { local keep `keep' `meanobs'_lci `meanobs'_uci } } sort `obsn' keep `keep' tempfile newvars qui save `"`newvars'"' restore merge 1:1 _n using `"`newvars'"', nogenerate noreport end mata: void gq(string scalar weightsname, string scalar nodesname) { 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)'] st_matrix(weightsname,weights) st_matrix(nodesname,nodes) } end