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

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

library(survival)
library(dplyr)
library(speff2trial)
library(devtools)
library(survRM2)
library(purrr)
# library(smim)
load_all()

Data Cleanning

The dataset is avaiable in speff2trial::ACTG175.

# Load dataset
data(ACTG175)

• 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

# Analysis data: antiretroviral naive and no history of intravenous drug use
db <- ACTG175 %>% subset(drugs == 0 & strat == 1) 
db <- db %>% mutate(time = days / 30.25, 
              age = (age - mean(age)) / sd(age), 
              status = cens, 
              group = arms, 
              pattern = ifelse(status == 1, 1, ifelse(time >= 24, 2, 3))) %>% 
         subset(group %in% c(0, 1)) %>% 
         mutate(group = ifelse(group == 0, 0, 1) )

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

time <- db$time
status <- db$status
group <- db$group
pattern <- db$pattern
x <- db %>% select(age, symptom)
tau <- 2 * 12

Summary statistics

nrow(db)

table(group)

table(status)

n_mi <- 50
n_b <- 100
seed <- 1234

Delta Adjusted Imputation

RMST within group

.delta <- seq(1, 5, by = 0.5) 
res <- map(.delta, function(delta0){
  delta   <- ifelse(group == 0, c(1,1,1)[pattern], c(1,1,delta0)[pattern])  # the third number control delta adjustment for MNAR 
  tmp <- rmst_delta(time, status, x = x, group, pattern, delta, tau, n_mi = n_mi, n_b = n_b, seed = seed)
  tmp
})
res0 <- map(res, function(x) as.data.frame(x$rmst)) 
names(res0) <- .delta
res0 <- bind_rows(res0, .id = "delta")

res0 %>% mutate(
  lower = rmst - 1.96 * sd, 
  upper = rmst + 1.96 * sd,
  lower_wb = rmst - 1.96 * wb_sd,
  upper_wb = rmst + 1.96 * wb_sd
) %>% mutate_if(is.numeric, formatC, digits = 2, format = "f")

RMST between group

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

res0 %>% group_by(delta) %>% summarise(
  diff_rmst = diff(rmst), 
  diff_sd = sqrt( sum(sd^2) ), 
  diff_wb_sd = sqrt(sum(wb_sd^2)),
  lower = diff_rmst - 1.96 * diff_sd,
  upper = diff_rmst + 1.96 * diff_sd,
  p_val = 2* (1 - pnorm( abs(diff_rmst/diff_sd) )),
  lower_wb = diff_rmst - 1.96 * diff_wb_sd, 
  upper_wb = diff_rmst + 1.96 * diff_wb_sd,
  p_val_wb = 2* (1 - pnorm( abs(diff_rmst/diff_wb_sd) )),
) %>% mutate_if(is.numeric, formatC, digits = 3, format = "f") %>% data.frame()

Control Based Imputation

delta   <- c(1,1,1)[pattern]  # the third number control delta adjustment for MNAR 
tmp <- rmst_control(time, status, x = x, 
                    group = group, 
                    ref_grp = 0, pattern = pattern, delta = delta, tau = tau, 
                    n_mi = n_mi, n_b = n_b, seed = seed)
tmp$rmst

RMST between group

• first row: Rubin's rule
• second row: Wild bootstrap

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

x <- rbind( diff_rmst(tmp$rmst[,"rmst"], tmp$rmst[, "sd"]), 
       diff_rmst(tmp$rmst[,"rmst"], tmp$rmst[, "wb_sd"]))
x <- as.data.frame(x)
x$lower <- x$diff - 1.96 * x$sd
x$upper <- x$diff + 1.96 * x$sd
round(x, digits = 3)

Direct Estimator (Tian et.al 2014)

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

Model Diagnosis

Placebo group

fit0 <- coxph(Surv(time, status) ~ age + symptom, data = subset(db, group == 0))
cox.zph(fit0)

fit_cen0 <- coxph(Surv(time, 1 - status) ~ age + symptom, data = subset(db, group == 0))
cox.zph(fit_cen0)

Treatment group

fit1 <- coxph(Surv(time, status) ~ age + symptom, data = subset(db, group == 1))
cox.zph(fit1)

fit_cen1 <- coxph(Surv(time, 1 - status) ~ age + symptom, data = subset(db, group == 1))
cox.zph(fit_cen1)

elong0527/smim documentation built on May 5, 2021, 1:03 p.m.