R中的coxph，beta受因子值的影响？答案

【问题标题】：coxph in R, beta affected by value of factor?R中的coxph，beta受因子值的影响？
【发布时间】：2018-12-07 16:50:37
【问题描述】：

我现在正在经营 coxph。我的设置：我有一个参考（无处理），然后是三种不同的处理（A、B 和 C）。我也有 A、B 和 C 的相互作用（例如，用 A 和 B 处理或 A 和 C 处理的样本等......）。我为这些治疗创建了虚拟变量，编码为 1 或 2（1 = 接受治疗，2 = 未接受治疗）。我使用as.factor() 来加载这些变量。

example:
A<-as.factor(Data$A)

我可以按如下方式运行，得到的结果表明接受治疗 B（又名 B = 1）对寿命有益（coef 为正）。这三者在某种程度上都很重要：

> coxph1<-coxph(Surv(Lifespan,Status)~A+B+C
> summary(coxph1)
Call:
coxph(formula = Surv(Life, Status) ~ A + B + C, data = FlyData, 
    method = "efron")

  n= 162, number of events= 140 

     coef exp(coef) se(coef)      z Pr(>|z|)    
A -0.3486    0.7057   0.1761 -1.980 0.047753 *  
B  0.5911    1.8059   0.1787  3.307 0.000944 ***
C -0.6956    0.4988   0.1815 -3.832 0.000127 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

  exp(coef) exp(-coef) lower .95 upper .95
A    0.7057     1.4170    0.4997    0.9966
B    1.8059     0.5537    1.2722    2.5635
C    0.4988     2.0050    0.3494    0.7119

Concordance= 0.822  (se = 0.095 )
Rsquare= 0.227   (max possible= 1 )
Likelihood ratio test= 41.75  on 3 df,   p=5e-09
Wald test            = 41.35  on 3 df,   p=6e-09
Score (logrank) test = 43.6  on 3 df,   p=2e-09

但是，当我运行带有交互项的 coxph 时，我想知道 A:B 或 A:C 等...是否有一些与 A 或 B 不同的交互，我得到以下信息：

> int.coxph <- coxph(Surv(Life, Status)~A*B*C, data=FlyData, method='efron')

警告信息：在 fitter(X, Y, strats, offset, init, control, weights = weights, : Loglik 在变量 1,2,3,4,5,6,7 之前收敛； beta 可能是无限的。

> summary(int.coxph)
Call:
coxph(formula = Surv(Life, Status) ~ A * B * C, data = FlyData, 
    method = "efron")

  n= 162, number of events= 140 

            coef  exp(coef)   se(coef)      z Pr(>|z|)
A      3.987e+01  2.066e+17  4.945e+03  0.008    0.994
B      1.856e+01  1.148e+08  2.472e+03  0.008    0.994
C      3.799e+01  3.144e+16  4.945e+03  0.008    0.994
A:B   -1.964e+01  2.967e-09  2.472e+03 -0.008    0.994
A:C   -3.954e+01  6.737e-18  4.945e+03 -0.008    0.994
B:C   -1.874e+01  7.241e-09  2.472e+03 -0.008    0.994
A:B:C  1.962e+01  3.318e+08  2.472e+03  0.008    0.994

      exp(coef) exp(-coef) lower .95 upper .95
A     2.066e+17  4.841e-18         0       Inf
B     1.148e+08  8.714e-09         0       Inf
C     3.144e+16  3.180e-17         0       Inf
A:B   2.967e-09  3.370e+08         0       Inf
A:C   6.737e-18  1.484e+17         0       Inf
B:C   7.241e-09  1.381e+08         0       Inf
A:B:C 3.318e+08  3.014e-09         0       Inf

Concordance= 0.869  (se = 0.095 )
Rsquare= 0.51   (max possible= 1 )
Likelihood ratio test= 115.6  on 7 df,   p=<2e-16
Wald test            = 9.24  on 7 df,   p=0.2
Score (logrank) test = 73.69  on 7 df,   p=3e-13

所以...这与其他一些问题相似...但是为什么 beta 接近无限？我对这个问题的额外扭曲是，如果我将变量重新编码为 0 或 1（而不是 1 和 2），那么我可以更改交互 coxph() 中的输出。为 coxph 重新编码：

coxph2<-coxph(Surv(Lifespan, Status)~A2+B2+C2))
summary(coxph2)
Call:
coxph(formula = Surv(Life, Status) ~ A2 + B2 + C2, data = FlyData, 
    method = "efron")

  n= 162, number of events= 140 

      coef exp(coef) se(coef)      z Pr(>|z|)    
A2  0.3486    1.4170   0.1761  1.980 0.047753 *  
B2 -0.5911    0.5537   0.1787 -3.307 0.000944 ***
C2  0.6956    2.0050   0.1815  3.832 0.000127 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

   exp(coef) exp(-coef) lower .95 upper .95
A2    1.4170     0.7057    1.0035     2.001
B2    0.5537     1.8059    0.3901     0.786
C2    2.0050     0.4988    1.4048     2.862

Concordance= 0.822  (se = 0.095 )
Rsquare= 0.227   (max possible= 1 )
Likelihood ratio test= 41.75  on 3 df,   p=5e-09
Wald test            = 41.35  on 3 df,   p=6e-09
Score (logrank) test = 43.6  on 3 df,   p=2e-09

正好相反，但交互 coxph 不同...

> full.coxph <- coxph(Surv(Life, Status)~A2*B2*C2, data=FlyData, method='efron')
Warning message:
In fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Loglik converged before variable  2,4,6,7 ; beta may be infinite. 
> summary(full.coxph)
Call:
coxph(formula = Surv(Life, Status) ~ A2 * B2 * C2, data = FlyData, 
    method = "efron")

  n= 162, number of events= 140 

               coef  exp(coef)   se(coef)      z Pr(>|z|)
A2       -7.067e-15  1.000e+00  3.204e-01  0.000    1.000
B2       -2.028e+01  1.558e-09  2.472e+03 -0.008    0.993
C2        9.821e-02  1.103e+00  3.204e-01  0.307    0.759
A2:B2     1.960e+01  3.266e+08  2.472e+03  0.008    0.994
A2:C2    -2.991e-01  7.415e-01  4.475e-01 -0.668    0.504
B2:C2     2.050e+01  7.970e+08  2.472e+03  0.008    0.993
A2:B2:C2 -1.962e+01  3.014e-09  2.472e+03 -0.008    0.994

         exp(coef) exp(-coef) lower .95 upper .95
A2       1.000e+00  1.000e+00    0.5337     1.874
B2       1.558e-09  6.417e+08    0.0000       Inf
C2       1.103e+00  9.065e-01    0.5888     2.067
A2:B2    3.266e+08  3.062e-09    0.0000       Inf
A2:C2    7.415e-01  1.349e+00    0.3085     1.782
B2:C2    7.970e+08  1.255e-09    0.0000       Inf
A2:B2:C2 3.014e-09  3.318e+08    0.0000       Inf

Concordance= 0.869  (se = 0.095 )
Rsquare= 0.51   (max possible= 1 )
Likelihood ratio test= 115.6  on 7 df,   p=<2e-16
Wald test            = 9.24  on 7 df,   p=0.2
Score (logrank) test = 73.69  on 7 df,   p=3e-13

为什么要更改分类变量的数值？ :S 我在这里错过了什么...用非数字变量（“否”和“是”）重新尝试这个结果与使用 0 和 1 相同。例如A 的上 0.95 是“1.874”，B 是“inf”。同样，coxph(Surv()~A+B+C) 为 B 提供负系数，就像上面一样。

【问题讨论】：

可以显示A、B、C的交叉表吗？
它有 163 行长。有没有好的方法来做到这一点？ dropbox.com/s/zb7sxra40virai0/David_FemaleManis-example.R?dl=0 dropbox.com/s/wnxcln78a6motre/… 现在我已经将数据和我的（写得很糟糕的）代码上传到了保管箱。 excel 文件具有以三种不同方式重新编码的 A、B 和 C 虚拟变量。感谢您的关注:)
三个虚拟变量的交叉表应该是3x3
恐怕我听不懂。虚拟变量表示是否对样本进行了处理 A、B 和/或 C。所以列的长度是行数？
嗯...有没有更好的方法来处理或组织这种数据？诚然，我不是统计学家，所以我认为“退化的帽子矩阵”意味着观察到的响应值与预测值几乎无限不同，这是否正确？因此，为什么置信区间接近无穷大？这是因为我将 A、B 和 C（二分结果）作为协变量输入吗？知道如果我使用 0 和 1 与 1 和 2 进行编码，为什么拟合会改变吗？

