Marginal effect for cox model

marginal_effect.R

library(survival)
library(tidyverse)
library(broom)
library(marginaleffects)


# Make a toy model
gr <- sample(letters[1:2], size = 10000, replace=T)
g <- as.numeric(gr=='a')
etime <- rexp(length(gr), rate = (5+1*g)/10)
d <- tibble(gr, etime, event=1)

# Fit the cox model
fit <- coxph(Surv(etime, event) ~ gr, data=d)

# Estimate surival times myself
b<-map_dfr(c('a','b'), ~{
  survfit(fit, 
          newdata = tibble(gr=.x)) %>% 
    tidy %>% 
    mutate(gr=.x)
})


# Difference in survival times by hand
bb <-b %>% 
  select(time, gr, estimate) %>% 
  pivot_wider(id_cols=time, names_from = gr, values_from = estimate) %>% 
  mutate(delta = b-a)
  
bb %>% 
  ggplot(aes(time, delta)) + 
  geom_line()

# Looks pretty close!
plot_comparisons(fit, 
                 variables = 'gr',  by='etime', type='survival')