# R-commands on frailty models 
# ============================

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

# We first look at the rat data (cf. lectures to weeks 46 and 47)

# 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)

# We compute empirical Bayes estimates of the frailties for each litter:
predict(fit.rats)

# We plot the empirical Bayes estimates:
plot(predict(fit.rats))



# We then consider the eye data from practical exercise 12

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

# For the sake of completeness, we start by repeating the command from question b in practical exercise 12

fit.b=parfm(Surv(time,status)~factor(trt),cluster="id", frailty="gamma", data=eye)
print(fit.b)

# We then compute and plot the empirical Bayes estimates of the frailties 
plot(predict(fit.b))