// Table 2.1
use https://statisticalhorizons.com/wp-content/uploads/rank.dta, clear
ltable dur promo, failure hazard noadjust
// Table 2.2
use https://statisticalhorizons.com/wp-content/uploads/rank.dta, clear
generate id = _n
reshape long art cit, i(id) j(year)
drop if year > dur // Remove empty observation created during reshape
replace promo = 0 if year < dur // Create time-varying failure variable
generate jobpres = prest1 // Create-time varying prestige varianle
replace jobpres = prest2 if year >= jobtime
eststo clear
eststo: logit promo undgrad phdmed phdprest jobpres art cit
estadd expb
eststo: logit promo undgrad phdmed phdprest jobpres art cit year c.year#c.year
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(2))") varwidth(17) label nonumber ///
mtitle("Model 1" "Model 2") interaction( X ) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Logistic Models Predicting the Probability of Promotion") ///
legend
// Likelihood ratio test on p. 13
lrtest est1 est2
// Additional model on p. 13
logit promo undgrad phdmed phdprest jobpres art cit year c.year#c.year c.phdprest#c.year
// Table 2.3
// Scenario A: All censored cases are in fact promoted
use https://statisticalhorizons.com/wp-content/uploads/rank.dta, clear
generate id = _n
reshape long art cit, i(id) j(year)
drop if year > dur // Remove empty observation created during reshape
replace promo = 0 if year < dur // Create time-varying failure variable
replace promo = 1 if year == dur & dur < 10 // All censored cases are promoted
generate jobpres = prest1 // Create-time varying prestige varianle
replace jobpres = prest2 if year >= jobtime
eststo: logistic promo undgrad phdmed phdprest jobpres art cit year c.year#c.year
// Scenario B: All censored cases did not experience event until the end of the observation period
use https://statisticalhorizons.com/wp-content/uploads/rank.dta, clear
generate id = _n
reshape long art cit, i(id) j(year)
drop if year > dur & promo == 1 // Remove empty observation created during reshape
// if they are promoted
replace promo = 0 if year < dur // Create time-varying failure variable
replace art = art[_n-1] if art == . // Carry observations forward
replace cit = cit[_n-1] if cit == . // Carry observations forward
generate jobpres = prest1 // Create-time varying prestige varianle
replace jobpres = prest2 if year >= jobtime
eststo: logistic promo undgrad phdmed phdprest jobpres art cit year c.year#c.year
esttab est2 est3 est4, eform b(2) z wide staraux varwidth(17) label nonumber ///
mtitle("Standard analysis" "Scenario A" "Scenario B") interaction( X ) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Extreme Case Scenarios for Informative Censoring") ///
legend modelwidth(18)
// Table 3.1
use https://statisticalhorizons.com/wp-content/uploads/recid.dta, clear
stset week, failure(arrest==1)
eststo clear
eststo: streg fin age race wexp mar paro prio, dist(exponential)
estadd expb
eststo: streg fin age race wexp mar paro prio, dist(weibull)
estadd expb
eststo: streg fin age race wexp mar paro prio, dist(ggamma)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(17) label nonumber ///
mtitle("Exponential" "Weibull" "Gamma") ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Estimates for Three Models of Recidivism") ///
legend keep(_t:*) eqlabels("", none) /// // Removes equation label
coeflabel(fin "Financial aid" age "Age at release" race "Black" ///
wexp "Work experience" mar "Married" paro "Paroled" ///
prio "Prior convictions")
// Table 3.2
use https://statisticalhorizons.com/wp-content/uploads/recid.dta, clear
stset week, failure(arrest==1)
eststo clear
qui eststo: streg fin age race wexp mar paro prio, dist(ggamma)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(lognormal)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(llogistic)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(weibull)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(gompertz)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(exponential)
qui estadd scalar dev = -2*e(ll)
qui esttab, cells(none) scalars("dev Deviance" "aic AIC" "bic BIC") sfmt(3) nomtitles noobs
esttab r(stats, transpose), coeflabel(est1 "Gamma" est2 "Log-normal" ///
est3 "Log-logistic" est4 "Weibull" ///
est5 "Gompertz" est6 "Exponential") ///
title("Goodness of Fit for Recidivism Models") ///
nomtitle collabels("Deviance" "AIC" "BIC")
// Figure 3.1: Hazard Function for Weibull Regression Model
qui streg fin age race wexp mar paro prio, dist(weibull)
stcurve, hazard xtitle("Weeks since release")
// Figure 3.2: Hazard Function for Log-Logistic Regression Model
qui streg fin age race wexp mar paro prio, dist(llogistic)
stcurve, hazard xtitle("Weeks since release")
// Table 4.1
use https://statisticalhorizons.com/wp-content/uploads/recid.dta, clear
stset week, failure(arrest==1)
eststo clear
eststo: stcox fin age race wexp mar paro prio
estadd expb
generate id = _n
reshape long work, i(id) j(stop)
generate start = stop - 1
drop if stop > week // Drop empty observations
replace arrest = 0 if week != stop // Create time-varying failure variable
stset stop, failure(arrest == 1) id(id) origin(start)
eststo: stcox fin age race wexp mar paro prio work
estadd expb
generate worklag = work[_n-1] if start > 0
eststo: stcox fin age race wexp mar paro prio worklag
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(18) label nonumber ///
mtitle("Basic" "Time-varying X" "Lagged X") modelwidth(15) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Cox Regression Estimates for Recidivism Data") legend ///
coeflabel(fin "Financial aid" age "Age at release" race "Black" ///
wexp "Work experience" mar "Married" paro "Paroled" ///
prio "Prior convictions" work "Employment") rename(worklag work)
// Table 4.2
list id stop start arrest work fin age if inlist(id, 339, 417), noobs sepby(id)
// Table 4.3
estimates restore est2
estat phtest, detail
// Table 4.4
eststo clear
eststo: stcox fin age race wexp mar paro prio work, tvc(age wexp)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("Interactions with time") modelwidth(22) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Cox Regression Estimates with Time Interactions") legend ///
eqlabels("Main effects" "Interactions with time") ///
coeflabel(fin "Financial aid" age "Age at release" race "Black" ///
wexp "Work experience" mar "Married" paro "Paroled" ///
prio "Prior convictions" work "Employment")
eststo clear
eststo: stcox fin age race wexp mar paro prio work, strata(wexp)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("Stratification") modelwidth(22) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Cox Regression Estimates with Stratification") legend ///
coeflabel(fin "Financial aid" age "Age at release" race "Black" ///
wexp "Work experience" mar "Married" paro "Paroled" ///
prio "Prior convictions" work "Employment") dropped(--)
// Table 4.5
qui stcox fin age race wexp mar paro prio work, tvc(age wexp)
foreach x of numlist 0 10 20 30 40 50 {
local b1 = _b[age] + `x' * _b[tvc:age]
local expb1 = exp(`b1')
local b2 = _b[wexp] + `x' * _b[tvc:wexp]
local expb2 = exp(`b2')
matrix results = (`x' , `b1' , `expb1' , `b2', `expb2')
if `x' == 0 matrix table45 = results
else matrix table45 = (table45\results)
}
esttab matrix(table45, fmt(0 3 3 3 3)), nomtitle ///
collabel("Weeks" "b" "Exp(b)" "b" "Exp(b)") varwidth(0) ///
title("Effects of Age and Work Experience at Different Times")
// Table 4.6
use https://statisticalhorizons.com/wp-content/uploads/recid.dta, clear
stset week, failure(arrest==1)
stcox fin age race wexp mar paro prio
stcurve, survival at(fin = 1 age = 21 race = 1 wexp = 1 mar = 0 paro = 1 prio = 4) ///
outfile(surv, replace)
use surv, clear
egen pickone = tag(_t)
list _t surv1 if inlist(_t, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 52) & pickone, noobs sep(0)
// Table 5.1
use https://statisticalhorizons.com/wp-content/uploads/tarp.dta, clear
eststo clear
stset arrstday, failure(type == 1 2)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
stset arrstday, failure(type == 1)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
stset arrstday, failure(type == 2)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("All arrests" "Property arrests" "Non-property arrests") modelwidth(18) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Estimates of Proportional Hazards Models for Different Arrest Types") legend ///
coeflabel(fin "Financial aid" age "Age at release" white "White" ///
male "Male" married "Married" paro "Paroled" ///
numprop "No. of property convictions" ///
crimprop "Imprisoned for property crime" ///
numarst "No. of arrests" ///
edcomb "Education")
// Figure 5.1
stset arrstday, failure(type == 1)
stcompet cumin = ci, compet1(2)
sort _t
twoway (line cumin _t if type == 1) ///
(line cumin _t if type == 2) ///
(scatteri .2 375 "Property", msymbol(none) mlabpos(0)) ///
(scatteri .15 370 "Non-property", msymbol(none) mlabpos(0)) ///
, legend(off) ///
xtitle("Days since release") ytitle("Probability of arrest") ///
title("Cumulative incidence")
// Table 5.2
eststo clear
stset arrstday, failure(type == 1 2)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
stset arrstday, failure(type == 2)
eststo: stcrreg fin age white male married paro numprop crimprop numarst edcomb, compete(type == 1)
estadd expb
stset arrstday, failure(type==1)
eststo: stcrreg fin age white male married paro numprop crimprop numarst edcomb, compete(type==2)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("All arrests" "Property arrests" "Non-property arrests") modelwidth(18) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Estimates of Subdistribution Hazards Models for Different Arrest Types") legend ///
coeflabel(fin "Financial aid" age "Age at release" white "White" ///
male "Male" married "Married" paro "Paroled" ///
numprop "No. of property convictions" ///
crimprop "Imprisoned for property crime" ///
numarst "No. of arrests" ///
edcomb "Education")
// Figure 5.2
qui stcrreg fin age white male married paro numprop crimprop numarst edcomb, compete(type==2)
stcurve, cif xtitle("Days since release") ///
title("Cumulative incidence of non-property arrests") ///
ylab(, format(%6.2f))
// Table 6.1
use https://statisticalhorizons.com/wp-content/uploads/tarp.dta, clear
eststo clear
estpost tabulate arrstcount
esttab, cell(b) nonumber collabel("Number of persons", lhs("Number of arrests")) nomtitle noobs ///
modelwidth(20) varwidth(20) varlabels(, blist(Total "{hline @width}{break}")) ///
title(Frequency Distribution for Number of Arrests)
// Table 6.2
eststo clear
eststo: nbreg arrstcount fin age white male married paro numprop crimprop numarst edcomb
estadd expb
use https://statisticalhorizons.com/wp-content/uploads/arrests.dta, clear
stset length, failure(arrind == 1)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
stcox fin age white male married paro numprop crimprop numarst edcomb, cluster(id)
set matsize 942
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb, shared(id)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("Negstive binomial count model" "Cox regression, gap time" "Cox regression, shared frailty") modelwidth(18) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Regression Models for Repeated Arrests") legend ///
coeflabel(fin "Financial aid" age "Age at release" white "White" ///
male "Male" married "Married" paro "Paroled" ///
numprop "No. of property convictions" ///
crimprop "Imprisoned for property crime" ///
numarst "No. of arrests" ///
edcomb "Education") keep(main:*) drop(_cons)
// Table 6.3
use https://statisticalhorizons.com/wp-content/uploads/arrests.dta, clear
stset length, failure(arrind == 1)
eststo clear
eststo: streg spellnum fin age white male married paro numprop crimprop numarst edcomb, cluster(id) dist(weibull)
estadd expb
eststo: streg spellnum fin age white male married paro numprop crimprop numarst edcomb, shared(id) dist(weibull)
estadd expb
stset end, failure(arrind==1) origin(begin)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb, cluster(id)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("Weibull, robust z" "Weibull, shared frailty" "Cox regression, origin times") modelwidth(18) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Regression Models for Repeated Arrests") legend ///
coeflabel(fin "Financial aid" age "Age at release" white "White" ///
male "Male" married "Married" paro "Paroled" ///
numprop "No. of property convictions" ///
crimprop "Imprisoned for property crime" ///
numarst "No. of arrests" ///
edcomb "Education") keep(main:)
// Analyses on p. 74
stcox fin age white male married paro numprop crimprop numarst edcomb, cluster(id) tvc(numarst) texp(_t/30.4)
Showing posts with label stset. Show all posts
Showing posts with label stset. Show all posts
Mar 11, 2018
Allison's (2014) book on event history and survival analysis
This replicates the analyses in Allison (2014). There is another approach here, but I like mine better.
Jul 4, 2017
Chapters 9 and 10 of Singer and Willett's (2003) book on longitudinal data analysis
// Figure 9.1use "C:\singer willett (2003)\teachers.dta", clear label define censor 0 "Not censored" 1 "Censored" label val censor censor histogram t, by(censor, legend(off) note("")) freq xlabel(1 (1) 12) ylabel(0 (100) 500) /// addlabels discrete xtitle("Years of teaching") name(figure91, replace) // Table 10.1generate event = !censor qui ltable t event, noadjust saving(lifetable, replace) preserve use lifetable, clear list t0 t1 start deaths lost hazard survival, sep(0) noobs erase lifetable.dta restore // Figure 10.1stset t, failure(event) sts generate h = h twoway scatter h t, msymbol(i) connect(l) ylabel(0 (.05) .15) xlabel(0 (1) 13) /// sort xtitle("Years in teaching") /// ytitle("Estimated hazard probability") /// name(figure101a, replace) nodraw ltable t event, noadjust notab graph noconf xlabel(0 (1) 13) ylabel(0 (.5) 1) /// yline(.5) xline(7.6) /// xtitle("Years in teaching") /// ytitle("Proportion surviving") /// name(figure101b, replace) nodraw graph combine figure101a figure101b, col(1) name(figure101, replace) ysize(8) // Figure 10.2use "C:\singer willett (2003)\relapse_days.dta", clear generate weeks = int(days / 7) + 1 generate event = !censor stset weeks, failure(event) sts generate h = h twoway (scatter h weeks, connect(l) msymbol(i)), xlabel(0 (1) 12) ylabel(, format(%6.2f)) /// ytitle("Estimated hazard probability") /// xtitle("Weeks after release") /// name(figure102a1, replace) nodraw ltable weeks event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) /// xlabel(0(1)12) yline(0.5) /// ytitle("Estimated survival probability") /// xtitle("Weeks after release") /// name(figure102a2, replace) nodraw graph combine figure102a1 figure102a2, col(2) /// title("{bf:A} Time to cocaine relapse", pos(11)) name(figure102a, replace) nodraw use "C:\singer willett (2003)\firstsex.dta", clear generate event = !censor stset time, failure(event) sts generate h = h twoway (scatter h time, connect(l) msymbol(i) sort), xlabel(6 (1) 12) ylabel(, format(%6.2f)) /// ytitle("Estimated hazard probability") /// xtitle("Grade") /// name(figure102b1, replace) nodraw ltable time event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) /// xlabel(6 (1) 12) yline(0.5) /// ytitle("Estimated survival probability") /// xtitle("Grade") /// name(figure102b2, replace) nodraw graph combine figure102b1 figure102b2, col(2) /// title("{bf:B} Age at first intercourse for males", pos(11)) name(figure102b, replace) nodraw use "C:\singer willett (2003)\suicide_orig.dta", clear generate event = !censor stset time, failure(event) sts generate h = h twoway (scatter h time, connect(l) msymbol(i) sort), xlabel(5 (2) 21) ylabel(, format(%6.2f)) /// ytitle("Estimated hazard probability") /// xtitle("Age") /// name(figure102c1, replace) nodraw ltable time event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) /// xlabel(5 (2) 21) yline(0.5) /// ytitle("Estimated survival probability") /// xtitle("Age") /// name(figure102c2, replace) nodraw graph combine figure102c1 figure102c2, col(2) /// title("{bf:C} Age at first suicide ideation", pos(11)) name(figure102c, replace) nodraw use "C:\singer willett (2003)\congress_orig.dta", clear generate event = !censor stset time, failure(event) sts generate h = h twoway (scatter h time, connect(l) msymbol(i) sort), xlabel(0 (1) 8) ylabel(0 (.1) .3, format(%6.2f)) /// ytitle("Estimated hazard probability") /// xtitle("Terms in office") /// name(figure102d1, replace) nodraw ltable time event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) /// xlabel(0 (1) 8) yline(0.5) /// ytitle("Estimated survival probability") /// xtitle("Terms in office") /// name(figure102d2, replace) nodraw graph combine figure102d1 figure102d2, col(2) /// title("{bf:D} Duration of congressional careers for females", pos(11)) name(figure102d, replace) nodraw graph combine figure102a figure102b figure102c figure102d, row(4) ysize(11) name(figure102, replace) // Table 10.2 use "C:\singer willett (2003)\teachers.dta", clear generate event = !censor qui ltable t event, noadjust saving(lifetable, replace) preserve use lifetable, clear list t0 hazard sehazard survival sesurvival, noobs sep(0) erase lifetable.dta restore // Figure 10.4 use "C:\singer willett (2003)\teachers.dta", clear list id t censor if inlist(id, 20, 126, 129), noobs sep(0) expand t bysort id: generate period = _n // Period identifier bysort id: generate event = !censor & _n == _N // Calculate outcome for // discrete-time analysis list id period event if inlist(id, 20, 126, 129), noobs sepby(id) // Table 10.3 table period, c(sum event n event mean event)
Reference
Singer, Judith D., and John B. Willett. 2003. Applied Longitudinal Data Analysis. Modeling Change and Event Occurrence. Oxford University Press. doi: 10.1093/acprof:oso/9780195152968.001.0001
Subscribe to:
Posts (Atom)