随机斜率 Cox 比例危险答案

【问题标题】：Random slopes Cox Proportional Hazards随机斜率 Cox 比例危险
【发布时间】：2018-12-17 04:09:58
【问题描述】：

我一直在尝试使用coxme 为我的模型中的每个协变量提取随机斜率。

library (coxme)
Start <- runif(5000, 1985, 2015)
Stop <- Start + runif(5000, 2, 10)
S <- data.frame (
  X1 <- runif(5000, 5.0, 7.5),
  X2 <- runif(5000, 5.0, 7.5),
  D <- rbinom(5000, 1, 0.8), 
  Letters <- sample(LETTERS, 5000, replace = TRUE), 
  Start <- Start,
  Stop <- Stop
  )

S_ind1 <- Surv (time = S$Start, time2 = S$Stop, event = S$D)
a <- coxme (S_ind1 ~ X1 + X2 + (X1 + X2|Letters), data = S)

我得到的只是：

Error in gchol(kfun(theta, varlist, vparm, ntheta, ncoef)) : 
NA/NaN/Inf in foreign function call (arg 5)
In addition: Warning messages:
1: In sqrt(xvar * zvar) : NaNs produced
2: In sqrt(xvar * zvar) : NaNs produced

当使用我自己的数据时，我经常得到：

Error in coxme.fit(X, Y, strats, offset, init, control, weights = weights,  : 
  'Calloc' could not allocate memory (56076596 of 8 bytes)

是否可以使用coxme 包含随机斜率？

如果没有，是否有其他使用其他包的替代方法？

【问题讨论】：

标签： cox-regression random-effects

【解决方案1】：

coxme 包的作者 Terry Therneau 通过电子邮件回答 - 他让我在这里发布。

下面是我对您的示例的重写，删除了 Surv 间接并在 data.frame 调用中使用“=”（我有点惊讶

library (coxme)

set.seed(1953)
time1 <- runif(5000, 1985, 2015)
time2 <- time1 + runif(5000, 2, 10)
test <- data.frame (
  x1 = runif(5000, 5.0, 7.5),
  x2 = runif(5000, 5.0, 7.5),
  death = rbinom(5000, 1, 0.8), 
  letters = sample(LETTERS, 5000, replace = TRUE), 
  time1 = time1,
  time2 = time2)

fit1 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1|letters), data=test)
fit2 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1+x1 | letters), test)
fit3 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1+x2 | letters), test)
fit4 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1+ x1 + x2 | letters),
              data=test, vinit= c(1e-6, 1e-8, 1e-8))

*1。所有模型都可以工作到 fit4。

我发现您的模型令人担忧，因为它有一个随机斜率但没有随机截距，就像所有通过原点的回归都让我担心一样：我很难对结果进行解释。尽管 lme 默认情况下将截距项放入，但 coxme 没有。
我希望 fit4 会起作用，并且也许有更好的初始估计会起作用。 coxme 的底层代码是我在所有生存工作中遇到的最难的最大化问题，从某种意义上说，最大化器很容易丢失并且永远找不到它的方式。这是一个有时需要手动操作的功能，通过有限的迭代次数和/或开始估计。我希望不是这样，我有一些长期计划通过添加一个基于 MCMC 的最大化器来改进这一点，理论上它永远不会丢失，但会以更长的计算时间为代价。
如果任何方差过于接近零，则 sqrt() 消息往往会作为舍入误差的函数出现。当然，在您的测试用例中，实际 MLE 的方差为 0。发生这种情况时，我通常会通过对一系列固定方差（vfixed 参数）进行拟合来直接检查方差是否为零。如果随着方差达到 1e-6 或更小的值，可能性是恒定的或增加的，那么我假设 MLE 为零并从模型中删除该随机项。特里 T.*

【讨论】：