标签： r interaction survival-analysis cox-regression

【解决方案1】：

您可能（实际上几乎可以肯定）有一个几乎退化的“帽子矩阵”，它是由模型矩阵与该交互作用形成的。你有所有的二阶交互以及三阶相互作用。根据因子中的级别数，完全填充模型矩阵所需的项数可能非常大。我接下来要尝试的是模型中项数稍少的模型。您可以使用 R 的公式接口来删除三阶项，并仅通过以下两种方式之一保留第一项和第二项：

int.coxph <- coxph(Surv(Life, Status)~( A+B+C)^2, data=FlyData, method='efron')

或者：

int.coxph <- coxph(Surv(Life, Status)~ A*B*C - A:B:C, data=FlyData, method='efron')

不确定您是否会通过这种方式获得满足感。您可能没有足够的数据来避免构建 XX^t 矩阵的退化，但如果您的结果没有像上面看到的那样明显地爆炸，那么结果可能是有意义的。另一种更安全的方法是先查看简化模型，然后再添加特定的交互：

 int.coxph.base <- coxph(Surv(Life, Status)~A+B+C,      data=FlyData, method='efron')
int.coxph.intAB <- coxph(Surv(Life, Status)~A+B+C +A:B, data=FlyData, method='efron')

第二个选项还有一个额外的优势，即允许您根据对数似然的变化轻松构建测试，而不是依赖于您在 print.coxph 或 @ 的默认打印输出中看到的不太可靠的 Wald 类型测试987654325@.

