r - 在R(coxph)中拟合 Gamma 脆弱模型

我有一个带有6个簇的data set，每个簇包含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

有人对调试有想法吗？

最佳答案

替代解决方案:更改方差因子
改变随机效应方差的方法似乎可以解决问题。
例如:

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)