cd "/Users/dave/UCSF/Projects/Teaching/Teaching 12/Lab 1/"
use "lab1.dta", clear
xtile median_cut=age, nq(2)
logistic died median_cut
predict out_median
recode age min/4=1 4/max=0, generate(lessthan4)
logistic died lessthan4
predict out_lt4
twoway (line out_median age, sort connect(stairstep)) (line out_lt4 age, sort connect(stairstep)), ytitle(Probability of Dying in One Year) xtitle(Age at Transplant) scheme(s1color) legend(order(1 "Cut at 14" 2 "Cut at 4"))
graph export "fig1.pdf", as(pdf) replace
xtile age4=age, nq(4)
xi: logistic died i.age4
predict p_age4
twoway (line out_median age, sort connect(stairstep)) (line out_lt4 age, sort connect(stairstep)) (line p_age4 age, sort connect(stairstep)), ytitle(Probability of Dying in One Year) xtitle(Age at Transplant) scheme(s1color) legend(order(1 "Cut at 14" 2 "Cut at 4" 3 "Quartiles"))
graph export "fig2.pdf", as(pdf) replace
logistic died age
predict plin
twoway (line plin age, sort) (line p_age4 age, sort connect(stairstep)), ytitle(Probability of Dying in One Year) xtitle(Age at Transplant) scheme(s1color) legend(order(1 "Linear" 2 "Quartiles"))
graph export "fig3.pdf", as(pdf) replace
*
mkspline cube_age = age, cubic nknots(3) displayknots
logistic died cube_age*
predict cube3_fit
est store A
quietly logistic died age
lrtest A

mkspline cube5_age = age, cubic nknots(5) displayknots
logistic died cube5_age*
predict cube5_fit
twoway (line plin age, sort) (line cube3_fit age, sort) (line cube5_fit age, sort), ytitle(Probability of Dying in One Year) xtitle(Age at Transplant) scheme(s1color) legend(order(1 "Linear Fit" 2 "Spline 3 Knot Fit" 3 "Spline 5 Knot Fit"))
graph export "fig4.pdf", as(pdf) replace
*
logistic died prevtx txtype cube_age*
est store B
gen out=_b[_cons]+_b[prevtx]*0.14+_b[txtype]*0.48+_b[cube_age1]*cube_age1+_b[cube_age2]*cube_age2
gen adjpr = exp(out)/(1+exp(out))
adjust prevtx txtype, gen(adjpr2) pr
quietly logistic died  prevtx txtype age
lrtest B
twoway (line adjpr age, sort) (line adjpr2 age, sort), ytitle(Probability of Dying in One Year) xtitle(Age at Transplant) scheme(s1color) legend(order( 1 "Adjusted #1" 2 "Adjusted #2"))
graph export "fig5.pdf", as(pdf) replace
*
mkspline cube7_age = age, cubic nknots(7) displayknots