*! version 2.0.0 4dec2011 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, 2011 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, 2011 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 inp el(_Spec,1,`x')
macro shift
}
}
if `tpar'~=-1 {
local newarg = `argnum'+1
forvalues x=2/`newarg' {
di as txt "`1'" " = " as inp 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, 2011 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