【讨论】：

再次感谢您的意见。删除三阶项有助于避免接近无穷大，删除所有交互项并将它们添加回来也会产生避免无穷大的输出。不幸的是，交互项对模型的贡献似乎完全奇怪……即使 A 和 C 明显有害而 B 明显有帮助，它输出的模型声称 A、B 和 C 都对生存有积极贡献。我将样本量增加了 8 倍，以尝试避免违反内在假设：事件数量等……并不是说我完全了解背景统计数据……
从系数中读取含义可能（几乎总是......）误导交互模型。更好的做法是使用 predict 函数并在 newdata 中放入您想要的对比度。

【解决方案2】：

我已经意识到导致我的问题的一个问题：我的生存数据根本没有足够的分辨率。我无法区分交互项的影响。如果我设计我的数据以产生答案，那么我可以获得合理的模型加载输出和有意义的交互项。归根结底，我计划使用所有三种模型类型的组合方法。即：

coxph(Surv(Time, Status)~A+B+C, data=data) #Additive effects
coxph(Surv(Time, Status)~Treatment, data=data) #Base treatment effects
coxph(Surv(Time, Status)~A+B+A:B, data=data) #Test interactions of interest

对加性效应的基本了解可以让您了解协变量在全球范围内对生存的贡献。分析治疗效果（即感兴趣的基本变量）可以让您了解各组是否不同，并从中可以使用加性效应和感兴趣的变量推断模式。

使用 42- 仅调查感兴趣的术语的方法在分析数据时也非常有用。无论我如何处理数据，当您在三方模型中包含所有交互项时，即使我设计为提供信息的数据也会遇到麻烦。但是只使用感兴趣的交互可以增加理解。

我认为这种事后分析需要通过第二次实验进行独立验证，该实验侧重于感兴趣的术语。

【讨论】：