colMeans(hazard) 中的错误：“x”必须是至少二维的数组

Question

我正在尝试从 R 中的 Cox PH 模型中开发估计生存曲线。由于它是一个多状态模型，我添加了一个新的数据参数，其中我保持协变量不变，除了我正在制作的 1 之外生存曲线.

newdata_LNI <- data.frame( PathGGS = rep("7",3),
                           PreTxPSA = rep(median(Cox_PH_Data$PreTxPSA)),
                           SVI_Binary = rep("0", 3),
                           LNI_Binary = rep(levels(Cox_PH_Data$LNI_Binary)))

我的考克斯模型编码如下，

Final_Cox_Model <- coxph(SurvObj ~   PreTxPSA + PathGGS  + LNI_Binary + SVI_Binary, 
                         data = Cox_PH_Data, id = Cox_PH_Data$`Sample ID`)

我的“SurvObj”的格式与时间和事件正确，并且我的变量都是相同的类型和级别。

我已经尝试了一切方法，但此错误仍然存在。上周更新 Rstudio 后我可以运行它，但现在这个错误不断发生。 这是

的输出

dput(head(Cox_PH_Data, 20))

structure(list(`Sample ID` = c("PCA0001", "PCA0002", "PCA0003", 
"PCA0004", "PCA0005", "PCA0006", "PCA0007", "PCA0008", "PCA0009", 
"PCA0010", "PCA0011", "PCA0012", "PCA0013", "PCA0014", "PCA0015", 
"PCA0016", "PCA0017", "PCA0018", "PCA0019", "PCA0020"), Race = structure(c(2L, 
2L, 2L, 4L, 2L, 4L, 2L, 2L, 2L, 4L, 2L, 2L, 2L, 4L, 2L, 4L, 2L, 
4L, 4L, 4L), levels = c("Asian", "Black", "Unknown", "White"), class = "factor"), 
    Type = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("MET", 
    "PRIMARY"), class = "factor"), DxAge = c(72.77089, 58.09296, 
    67.93575, 68.92686, 64.69953, 56.63913, 56.76507, 64.23682, 
    56.5789, 49.86554, 52.1161, 66.17527, 54.20513, 60.349, 58.71446, 
    67.35531, 55.91632, 55.23458, 68.11918, 59.88629), PreDxBxPSA = c(132, 
    8.1, 4.6, 27.5, 8.7, 15.7, 14, 10.4, 14.6, 5.2, 2, 3.3, 9.4, 
    4.2, 2.9, 12, 9.32, 8.5, 6.6, 3.8), BxGGS = structure(c(2L, 
    3L, 3L, 1L, 2L, 2L, 3L, 3L, 2L, 2L, 2L, 2L, 3L, 2L, 2L, 2L, 
    3L, 2L, 3L, 3L), levels = c("5", "6", "7", "8", "9"), class = "factor"), 
    PreTxPSA = c(43.9, 43.2, 11.3, 11.8, 5.9, 8.2, 3.8, 10.4, 
    12.9, 5.2, 6.7, 3.26, 9.35, 2.91, 2.9, 12, 7.65, 9.23, 6.6, 
    3.8), ClinT_Stage = structure(c(4L, 2L, 5L, 5L, 2L, 5L, 2L, 
    6L, 4L, 2L, 2L, 6L, 2L, 2L, 6L, 5L, 4L, 4L, 2L, 5L), levels = c("NA", 
    "T1C", "T2", "T2A", "T2B", "T2C", "T3", "T3A", "T3B", "T3C"
    ), class = "factor"), SMS = structure(c(1L, 2L, 2L, 1L, 2L, 
    1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L
    ), levels = c("Negative", "Positive"), class = "factor"), 
    ECE = structure(c(2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 
    1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L), levels = c("Absent", 
    "Present"), class = "factor"), SVI = structure(c(2L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L), levels = c("Negative", "Positive"), class = "factor"), 
    LNI = structure(c(2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 3L, 3L, 2L, 2L, 3L, 2L, 2L), levels = c("Abnormal_N1", 
    "Normal_N0", "Not Done_NX"), class = "factor"), Num_Nodes_Removed = c(12, 
    13, 7, 13, 5, 4, 8, 7, 20, 3, 10, 10, 3, 0, 0, 2, 8, 0, 1, 
    3), Num_Nodes_Positive = c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0), PathStage = structure(c(6L, 
    1L, 4L, 4L, 3L, 3L, 3L, 4L, 3L, 7L, 3L, 3L, 4L, 3L, 4L, 3L, 
    4L, 4L, 3L, 4L), levels = c("T2A", "T2B", "T2C", "T3A", "T3B", 
    "T3C", "T4"), class = "factor"), PathGGS = structure(c(2L, 
    3L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 
    2L, 2L, 1L, 2L), levels = c("6", "7", "8", "9"), class = "factor"), 
    `Nomogram PFP_PostRP` = c("NA", "NA", "NA", "NA", "NA", "NA", 
    "NA", "NA", "NA", "NA", "NA", "NA", "NA", "98.092037304164904", 
    "95.644743112597098", "99", "99", "69.729356989036404", "NA", 
    "97.110154646526297"), BCR_FreeTime = c(18.49731, 58.02177, 
    93.14367, 152.5453, 126.0971, 160.9562, 98.59759, 149.1942, 
    64.75703, 35.05619, 82.17014, 128.4298, 10.38215, 76.45338, 
    22.70274, 74.21925, 104.3801, 18.95728, 110.3268, 56.93756
    ), BCR_Event_New = structure(c(2L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L), levels = c("0", 
    "1"), class = "factor"), SVI_Binary = structure(c(2L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L), levels = c("0", "1"), class = "factor"), LNI_Binary = structure(c(1L, 
    1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 
    1L, 3L, 1L, 1L), levels = c("0", "1", "2"), class = "factor")), row.names = c(NA, 
-20L), class = c("tbl_df", "tbl", "data.frame"))

