在R(coxph)中拟合伽玛脆弱模型

Question

我有一个包含6个集群的数据集,每个集群包含48个(在这种情况下可能被审查event = 0)生存时间.该x列包含特定于集群的解释变量.我尝试使用伽玛脆弱模型描述该数据,如下所示

 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

有没有人知道如何调试？

Answer 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. 集群1,3,5:事件均匀分布
2. 集群2,4:仅适用于小时间的事件.
3. 集群6:惊人的一个(只有事件1)

Answer 2

这是 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”来强制非稀疏计算。