*****************************************************************************
* mayo-1-do.txt                                                             *
*                                                                           *
* PURPOSE:  analysis of the PBC data                                        *
*                                                                           *
* DATE:  01/05/05                                                           *
*                                                                           *
*****************************************************************************
***
*** Mayo Model Variables
***
***	 1. age					
***	 2. albumin				
***	 5. serum bilirubin			
***	 8. edema treatment			
***	12. prothrombin time			
***	16. stage				
***
***	10. time				
***	17. censoring				
***						
*****************************************************************************

***
*** input the data and recode -9 to missing
***

#delimit;
infile 
         age
         albu
         alkphos
         ascites
         bili
         chol
         edema
         edemaTx
         hepmeg
         time
         plate
         proth
         sex
         sgot
         spiders
         stage
         status
         treat
         trigly
         copper        using mayo.dat;
#delimit cr

***
*** recode
***
generate logalb = log(albu)
generate logbil = log(bili)
generate logpro = log(proth)


recode alkphos  -9=.
recode ascites  -9=.
recode hepmeg   -9=.
recode sex      -9=.
recode sgot     -9=.
recode spiders  -9=.
recode stage    -9=.
recode treat    -9=.
recode trigly   -9=.
recode copper   -9=.


***
*** labels
***
label variable time "time (days)"
label variable status "status"

***
*** EDA for predictors
***
summarize logalb
centile logalb, centile( 10 25 50 75 90 )

summarize logbil
centile logbil, centile( 10 25 50 75 90 )

summarize age
centile age, centile( 10 25 50 75 90 )

summarize logpro
centile logpro, centile( 10 25 50 75 90 )

tabulate stage

tabulate edema

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

***
*** overall survival
***

sts graph, gwood saving("mayo-1-0.plot")
gphpen -deps "mayo-1-0.plot"

sts list


***
*** look at logpro groups
***

graph time logpro, symbol([status]) saving("mayo-1-1a.plot") xline( 2.3, 2.4 )
gphpen -deps "mayo-1-1a.plot"

generate procat = logpro
recode procat min/2.2999=1 2.30/2.3999=2 2.40/max=3
label variable procat "Categories of log(prothombin)"
label define plab 1 "log(pro)<2.30" 2 "log(pro) 2.30-2.39" 3 "log(pro)>=2.40"
label values procat plab

tabulate procat

sts graph, by(procat) saving("mayo-1-1b.plot")
gphpen -deps "mayo-1-1b.plot"

sts list, by(procat)

***
*** test for equality of survival curves
***

sts test procat, logrank 

sts test procat, wilcoxon