在 Cox 回归中获得零的 P 值：R
Getting P-Values of Zero in Cox Regression: R
我是一名在 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
除了前两个系数（EFST和Status），其他所有基因的系数要么极小要么极大，导致非常大negative/positive t-statistics，它解释了您看到的 p 值。
我不太确定我理解你在做什么。对 249 名患者数据中的 6,000 个基因进行回归不会意味着您拥有的参数比观察值多得多吗？
在这种情况下，您 运行 陷入了可以解释参数估计的过度拟合问题。
不要在数据框中包含 Surv() 对象。
firstFinalData <- firstFinalData[,-1]
firstCox <- coxph(Surv(EFST, Status) ~ ., data = firstFinalData)
它应该有效（编辑：在较少数量的变量上）。
正如 Maurits Evers 所说，运行 一个只有 249 个受试者的具有 6,000 个预测变量（基因）的模型将导致收敛问题。考虑减少基因数量（或获得更多患者！）
如果您想使用单个预测变量执行多个 Cox 回归模型，您可以使用您发布的示例数据使用以下代码。首先，我删除了第一列中的生存对象。
myData <- finalData[,-1]
library(survival)
firstCox <- co
coxph(Surv(EFST, Status) ~ ., data = myData)
这 returns 如前所述的警告（预测变量太多）
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
...
我是一名在 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
除了前两个系数（EFST和Status），其他所有基因的系数要么极小要么极大，导致非常大negative/positive t-statistics，它解释了您看到的 p 值。
我不太确定我理解你在做什么。对 249 名患者数据中的 6,000 个基因进行回归不会意味着您拥有的参数比观察值多得多吗？
在这种情况下，您 运行 陷入了可以解释参数估计的过度拟合问题。
不要在数据框中包含 Surv() 对象。
firstFinalData <- firstFinalData[,-1]
firstCox <- coxph(Surv(EFST, Status) ~ ., data = firstFinalData)
它应该有效（编辑：在较少数量的变量上）。
正如 Maurits Evers 所说，运行 一个只有 249 个受试者的具有 6,000 个预测变量（基因）的模型将导致收敛问题。考虑减少基因数量（或获得更多患者！）
如果您想使用单个预测变量执行多个 Cox 回归模型，您可以使用您发布的示例数据使用以下代码。首先，我删除了第一列中的生存对象。
myData <- finalData[,-1]
library(survival)
firstCox <- co
coxph(Surv(EFST, Status) ~ ., data = myData)
这 returns 如前所述的警告（预测变量太多）
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
...