在 R 中指定混合效应 Cox 模型中的方差结构答案

【问题标题】：Specifying variance structure in mixed effects Cox model in R在 R 中指定混合效应 Cox 模型中的方差结构
【发布时间】：2014-03-27 01:38:56
【问题描述】：

我正在使用 coxme 包中的函数 coxme() 在 R 中拟合混合效果 Cox 模型。在我的模型中，我有一个审查的生存时间 $X$、一个协变量 $Z$ 和一个分组变量 $Group$。有一个随机截距和一个随机斜率，即我正在拟合模型 $\lambda (t|Z,b_0,b_1) = \lambda_0(t) e^{\beta Z + b_0 + b_1 Z}$。

我想指定随机效应的协方差矩阵的结构；特别是，我想指定 $b_0$ 和 $b_1$ 是不相关的，因此协方差矩阵是对角线。似乎coxme() 在指定协方差结构方面非常灵活。我认为我应该将矩阵列表传递给此函数的 varlist 选项，但到目前为止我的尝试都失败了，我认为我并不完全理解它应该如何工作。

也可以将自定义方差函数传递给此选项，事实上，其中一个包 vignettes 提供了一些如何完成此操作的示例。但是，在这种情况下，这个过程似乎很乏味并且（我希望）是不必要的。所以我的问题是，如何在coxme() 函数中轻松指定对角协方差矩阵结构？

以下是一些模拟示例数据和指定协方差结构的第一次尝试。我希望我告诉coxme() 使用线性组合 $V = \sigma^2_1 A + \sigma^2_2 B$，$A$ 和 $B$ 定义如下，并且这将有效地适合具有任意对角元素的对角协方差矩阵。

> 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

【问题讨论】：

我不知道该怎么做，但几个月前我使用了coxmepackage 并且有一个无法在这里回答的问题。因此，我向包的开发者 Terry Therneau 发送了一封电子邮件，他回复我的速度非常快而且很有帮助，所以不妨问问他！

标签： r mixed-models cox-regression

【解决方案1】：

经过大量的谷歌搜索和试验，我再次查看了 one of the vignettes 的 coxme 包。我在第 3 节的最后一个要点中找到了该问题的一个可能答案，它说

"- 默认情况下，假定一个完整的协方差矩阵。模型 2 表明，指定独立性的一种简单方法是将效果放在单独的项中。"

所以在数据示例中，可以使用以下命令来获得独立的随机效果：

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

$Group
        Z 
0.2822648

仅显示随机截距和随机斜率的方差，因为假设相关性为零。

一种可能的替代方法是使用phmm 包中的phmm() 函数。由于估计是基于完全似然的，因此这种方法会给出稍微不同的拟合。

【讨论】：