Specifying variance structure in mixed effects Cox model in R

Question

I am fitting a mixed effects Cox model in R using the function coxme() in the coxme package. In my model I have a censored survival time $X$, a single covariate $Z$, and a grouping variable $Group$. There is a random intercept and a random slope, i.e. I am fitting the model $\lambda (t|Z,b_0,b_1) = \lambda_0(t) e^{\beta Z + b_0 + b_1 Z}$.

I would like to specify the structure of the covariance matrix of random effects; in particular, I would like to specify that $b_0$ and $b_1$ are uncorrelated, so that the covariance matrix is diagonal. It seems that coxme() is quite flexible in specifying the covariance structure. I think that I should be passing a list of matrices to the varlist option of this function, but so far my attempts at this have failed and I don't think I fully understand how it is supposed to work.

It is also possible to pass a custom variance function to this option, and in fact one of the package vignettes gives some examples of how this can be done. However, the process seems tedious and (I'm hoping) unnecessary in this case. So my question is, how can I easily specify a diagonal covariance matrix structure in the coxme() function?

Below is some simulated example data and a first attempt at specifying the covariance structure. My hope was that I was telling coxme() to use the linear combination $V = \sigma^2_1 A + \sigma^2_2 B$, with $A$ and $B$ as defined below, and that this would effectively fit a diagonal covariance matrix with arbitrary diagonal elements.

> n = 25  # Size of each cluster
> K = 25  # Number of clusters
> N = n*K  # Total number of observations
> 
> Z = rnorm(n=N, mean=0.5, sd=0.5)  # Covariate
> b0 = rep(rnorm(n=K, mean=0, sd=0.5), each=n)  # Random intercept
> b1 = rep(rnorm(n=K, mean=0, sd=0.5), each=n)  # Random slope
> Group = factor(x=rep(1:K, each=n))
> 
> beta = 2
> eta = beta*Z + b0 + b1*Z
> T = rexp(n=N, rate=exp(eta))  # Exponential failure time, conditional on Z, b0, and b1
> C = runif(n=N, min=0, max=2.5)  # Uniform censoring time to get about 20% censoring
> 
> time = pmin(T,C)  # Censored observation time
> delta = T < C  # Event indicator
> 
> A = matrix(c(1, 0, 0, 0), nrow=2)
> B = matrix(c(0, 0, 0, 1), nrow=2)
> my.covariance = list(A, B)
> fit = coxme(Surv(time, delta) ~ Z + (1 + Z | Group), varlist = my.covariance)
Error in coxme(Surv(time, delta) ~ Z + (1 + Z | Group), varlist = my.covariance) : 
  In random term 1: Mlist cannot have both covariates and grouping

I don't know how to do that but I used the coxmepackage a few months ago and had a question that couldn't be answered here. Therefore, I sent an email to Terry Therneau, the package's developper and he answered me quite fast and usefully so maybe ask him ! — Vincent

Breaking Waves Breaking Waves · Accepted Answer · 2014-04-04T19:10:51

After much googling and experimenting, I took a second look at one of the vignettes for the coxme package. I found one possible answer to the problem in the last bullet point of Section 3, which says

"- By defaut, a full covariance matrix is assumed. Model 2 shows that a simple way to specify independence is to place the effects in separate terms."

So in the data example, the following command can be used to get independent random effects:

> fit = coxme(Surv(time, delta) ~ Z + (1 | Group) + (Z | Group))
> fit$vcoef
$Group
Intercept 
0.1181417 

$Group
        Z 
0.2822648

Only the variances of the random intercept and random slope are shown since the correlation is assumed to be zero.

A possible alternative is to use the phmm() function in the phmm package. This method will give a slightly different fit since estimation is based on the full likelihood.

Specifying variance structure in mixed effects Cox model in R

1 Answers