*****************************************************************
***  leuk2.do                                                 ***
*****************************************************************  
infile time status tx logwbc using leuk2.dat

label variable time "time (weeks)"
label variable status "status"
label variable tx "treatment"
label variable logwbc "log(white blood cell count)"

list

***
*** recode tx
***
recode tx 1=0 2=1

label define tlab 0 "treatment" 1 "placebo"
label values tx tlab


***
*** summarize wbc by tx
***
sort tx
by tx:  summarize logwbc

***
*** center logwbc = important for survival!
***
generate newlwbc = logwbc-3.00


***
*** this defines the failure outcome 
***
stset time, failure(status)
stset, noshow


***
*** Univariate analysis with treatment only
***
sts graph, by(tx) clwidth( medthick medthick ) clpattern( solid dash )
graph export "c:\courses\va-survival\leuk2-0.eps", as(eps) replace

***
*** Cox regression with TX
***
stcox tx, nohr basesurv(shat)
stcox tx

graph twoway (scatter shat time)
graph export "c:\courses\va-survival\leuk2-0a.eps", as(eps) replace

stcoxkm, by(tx)
graph export "c:\courses\va-survival\leuk2-0b.eps", as(eps) replace

***
*** let's look at KM curves for levels of WBC
***
generate wbccat = logwbc
recode wbccat min/1.99=1 2.00/2.99=2 3.00/3.99=3 4.00/max=4
label define wlab 1 "log(wbc) < 2.00" 2 "log(wbc) 2.00-2.99" 3 "log(wbc) 3.00-3.99" 4 "log(wbc) >= 4.00"
label values wbccat wlab

table wbccat

sts graph, by(wbccat) clwidth( medthick medthick medthick medthick ) ///
           clpattern( solid dash solid dash )
graph export "c:\courses\va-survival\leuk2-1.eps", as(eps) replace

***
*** and log-rank test
***
sts test wbccat, logrank


***
*** Cox regression
***
stcox tx, nohr
stcox tx
*lrtest, saving(1)
est store mod1

stcox tx newlwbc, nohr
stcox tx newlwbc
*lrtest, saving(2)
est store mod2

xi: stcox i.tx*newlwbc, nohr
xi: stcox i.tx*newlwbc
*lrtest, saving(3)
est store mod3

*lrtest, using(3) model(2)
lrtest mod3 mod2, stats

*lrtest, using(2) model(1)
lrtest mod2 mod1, stats


***
*** use model 2
***
stcox tx newlwbc, nohr basesurv( s0hat ) 


***
*** KM and adjusted KM
***
sts graph, by(tx) clwidth( medthick medthick ) clpattern( solid dash )
graph export "c:\courses\va-survival\leuk2-2.eps", as(eps) replace

sts graph, by(tx) adjustfor(newlwbc) clwidth( medthick medthick ) ///
           clpattern( solid dash )
graph export "c:\courses\va-survival\leuk2-3.eps", as(eps) replace

sts graph, strata(tx) adjustfor(newlwbc) clwidth( medthick medthick ) ///
           clpattern( solid dash )
graph export "c:\courses\va-survival\leuk2-3b.eps", as(eps) replace


* stcoxkm, strata(tx) adjustfor(newlwbc)