随机斜率 Cox 比例风险

Random slopes Cox Proportional Hazards

我一直在尝试使用 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 来包含随机坡度?

如果没有,是否有使用其他包的其他替代方案?

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

下面是我对您的示例的重写,删除了 Surv 间接寻址并在 data.frame 调用中使用了“=”(我有点惊讶 <- 在该上下文中有效) ,并添加 set.seed 以便示例可重现。

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.

  1. 我发现你的模型令人担忧,因为它有一个随机斜率但没有随机截距,就像所有通过原点的回归让我担心的一样:我很难解释结果。尽管 lme 默认包含截距项,但 coxme 没有。

  2. 我曾希望 fit4 能奏效,也许有了更好的初始估计它就会奏效。 coxme 的底层代码是我在所有生存工作中遇到的最难的最大化问题,从最大化器很容易迷路并且永远找不到它的方式的意义上来说很难。这是一个有时需要手动操作的函数,通过有限的迭代计数 and/or 开始估计。我希望不是这样,我有一些长期计划通过添加一个替代的基于 MCMC 的最大化器来改进这一点,理论上它永远不会丢失但以更长的计算时间为代价。

  3. 如果任何方差太接近于零,则 sqrt() 消息往往会作为舍入误差的函数出现。当然,在您的测试用例中,实际 MLE 的方差为 0。发生这种情况时,我通常会通过使用一系列固定方差(vfixed 参数)进行拟合来直接检查零方差。如果可能性不变或随着方差值达到 1e-6 或更小而增加,那么我假设 MLE 为零并从模型中删除该随机项。 特里 T.*