R 中的 coxph，受因子值影响的 beta？
coxph in R, beta affected by value of factor?
我现在正在 运行宁 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）对寿命有益（系数为正）。这三者在某种程度上都很重要：
> 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 我在这里遗漏了什么......用非数字变量（"no" 和 "yes"）重新尝试给出与使用 0 和 1 相同的结果。例如A 的上限 .95 是“1.874”，B 的上限是 "inf"。类似地，coxph(Surv()~A+B+C) 给出 B 的负系数，就像上面的那样。
您可能（事实上几乎可以肯定）有一个近乎退化的 "hat matrix" 这是模型矩阵与该相互作用形成的。你有所有的二阶相互作用以及
三阶相互作用。根据因子中的水平数，完全填充模型矩阵所需的项数可能非常大。接下来我要尝试的是模型中的项略少的模型。您可以使用 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')
第二个选项有一个额外的优势，即允许您根据对数似然的变化轻松构建测试，而不是依赖于您在 [=13 的默认打印输出中看到的不太可靠的 Wald 型测试=] 或 summary.coxph.
我已经意识到导致我出现问题的原因之一：我的生存数据根本没有足够的分辨率。我无法区分交互项的影响。如果我设计我的数据来产生答案，那么我可以获得合理的模型加载输出和有意义的交互项。最后，我计划使用所有三种模型类型的组合方法。即：
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- 仅调查感兴趣的术语的方法在分析数据时也非常有用。无论我如何处理数据，当您将所有交互项都包含在三方模型中时，即使是我设计为提供信息的数据也会遇到麻烦。但是只使用感兴趣的交互可以增加理解。
我想这种 post-hoc 分析需要来自关注感兴趣术语的第二个实验的独立验证。
我现在正在 运行宁 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）对寿命有益（系数为正）。这三者在某种程度上都很重要：
> 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 我在这里遗漏了什么......用非数字变量（"no" 和 "yes"）重新尝试给出与使用 0 和 1 相同的结果。例如A 的上限 .95 是“1.874”，B 的上限是 "inf"。类似地，coxph(Surv()~A+B+C) 给出 B 的负系数，就像上面的那样。
您可能（事实上几乎可以肯定）有一个近乎退化的 "hat matrix" 这是模型矩阵与该相互作用形成的。你有所有的二阶相互作用以及 三阶相互作用。根据因子中的水平数，完全填充模型矩阵所需的项数可能非常大。接下来我要尝试的是模型中的项略少的模型。您可以使用 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')
第二个选项有一个额外的优势，即允许您根据对数似然的变化轻松构建测试，而不是依赖于您在 [=13 的默认打印输出中看到的不太可靠的 Wald 型测试=] 或 summary.coxph.
我已经意识到导致我出现问题的原因之一：我的生存数据根本没有足够的分辨率。我无法区分交互项的影响。如果我设计我的数据来产生答案，那么我可以获得合理的模型加载输出和有意义的交互项。最后，我计划使用所有三种模型类型的组合方法。即：
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- 仅调查感兴趣的术语的方法在分析数据时也非常有用。无论我如何处理数据，当您将所有交互项都包含在三方模型中时，即使是我设计为提供信息的数据也会遇到麻烦。但是只使用感兴趣的交互可以增加理解。
我想这种 post-hoc 分析需要来自关注感兴趣术语的第二个实验的独立验证。