Support for 'response' for Cox survival models ?

[katossky]

Would it make sense to add support for type "response" for Cox models, that would predict the average survival time ? I've seen at least one person doing something similar for discrete models here. In theory, the expected survival time is computable from the hazard by integration, which should be numerically tractable, shouldn't it ? The benefit would be to be able to compare Cox models with accelerated time models, which at the moment one can not. (Disclaimer : I am by no way an expert in survival models, just a generalist statistician.)

If it were to make sense, what would it take to implement it ? From the code it looks like type is passed down until it reaches stat::predict so this is not a marginaleffects issue ?

[vincentarelbundock]

This should be very easy for users to do. See here:

https://marginaleffects.com/bonus/extensions.html#modify-or-extend-supported-models

That said, adding new types is out of scope for marginaleffects, so this would be a user-side project.

[katossky]

I tried the suggested approach. I got part of the package to work, but I also run into some trouble. I take any feedbacks on the culprits for the errors.

library(survival)
library(lubridate) 
library(dplyr)
library(marginaleffects)

# simulated toy data ---------------------------------------------------------------------

set.seed(123)

n  <- 500
id <- 1:n

birth_year <- sample(1940:1970, n, replace = TRUE)
gender     <- factor(sample(c("Male", "Female"), n, replace = TRUE))

# True retirement age follows a normal distribution,
# with women retiring ~2 years earlier on average
true_retire_age <- rnorm(n,
                         mean = ifelse(gender == "Female", 63, 65),
                         sd   = 3)

# Current date (analysis time) and current age
today      <- ymd("2025-06-17")        # keep script reproducible
age_today  <- year(today) - birth_year

# Event indicator: retired (1) or right-censored (0)
retired <- true_retire_age <= age_today

# Realised follow-up time (in years) -------------------
# – if retired: age at retirement
# – if censored: age at time-of-analysis
time <- ifelse(retired, true_retire_age, age_today)

# Optional: a retirement-date column (NA for censored)
date_retirement <- rep(NA_Date_, n)
date_retirement[retired] <-
  ymd(paste0(birth_year[retired] + round(true_retire_age[retired]), "-07-01"))

# Assemble dataset
df <- tibble(
  id,
  birth_year,
  gender,
  date_retirement,
  time,          # time-to-event or censoring
  retired        # event indicator
)


# Cox proportional hazards ---------------------------------------------------------------

options("marginaleffects_model_classes" = "mycoxph")

m_cox1  <- coxph(Surv(time, retired) ~ gender + birth_year, data = df, ties = "efron")
class(m_cox1) # "coxph"
class(m_cox1) <- c("mycoxph", class(m_cox1))

# prediction functions

get_predict.mycoxph <- function(model, newdata, type="response", ...) {
  if(type != "response") get_predict.coxph(model, newdata, type, ...) else {

    sf <- survfit(model, newdata = newdata)

    fs <-  rbind(1, sf$surv) # distribution de la durée de vie
    t  <- c(0, sf$time) # temps où la durée de vie est évaluée
    dt <- diff(t) # intervalles

    return(data.frame(
      rowid = 1:ncol(fs), # --- edited used to be id ---
      estimate = colSums(((fs[-nrow(fs), ] + fs[-1, ]) / 2) * dt)
    ))
  }
}

# tests that suceeded --------------------------------------------------------------------

datagrid(
  model = m_cox1,
  gender = unique,
  birth_year = unique
)

get_predict(
  model = m_cox1,
  newdata = datagrid(
    model = m_cox1,
    gender = unique,
    birth_year = unique
  )
)

predictions(
  m_cox1,
  newdata = datagrid(gender = unique, birth_year = unique)
)

avg_predictions(m_cox1)

comparisons(
  m_cox1, variables = "gender", type = "response",
  newdata =  datagrid(birth_year = unique)
) # OK but does not show birth_year in the output contrary to in doc

plot_predictions(m_cox1, condition = "birth_year")
plot_predictions(m_cox1, condition = "gender")

slopes(m_cox1, variables = "birth_year", newdata =  datagrid(
  gender = unique,
  birth_year = unique
))

avg_slopes(m_cox1, variables = "birth_year")  # --- edited now works ---
plot_slopes(m_cox1, variables = "birth_year", condition = "birth_year") # --- edited now works ---

# tests that failed ----------------------------------------------------------------------

predictions(m_cox1, newdata = "mean")

# Unable to compute predicted values with this model. You can try to
#   supply a different dataset to the `newdata` argument.

predictions(m_cox1, vcov = "boot") # does not work but help says it should
plot_predictions(m_cox1, condition = "birth_year", vcov = "boot")
# Must be element of set
#  * {'stata','robust','HC','HC0','HC1','HC2','HC3','HC4','HC4m','HC5',
#    'HAC','NeweyWest','kernHAC','OPG','satterthwaite','kenward-roger'},
#  * but is 'boot'

# --- editted did not use with id instead of rowid ----
# predictions(m_cox1) |> inferences(method = "boot", R = 1000)
# Assertion on 'type' failed: Must be element of set {'survival','expected','lp','risk'},
# but is 'response'.
 # does not work
# Error in `ggplot2::geom_ribbon()`:
# ! Problem while computing aesthetics.
# ℹ Error occurred in the 1st layer.
# Caused by error in `.data[["birth_year"]]`:
# ! Column `birth_year` not found in `.data`.

[katossky]

I was using version 0.25.1 with R 4.3.2 or R 4.4.2

[vincentarelbundock]

the first thing i notice is that the column should be called rowid and not id

[katossky]

Thank you for pointing that out. This fixes warning / error for avg_slopes and plot_slopes. Errors with vcov="boot", newdata = 'mean' and type = 'response' remain. Should I open issues ?

[vincentarelbundock]

plot_*() does not support vcov="boot". I investigated this a while back and it is very difficult. I'd be happy to merge a PR if someone submits one, but am unlikely to work on it myself. #1460

[vincentarelbundock]

I'm not sure about the other issues. Yes, It would be worth opening issues for each problem separately. Thanks!

[vincentarelbundock]

predictions(m_cox1, newdata = "mean") fails because your custom function fails when supplied a single row data frame.

[vincentarelbundock]

Unfortunately, we cannot support custom classes in the inferences() function, because there is no guarantee that these come with an appropriate update() function. Enabling them would be dangerous.

This is now documented more clearly (in the dev version) and we raise an error immediately when called.

[vincentarelbundock]

Of course, it's quite easy to bootstrap it "manually" with rsample or boot. That's what I would recommend.