SurvObj 看起来像这样

SurvObj <- Surv(time = Cox_PH_Data$BCR_FreeTime, event = Cox_PH_Data$BCR_Event_New)

导致错误的代码行是

adjusted_survival_curve_LNI <- survfit(Final_Cox_Model, newdata = newdata_LNI)

产生的错误是， colMeans(危险)中的错误： 'x' 必须是至少二维的数组

Answer 1

如果我正确评估问题，则您没有正确指定多状态 Cox-PH 模型，因为您的事件变量只有两个级别。从

coxph()

函数获取文档，事件必须指定如下：

The status indicator, normally 0=alive, 1=dead. Other choices are TRUE/FALSE 
(TRUE = death) or 1/2 (2=death). For interval censored data, the status indicator
 is 0=right censored, 1=event at time, 2=left censored, 3=interval censored. 
For multiple endpoint data the event variable will be a factor, whose first 
level is treated as censoring. Although unusual, the event indicator can be 
omitted, in which case all subjects are assumed to have an event.

事件的第一级指定审查。对于正确的多状态 Cox-PH 模型，您的事件中应该至少有三个级别（另请参阅 https://cran.r-project.org/web/packages/survival/vignettes/compete.pdf（结束）第 6 页）

> mgus2$etime <- with(mgus2, ifelse(pstat==0, futime, ptime))
> event <- with(mgus2, ifelse(pstat==0, 2*death, 1))
> mgus2$event <- factor(event, 0:2, labels=c("censor", "pcm", "death"))
> table(mgus2$event)
censor pcm death
409 115 860
> mfit2 <- survfit(Surv(etime, event) ~ sex, data=mgus2)
> print(mfit2, rmean=240, scale=12)

您是否有机会以其他方式指定模型，或者您确定您拥有多状态 Cox-PH？

multihaz()

包中的

survival

函数是通过一维风险调用的，因此不是一个数组（正如

colMeans()

所期望的那样）。这就是错误发生的原因。