9

我有一个包含 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

有谁知道如何调试？

4

3 回答 3

8

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

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

例如：

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）
于 2012-12-06T19:20:57.593 回答
4

问题出在数据上；x如果x每个cluster. _

通过集群查看数据的分布，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如果除了and之外没有其他变量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：由于数据的结构方式，您无法区分差异。

于 2012-12-05T16:53:07.603 回答
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”强制进行非稀疏计算。

于 2012-12-08T06:11:19.430 回答