### PubH 7450
### Example: using the larynx cancer data as discussed in Example 12.2,
###   1. fit a semi-parametric AFT model using the Buckley-James method;
###   2. compare with a parametric (e.g., Weibull) AFT;
###   3. compare a semi-parametric PHM with a parametric (e.g., Weibull) PHM
###

dat<-read.table('/home/merganser/weip/public_html/course/7450/data/larynx.txt',
                skip=18,
                col.names=c("stage","time", "age", "year", "status"))
library(rms)

fit1<-bj(Surv(time, status)~ factor(stage)+age, data=dat)

fit1
#Buckley-James Censored Data Regression
#
#bj(formula = Surv(time, status) ~ factor(stage) + age, data = dat)
#
#       Obs     Events       d.f. error d.f.      sigma 
#        90         50          4         45      1.052 
#
#             Value Std. Error       Z  Pr(>|Z|)
#Intercept  3.19709    0.92298  3.4639 5.324e-04
#stage=2   -0.20702    0.48840 -0.4239 6.717e-01
#stage=3   -0.92008    0.37272 -2.4685 1.357e-02
#stage=4   -1.83530    0.41932 -4.3768 1.204e-05
#age       -0.01695    0.01349 -1.2566 2.089e-01

###compared with a Weibull AFT model (see Table 12.1):

library(survival)

fit2<-survreg(Surv(time, status)~ factor(stage)+age, data=dat, dist="weibull")

summary(fit2)
#Call:
#survreg(formula = Surv(time, status) ~ factor(stage) + age, data = dat, 
#    dist = "weibull")
#                 Value Std. Error      z        p
#(Intercept)     3.5288     0.9041  3.903 9.50e-05
#factor(stage)2 -0.1477     0.4076 -0.362 7.17e-01
#factor(stage)3 -0.5866     0.3199 -1.833 6.68e-02
#factor(stage)4 -1.5441     0.3633 -4.251 2.13e-05
#age            -0.0175     0.0128 -1.367 1.72e-01
#Log(scale)     -0.1223     0.1225 -0.999 3.18e-01
#
#Scale= 0.885 
#
#Weibull distribution
#Loglik(model)= -141.4   Loglik(intercept only)= -151.1
#        Chisq= 19.37 on 4 degrees of freedom, p= 0.00066 
#Number of Newton-Raphson Iterations: 5 
#n= 90 

######################################
# correspondongly, one can compare a semi-paramrtic PHM with
# a parametric (e.g. Weibull) PHM:
#semi-parametric PHM:
fit3<-coxph(Surv(time, status)~ factor(stage)+age, data=dat)
fit3
#                 coef exp(coef) se(coef)     z       p
#factor(stage)2 0.1400      1.15   0.4625 0.303 7.6e-01
#factor(stage)3 0.6424      1.90   0.3561 1.804 7.1e-02
#factor(stage)4 1.7060      5.51   0.4219 4.043 5.3e-05
#age            0.0190      1.02   0.0143 1.335 1.8e-01
#
#Likelihood ratio test=18.3  on 4 df, p=0.00107  n= 90 


library(eha)
fit3b<-coxreg(Surv(time, status)~ factor(stage)+age, data=dat)
fit3b
#Covariate           Mean       Coef     Rel.Risk   S.E.    Wald p
#factor(stage) 
#               1    0.459     0         1 (reference)
#               2    0.197     0.140     1.150     0.462     0.762 
#               3    0.281     0.642     1.901     0.356     0.071 
#               4    0.063     1.706     5.507     0.422     0.000 
#age                64.388     0.019     1.019     0.014     0.182 
#
#Events                    50 
#Total time at risk         377.8 
#Max. log. likelihood      -187.71 
#LR test statistic         18.3 
#Degrees of freedom        4 
#Overall p-value           0.00107220


####parametric PHM: compare to Table 12.2
fit4<-phreg(Surv(time, status)~ factor(stage)+age, data=dat, dist="weibull")
fit4
#Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
#factor(stage) 
#               1    0.459     0         1           (reference)
#               2    0.197     0.167     1.182     0.461     0.717 
#               3    0.281     0.663     1.940     0.356     0.062 
#               4    0.063     1.745     5.726     0.415     0.000 
#age                64.388     0.020     1.020     0.014     0.166 
#
#log(scale)                    1.974     7.196     0.138     0.000 
#log(shape)                    0.122     1.130     0.123     0.318 
#
#Events                    50 
#Total time at risk         377.8 
#Max. log. likelihood      -141.42 
#LR test statistic         19.4 
#Degrees of freedom        4 
#Overall p-value           0.000663697