* do "~/Documents/teaching/c2017/biostat208/labs/lab8/lab8.do"

cd "~/Documents/teaching/c2016/biostat208/labs/lab8/"
use lab8, clear

log using "lab8", replace
* summaries
des
sum

* assoc. between CHD & arcus

* using cs
cs chd69 arcus, or

* using logistic
logistic chd69 arcus

* descriptive summaries of age by outcome status
graph box age, over(chd69) name(age_box)
twoway (kdensity age if chd69==0, bw(3) area(1) lpattern(solid)) (kdensity age if chd69==1, bw(3) area(1) lpattern(longdash)), ytitle("Density") xtitle("age") legend(order(1 "no CHD" 2 "CHD")) name(age_den)

* categorical version of age
recode age 36/40=0 41/45=1 46/50=2 51/55=3 56/60=4, gen(agec)
label define agelab 0 "36-40" 1 "41-45" 2 "46-50" 3 "51-55" 4 "56-60" 
label val agec agelab

* CHD by categorical age - frequency tables
tabulate chd69 agec, col
tabodds chd69 agec, or

* CHD by categorical age - logistic model

logistic chd69 i.agec
testparm i.agec
contrast q(1).agec, noeffects
test -1.agec + 3.agec + 2*4.agec=0
* evaluate evidence for departue from linearity
contrast q(2/4).agec, noeffects

lincom 3.agec-2.agec
logistic chd69 ib2.agec

* lowess fit & linear logistic model for age as continuous

lowess chd69 age, bw(0.5) logit generate(chdsm) name(age_lws)
logistic chd69 age
logistic chd69 age, coef
predict lp, xb

* generate spline basis and fit model
mkspline agesp=age, cubic
logistic chd69 agesp*
predict lps, xb

* test for nonlinearity of log outcome odds in age
test agesp2 agesp3 agesp4

twoway (connected chdsm age, sort msymbol(none)) (connected lp age, sort msymbol(none)) (connected lps age, sort msymbol(none)), legend(order(1 2 3) label(1 "lowess") label(2 "linear") label(3 "spline") pos(6) ring(0)) name(chdage_pr)

tabodds chd69 age, ciplot


* generate distinct levels of age
levelsof age

* generate log-odds estimates and pointwise 95% CI for unique ages based on the previously fitted spline model, save results as new variables for graphing ("x" represents age)
xblc agesp*, covname(age) at(`r(levels)') generate(x losp lolb loub)

* plot, with style options suitable for printing/presentation
twoway (rarea lolb loub x, color(gray)) (connected losp x, sort msymbol(none) color(black) lpattern(solid)) , xtitle("Age") ytitle("Log Odds of CHD") legend(order(2 1) label(1 "95% CI") label(2 "Log odds") pos(5) ring(0) region(lwidth(none))) plotregion(style(none)) scheme(s1mono) name(chdage_lsp)

* convert to probability scale from log odds
gen prsp = exp(losp)/(1 + exp(losp))
gen plb = exp(lolb)/(1 + exp(lolb))
gen pub = exp(loub)/(1 + exp(loub))
twoway (rarea plb pub x, color(gray)) (connected prsp x, sort msymbol(none) color(black) lpattern(solid)) , xtitle("Age") ytitle("Probability CHD") legend(order(2 1) label(1 "95% CI") label(2 "Probability") pos(6) ring(0) region(lwidth(none))) plotregion(style(none)) scheme(s1mono) name(chdage_psp)
graph export pr_age_sp.pdf, replace

log close