quietly: infile age albumin alkphos ascites bili cholest edema edmadj hepmeg   ///
	obstime platelet protime sex sgot spiders stage status treatmnt   ///
	triglyc urinecu using "http://www.emersonstatistics.com/datasets/liver.txt

drop in 1

g male = 1 - sex
recode bili 0/1=1 1/2=2 2/4=4 4/8=8 8/max= 16, gen(bilictg)
bysort bilictg: summ bili
g cbili= bili-1
g logbili= log(bili) / log(2)
replace obstime = obstime / 365.25

****
**** Problem 1: Descriptive statistics
****

*** Censoring distribution
g cens= 1 -status
stset obstime cens
sts graph, t1title("Distribution of Time to Censoring") xtitle("Time on Study (years)")
sts graph, by(bilictg)
sts list, at(2 4 6 8)
stci, rmean
g obssqr= obstime^2
stset obssqr cens
stci, rmean
di sqrt(56.00238 - 6.879207^2)

*** Baseline variables
tabstat male age bili, by(bilictg) stat(n mean sd min q max) col(stat) m long

*** Survival distribution
stset obstime status
sts list, by(bilictg) at(2 4 6 8)
sts list, at(2 4 6 8)
sts graph, risktable(0 2 4 6 8 10) cen(s) by(bilictg)

****
**** Problem 2: Minimum censoring time
****

bysort status: summ obstime
table status if obstime <= 1.459274

****
**** Problem 3: Linear continuous (untransformed) bilirubin
****

stset obstime status
stcox cbili, robust
predict fithr3

****
**** Problem 4: Logarithmically transformed bilirubin
****

stcox logbili, robust
predict fithr4

****
**** Problem 5: Testing linearity of association
****

stcox cbili logbili, robust
predict fithr5

****
**** Problem 6: Fitted values
****

scatter fithr3 fithr4 fithr5 bili, ytitle("Fitted Hazard Ratio Estimates") ///
	xtitle("Serum Bilirubin (mg/dL)") t1title("Fitted Hazard Ratio Estimates") ///
	legend(label(1 "Untransf bili") label(2 "Log bili") label(3 "Both"))
	
****
**** Problem 7: Adjustment for sex, age
****

stcox logbili age male, robust
stcox logbili age if male!=., robust
stcox logbili male, robust
stcox logbili if male!=., robust