*! version 2.0.2 16dec2015 by alexis dot dinno at pdx dot edu *! discrete-time event history analysis command with parsimonious *! smoothing of time, different estimation strategies, graphing *! and more * Copyright Notice * dthaz and dthaz.ado are Copryright (c) 2001, 2015 Alexis Dinno * * This file is part of dthaz. * * dthaz is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) at any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program (dthaz.copying); if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA * Syntax: dthaz [varlist] [if] [in] [fweight pweight iweight/] [, * specify(#[,#,#,...]) tpar(#) truncate(#) pretrunc(#) display(#) level(#) * model Link(string) cluster(varlist) reuse suppress graph(#) XLAbel(string) * YLAbel(string) XTick(numlist) YTick(numlist) graph_options copyleft] * Check for version compatibility and notify user of version incompatibility * and let them know I am ammennable to making back-compatible revisions. program define dthaz if int(_caller())<7 { di in r "dthaz- does not support this version of Stata." _newline di as txt "Requests for a v6 compatible version may be challenging to honor." di as txt "Requests for a version compatible with versions of STATA earlier than v6 are " di as txt "untenable since I do not have access to the software." _newline di as txt "All requests are welcome and will be considered." exit } if int(_caller()) >= 8 { dthaz8 `0' } if int(_caller()) == 7 { dthaz7 `0' } end ******************************************************************************* ******************************************************************************* * dthaz for STATA Versions 8+ * ******************************************************************************* ******************************************************************************* program define dthaz8, eclass version 8.0 syntax [varlist] [if] [in] [fweight pweight iweight] [, SPecify(numlist) /* */ TPar(integer -1) Truncate(integer 0) Pretrunc(integer 0) /* */ Display(integer 0) Model Link(string) level(cilevel) /* */ CLUSter(varlist numeric min=1 max=1) reuse suppress GRaph(integer 0) /* */ YLAbel(string) XLAbel(string) XTick(numlist ascending) /* */ YTick(numlist ascending) copyleft * ] *NOTE: The reuse switch is a programming aid for use by msdthaz. This option *tells dthaz to calculate the specified probabilities, but use whatever the *most recent estimate was. In this way, computing time is reduced by avoiding *repeadted estimation of the model in msdthaz. This is especially helpful when *using the complementary log-log link. Reuse is NOT intended for user *application of dthaz. quietly { preserve ******************************************************************************* *Set up shop, prepare variables, confirm truncate, pretrunc, tpar, etc. * ******************************************************************************* * display the copyleft information if requested if "`copyleft'" == "copyleft" { noisily { di _newline "Copyright Notice" di "dthaz and dthaz.ado are Copyright (c) 2001, 2015 alexis dinno" _newline di "This file is part of dthaz." _newline di "dthaz is free software; you can redistribute it and/or modify" di "it under the terms of the GNU General Public License as published by" di "the Free Software Foundation; either version 2 of the License, or" di "(at your option) at any later version." _newline di "This program is distributed in the hope that it will be useful," di "but WITHOUT ANY WARRANTY; without even the implied warranty of" di "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" di "GNU General Public License for more details." _newline di "You should have received a copy of the GNU General Public License" di "along with this program (dthaz.copying); if not, write to the Free Software" di "Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA" _newline } } *Validate link function if "`link'"=="" { local link = "logit" } if (!("`link'"=="logit" | "`link'"=="cloglog" | "`link'"=="probit")) { noisily { di as err "invalid link(): `link'." _newline "valid link() options are logit, cloglog, or probit." } error (199) } *Validate for use of both truncate and pretrunc. truncate+pretrunc must *be < maxl (of the original dataset) minus 2 or there will be no periods *to analyse change for! sum _period if (`truncate' ~= 0 & `pretrunc' ~=0) { if ( (`truncate'-`pretrunc') <= 2 ) { noisily { di as err "truncate() and pretrunc() values are incompatible." di as err "not enough periods left to analyze: (`truncate'-`pretrunc') <= 2." } exit } } *Set up truncate appropriately; maxl is the MAXimum Length; mxm1 is the *MaXimum length Minus 1 sum _period if (`truncate'<1 | `truncate'>r(max)) { local maxl=r(max) if (`truncate' ~= 0) { noisily { di as err "Truncation value out of range." di as white "No truncation used." _newline } } } else { local maxl=`truncate' } local mxm1=`maxl'-1 local newobs=`maxl'+1 *Evaluate whether any of the time periods are entirely absent event occurence. *Terminate program if this is the case and notify user. forvalues x=1/`maxl' { count if _period==`x' & _Y==1 local events=r(N) count if _period==`x' & _Y==0 local noevents=r(N) if `events'==0 | `noevents'==0 { noisily { di as err "No variation in event occurence for period " `x' "!" di as err "(no events occured at this time period)" } error `events'==0 exit } } *Deal with pretrunc... make sure that pretrunc is either greater than or *equal to zero AND less than maxl-l, or do pre-truncation. if (`pretrunc'<0 | `pretrunc'>=`mxm1') { noisily { di as err "Pre-truncation value out of range." di as white "No pre-truncation used." _newline } } else { if (`pretrunc' ~= 0) { * Dropping the first `pretrunc' number of period indicator variables forvalues x=1/`pretrunc' { drop _d`x' } * Sequentially renaming the remaining period indicators. local maxl = `maxl'-`pretrunc' local mxm1 = `mxm1'-`pretrunc' local newobs = `newobs'-`pretrunc' forvalues x=1/`maxl' { local curper = `pretrunc'+`x' rename _d`curper' _d`x' } * Dropping observations in the pre-truncated periods drop if _period<=`pretrunc' * Revaluing _period so that the earliest time value is 1 (i.e. * subtracting the pre-truncation number from _period. replace _period = _period-`pretrunc' } } *Deal with tpar values. For those too low, slap on wrist and set to default. *For those overspecified, alert and lower to maximum order polynomial *representation of time. *Here's the slap on wrist if `tpar'<-2 { noisily { di as error "The specified value for time parameterization out of range." di as white "Time modeled as fully discreet." _newline } local tpar = -1 } *Here's the alert if `tpar'>=`maxl' { noisily { di as error "Polynomial parameterization of time overspecified. " di as error `mxm1' " is the highest order polynomial parameterization allowed for this dataset." di as white "Time parameterized as a polynomial of order " `mxm1' "." _newline } local tpar = `mxm1' } *Create an explicit constant term for later funky parameterizations of time if `tpar'~=-1 { generate _one=1 } *This line is here to facilitate the "macro shift" style yumminess later tokenize `varlist' ******************************************************************************* *Specify which time range to estimate for different parameterizations of time * ******************************************************************************* *For fully discrete time if `tpar'==-1 { local tvars = "_d1-_d"+"`maxl'" } *For constant effect of time if `tpar'==0 { local tvars="_one" } *Take care of n-polynomial time variables for n>=2 if `tpar'>=1 { *The local variable npoly is the order of the polynomial local npoly=`tpar' *_period_1 through _periodX are the predictors representing npoly time for num 1/`npoly': generate _period_X = _period^X *These predictors are arranged from highest to lowest order for num 1/`npoly': order _period_X *And the local variable tvars which describes which time-predictors to use *in the regression model is specified. local tvars="_period_"+"`npoly'"+"-_period_1 _one" local tvars = "`tvars'" } *And let's not forget root parameterization of time... if `tpar'==-2 { generate _periodrt=_period^.5 local tvars="_periodrt _period _one" } ******************************************************************************* *Prepare to tortuously create the specification matrix for predictors * *The matrix _Zero will be prepended to _Spec to create an NxN matrix * *containing the estimate specification values from specify(). * ******************************************************************************* if `tpar'==-1 { matrix _Zero=J(1,`maxl',0) } if `tpar'==0 { matrix _Zero=J(1,`maxl',0) } if `tpar'==-2 { local maxlp=2*`maxl' matrix _Zero = J(1,`maxlp',0) } *_Spec contains the variable values for the specified estimate matrix _One = (1) matrix input _Spec = (`specify') local argnum = colsof(_Spec) if `tpar'~=-1 { matrix _Spec = _One,_Spec } *End the quietly block so the following output is visible } ******************************************************************************* *Initial output. Tell the user what's in store, reiterate specification for * *estimate. At this point _Spec holds values corresponding to the tokenized * *variables in varlist. * ******************************************************************************* if "`suppress'"=="" { noisily { di _newline di as txt "Discrete-Time Estimation of Conditional Hazard and Survival Probabilities" di as txt "------------------------------------------------------------------------------" if `tpar'==-2 { di as txt "Time Parameterization: Square Root" } if `tpar'==-1 { di as txt "Time Parameterization: Fully Discrete" } if `tpar'==0 { di as txt "Time Parameterization: Constant Effect polynomial of order 0)" } if `tpar'==1 { di as txt "Time Parameterization: Linear (polynomial of order 1)" } if `tpar'==2 { di as txt "Time Parameterization: Quadratic (polynomial of order 2)" } if `tpar'==3 { di as txt "Time Parameterization: Cubic (polynomial of order 3)" } if `tpar'>=4 { di as txt "Polynomial Time Parameterization of Order " `tpar' } if "`specify'"=="" { di as txt "Baseline model (no additional predictors)" } if "`specify'"~="" { di _newline as txt "Additional predictors specified as:" } if `tpar'==-1 { forvalues x=1/`argnum' { di as txt "`1' = " as res el(_Spec,1,`x') macro shift } } if `tpar'~=-1 { local newarg = `argnum'+1 forvalues x=2/`newarg' { di as txt "`1' = " as res el(_Spec,1,`x') macro shift } } } } quietly { ******************************************************************************* *Now _Spec can be transmogrified into a form applicable to later estimate. * ******************************************************************************* if `tpar'==-1 | `tpar'==-2 { * if `tpar'<=-1 { matrix _Spec = _Zero,_Spec } *Make _Spec diagonal so it can multiply with _Q2 for hazard calculation matrix _Spec = diag(_Spec) ******************************************************************************* *Whip out the estimation... ensuring the appropriate predictors (if `varlist' * *is empty, then the numbered macros `1', `2', etc. contain the variable names * *from the dataset in order of appearance. Please see the note about the intent* *of the reuse switch at the start of this code. * ******************************************************************************* if "`reuse'"=="" { if "`link'"=="logit" { if "`specify'"=="" { logit _Y `tvars' `if' `in' [`weight' `exp'], nocons cluster(`cluster') } else { logit _Y `tvars' `varlist' `if' `in' [`weight' `exp'], nocons cluster(`cluster') } if "`model'"~="" { noisily logit } if "`suppress'"=="suppress" { exit } } if "`link'"=="cloglog" { if "`specify'"=="" { cloglog _Y `tvars' `if' `in' [`weight' `exp'], nocons cluster(`cluster') } else { cloglog _Y `tvars' `varlist' `if' `in' [`weight' `exp'], nocons cluster(`cluster') } if "`model'"~="" { noisily cloglog } if "`suppress'"=="suppress" { exit } } if "`link'"=="probit" { if "`specify'"=="" { probit _Y `tvars' `if' `in' [`weight' `exp'], nocons cluster(`cluster') } else { probit _Y `tvars' `varlist' `if' `in' [`weight' `exp'], nocons cluster(`cluster') } if "`model'"~="" { noisily probit } if "`suppress'"=="suppress" { exit } } *Clean up for constant effect time if `tpar'==0 { drop _one } *Clean up for n-polynomial time if `tpar'>=1 { drop `tvars' } *Clean up for root time if `tpar'==-2 { drop _one _periodrt } } ******************************************************************************* *At this point this program departs to highly convoluted realms. The basic * *strategy is to obtain a representation of estimates for the effect of time at* *each period of the dataset. How this is done varies depending on the * *parameterization of time, but usually involves creating a vector "_Q" holding* *estimates for each period. The second part of the strategy is the creation of* *a constant estimate term called "est" which contains the parameter estimates * *for the constant, pre-specified portion of the model. * ******************************************************************************* *It all starts with the estimates matrix _Q = e(b) ******************************************************************************* *Discrete-time estimates of time's effect: copy estimates from coeficient * *matrix. _Q will be used to describe the period estimates, while _Q2 will * *provide the estimates from the estimate's specification. * ******************************************************************************* if `tpar'==-1 { matrix _Q2 = diag(e(b)) } *for constant effect of time... if `tpar'==0 { matrix _Q = e(b) matrix _Q2 = diag(_Q) } ******************************************************************************* *The following code generates a number of matrices equal to the order of the * *polynomial time parameterization. Each matrix is named PolyTimeXO where X is * *the power that the value of period is raised to. These values are then * *aggregated for each time period in the matrix _ParameterizedTime which will * *be used in the conditional hazard probability calculation. * ******************************************************************************* if `tpar'>=1 { *This loop creates the separate matrix for each exponentiation of time forvalues npolynum=1/`npoly' { local no = 1+`npoly'-`npolynum' for num `no': matrix _PolyTimeXO = J(1,`maxl',0) *This loop creates the values within the exponentiation for each period forvalues per=1/`maxl' { *Replace the placeholder in each period position of PolyTimeXO with the *appropriate exponentiation of period times the parameter estimate. matrix _PolyTime`no'O[1,`per']=(`per'^`no')*el(_Q,1,`npolynum') } } *Get collapse the previously arranged matrices into one and loose 'em... matrix _ParameterizedTime = J(1,`maxl',0) forvalues order=1/`npoly' { matrix _ParameterizedTime = _ParameterizedTime + _PolyTime`order'O matrix drop _PolyTime`order'O } *Create the matrix _Q2 which holds the model specifications for the constant * *term _one plus any predictors in the model. local specs = `npoly'+1 *Do so by copying the estimates _after_ the terms representing time matrix _Q2 = _Q[1,`specs'...] matrix _Q2 = diag(_Q2) } ******************************************************************************* *For root parameterizations of time Matrix _Root holds the estimates for * *linearly parameterized time for each period. _Toor is _Root reversed... If * *this looks wierd, you ought to have seen the code before I wrote the above * *bit of prettiness for n-polynomial representation... * ******************************************************************************* if `tpar'==-2 { *Create linear time vector forvalues x=2/`maxl' { local place = (`x')*el(_Q,1,2) if `x'==2 { matrix _Linear = (`place') matrix _Raenil = _Linear } else { matrix _Linear = (`place'),_Linear matrix _Raenil = _Raenil,(`place') } } local place = el(_Q,1,2) matrix _Raenil = (`place'),_Raenil *Create root time vector forvalues x=1/`maxl' { local place = ((`x')^.5)*el(_Q,1,1) if `x'==1 { matrix _Root = (`place') matrix _Toor = _Root } else { matrix _Root = (`place'),_Root matrix _Toor = _Toor,(`place') } } *At this point _Q is B_psq B_p B_one B_predictors... We need to get rid of the *B_sq term, so: local bark = (`argnum'+3) forvalues x=2/`bark' { local sploof = el(_Q,1,`x') if `x'==2 { matrix _Blar = (`sploof') } else { matrix _Blar = _Blar,(`sploof') } } matrix _Q = _Blar matrix _Q2 = _Root,_Linear,_Q matrix _Q2 = diag(_Q2) matrix drop _Blar _Linear _Root } ******************************************************************************* *_Q2 contains the parameter*specified variable values for the estimate at hand* ******************************************************************************* matrix _Q2 = _Q2*_Spec ******************************************************************************* *est contains the actual quantity employed in adjusting the calculation of * *hazard. The wierdness in producing the estimates for root-represented time is* *apparent in the need to clean up for it especially here. * ******************************************************************************* if `tpar'>=-1 { local est = trace(_Q2) } if `tpar'==-2 { local est = trace(_Q2) matrix _Q = _Toor,_Raenil matrix drop _Toor _Raenil } ******************************************************************************* *Into the home stretch. Matrix _Hazard will be created and populated using * *some variation of the basic algorithm. * ******************************************************************************* *Generate hazard matrix, with a non-calculated 0th period matrix _Hazard = (0, 0) *The following bit of flow-control limits the number of periods displayed if *the user asked for it with the display option if (`display'>0 & `display'<`maxl') { local maxl=`display' local newobs=`maxl'+1 } *Set the number of coeficients from the estimate matrix and specification matrix *Hazard probabilities forvalues num=1/`maxl' { *for fully discrete and _Linear time if `tpar'==-1 { local specest = el(_Q,1,`num')+`est' } *for constant effect of time if `tpar'==0 { local specest = `est' } *for root time if `tpar'==-2 { local lin = `num'+`maxl' local specest = el(_Q,1,`num')+el(_Q,1,`lin')+`est' } *for polynomial time if `tpar'>=1 { local specest = el(_ParameterizedTime,1,`num')+`est' } if "`link'"=="logit" { local haz = 1/(1+(exp(-1*(`specest')))) } if "`link'"=="cloglog" { local haz = 1-(exp(-1*exp(`specest'))) } if "`link'"=="probit" { local haz = normal(`specest') } matrix _Period = (`num',`haz') matrix _Hazard = _Hazard\_Period } *Produce survival probabilities and append to _Hazard forvalues num=1/`newobs' { if `num'==1 { local lastsur = 1 matrix _Survival = (1) } else { local back = `num'-1 local lastsur = el(_Survival,`back',1) } local haz = el(_Hazard,`num',2) local sur = (1-`haz')*(`lastsur') matrix _Period = (`sur') if `num'~=1 { matrix _Survival = _Survival\_Period } } *Produce standard errors for h_t and S_t local q = e(df_m) matrix _sigma_h = vecdiag(I(`maxl')) matrix _sigma_Sa = vecdiag(I(`maxl')) matrix _sigma_Sb = vecdiag(I(`maxl')) matrix _sigma_S = vecdiag(I(`maxl')) matrix betas = e(b) matrix V = e(V) forvalues t=1/`maxl' { if `tpar' == -2 { matrix Z = (`t'^.5, `t', 1) } if `tpar' == -1 { matrix Z = I(`maxl') matrix Z = Z[`t',1..`maxl'] } if `tpar' == 0 { matrix Z = (1) } if `tpar' >0 { matrix Z = (1) forvalues polynomial=1/`tpar' { matrix Z = `t'^`polynomial',Z } } if "`specify'" ~= "" { tokenize `specify' foreach specification in `*' { matrix Z = Z,(`specification') } } local lZ = colsof(betas) matrix bZ = Z forvalues i=1/`lZ' { matrix bZ[1,`i'] = betas[1,`i']*Z[1,`i'] } if "`link'"=="logit" { matrix _sigma_h[1,`t'] = ( trace(( ((exp( trace(diag(bZ)) ))/((1+exp( trace(diag(bZ)) ))^2))^2 * Z*V*(Z')) )^.5) matrix _sigma_Sa[1,`t'] = ( trace( ((exp( trace(diag(bZ)) ))/(1+exp( trace(diag(bZ)) )))^2 * Z*V*(Z')) ) matrix _sigma_Sb[1,`t'] = (1/(1 + exp( trace(diag(bZ)) ))) } if "`link'"=="cloglog" { matrix _sigma_h[1,`t'] = ( trace(( ( exp( trace(diag(bZ)) - exp(trace(diag(bZ))) ) )^2 * Z*V*(Z')) )^.5) matrix _sigma_Sa[1,`t'] = ( trace( ( exp(trace(diag(bZ)) ) )^2 * Z*V*(Z')) ) matrix _sigma_Sb[1,`t'] = exp( -exp( trace(diag(bZ)) ) ) } if "`link'"=="probit" { matrix _sigma_h[1,`t'] = ( trace(( ( (1/sqrt(2*_pi))*exp((-(trace(diag(bZ))^2))/(2)) )^2 * Z*V*(Z')) )^.5) matrix _sigma_Sa[1,`t'] = ( trace( ( ((1/sqrt(2*_pi))*exp((-(trace(diag(bZ))^2))/(2))) / (1-normal(trace(diag(bZ)))) )^2 * Z*V*(Z')) ) matrix _sigma_Sb[1,`t'] = 1 - normal(trace(diag(bZ))) } if `t' > 1 { matrix _sigma_Sa[1,`t'] = (_sigma_Sa[1,`t'] + _sigma_Sa[1,(`t'-1)]) matrix _sigma_Sb[1,`t'] = (_sigma_Sb[1,`t'] * _sigma_Sb[1,(`t'-1)]) } matrix _sigma_S[1,`t'] = ( (_sigma_Sa[1,`t']*((_sigma_Sb[1,`t'])^2))^.5 ) } matrix _sigma_h = (0),_sigma_h matrix _sigma_S = (0),_sigma_S matrix _Hazard = _Hazard,_Survival,(_sigma_h'),(_sigma_S') matname _Hazard Period Hazard Survival sigma_h sigma_S, columns(1...) explicit } ******************************************************************************* *Final bit of output: display results to user * ******************************************************************************* if "`suppress'"=="" { noisily { di _newline di as txt "-------------------------------------------------------------" di as txt _col(4) "Period" _col(16) "p(Hazard)" _col(28) "Std. Err. " _col(39) "p(Survival)" _col(52) "Std. Err. " di as txt _col(4) " (t)" _col(16) " ^h(t)" _col(28) " ^h(t)" _col(39) " ^S(t)" _col(52) " ^S(t)" di as txt "-------------------------------------------------------------" di as result _col(6) "0" _col(19) "--" _col(31) "--" _col(39) " 1" _col(52) " 0" as result forvalues i=2/`newobs' { display as result _col(6) el(_Hazard,`i',1) _col(16) %10.7f el(_Hazard,`i',2) _col(28) %10.7f el(_Hazard,`i',4) _col(39) %10.7f el(_Hazard,`i',3) _col(52) %10.7f el(_Hazard,`i',5) as result } di as txt "-------------------------------------------------------------" *Note estimation assumptions if "`link'"=="logit" { di as txt "Logit Link (assumes proportional odds)" } if "`link'"=="cloglog" { di as txt "Complementary Log-Log Link (assumes proportional hazards)" } if "`link'"=="probit" { di as txt "Probit Link (assumes proportional probits)" } } } quietly { ******************************************************************************* *Clean this mess on up! * ******************************************************************************* *Drop that which must be dropped matrix drop _One _Period _Q _Q2 _Spec _Survival Z bZ if `tpar'==-2 | `tpar'==-1 { matrix drop _Zero } if `tpar'>=1 { matrix drop _ParameterizedTime } ******************************************************************************* *Graph some output! * ******************************************************************************* restore, preserve svmat _Hazard, names(col) *Generate cumulative incidence gen InvSurvival = 1 - Survival *Generate upper and lower bounds for Hazard gen Hub = Hazard+(invnormal(1-((1-(`level'/100))/2))*sigma_h) if Period > 0 replace Hub = 1 if Hub > 1 & Period > 0 gen Hlb = Hazard-(invnormal(1-((1-(`level'/100))/2))*sigma_h) if Period > 0 replace Hlb = 0 if Hlb < 0 & Period > 0 *Generate upper and lower bounds for Survival gen Sub = Survival+(invnormal(1-((1-(`level'/100))/2))*sigma_S) if Period > 0 replace Sub = 1 if Sub > 1 & Period > 0 gen Slb = Survival-(invnormal(1-((1-(`level'/100))/2))*sigma_S) if Period > 0 replace Slb = 0 if Slb < 0 & Period > 0 *Generate upper and lower bounds for cumulative incidence gen Cub = InvSurvival+(invnormal(1-((1-(`level'/100))/2))*sigma_S) if Period > 0 replace Cub = 1 if Cub > 1 & Period > 0 gen Clb = InvSurvival-(invnormal(1-((1-(`level'/100))/2))*sigma_S) if Period > 0 replace Clb = 0 if Clb < 0 & Period > 0 set more on *Manage graphing options input for all graphs if `graph'> 4 | `graph'<1 { local `graph'=0 } if `"`xtitle'"'==`""' { local b1title = `"b1title("Period")"' local xtitle = `"xtitle("Period")"' } else { local b1title = `"b1title(`xtitle')"' local xtitle = `"xtitle(`xtitle')"' } if "`ytick'"=="" { local ytick "yt(0(.2)1)" } else { local ytick "yt(`ytick')" } if "`ylabel'"=="" { local ylabel "yla(0(.2)1)" } else { local ylabel "ylabel(`ylabel')" } *Graph conditional hazard probabilities if `graph'==1 { *Deal with axes specific to conditional hazard curves if "`xtick'"=="" { local xtick "xt(1(1)`maxl')" } else { local xtick "xs(`xtick')" } if "`xlabel'"=="" { local xlabel "xla(1(1)`maxl')" } else { local xlabel "xlabel(`xlabel')" } if `"`ytitle'"'==`""' { local ytitle = `"ytitle("Estimated Conditional Hazard Probability")"' } else { local ytitle = `"ytitle(`ytitle')"' } graph twoway line Hazard Period if Period~=0, `ytick' `ylabel' `ytitle' `xtick' `xlabel' `xtitle' title("Hazard Curve") subtitle("Conditional hazard probability vs. period") `options' legend(off) || line Hub Period, lwidth(vvthin) lcolor(gs8) || line Hlb Period, lwidth(vvthin) lcolor(gs8) } *Graph survival probabilities if `graph'==2 | `graph'==4 { *Deal with axes specific to survival curves if "`xtick'"=="" { local xtick = "xt(0(1)`maxl')" } else { local xtick = "xt(`xtick')" } if "`xlabel'"=="" { local xlabel = "xla(0(1)`maxl')" } else { local xlabel = "xlabel(`xlabel')" } if `"`ytitle'"'==`""' { if `graph' == 2 { local ytitle = `"ytitle("Estimated Survival Probability")"' } if `graph' == 4 { local ytitle = `"ytitle("Estimated Cumulative Incidence Probability")"' } } else { local ytitle = `"ytitle(`ytitle')"' } if `graph' == 2 { line Survival Period, `ytick' `ylabel' `ytitle' `xtick' `xlabel' `xtitle' title("Survival Curve") subtitle("survival probability vs. period") `options' legend(off) || line Sub Period, lwidth(vvthin) lcolor(gs8) || line Slb Period, lwidth(vvthin) lcolor(gs8) } if `graph' == 4 { line InvSurvival Period, `ytick' `ylabel' `ytitle' `xtick' `xlabel' `xtitle' title("Cumulative Incidence Curve") subtitle("cumulative incidence probability vs. period") `options' legend(off) || line Cub Period, lwidth(vvthin) lcolor(gs8) || line Clb Period, lwidth(vvthin) lcolor(gs8) } } *Graph both conditional hazard and survival probabilities if `graph'==3 { *Prepare for different x-axis labeling needs for hazard and survival local tempxt = "`xtick'" local tempxla = "`xlabel'" local tempytitle = "`ytitle'" if "`xtick'"=="" { local xtick = "xt(1(1)`maxl')" } else { local xtick = "xt(`xtick')" } if "`xlabel'"=="" { local xlabel = "xla(1(1)`maxl')" } else { local xlabel = "xlabel(`xlabel')" } if "`ytitle'"=="" { local ytitle = `"ytitle("Estimated Conditional Hazard Probability")"' } else { local ytitle = "ytitle(`ytitle')" } line Hazard Period if Period~=0, `ytick' `ylabel' `ytitle' `xtick' `xlabel' `xtitle' title("Hazard Curve") subtitle("conditional hazard probability vs. period") `options' legend(off) || line Hub Period, lwidth(vvthin) lcolor(gs8) || line Hlb Period, lwidth(vvthin) lcolor(gs8) more local xtick = "`tempxt'" local xlabel = "`tempxla'" local ytitle = "`tempytitle'" if "`xtick'"=="" { local xtick = "xt(0(1)`maxl')" } else { local xtick = "xt(`xtick')" } if "`xlabel'"=="" { local xlabel = "xla(0(1)`maxl')" } else { local xlabel = "xlabel(`xlabel')" } if `"`ytitle'"'==`""' { local l1title = `"l1title("Estimated Survival Probability")"' local ytitle = `"ytitle("Estimated Survival Probability")"' } else { local l1title = `"l1title(`ytitle')"' local ytitle = `"ytitle(`ytitle')"' } line Survival Period, `ytick' `ylabel' `ytitle' `xtick' `xlabel' `xtitle' title("Survival Curve") subtitle("estimated survival probability vs. period") `options' legend(off) || line Sub Period, lwidth(vvthin) lcolor(gs8) || line Slb Period, lwidth(vvthin) lcolor(gs8) } restore, preserve ******************************************************************************* *Say goodbye to the nice user! * ******************************************************************************* mata: st_matrix("Hazard", st_matrix("_Hazard")[.,2]') mata: st_matrix("HazardSE", st_matrix("_Hazard")[.,4]') mata: st_matrix("Survival", st_matrix("_Hazard")[.,3]') mata: st_matrix("SurvivalSE", st_matrix("_Hazard")[.,5]') * clean up matrix drop _Hazard _sigma_S _sigma_h V betas _sigma_Sb _sigma_Sa *Retain _Hazard for user... ereturn matrix SurvivalSE = SurvivalSE ereturn matrix Survival = Survival ereturn matrix HazardSE = HazardSE ereturn matrix Hazard = Hazard return clear } end ******************************************************************************* ******************************************************************************* * dthaz for STATA v7 ******************************************************************************* ******************************************************************************* program define dthaz7, rclass version 7.0 syntax [varlist] [if] [in] [fweight pweight iweight] [, SPecify(numlist) /* */ TPar(integer -1) Pretrunc(integer 0) Display(integer 0) /* */ Model Link(string) level(cilevel) Truncate(integer 0) /* */ CLUSter(varlist numeric min=1 max=1) reuse YLAbel(string) /* */ XLAbel(string) XTick(numlist ascending) suppress GRaph(integer 0) /* */ YTick(numlist ascending) copyleft * ] *NOTE: The reuse switch is a programming aid for use by msdthaz. This option *tells dthaz to calculate the specified probabilities, but use whatever the *most recent estimate was. In this way, computing time is reduced by avoiding *repeadted estimation of the model in msdthaz. This is especially helpful when *using the complementary log-log link. Reuse is NOT intended for user *application of dthaz. quietly { preserve ******************************************************************************* *Set up shop, prepare variables, confirm truncate, pretrunc, tpar, etc. * ******************************************************************************* * display the copyleft information if requested if "`copyleft'" == "copyleft" { noisily { di _newline "Copyright Notice" di "dthaz and dthaz.ado are Copyright (c) 2001, 2015 alexis dinno" _newline di "This file is part of dthaz." _newline di "dthaz is free software; you can redistribute it and/or modify" di "it under the terms of the GNU General Public License as published by" di "the Free Software Foundation; either version 2 of the License, or" di "(at your option) at any later version." _newline di "This program is distributed in the hope that it will be useful," di "but WITHOUT ANY WARRANTY; without even the implied warranty of" di "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" di "GNU General Public License for more details." _newline di "You should have received a copy of the GNU General Public License" di "along with this program (dthaz.copying); if not, write to the Free Software" di "Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA" _newline } } *Confirm the validity of cluster (for v7) if "`cluster'" ~= "" { local cluster = "cluster(`cluster')" } *Validate link function if "`link'"=="" { local link = "logit" } if (!("`link'"=="logit" | "`link'"=="cloglog" | "`link'"=="probit")) { noisily { di as err "invalid link(): `link'." _newline "valid link() options are logit, cloglog, or probit." } error (199) } *Validate for use of both truncate and pretrunc. truncate+pretrunc must *be < maxl (of the original dataset) minus 2 or there will be no periods *to analyse change for! sum _period if (`truncate' ~= 0 & `pretrunc' ~=0) { if ( (`truncate'-`pretrunc') <= 2 ) { noisily { di as err "truncate() and pretrunc() values are incompatible." di as err "not enough periods left to analyze: (`truncate'-`pretrunc') <= 2." } exit } } *Set up truncate appropriately; maxl is the MAXimum Length; mxm1 is the *MaXimum length Minus 1 sum _period if (`truncate'<1 | `truncate'>r(max)) { local maxl=r(max) if (`truncate' ~= 0) { noisily { di as err "Truncation value out of range." di as white "No truncation used." _newline } } } else { local maxl=`truncate' } local mxm1=`maxl'-1 local newobs=`maxl'+1 *Evaluate whether any of the time periods are entirely absent event occurence. *Terminate program if this is the case and notify user. forvalues x=1/`maxl' { count if _period==`x' & _Y==1 local events=r(N) count if _period==`x' & _Y==0 local noevents=r(N) if `events'==0 | `noevents'==0 { noisily { di as err "No variation in event occurence for period " `x' "!" di as err "(no events occured at this time period)" } error `events'==0 exit } } *Deal with pretrunc... make sure that pretrunc is either greater than or *equal to zero AND less than maxl-l, or do pre-truncation. if (`pretrunc'<0 | `pretrunc'>=`mxm1') { noisily { di as err "Pre-truncation value out of range." di as white "No pre-truncation used." _newline } } else { if (`pretrunc' ~= 0) { * Dropping the first `pretrunc' number of period indicator variables forvalues x=1/`pretrunc' { drop _d`x' } * Sequentially renaming the remaining period indicators. local maxl = `maxl'-`pretrunc' local mxm1 = `mxm1'-`pretrunc' local newobs = `newobs'-`pretrunc' forvalues x=1/`maxl' { local curper = `pretrunc'+`x' rename _d`curper' _d`x' } * Dropping observations in the pre-truncated periods drop if _period<=`pretrunc' * Revaluing _period so that the earliest time value is 1 (i.e. * subtracting the pre-truncation number from _period. replace _period = _period-`pretrunc' } } *Deal with tpar values. For those too low, slap on wrist and set to default. *For those overspecified, alert and lower to maximum order polynomial *representation of time. *Here's the slap on wrist if `tpar'<-2 { noisily { di as error "The specified value for time parameterization out of range." di as white "Time modeled as fully discreet." _newline } local tpar = -1 } *Here's the alert if `tpar'>=`maxl' { noisily { di as error "Polynomial parameterization of time overspecified. " di as error `mxm1' " is the highest order polynomial parameterization allowed for this dataset." di as white "Time parameterized as a polynomial of order " `mxm1' "." _newline } local tpar = `mxm1' } *Create an explicit constant term for later funky parameterizations of time if `tpar'~=-1 { generate _one=1 } *This line is here to facilitate the "macro shift" style yumminess later tokenize `varlist' ******************************************************************************* *Specify which time range to estimate for different parameterizations of time * ******************************************************************************* *For fully discrete time if `tpar'==-1 { local tvars = "_d1-_d"+"`maxl'" } *For constant effect of time if `tpar'==0 { local tvars="_one" } *Take care of n-polynomial time variables for n>=2 if `tpar'>=1 { *The local variable npoly is the order of the polynomial local npoly=`tpar' *_period_1 through _periodX are the predictors representing npoly time for num 1/`npoly': generate _period_X = _period^X *These predictors are arranged from highest to lowest order for num 1/`npoly': order _period_X *And the local variable tvars which describes which time-predictors to use *in the regression model is specified. local tvars="_period_"+"`npoly'"+"-_period_1 _one" local tvars = "`tvars'" } *And let's not forget root parameterization of time... if `tpar'==-2 { generate _periodrt=_period^.5 local tvars="_periodrt _period _one" } ******************************************************************************* *Prepare to tortuously create the specification matrix for predictors * *The matrix _Zero will be prepended to _Spec to create an NxN matrix * *containing the estimate specification values from specify(). * ******************************************************************************* if `tpar'==-1 { matrix _Zero=J(1,`maxl',0) } if `tpar'==0 { matrix _Zero=J(1,`maxl',0) } if `tpar'==-2 { local maxlp=2*`maxl' matrix _Zero = J(1,`maxlp',0) } *_Spec contains the variable values for the specified estimate matrix _One = (1) matrix input _Spec = (`specify') local argnum = colsof(_Spec) if `tpar'~=-1 { matrix _Spec = _One,_Spec } *End the quietly block so the following output is visible } ******************************************************************************* *Initial output. Tell the user what's in store, reiterate specification for * *estimate. At this point _Spec holds values corresponding to the tokenized * *variables in varlist. * ******************************************************************************* if "`suppress'"=="" { noisily { di _newline di as txt "Discrete-Time Estimation of Conditional Hazard and Survival Probabilities" di as txt "------------------------------------------------------------------------------" if `tpar'==-2 { di as txt "Time Parameterization: Square Root" } if `tpar'==-1 { di as txt "Time Parameterization: Fully Discrete" } if `tpar'==0 { di as txt "Time Parameterization: Constant Effect polynomial of order 0)" } if `tpar'==1 { di as txt "Time Parameterization: Linear (polynomial of order 1)" } if `tpar'==2 { di as txt "Time Parameterization: Quadratic (polynomial of order 2)" } if `tpar'==3 { di as txt "Time Parameterization: Cubic (polynomial of order 3)" } if `tpar'>=4 { di as txt "Polynomial Time Parameterization of Order " `tpar' } if "`specify'"=="" { di as txt "Baseline model (no additional predictors)" } if "`specify'"~="" { di _newline as txt "Additional predictors specified as:" } if `tpar'==-1 { forvalues x=1/`argnum' { di as txt "`1' = " as result el(_Spec,1,`x') macro shift } } if `tpar'~=-1 { local newarg = `argnum'+1 forvalues x=2/`newarg' { di as txt "`1' = " as result el(_Spec,1,`x') macro shift } } } } quietly { ******************************************************************************* *Now _Spec can be transmogrified into a form applicable to later estimate. * ******************************************************************************* if `tpar'==-1 | `tpar'==-2 { * if `tpar'<=-1 { matrix _Spec = _Zero,_Spec } *Make _Spec diagonal so it can multiply with _Q2 for hazard calculation matrix _Spec = diag(_Spec) ******************************************************************************* *Whip out the estimation... ensuring the appropriate predictors (if `varlist' * *is empty, then the numbered macros `1', `2', etc. contain the variable names * *from the dataset in order of appearance. Please see the note about the intent* *of the reuse switch at the start of this code. * ******************************************************************************* if "`reuse'"=="" { if "`link'"=="cloglog" { if "`link'"=="logit" { if "`specify'"=="" { logit _Y `tvars' `if' `in' [`weight' `exp'], nocons `cluster' } else { logit _Y `tvars' `varlist' `if' `in' [`weight' `exp'], nocons `cluster' } if "`model'"~="" { noisily logit } if "`suppress'"=="suppress" { exit } } if "`specify'"=="" { cloglog _Y `tvars' `if' `in' [`weight' `exp'], nocons `cluster' } else { cloglog _Y `tvars' `varlist' `if' `in' [`weight' `exp'], nocons `cluster' } if "`model'"~="" { noisily cloglog } if "`suppress'"=="suppress" { exit } } if "`link'"=="probit" { if "`specify'"=="" { probit _Y `tvars' `if' `in' [`weight' `exp'], nocons `cluster' } else { probit _Y `tvars' `varlist' `if' `in' [`weight' `exp'], nocons `cluster' } if "`model'"~="" { noisily probit } if "`suppress'"=="suppress" { exit } } *Clean up for constant effect time if `tpar'==0 { drop _one } *Clean up for n-polynomial time if `tpar'>=1 { drop `tvars' } *Clean up for root time if `tpar'==-2 { drop _one _periodrt } } ******************************************************************************* *At this point this program departs to highly convoluted realms. The basic * *strategy is to obtain a representation of estimates for the effect of time at* *each period of the dataset. How this is done varies depending on the * *parameterization of time, but usually involves creating a vector "_Q" holding* *estimates for each period. The second part of the strategy is the creation of* *a constant estimate term called "est" which contains the parameter estimates * *for the constant, pre-specified portion of the model. * ******************************************************************************* *It all starts with the estimates matrix _Q = e(b) ******************************************************************************* *Discrete-time estimates of time's effect: copy estimates from coeficient * *matrix. _Q will be used to describe the period estimates, while _Q2 will * *provide the estimates from the estimate's specification. * ******************************************************************************* if `tpar'==-1 { matrix _Q2 = diag(e(b)) } *for constant effect of time... if `tpar'==0 { matrix _Q = e(b) matrix _Q2 = diag(_Q) } ******************************************************************************* *The following code generates a number of matrices equal to the order of the * *polynomial time parameterization. Each matrix is named PolyTimeXO where X is * *the power that the value of period is raised to. These values are then * *aggregated for each time period in the matrix _ParameterizedTime which will * *be used in the conditional hazard probability calculation. * ******************************************************************************* if `tpar'>=1 { *This loop creates the separate matrix for each exponentiation of time forvalues npolynum=1/`npoly' { local no = 1+`npoly'-`npolynum' for num `no': matrix _PolyTimeXO = J(1,`maxl',0) *This loop creates the values within the exponentiation for each period forvalues per=1/`maxl' { *Replace the placeholder in each period position of PolyTimeXO with the *appropriate exponentiation of period times the parameter estimate. matrix _PolyTime`no'O[1,`per']=(`per'^`no')*el(_Q,1,`npolynum') } } *Get collapse the previously arranged matrices into one and loose 'em... matrix _ParameterizedTime = J(1,`maxl',0) forvalues order=1/`npoly' { matrix _ParameterizedTime = _ParameterizedTime + _PolyTime`order'O matrix drop _PolyTime`order'O } *Create the matrix _Q2 which holds the model specifications for the constant * *term _one plus any predictors in the model. local specs = `npoly'+1 *Do so by copying the estimates _after_ the terms representing time matrix _Q2 = _Q[1,`specs'...] matrix _Q2 = diag(_Q2) } ******************************************************************************* *For root parameterizations of time Matrix _Root holds the estimates for * *linearly parameterized time for each period. _Toor is _Root reversed... If * *this looks wierd, you ought to have seen the code before I wrote the above * *bit of prettiness for n-polynomial representation... * ******************************************************************************* if `tpar'==-2 { *Create linear time vector forvalues x=2/`maxl' { local place = (`x')*el(_Q,1,2) if `x'==2 { matrix _Linear = (`place') matrix _Raenil = _Linear } else { matrix _Linear = (`place'),_Linear matrix _Raenil = _Raenil,(`place') } } local place = el(_Q,1,2) matrix _Raenil = (`place'),_Raenil *Create root time vector forvalues x=1/`maxl' { local place = ((`x')^.5)*el(_Q,1,1) if `x'==1 { matrix _Root = (`place') matrix _Toor = _Root } else { matrix _Root = (`place'),_Root matrix _Toor = _Toor,(`place') } } *At this point _Q is B_psq B_p B_one B_predictors... We need to get rid of the *B_sq term, so: local bark = (`argnum'+3) forvalues x=2/`bark' { local sploof = el(_Q,1,`x') if `x'==2 { matrix _Blar = (`sploof') } else { matrix _Blar = _Blar,(`sploof') } } matrix _Q = _Blar matrix _Q2 = _Root,_Linear,_Q matrix _Q2 = diag(_Q2) matrix drop _Blar _Linear _Root } ******************************************************************************* *_Q2 contains the parameter*specified variable values for the estimate at hand* ******************************************************************************* matrix _Q2 = _Q2*_Spec ******************************************************************************* *est contains the actual quantity employed in adjusting the calculation of * *hazard. The wierdness in producing the estimates for root-represented time is* *apparent in the need to clean up for it especially here. * ******************************************************************************* if `tpar'>=-1 { local est = trace(_Q2) } if `tpar'==-2 { local est = trace(_Q2) matrix _Q = _Toor,_Raenil matrix drop _Toor _Raenil } ******************************************************************************* *Into the home stretch. Matrix _Hazard will be created and populated using * *some variation of the basic algorithm. * ******************************************************************************* *Generate hazard matrix, with a non-calculated 0th period matrix _Hazard = (0, 0) *The following bit of flow-control limits the number of periods displayed if *the user asked for it with the display option if (`display'>0 & `display'<`maxl') { local maxl=`display' local newobs=`maxl'+1 } *Set the number of coeficients from the estimate matrix and specification matrix *Hazard probabilities forvalues num=1/`maxl' { *for fully discrete and _Linear time if `tpar'==-1 { local specest = el(_Q,1,`num')+`est' } *for constant effect of time if `tpar'==0 { local specest = `est' } *for root time if `tpar'==-2 { local lin = `num'+`maxl' local specest = el(_Q,1,`num')+el(_Q,1,`lin')+`est' } *for polynomial time if `tpar'>=1 { local specest = el(_ParameterizedTime,1,`num')+`est' } if "`link'"=="logit" { local haz = 1/(1+(exp(-1*(`specest')))) } if "`link'"=="cloglog" { local haz = 1-(exp(-1*exp(`specest'))) } if "`link'"=="probit" { local haz = normal(`specest') } matrix _Period = (`num',`haz') matrix _Hazard = _Hazard\_Period } *Produce survival probabilities and append to _Hazard forvalues num=1/`newobs' { if `num'==1 { local lastsur = 1 matrix _Survival = (1) } else { local back = `num'-1 local lastsur = el(_Survival,`back',1) } local haz = el(_Hazard,`num',2) local sur = (1-`haz')*(`lastsur') matrix _Period = (`sur') if `num'==1 { } else { matrix _Survival = _Survival\_Period } } *Produce standard errors for h_t and S_t local q = e(df_m) matrix _sigma_h = vecdiag(I(`maxl')) matrix _sigma_Sa = vecdiag(I(`maxl')) matrix _sigma_Sb = vecdiag(I(`maxl')) matrix _sigma_S = vecdiag(I(`maxl')) matrix betas = e(b) matrix V = e(V) forvalues t=1/`maxl' { if `tpar' == -2 { matrix Z = (`t'^.5, `t', 1) } if `tpar' == -1 { matrix Z = I(`maxl') matrix Z = Z[`t',1..`maxl'] } if `tpar' == 0 { matrix Z = (1) } if `tpar' >0 { matrix Z = (1) forvalues polynomial=1/`tpar' { matrix Z = `t'^`polynomial',Z } } if "`specify'" ~= "" { tokenize `specify' foreach specification in `*' { matrix Z = Z,(`specification') } } local lZ = colsof(betas) matrix bZ = Z forvalues i=1/`lZ' { matrix bZ[1,`i'] = betas[1,`i']*Z[1,`i'] } if "`link'"=="logit" { matrix _sigma_h[1,`t'] = ( trace(( ((exp( trace(diag(bZ)) ))/((1+exp( trace(diag(bZ)) ))^2))^2 * Z*V*(Z')) )^.5) matrix _sigma_Sa[1,`t'] = ( trace( ((exp( trace(diag(bZ)) ))/(1+exp( trace(diag(bZ)) )))^2 * Z*V*(Z')) ) matrix _sigma_Sb[1,`t'] = (1/(1 + exp( trace(diag(bZ)) ))) } if "`link'"=="cloglog" { matrix _sigma_h[1,`t'] = ( trace(( ( exp( trace(diag(bZ)) - exp(trace(diag(bZ))) ) )^2 * Z*V*(Z')) )^.5) matrix _sigma_Sa[1,`t'] = ( trace( ( exp(trace(diag(bZ)) ) )^2 * Z*V*(Z')) ) matrix _sigma_Sb[1,`t'] = exp( -exp( trace(diag(bZ)) ) ) } if "`link'"=="probit" { matrix _sigma_h[1,`t'] = ( trace(( ( (1/sqrt(2*_pi))*exp((-(trace(diag(bZ))^2))/(2)) )^2 * Z*V*(Z')) )^.5) matrix _sigma_Sa[1,`t'] = ( trace( ( ((1/sqrt(2*_pi))*exp((-(trace(diag(bZ))^2))/(2))) / (1-normal(trace(diag(bZ)))) )^2 * Z*V*(Z')) ) matrix _sigma_Sb[1,`t'] = 1 - normal(trace(diag(bZ))) } if `t' > 1 { matrix _sigma_Sa[1,`t'] = (_sigma_Sa[1,`t'] + _sigma_Sa[1,(`t'-1)]) matrix _sigma_Sb[1,`t'] = (_sigma_Sb[1,`t'] * _sigma_Sb[1,(`t'-1)]) } matrix _sigma_S[1,`t'] = ( (_sigma_Sa[1,`t']*((_sigma_Sb[1,`t'])^2))^.5 ) } matrix _sigma_h = (0),_sigma_h matrix _sigma_S = (0),_sigma_S matrix _Hazard = _Hazard,_Survival,(_sigma_h'),(_sigma_S') matname _Hazard Period Hazard Survival sigma_h sigma_S, columns(1...) explicit } ******************************************************************************* *Final bit of output: display results to user * ******************************************************************************* if "`suppress'"=="" { noisily { di _newline di as txt "-------------------------------------------------------------" di as txt _col(4) "Period" _col(16) "p(Hazard)" _col(28) "Std. Err. " _col(39) "p(Survival)" _col(52) "Std. Err. " di as txt _col(4) " (t)" _col(16) " ^h(t)" _col(28) " ^h(t)" _col(39) " ^S(t)" _col(52) " ^S(t)" di as txt "-------------------------------------------------------------" di as result _col(6) "0" _col(19) "--" _col(31) "--" _col(39) " 1" _col(52) " 0" as result forvalues i=2/`newobs' { display as result _col(6) el(_Hazard,`i',1) _col(16) %10.7f el(_Hazard,`i',2) _col(28) %10.7f el(_Hazard,`i',4) _col(39) %10.7f el(_Hazard,`i',3) _col(52) %10.7f el(_Hazard,`i',5) as result } di as txt "-------------------------------------------------------------" *Note estimation assumptions if "`link'"=="logit" { di as txt "Logit Link (assumes proportional odds)" } if "`link'"=="cloglog" { di as txt "Complementary Log-Log Link (assumes proportional hazards)" } if "`link'"=="probit" { di as txt "Probit Link (assumes proportional probits)" } } } quietly { ******************************************************************************* *Clean this mess on up! * ******************************************************************************* *Drop that which must be dropped matrix drop _One _Period _Q _Q2 _Spec _Survival if `tpar'==-2 | `tpar'==-1 { matrix drop _Zero } if `tpar'>=1 { matrix drop _ParameterizedTime } ******************************************************************************* *Graph some output! * ******************************************************************************* restore, preserve svmat _Hazard, names(col) *Generate cumulative incidence gen InvSurvival = 1 - Survival *Generate upper and lower bounds for Hazard gen Hub = Hazard+(invnormal(1-((1-(`level'/100))/2))*sigma_h) if Period > 0 replace Hub = 1 if Hub > 1 & Period > 0 gen Hlb = Hazard-(invnormal(1-((1-(`level'/100))/2))*sigma_h) if Period > 0 replace Hlb = 0 if Hlb < 0 & Period > 0 *Generate upper and lower bounds for Survival gen Sub = Survival+(invnormal(1-((1-(`level'/100))/2))*sigma_S) if Period > 0 replace Sub = 1 if Sub > 1 & Period > 0 gen Slb = Survival-(invnormal(1-((1-(`level'/100))/2))*sigma_S) if Period > 0 replace Slb = 0 if Slb < 0 & Period > 0 *Generate upper and lower bounds for cumulative incidence gen Cub = InvSurvival+(invnormal(1-((1-(`level'/100))/2))*sigma_S) if Period > 0 replace Cub = 1 if Cub > 1 & Period > 0 gen Clb = InvSurvival-(invnormal(1-((1-(`level'/100))/2))*sigma_S) if Period > 0 replace Clb = 0 if Clb < 0 & Period > 0 set more on *Manage graphing options input for all graphs if `graph'>4 | `graph'<1 { local `graph'=0 } if `"`b1title'"'==`""' { local b1title = `"b1title("Period")"' } else { local b1title = `"b1title(`b1title')"' } if "`ytick'"=="" { local ytick "yt(0(.2)1)" } else { local ytick "yt(`ytick')" } if "`ylabel'"=="" { local ylabel "yla(0(.2)1)" } else { local ylabel "ylabel(`ylabel')" } *Graph conditional hazard probabilities if `graph'==1 { *Deal with axes specific to conditional hazard curves if "`xtick'"=="" { local xtick "xt(1(1)`maxl')" } else { local xtick "xs(`xtick')" } if "`xlabel'"=="" { local xlabel "xla(1(1)`maxl')" } else { local xlabel "xlabel(`xlabel')" } if `"`l1title'"'==`""' { local l1title = `"l1title("Estimated Conditional Hazard Probability")"' } else { local l1title = `"l1title(`l1title')"' } gr Hazard Period if Period~=0, `ytick' `ylabel' `l1title' `xtick' `xlabel' `b1title' title("Hazard Curve: conditional hazard probability vs. period") connect(l) `options' } *Graph survival probabilities if `graph'==2 | `graph'==4 { *Deal with axes specific to survival curves if "`xtick'"=="" { local xtick = "xt(0(1)`maxl')" } else { local xtick = "xt(`xtick')" } if "`xlabel'"=="" { local xlabel = "xla(0(1)`maxl')" } else { local xlabel = "xlabel(`xlabel')" } if `"`l1title'"'==`""' { if `graph' == 2 { local l1title = `"l1title("Estimated Survival Probability")"' } if `graph' == 4 { local l1title = `"l1title("Estimated Cumulative Event Probability")"' } } else { local l1title = `"l1title(`l1title')"' } gen InvSurvival = 1 - Survival if `graph' == 2 { gr Survival Period, `ytick' `ylabel' `l1title' `xtick' `xlabel' `b1title' title("Survival Curve: survival probability vs. period") connect(l) `options' } if `graph' == 4 { gr InvSurvival Period, `ytick' `ylabel' `l1title' `xtick' `xlabel' `b1title' title("Cumulative Event Curve: cumulative event probability vs. period") connect(l) `options' } } *Graph both conditional hazard and survival probabilities if `graph'==3 { *Prepare for different x-axis labeling needs for hazard and survival local tempxt = "`xtick'" local tempxla = "`xlabel'" local tempytitle = "`ytitle'" if "`xtick'"=="" { local xtick = "xt(1(1)`maxl')" } else { local xtick = "xt(`xtick')" } if "`xlabel'"=="" { local xlabel = "xla(1(1)`maxl')" } else { local xlabel = "xlabel(`xlabel')" } if "`l1title'"=="" { local l1title = `"l1title("Estimated Conditional Hazard Probability")"' } else { local l1title = "l1title(`l1title')" } gr Hazard Period if Period~=0, `ytick' `ylabel' `l1title' `xtick' `xlabel' `b1title' title("Hazard Curve: conditional hazard probability vs. period") connect(l) `options' local xtick = "`tempxt'" local xlabel = "`tempxla'" local ytitle = "`tempytitle'" if "`xtick'"=="" { local xtick = "xt(0(1)`maxl')" } else { local xtick = "xt(`xtick')" } if "`xlabel'"=="" { local xlabel = "xla(0(1)`maxl')" } else { local xlabel = "xlabel(`xlabel')" } if `"`l1title'"'==`""' { local l1title = `"l1title("Estimated Survival Probability")"' } else { local l1title = `"l1title(`l1title')"' } gr Survival Period, `ytick' `ylabel' `t1title' `xtick' `xlabel' `b1title' title("Survival Curve: estimated survival probability vs. period") connect(l) `options' } restore, preserve ******************************************************************************* *Say goodbye to the nice user! * ******************************************************************************* *Retain _Hazard for user... ereturn matrix SurvivalSE = _Hazard[1..`maxl',2] ereturn matrix Survival = _Hazard[1..`maxl',4] ereturn matrix HazardSE = _Hazard[1..`maxl',3] ereturn matrix Hazard = _Hazard[1..`maxl',5] return clear * clean up matrix drop _Hazard _sigma_S _sigma_h V betas _sigma_Sb _sigma_Sa } end