Hello,
I have a question concerning fitting a cox model with a random intercept, also known as a frailty model. I am using both the coxme package, and the frailty statement in coxph. Often 'shared' frailty models are implemented in practice, to group people who are from a cluster to account for homogeneity in outcomes for people from the same cluster. I am more interested in the classic frailty model, 'Early frailty models incorporated subject-specific random effects to account for unmeasured subject characteristics that influenced the hazard of the occurrence of the outcome'. This is because I have data where I would like to estimate the amount of variation between patients risks that has not been accounted for by adjusting for the variables that I have. I initially ran some simulations to check that I could estimate the variance of such random effects consistently with no bias. When the random effect is applied on the group level (shared frailty model), the variance of the random effect can be calculated. When the random effect is on an observation level, it is consistently underestimated. This happens when using both the coxme<https://cran.r-project.org/web/packages/coxme/coxme.pdf> package, and the and frailty<https://www.rdocumentation.org/packages/survival/versions/2.11-4/topics/frailty> statement in coxph. Questions: 1) Why is variance of the random effect underestimating when the random effect is on an individual level, is this a mathematical problem or a coding issue? 2) Is there any literature/R packages that look specifically at random effects applied on the observation level? Many thanks for any solutions/or references which may be helpful. Reproducible example here (using coxme): setwd("/mnt/ja01-home01/mbrxsap3/phd_risk/R/p4_run_analysis/") library(coxme) library(survival) ### Create data with a group level random effect simulWeib.group <- function(N, lambda, rho, beta1, beta2, beta3, beta4, rateC, sigma, M) { # covariate --> N Bernoulli trials x1 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5)) x2 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5)) x3 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5)) x4 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5)) # Now create random effect stuff # Create vector of groups re.group <- sample(x=1:M, size=N, replace=TRUE, prob=rep(1/M, M)) # Create vector of r.e effects re.effect <- rnorm(M,0,sigma) # Now create a vector which assigns the re.effect depending on the group re.group.effect <- re.effect[re.group] # Weibull latent event times v <- runif(n=N) Tlat <- round((- log(v) / (lambda * exp(x1 * beta1 + x2 * beta2 + x3 * beta3 + x4 * beta4 + re.group.effect)))^(1 / rho)) # censoring times #C <-rep(100000,N) C <- rexp(n=N, rate=rateC) # follow-up times and event indicators time <- pmin(Tlat, C) #status <- as.numeric(rep(1,N)) status <- as.numeric(Tlat <= C) # data set data.frame(id=1:N, time=time, status=status, x1 = x1, x2 = x2, x3 = x3, x4 = x4, group=re.group) } ### Create data with an individual level random effect simulWeib <- function(N, lambda, rho, beta1, beta2, beta3, beta4, rateC, sigma) { # covariate --> N Bernoulli trials x1 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5)) x2 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5)) x3 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5)) x4 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5)) # Now create random effect stuff # Create one vector of length N, all drawn from same normal distribution rand.effect <- rnorm(N,0,sigma) # Weibull latent event times v <- runif(n=N) Tlat <- round((- log(v) / (lambda * exp(x1 * beta1 + x2 * beta2 + x3 * beta3 + x4 * beta4 + rand.effect)))^(1 / rho)) # censoring times #C <-rep(100000,N) C <- rexp(n=N, rate=rateC) # follow-up times and event indicators time <- pmin(Tlat, C) #status <- as.numeric(rep(1,N)) status <- as.numeric(Tlat <= C) # data set data.frame(id=1:N, time=time, status=status, x1 = x1, x2 = x2, x3 = x3, x4 = x4) } set.seed(101) ### Individual random effects (frailty). sd.2<-rep(0,50) for (i in 1:50){ data2<-simulWeib(25000,lambda=0.0001,rho=2,beta1=0.33,beta2=5,beta3=0.25,beta4=0,rateC=0.0000000001, sigma = 0.25) data2$id<-as.factor(data2$id) fit.cox2<-coxme(Surv(time,status) ~ x1 + x2 + x3 + (1 | id), data=data2) sd.2[i]<-sqrt(as.numeric(fit.cox2$vcoef)) print(i) } print("model 2 done") ### Same as previous example, but patients are grouped sd.10<-rep(0,50) for (i in 1:50){ data10<-simulWeib.group(25000,lambda=0.0001,rho=2,beta1=0.33,beta2=5,beta3=0.25,beta4=0,rateC=0.0000000001, sigma = 0.25, M=40) data10$group<-as.factor(data10$group) fit.cox10<-coxme(Surv(time,status) ~ x1 + x2 + x3 + (1 | group), data=data10) sd.10[i]<-sqrt(as.numeric(fit.cox10$vcoef)) print(i) } print("model 10 done") PhD student Health e-Research Centre � Farr Institute Rm 2.006 � Vaughan House � Portsmouth Street � University of Manchester � M13 9GB [[alternative HTML version deleted]]
______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
