【问题标题】：Fitting a gamma frailty model in R (coxph)在 R (coxph) 中拟合伽马脆弱模型
【发布时间】：2012-11-20 15:05:27
【问题描述】：

我有一个 data set 有 6 个集群，每个集群包含 48 个（可能被审查，在这种情况下 event = 0）生存时间。 x 列包含一个特定于集群的解释变量。我尝试使用 gamma 脆弱模型来描述该数据，如下所示

 library(survival)

 mod <- coxph(Surv(time, event) ~ 
   x + frailty.gamma(cluster, eps=1e-10, method="em", sparse=0),
              outer.max=1000, iter.max=10000,
              data=data)

这是错误信息：

Error in if (history[2, 3] < (history[1, 3] + 1)) theta <- mean(history[1:2,  : 
  missing value where TRUE/FALSE needed

有人知道如何调试吗？

【问题讨论】：

frailty 的帮助页面显示：“对于 Cox 模型，coxme 包已取代此方法。”
感谢您的共同关注@DWin。然而，coxme 适合高斯弱点，而不是伽马弱点。
没时间进一步研究这个，但要看的地方在"cfun"，由fraily.gamma():FG <- frailty.gamma(data$cluster, eps=1e-10, method="em", sparse=TRUE); FG["cfun"]制作。
@tim riffe：如果你有时间，我很想知道更多；-)
投了赞成票，这样你就可以收回一些巨额赏金......

标签： r error-handling

【解决方案1】：

替代解决方案：更改方差因子

改变随机效应方差的方法，似乎可以解决问题。

例如：

mod.aic <- coxph(Surv(time, event) ~ 
               x + frailty.gamma(cluster, eps=1e-10, method="aic", sparse=0),
             outer.max=1000, iter.max=10000,
             data=dat)

plot(survfit(mod.aic), col=4)

ddd hoc 解决方案：如果我们删除一个集群就可以工作

也许这不能完全回答你的问题，但是当我删除任何集群时，例如：

par(mfrow=c(2,3))
res <- sapply( 1:6 , function(x) {
                      mod <- 
                        coxph(Surv(time, event) ~ 
                        x + frailty.gamma(cluster, eps=1e-10, method="em", sparse=0),
                        outer.max=1000, iter.max=10000,
                        data=subset(dat,cluster != x)
                        )
                     plot(survfit(mod), col=4,main= paste ('cluster', x, 'is removed'))
                     legend(10,1,mod$iter)
              })

coxph 收敛，我对所有样本都有相同的结果。

我没有足够的关于您的数据的信息来进行进一步分析，但我尝试在不同的集群之间进行一些比较。

library(ggplot2
qplot(data = dat, x=time , y = x , facets= event~cluster)

我注意到 3 组：

集群 1,3,5：事件均匀分布
集群 2 ,4 ：仅在小时候发生的事件。
第 6 组：令人惊叹的一个（唯一事件 1）

【讨论】：

感谢您指出这一点。然而，“method=aic”不是标准的，与“method=em”相比，我不知道它的性能如何，“method=em”已被证明与 EM 算法给出相同的估计......
所以即使有更多解释也没有办法使用aic来选择方差？
对不起，我没有得到这个问题...更多解释什么？顺便说一句，无论如何，为你的答案 +1
如果我看到模拟结果显示它表现良好，我会相信 'method=aic'。但是，对我来说，它看起来像是一种临时方法……我从未见过有人在这种情况下使用这种方法。
谢谢！删除一个集群确实是一个好主意。如果我可以投票两次，我会的！我会尽快发布 Terry Therneau（coxph 的作者）的答案。

【解决方案2】：

问题出在数据上；如果每个cluster 中的所有x 都相同，则您无法将特定于集群的效果与x 分开。

通过集群查看x 在您的数据中的分布，我们可以看到：

table(data$x,data$cluster)
     1  2  3  4  5  6
  0  0 48  0 48 48  0
  1 48  0 48  0  0 48

我认为集群特定的解释变量是什么意思。这在任何模型中都将是一个问题，因为x 与cluster 共线（我认为是这个词）。甚至尝试最基本的模型：

data$cluster<-as.factor(data$cluster)
mod <- coxph(Surv(time, event) ~ x + cluster, data=data)

Warning message:
In coxph(Surv(time, event) ~ x + cluster, data = data) :
  X matrix deemed to be singular; variable 5

矩阵是奇异的，因为无法区分cluster 和x 的效果。

如果你除了cluster和x之外没有其他变量，那么你真正能做的就是单独运行集群的效果：

data$cluster<-as.factor(data$cluster)
coxph(Surv(time, event) ~ cluster,data=data)

Call:
coxph(formula = Surv(time, event) ~ cluster, data = data)


          coef exp(coef) se(coef)     z       p
cluster2 1.070      2.92    0.382  2.80 5.1e-03
cluster3 0.499      1.65    0.384  1.30 1.9e-01
cluster4 1.705      5.50    0.365  4.68 2.9e-06
cluster5 2.058      7.83    0.370  5.56 2.7e-08
cluster6 4.415     82.69    0.399 11.06 0.0e+00

考虑cluster1和cluster6具有相同的值x，它们之间的危险比是83。也许cluster6不同，也许x在cluster6中的行为不同：你由于数据的结构方式，无法区分。

【讨论】：

非常感谢您的回答。但是，将“集群效应”视为随机应该允许拟合模型。顺便说一句，将“frailty.gamma(cluster, eps=1e-10, method="em", sparse=0)" 更改为 "frailty.gaussian(cluster, eps=1e-10, method="reml", sparse =0)" 给出一些东西。我真的认为问题出在我的 R 代码上。无论如何，谢谢你的洞察力。
好点。到目前为止，我已经做了一些调试，似乎在survival:::coxpenal.fit 内的一个拟合函数中，loglik 返回为NaN，表明在 C 代码中某些内容被零除。我怀疑除以零本身是由于数据结构造成的“完美拟合”，但这只是猜测。

【解决方案3】：

这是 Terry Therneau（coxph 的作者）给我的答案。

我查看了您的数据：

> table(x, cluster)
     1  2  3  4  5  6
  0  0 48  0 48 48  0
  1 48  0 48  0  0 48

您的协变量“x”完全由集群变量预测。如果您适合固定效应模型： coxph(Surv(time, event) ~ factor(cluster) +x)

然后“x”变量被声明为多余的。当随机效应的方差为足够大，当方差足够大时，伽马模型中也会发生同样的情况。您的模型接近此限制，并且解决方案失败。如手册页所述，现在首选 coxme 函数。

最后，您的特定错误消息是由“稀疏”的无效值引起的。我将在程序中添加一个检查。您可能希望“sparse=10”强制进行非稀疏计算。

【讨论】：