In elong0527/smim: Survival Sensitivity Analysis using Multiple Imputation and Martingale

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(survival)
library(dplyr)
library(ISwR)
library(devtools)
library(survRM2)
# library(smim)
load_all()

Data Cleanning

Pattern
1: observed event time
2: MAR imputation
3: MNAR imputation
Delta
1: no adjustment
larger than 1: hazard is larger than MAR
less than 1: hazard is smaller than MAR

data("melanom")

db <- melanom
db <- db %>% mutate(
  time = days/30.25,   # unit in Month
  status = ifelse(status %in% c(1,3), 0, 1 ),
  group  = ifelse(thick > median(thick), 1, 0),
  pattern = ifelse(status == 1, 1, ifelse(time < 8, 2, 3)),
  log_thick = log(thick)
) %>% arrange(group, status, pattern, time) %>% 
  mutate(id = 1:n())

time <- db$time
status <- db$status
group <- db$group
pattern <- db$pattern
x <- db %>% select(sex, ulc)
tau <- 10*12
# db0 <- db %>% select(time, status, sex, log_thick)
# fit <- coxph(Surv(time, status) ~ rep(1, nrow(db0)), data = db0, x = TRUE, y = TRUE)

MAR analysis

rmst2(time, status, group, tau= tau)

fit <- survfit(Surv(time, status) ~ group, data = db)
# Visualize with survminer
plot(fit, col = c(1,2))

Delta Adjusted Imputation (MAR)

RMST within group

delta   <- c(1,1,1)[pattern]  # the third number control delta adjustment for MNAR 
tmp <- rmst_delta(time, status, x = rep(1, length(time)), group, pattern, delta, tau, n_mi = 5, n_b = 100, seed = 1234)
tmp$rmst

RMST between group

diff_rmst <- function(rmst, sd){
  diff <- diff(rmst)
  diff_sd <- sqrt( sum(sd^2) )
  p_val <- 2* (1 - pnorm( abs(diff/diff_sd) ))
  c(diff = diff, sd = diff_sd, p = p_val)
}

rbind( diff_rmst(tmp$rmst[,"rmst"], tmp$rmst[, "sd"]), 
       diff_rmst(tmp$rmst[,"rmst"], tmp$rmst[, "wb_sd"]))

Survival Curve After Imputation

u_time <- sort(unique(time))
plot(tmp$surv[[1]][, "time"], tmp$surv[[1]][, "survival"], type = "l")
lines(tmp$surv[[2]][, "time"], tmp$surv[[2]][, "survival"], col = 2)

Delta Adjusted Imputation (MNAR)

RMST within group

delta   <- c(1,1,1.25)[pattern]  # the third number control delta adjustment for MNAR 
tmp <- rmst_delta(time, status, x = x, group, pattern, delta, tau, n_mi = 5, n_b = 100, seed = 1234)
tmp$rmst

RMST between group

diff_rmst <- function(rmst, sd){
  diff <- diff(rmst)
  diff_sd <- sqrt( sum(sd^2) )
  p_val <- 2* (1 - pnorm( abs(diff/diff_sd) ))
  c(diff = diff, sd = diff_sd, p = p_val)
}

rbind( diff_rmst(tmp$rmst[,"rmst"], tmp$rmst[, "sd"]), 
       diff_rmst(tmp$rmst[,"rmst"], tmp$rmst[, "wb_sd"]))

Survival Curve After Imputation

u_time <- sort(unique(time))
plot(tmp$surv[[1]][, "time"], tmp$surv[[1]][, "survival"], type = "l")
lines(tmp$surv[[2]][, "time"], tmp$surv[[2]][, "survival"], col = 2)