【问题标题】：Getting P-Values of Zero in Cox Regression: R在 Cox 回归中获得零 P 值：R
【发布时间】：2020-02-25 01:17:50
【问题描述】：

我是一名在 R 中进行基因表达生存分析的学生。我有 249 名患者的表达数据，我使用 6,000 个基因以及它们的无事件生存时间和生命状态作为响应变量。当我尝试在我的数据集上运行 Cox 回归时，我得到了非常奇怪的结果（p 值为 0.00 和奇怪的风险比）。我已经多次检查了我的代码，但我无法发现我的错误（当我早些时候尝试只使用一个基因时，它工作得很好，但是当我尝试使用'.'函数测试多个基因时，我不是得到更好的结果）。我非常感谢任何帮助，并附上了我的代码和输出！如果需要更多信息，请告诉我。

library(survival)
options(expressions = 5e5)
firstSplitData <- read.delim("/Users/menon/OneDrive/Desktop/csrsef files/FirstSplitDataFrame.txt")
firstInitialData <- data.frame(firstSplitData)
firstEventFreeTime <- firstInitialData[ , c("EFST")] 
firstVitalStatus <- firstInitialData[, c("Status")]
#create a temporary object to use in the final object in order to be able to use '.'
temporaryObj <- Surv(as.numeric(firstEventFreeTime), firstVitalStatus == 2)
firstFinalData <- data.frame(SurvObj = temporaryObj)
#bind the two together for the final data 
firstFinalData <- cbind(firstFinalData, firstInitialData[, 2:ncol(firstInitialData)])
#create final cox model
firstCox <- coxph(SurvObj ~ ., data =  firstFinalData)
summary(firstCox)$coefficients

这是我的（部分）输出：

> summary(firstCox)$coefficients
                     coef     exp(coef)     se(coef)             z      Pr(>|z|)
EFST         3.644083e-03  1.003651e+00 0.0001340611    27.1822581 1.052851e-162
Status      -2.926090e+00  5.360625e-02 0.3182658189    -9.1938542  3.790122e-20
AADACL3      1.502153e+02  1.728460e+65 0.3665374081   409.8224582  0.000000e+00
AADACL4      5.857192e+01  2.738174e+25 0.3681708023   159.0889828  0.000000e+00
ACADM        2.455978e+02 4.589695e+106 0.2175220391  1129.0710334  0.000000e+00
ACAP3        4.093913e+02 6.256964e+177 0.2756635268  1485.1121632  0.000000e+00
ACOT11       1.940976e+01  2.688751e+08 0.3251033140    59.7033512  0.000000e+00
ACOT7       -2.841794e+02 3.823403e-124 0.3139848504  -905.0736377  0.000000e+00
ACTB        -5.562202e+01  6.976896e-25 0.3173481100  -175.2713234  0.000000e+00
ACTL8       -4.017414e+02 3.356676e-175 0.3435128215 -1169.5093020  0.000000e+00
ACTRT2      -7.613568e+01  8.603881e-34 0.2861088372  -266.1074036  0.000000e+00
ADC         -1.244476e+02  8.976070e-55 0.3201452217  -388.7223972  0.000000e+00
ADPRHL2      4.887427e+01  1.681998e+21 0.2895110526   168.8165913  0.000000e+00
AGMAT        7.266946e+02           Inf 0.4295874196  1691.6104194  0.000000e+00
AGO1         3.352041e+02 3.778188e+145 0.2633158947  1273.0111995  0.000000e+00
...

这是dput(firstFinalData[1:10, 1:10]) 产生的结果：

structure(list(SurvObj = structure(c(444, 5553, 5296, 922, 205, 
47, 401, 245, 263, 5564, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0), .Dim = c(10L, 
2L), .Dimnames = list(NULL, c("time", "status")), type = "right", class = "Surv"), 
    EFST = c(444L, 5553L, 5296L, 922L, 205L, 47L, 401L, 245L, 
    263L, 5564L), Status = c(2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 1L), AADACL3 = c(5.52132, 5.64712, 5.45876, 5.71481, 
    5.1269, 5.88764, 5.08912, 4.91729, 5.65387, 5.59824), AADACL4 = c(5.17251, 
    5.41843, 5.10969, 5.23402, 4.60353, 5.70923, 5.02245, 5.1466, 
    4.8355, 4.83986), ACADM = c(7.47834, 7.43494, 7.91155, 7.86337, 
    8.39009, 6.16251, 7.83793, 7.71742, 6.98061, 7.78087), ACAP3 = c(7.80589, 
    8.00354, 7.75014, 7.61566, 7.55267, 7.9449, 7.20561, 7.99776, 
    7.72778, 7.43355), ACOT11 = c(6.75915, 6.30386, 6.38214, 
    6.54392, 6.64743, 6.78981, 6.42641, 6.58761, 6.66693, 6.53731
    ), ACOT7 = c(8.11807, 8.38011, 7.8349, 8.43645, 8.11502, 
    8.0109, 7.6866, 8.55327, 8.17004, 7.44455), ACTB = c(10.8227, 
    11.4556, 11.4216, 11.332, 10.9536, 9.83797, 11.2352, 11.5006, 
    11.1817, 10.895)), row.names = c(NA, 10L), class = "data.frame")

非常感谢！

编辑：

我在运行firstCox <- coxph(SurvObj ~ ., data = firstFinalData) 时也收到此警告消息：

