/// Reading in and preparing data

quietly:                                                                        ///
infile ptid mridate age male race weight height packyrs yrsquit alcoh physact   ///
       chf chd stroke diabetes genhlth ldl alb crt plt sbp aai fev dsst atrophy ///
	   whgrd numinf volinf obstime death                                        ///
	   using http://www.emersonstatistics.com/datasets/mri.txt

list in 1
drop in 1

replace obstime= obstime / 365.25
egen ldlCTG = cut(ldl), at(0 70 100 130 160 190 250)

/// Descriptive statistics for censoring distribution

bysort death: tabstat obstime, stat(n max) col(stat) by(ldlCTG)
g censor = 1 - death
stset obstime censor
sts graph, t1("Censoring Distribution") 
stci, p(10)
stci, p(25)
stci, p(50)
stci, p(75)
stci, p(90)
stci, rmean

/// Descriptive statistics for LDL

tabstat ldl, stat(n mean sd min q max) col(stat)
list obstime death if ldl==.

/// Descriptive statistics for survival by LDL

stset obstime death
sts graph, by(ldlCTG) plotopts(lwid(1.25)) ylabel(0.5(0.1)1.0) ///
	legend(label(1 "LDL 11 - 69 mg/dL") label(2 "LDL 70 - 99 mg/dL") ///
		label(3 "LDL 100 - 129 mg/dL") label(4 "LDL 130 - 159 mg/dL") ///
		label(5 "LDL 160 - 189 mg/dL") label(6 "LDL 190 - 247 mg/dL"))
sts list, by(ldlCTG) at(2 5)
sts list, at(2 5)
stci, p(10) by(ldlCTG)
stci, p(20) by(ldlCTG)

/// To compute restricted mean I restrict observations to 5.75 years based on censoring
/// distribution within the LDL strata

g obs575= obstime
replace obs575= 5.75 if obstime > 5.75
g death575= death
replace death575= 0 if obstime > 5.75
stset obs575 death575
stci, rmean by(ldlCTG)

/// Problem 1

stset obstime death
stcox ldl, robust
g cldl= ldl - 160
stcox cldl, robust
predict fithrA

/// Problem 2 (I use base 1.1 logarithm to be able to talk about 10% difference)

g logldl= log(ldl) / log(1.1)
stcox logldl, robust
g clogldl= log(ldl / 160) / log(1.1)
stcox clogldl, robust
predict fithrB

/// Problem 3

g ldlsqr= ldl^2
stcox ldl ldlsqr, robust
g cldlsqr= cldl^2
stcox cldl cldlsqr, robust
predict fithrC
summ fithrC
list ldl fithrC if fithrC < 0.9908

/// Problem 4

twoway (scatter fithrA fithrB fithrC ldl, color(black blue green)), ///
	legend(label(1 "Linear fit") label(2 "Log fit") label(3 "Quadratic fit"))