[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-help
Subject:    [R] coxme in R underestimates variance of random effect, when random effect is on observation level
From:       Alexander Pate <alexander.pate () manchester ! ac ! uk>
Date:       2018-03-28 17:27:16
Message-ID: 1FBFD3C651590248A8EE6C86B93BE4173318D9AA () MBXP01 ! ds ! man ! ac ! uk
[Download RAW message or body]

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.


[prev in list] [next in list] [prev in thread] [next in thread]