Error in colMeans(hazard) : 'x' must be an array of at least two dimensions

Question

I'm trying to develop estimated survival curves from my Cox PH Model in R. Since it is a multi-state model I have added a new data argument where I've kept covariates constant except for the 1 for which I am making the survival curve for.

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)))

My cox model is coded as follows,

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

My "SurvObj" is formatted correctly with time and event, and my variables are all the same type and level.

I've tried everything there is to try and still this error persists. I was able to run it last week after updating my Rstudio but now this error keeps occuring. Here is the output from

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"))

And the SurvObj looks like this

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

The line of code that results in the error is

adjusted_survival_curve_LNI <- survfit(Final_Cox_Model, newdata = newdata_LNI)

The resulting error is, Error in colMeans(hazard) : 'x' must be an array of at least two dimensions

Answer 1

If I am assessing the problem correctly, you do not specify the multi-state Cox-PH model correctly, as your event variable has only two levels. Taking the documentation from the coxph()-function, the event has to be specified as follows:

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.

The first level in the event specifies the censoring. For a proper multi-state Cox-PH model you should have at least three levels in your event (see also https://cran.r-project.org/web/packages/survival/vignettes/compete.pdf (end of page 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)

Do you have any chance of specifying the model in another way or are you sure you have a multi-state Cox-PH?

The multihaz() function inside the survival-package is called with a one-dimensional hazard and, therefore, not an array (as it is expected by colMeans()). This is why the error occurs.