[R] Any examples of a frailty model actually used for prediction ?

[Koffijberg, H.]

Hi everyone,

I'm looking for any examples of useful frailty models, in particular any 
situation in which a cox proportional hazards model with frailty outperforms a 
regular cox proportional hazards model with respect to prediction of the time 
to event (or the X-year risk of an event). I have defined my own gamma-frailty 
cox PH model in R but on my simulated data sample it does not predict any 
better than a regular cox model. In fact the predictions of the frailty model 
are often worse, *even* when I purposely add gamma-distributed frailty to the 
simulated sample (as I have previously posted). Can anyone help me, e.g. by 
pointing me to some R code in which the performance of models with and without 
frailty component are assessed ? Also any examples in which the actual 
estimated individual frailty values are used for predictions are very much 
welcome.

Thanks.
Erik Koffijberg.

=
Julius Center for Health Sciences and Primary Care
University Medical Center Utrecht Str. 6.131
P.O.Box 85500, 3508 GA Utrecht, The Netherlands
Tel:+31 30 250 3013,
Fax:+31 30 250 5480
E-mail: [EMAIL PROTECTED]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to interpret these results from a simple gamma-frailty model

[Koffijberg, H.]

Dear R users,

I'm trying to fit a gamma-frailty model on a simulated dataset, with 6 
covariates, and I'm running into some results I do not understand. I 
constructed an example from my simulation code, where I fit a coxph model 
without frailty (M1) and with frailty (M2) on a  number of data samples with a 
varying degree of heterogeneity (I'm running R 2.3.1, running takes ~1 min).

library(survival); set.seed(1)
lambda - 0.01  # Exp. hazard rate
# Beta coefficients for Age,TC,HDLC,SBP,Diab,Smok
beta - c(0.0483,0.0064,-0.0270,0.0037,0.4284,0.5234)   
n - 1000; Ngrp - 2;   # Nr patients, Nr frailty groups
# Thetas for gamma-frailty  
thetaset - c(1,2,10,100); Ntheta - length(thetaset);  

# Define the simulated population
age -rnorm(n,48.6,11.7);tc-rnorm(n,200,30)
hdlc-rnorm(n,47,6);sbp-rnorm(n,135,6) 
rtmp - runif(0,1,n=n); diab - rep(0, n); diab[rtmp  0.05223] - 1;
rtmp - runif(0,1,n=n); smok - rep(0, n); smok[rtmp  0.40458] - 1;
# Storage for the estimated coefficients in both models
cf - 0; length(cf) - 2*Ntheta*6; dim(cf) - c(Ntheta,2,6); 

for (i in 1:Ntheta)
{
l = thetaset[i]
v=rep(rgamma(n/2,shape=l,scale=1/l),2)  # Shared frailty values
print(paste(Theta = ,l,  variance in introduced frailty v   = 
,var(v)))
id=rep(seq(1:(n/2)),2); # Frailty id's, 
c - rep(1,n);  # censoring flags
r - runif(n,0,1);  # For random variate 
generation
t - (-1*log(r)) / 
(lambda*v*exp(age*beta[1]+tc*beta[2]+hdlc*beta[3]+sbp*beta[4]+diab*beta[5]+smok*beta[6]))
 
fitM1=coxph(Surv(t,c)~age+tc+hdlc+sbp+diab+smok);   # Model 1, no 
frailty
cf[i,1,] - round(fitM1$coef,6) # Model 1 
coefficients
fitM2=coxph(Surv(t,c)~age+tc+hdlc+sbp+diab+smok+
frailty(id,dist=gamma,sparse=T,method=em))  # Model 2, 
frailty
cf[i,2,] - round(fitM2$coef,6) # Model 2, 
coefficients

# store estimated Variance of random effect

varre=fitM2$history$frailty$history[length(fitM2$history$frailty$history[,1])]  

print(paste(Theta = ,l,  variance in estimated random effect= 
,varre))
print(paste(Theta = ,l,  variance in estimated ind. frailty = , 
var(exp(fitM2$frail
print(paste(Theta=,l, org beta ,beta, M1: ,cf[i,1,],  M2: , 
cf[i,2,]))
# Calculate the 'actual' 10-year risk and the risk estimated using M1 
and M2
estv - exp(rep(fitM2$frail,2)) # Individual frailty values 
from M2
# Original risk score
rs- 
age*beta[1]+tc*beta[2]+hdlc*beta[3]+sbp*beta[4]+diab*beta[5]+smok*beta[6]   
   
pevent - exp(-10*lambda)^(exp(rs)) # Original 
10-year risk
rsM1- 
(cf[i,1,1]*age+cf[i,1,2]*tc+cf[i,1,3]*hdlc+cf[i,1,4]*sbp+cf[i,1,5]*diab+cf[i,1,6]*smok)
peventM1 - exp(-10*lambda)^(exp(rsM1))
rsM2- 
(cf[i,1,1]*age+cf[i,1,2]*tc+cf[i,1,3]*hdlc+cf[i,1,4]*sbp+cf[i,1,5]*diab+cf[i,1,6]*smok)
peventM2 - exp(-10*lambda)^(estv*exp(rsM1))
# Proportion of more accurate predictions
pred - sum(abs(pevent-peventM2)  abs(pevent-peventM1))/n  
print(paste(Theta = ,l, proportion pred. M2 more accurate than M1 = 
,pred)); 
}


corresponding output:

[1] Theta =  1  variance in introduced frailty v   =  0.948105787326522
[1] Theta =  1  variance in estimated random effect=  1.04435029128552
[1] Theta =  1  variance in estimated ind. frailty =  0.520078219281493
[1] Theta= 1  org beta  0.0483  M1:  0.025015  M2:  0.045987  
[2] Theta= 1  org beta  0.0064  M1:  0.002274  M2:  0.005178  
[3] Theta= 1  org beta  -0.027  M1:  -0.014239  M2:  -0.028824
[4] Theta= 1  org beta  0.0037  M1:  0.00071  M2:  0.010826   
[5] Theta= 1  org beta  0.4284  M1:  0.421602  M2:  0.648227  
[6] Theta= 1  org beta  0.5234  M1:  0.360069  M2:  0.593551  
[1] Theta =  1  proportion pred. M2 more accurate than M1 =  0.448
[1] Theta =  2  variance in introduced frailty v   =  0.515806773271924
[1] Theta =  2  variance in estimated random effect=  0.644954876113229
[1] Theta =  2  variance in estimated ind. frailty =  0.277209957022078
[1] Theta= 2  org beta  0.0483  M1:  0.033213  M2:  0.046216  
[2] Theta= 2  org beta  0.0064  M1:  0.005633  M2:  0.009491  
[3] Theta= 2  org beta  -0.027  M1:  -0.010048  M2:  -0.016658
[4] Theta= 2  org beta  0.0037  M1:  -0.006032  M2:  -0.003639
[5] Theta= 2  org beta  0.4284  M1:  0.402992  M2:  0.403226  
[6] Theta= 2  org beta  0.5234  M1:  0.457665  M2:  0.691871  
[1] Theta =  2  proportion pred. M2 more accurate than M1 =  0.482
[1] Theta =  10  variance in introduced frailty v   =  0.104197974543906
[1] Theta =  10  variance in estimated random effect=  0.00331174766884151
[1] Theta =  10  variance in estimated ind. frailty =