# R-commands on frailty models for the rat data
# =============================================

# The data are avaiable in the survival library,
# but note that the documentation given there is wrong.


# We will use the R package "parfm", so this has to be installed

# Load the library (which automatically loads the survival library, and more):
library(parfm)

# Inspect the rat data:
rats

# Fit a frailty model with gamma frailty (with mean 1 and variance theta)
# and Weibull baseline. Treatment (rx) is a covariate, and we use years as time scale:
fit.rats=parfm(Surv(time/52,status)~rx,cluster="litter",frailty="gamma",data=rats)
print(fit.rats)


# Fit model with no frailty:
fit.rats0=parfm(Surv(time/52,status)~rx,data=rats)
print(fit.rats0)

# Compute likelihood ratio statistic (i.e. two times the difference in log-likelihoods)
LR=2*(attributes(fit.rats)$loglik- attributes(fit.rats0)$loglik)

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