# R commands for practical exercise 12

# We read the data:
eye=read.table("http://folk.uio.no/borgan/BGC1-2012/data/retinopathy.txt", header=T)

# We will use the parfm-library, so this must be installed and loaded

# a)
# We first fit a model with Weibull baseline, no frailty and treatment as covariate:
fit.a=parfm(Surv(time,status)~factor(trt),data=eye)
print(fit.a)

#b) 
# We then fit a model with Weibull baseline, gamma frailty and treatment as covariate:
fit.b=parfm(Surv(time,status)~factor(trt),cluster="id", frailty="gamma", data=eye)
print(fit.b)

# We compute likelihood ratio statistic (i.e. two times the difference in log-likelihoods)
LR=2*(attributes(fit.b)$loglik- attributes(fit.a)$loglik)

# We compute P-value for one-sided likelihood ratio test:
(1-pchisq(LR,1))/2

# The frailty is significant.

# c) 
# We fit a model with Weibull baseline, gamma frailty, and treatment, age at onset and type of diabetes as covariates. 
fit.c=parfm(Surv(time,status)~factor(trt)+factor(typediab)+ageonset, cluster="id", frailty="gamma", data=eye)
print(fit.c)

# Type of diabetes and age at onset do not have a significant effect