* do file for running Poisson regression on AML dataset

* create log file
log using "poisson reg.txt", text replace


* The analysis consists of the following steps:

* Part a.  Input data, define as survival dataset
* Part b.  Define survival variables:  stset
* Part c.  Fit the Cox PH model
* Part d.  Perform some Poisson regression models on the data



* Part a.  Input data, define as a survival dataset

*   id,  x(0=no maint 1=maint), t = time to relapse,  failed=(1=relaspsed 0=censored)

input id x t failed
1   1  9   0
2   1  13  0
3   1  13  1
4   1  18  0
5   1  23  0
6   1  28  1
7   1  31  0
8   1  34  0
9   1  45  1
10  1  48  0
11  1  161 1
12  0  5   0
13  0  5   0
14  0  8   0
15  0  8   0
16  0  12  0
17  0  16  1
18  0  23  0
19  0  27  0
20  0  30  1
21  0  33  0
22  0  43  0
23  0  45  0
end


* Part b.  Define survival variables:  stset

stset t , failure(failed==1) id(id)


* Save as Stata dataset

save cl10ex1.dta , replace



* Part c.  Fit the Cox PH model

use cl12ex1.dta , clear
stcox x

* Part d.  Perform some Poisson regression models on this data:

* bin the survival data

stsplit tbin , at( 10(10)50)
strate tbin x , output(binrates.dta,replace)

* use some Poisson models to estimate the treatment effect
use binrates.dta, clear
gen midp = tbin + 5

* Graph the observed log rates!
gen lograte = log(_D/_Y)
twoway (scatter lograte midp if x==1, ms(+) legend(off)) (scatter lograte midp if x==0, ms(x) xtitle(Midpoit of Time Interval) ytitle(Log Rate of Relapse) title("+ = Maintained, x = Not Maintained"))
graph export "amldata.wmf", replace

* Linear Time
poisson _D x midp, exposure(_Y) irr

* Non-Linear Time, change of slope at 25 weeks
gen midpsp = 0
replace midpsp = midp if midp>25
poisson _D x midp midpsp, exposure(_Y) irr

* Non-proportional hazards?
gen inter = x * midp
poisson _D x midp inter, exposure (_Y) irr

* Step function of time
xi: poisson _D x i.midp, exposure(_Y) irr

log close