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

From: Koffijberg, H. <H.Koffijberg_at_umcutrecht.nl>
Date: Tue 19 Sep 2006 - 12:04:27 GMT

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(10000)
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 = 2.30691702541282e-05"
[1] "Theta= 10 org beta 0.0483 M1: 0.048761 M2: 0.048885"
[2] "Theta= 10 org beta 0.0064 M1: 0.002817 M2: 0.002833"
[3] "Theta= 10 org beta -0.027 M1: -0.024016 M2: -0.02411"
[4] "Theta= 10 org beta 0.0037 M1: 0.002227 M2: 0.002206"
[5] "Theta= 10 org beta 0.4284 M1: 0.36471 M2: 0.36624"
[6] "Theta= 10 org beta 0.5234 M1: 0.531466 M2: 0.532811"
[1] "Theta = 10 proportion pred. M2 more accurate than M1 = 0.604"
[1] "Theta = 100 variance in introduced frailty v = 0.0101553285954112"
[1] "Theta = 100 variance in estimated random effect= 0.0245876188637748"
[1] "Theta = 100 variance in estimated ind. frailty = 0.00113624325597910"
[1] "Theta= 100 org beta 0.0483 M1: 0.051439 M2: 0.052775"
[2] "Theta= 100 org beta 0.0064 M1: 0.002654 M2: 0.002834"
[3] "Theta= 100 org beta -0.027 M1: -0.025281 M2: -0.025606"
[4] "Theta= 100 org beta 0.0037 M1: 0.005468 M2: 0.005903"
[5] "Theta= 100 org beta 0.4284 M1: 0.467735 M2: 0.48221"
[6] "Theta= 100 org beta 0.5234 M1: 0.673852 M2: 0.683147"
[1] "Theta = 100 proportion pred. M2 more accurate than M1 = 0.564"
Warning message:
Inner loop failed to coverge for iterations 2 3 in: coxpenal.fit(X, Y, strats, offset, init = init, control, weights = weights, >

I don't understand the following:
[1] The variance of the estimated individual frailty values is always much lower than both the original variance and the estimated variance of the random effect. Why is this ? Obviously the variance of the random effect is not calculated directly from the individual frailties after exponentiating them.

[2] Random effects that are not very large are not even picked up by M2, e.g. theta = 10, yields a variance of ~1/10 in the frailty distribution used to set up the data, however, the estimated variance of that random effect equals only 0.003311. Why is frailty not picked up in this case ? Is it just that the variance of the heterogeneity inserted into the data is too small ?

[3] When I compare the predictive abilities of M1 and M2, by calculating the differences, e.g. between predicted 10-year risks, and determine the proportion of more accurate predictions over the complete sample, when using M2 instead of M1, M2 does not perform relevantly better than M1. For theta = 1 or 2 (substantial heterogeneity) M2 performs even worst than M1, while it should be able to take this heterogeneity into account and model M1 cannot do that. Can anyone explain why using M2 does not lead to better predictions in this situation ?

[4] Intuitively it seems to me that in absence of heterogeneity (or with very little heterogeneity, e.g. theta=100) both models should be able to estimate the regression coefficients for the 6 covariates accurately, given enough events. With n=1000, the estimated coefficients are still very inaccurate (see e.g. theta=100, beta 2, 4 and 6) even though we have over 150 events/parameter. Even with lots of events (e.g. 50000) the estimates get more accurately but still are way from perfect. Is there an underlying reason for this slow convergence ?

Is there anyone who can clarify some of these results ? Many thanks in advance for your help. Erik.



Erik Koffijberg, MSc
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: H.Koffijberg@umcutrecht.nl

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. Received on Tue Sep 19 22:14:38 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 19 Sep 2006 - 12:30:05 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.