In fitter(X, Y, istrat, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge

【问题讨论】：

你的命令太多了。导入数据后，您可以运行模型。请看下面我的回答。将Surv() 函数直接放在对coxph 的调用中。 R 知道不要在. 中包含结果。
哦，我刚刚意识到您可能想要运行多个回归，每个基因一个。那是对的吗？你用coxph(SurvObj ~ ., data = firstFinalData) 命令把我扔了，因为. 表示你想运行一个多变量模型，这在这里是不可能的。

标签： r bioinformatics survival-analysis cox-regression

【解决方案1】：

如果您想使用单个预测变量执行多个 Cox 回归模型，您可以使用以下代码和您发布的示例数据。首先我删除第一列中的生存对象。

myData <- finalData[,-1]

library(survival)
firstCox <- co

coxph(Surv(EFST, Status) ~ ., data =  myData)

这会返回一个警告（预测变量过多）

Warning message:
In fitter(X, Y, istrat, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge

要运行多个单变量模型，首先创建单变量公式列表：

formulas <- sapply(names(myData)[3:9], function(x) as.formula(paste('Surv(EFST, Status) ~ ',x)))

使用coxph 函数创建模型列表：

models <- lapply(formulas, function(x) coxph(x, data=myData))

提取风险比 (exp(coef)) 和 95% 置信区间：

res <- lapply(models, function(x) return(cbind(HR=exp(coef(x)), exp(confint(x)), Pval=coef(summary(x))[5])))
res

$AADACL3
               HR       2.5 %   97.5 %     Pval
AADACL3 0.1858129 0.008579879 4.024119 0.283442

$AADACL4
               HR      2.5 %   97.5 %      Pval
AADACL4 0.8481017 0.02748128 26.17333 0.9249839
...

【讨论】：

抱歉所有问题 - 我还想使用此代码提取 p 值，并尝试将 summary() 应用于我的模型列表对象以获得 p 值，但我收到的是这样的东西： GIPC2 19 coxph list GJA9.MYCBP 19 coxph list GLTPD1 19 coxph list ... 如果这不是问这个问题的合适地方，我会很乐意移动它。非常感谢！
只需在exp(confint(x)) 之后添加：, Pval=coef(summary(x))[5])。
您好，爱德华，再次感谢您！

【解决方案2】：

除了前两个系数（EFST 和Status）外，所有其他基因的系数要么极小要么极大，导致非常大的负/正t-统计，这解释了您看到的 p 值。

我不确定我是否理解你在做什么。对 249 名患者数据中的 6,000 个基因进行回归难道不意味着您的参数比观察值多得多吗？

在这种情况下，您会遇到可以解释参数估计的过度拟合问题。

【讨论】：

您好，感谢您的回复 - 我有 249 名患者的 23,000 个基因的基因表达数据，我正在尝试对他们的 6,000/23,000 个基因进行 Cox 回归（测试每个基因单独）通过使用无事件生存时间和生命状态。我对自己想做的事情有根本的误解吗？我的假设是，由于我分别测试每个基因，我不会遇到任何问题。感谢您的帮助！
嗨安基塔。仍然不确定我是否完全理解你在做什么。除此之外，单独测试每个基因通常不是一个好主意。您会遇到多个测试问题，并且通常无法（有意义地）比较结果。毕竟，你工作的重点不就是评估不同基因对患者生存的影响吗？其次，在公式表达式SurvObj ~ . 中，右侧扩展为左侧未包含的所有变量/列。所以你实际上是在一个模型中回归所有 6000 个基因（不是单独的）
我正在尝试寻找与生存相关的生物标志物，因此我想针对无事件生存时间和生命状态分别测试每个基因（我已经用我拥有的基因做过一次）假设并且它有效，但我想对所有 23,000 个基因都这样做）。但是，除了使用“。”之外，我不确定如何在 R 中同时执行此操作。功能，但我现在意识到这需要一起测试基因，而不是单独测试。
不，我认为单独测试基因在统计上是不合理/有效的，因为您会忽略所有其他基因的所有混杂效应。这是一个非常滑的斜坡，最多可能会导致无法解释的结果，最坏的情况可能是采摘樱桃。
好的，我明白了。我的原始程序/项目基于我看到的这篇文章 - biostars.org/p/344233 - 但我意识到我现在可能不得不改变一些事情。感谢您的建议！

【解决方案3】：

不要在数据框中包含 Surv() 对象。

firstFinalData <- firstFinalData[,-1]
firstCox <- coxph(Surv(EFST, Status) ~ ., data =  firstFinalData)

它应该可以工作（编辑：在较少数量的变量上）。

正如 Maurits Evers 所说，仅对 249 个受试者运行一个包含 6,000 个预测变量（基因）的模型将导致收敛问题。考虑减少基因数量（或获得更多患者！）

【讨论】：

不确定这是否可行。我可能误解了，但 OP 似乎有 6,000 个回归变量（基因），但只有 249 个观察值（患者）。这对我来说听起来像是一个不适定的问题，如果没有某种形式的正则化（或贝叶斯框架中参数的一些信息/正则化先验），这将导致过度拟合。
我同意。我只是指出了多余的命令，特别是 Surv 对象不应包含在数据框中，因为 R 将尝试在 . 中包含所有其他变量，以及 OP dputted 包含的数据变量太多。看看输出。 EFST 和 Status 包含在预测变量中！
感谢您的回复。即使我单独测试基因，这是否仅仅因为我正在对 249 名患者测试 6,000 个基因而引起问题？我不确定我是否理解这一点，因为我能够在单独和一个一个地测试基因时成功获得 Cox 回归结果（而不是让程序同时完成所有这些）。
这不是编码问题。我的回答解决了您的编码问题。我建议您在另一篇文章或其他地方提出建模问题，或咨询合格的统计学家。
@Edward “看看输出。EFST 和 Status 作为预测变量包括在内！” 是的，你在修复 OP 语法方面是正确的。这似乎是一个多因素问题（编码语法和统计信息）。