Augmented Inverse Probability Weighted Complete Case Estimator for Survival Analysis

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

set.seed(123)
library(truncnorm)
library(survival)
library(adjustedCurves)
library(MASS)
library(knitr)
library(kableExtra)
library(dagitty)
library(parallel)
library(future.apply)
library(future)

n_cores <- detectCores() - 2

Background

We continue our previous discussion on doubly robust (DR) estimators for time-to-event outcomes using pseudo observations but focus on the situation where the loss to follow-up depends on pre-baseline and post-baseline characteristics, referred to as informative censoring. For example, patients who are sicker at baseline or who develop treatment-related adverse events may be more likely to drop out early. In such situations, the censoring mechanism must be explicitly modeled as a function of baseline and time-dependent covariates.

Augmented Inverse Probability Weighted Complete Case Estimator

We use similar notations as before and consider \(n\) independent and identically distributed individuals indexed by \(i\). Let \(D\) denote a binary treatment (0 or 1), \(T_d\) denote the potential survival time under treatment \(D=d\), \(C\) the censoring time, \(U_i=\mathrm{min}(T_i, C_i)\), and \(\boldsymbol{Z}\) denotes pre-baseline covariates.

Let \(S_d(t)\) be the survival function under treatment \(d\) at time \(t\). We define the average treatment effect (ATE) as the difference in survival probabilities between treatment groups at time \(t\), i.e., \(S_1(t)-S_0(t)\). Without loss of generality, we focus our illustration on estimating \(S_1(t)\), and \(S_0(t)\) can be estimated analogously.

Note that we only get to observe \(T_1\) when \(D=1\) and \(\Delta = 1(C \ge T_1)\); otherwise \(T_1\) is missing or coarsened. Further, this is monotone coarsening (Tsiatis 2006) and we assume coarsening at random, i.e., \((D,C) \perp T_1|Z\). Intuitively, coarsening in this case is induced by observed treatment and loss to follow-up, with treatment assignment depends only on pre-baseline characteristics and loss to follow-up relies on both pre- and post-baseline covariates, that is, informative censoring.

To handle bias induced by monotonic data coarsening, we apply the augmented inverse probability weighted Complete Case (AIPWCC) estimator described in Bai, Tsiatis, and O’Brien (2013). Specifically, the AIPWCC estimator

accounts for coarsening due to treatment assignment using the propensity score model: \(\hat{e}(\boldsymbol{Z})=P(D=1|Z)\)
handles informative censoring using the treatment-specific censoring model: \(\hat{K}_1^c(r, Z) = P(C \ge r|Z, D=1)\)
leverages the outcome regression model like in standard DR estimators: \(\hat{m}_1(\boldsymbol{Z})=P(T_1>r|Z)\)

The optimal AIPWCC estimator is \[ \hat{S}_1(u) = \frac{1}{n} \sum_{i=1}^{n} \left[ \frac{D_i\mathbb{1}(U_i\ge u)}{\hat{e}(Z_i)\hat{K}_1^c(u, Z_i)}-\left\{\frac{D_i-\hat{e}(Z_i)}{\hat{e}(Z_i)} \hat{m}_1(u, Z_i)\right\}+\int_0^u \frac{D_i}{\hat{e}(Z_i)}\frac{\mathrm{d}M_{1, i}^c(r, Z_i)}{\hat{K}_1^c(r, Z_i)}\frac{\hat{m}_1(u, Z_i)}{\hat{m}_1(r, Z_i)}\right]\] where \(\mathrm{d}M_{1, i}^c(r, Z_i) = \frac{\mathrm{d}N^c(r)-\lambda_1^c(r, Z_i)Y(r)}{\hat{K}_1^c(r, Z_i)}\) is the martingale increment for the censoring distribution with \(N^c(r)=\mathbb{1}(U\le r, \Delta=0)\), \(Y(r)=\mathbb{1}(U\ge r)\), and \(\lambda_1^c(r, Z_i)=\frac{-\mathrm{d} \mathrm{log}\hat{K}_1^c(r, Z_i)}{\mathrm{d}r}\).

The first term in the bracket in the equation above is the inverse probability weighted complete case term. The second term is the outcome-regression augmentation term in a standard DR estimator and has a mean of zero if the PS model is correct. The third term is the censoring augmentation term and also has a mean of zero if the censoring model is correct. The third term is essential to making the estimator robust to misspecification in the censoring model when the outcome model is correct.

The AIPWCC estimator is doubly robust:

If the coarsening models, including the propensity score model and censoring model, are correct
If the outcome model, that is, the conditional survival distribution is correct

The variance of the AIPWCC estimator for \(\hat{S}_1(u)\) can be estimated using Equation (5) in Bai, Tsiatis, and O’Brien (2013). If all three models (PS model, censoring model, and outcome model) are correctly specified, the estimator is locally semiparametric efficient, meaning it attains the smallest possible variance. One caveat is that if the outcome model or coarsening models are misspecified, the AIPWCC estimator for ATE remains unbiased, but the standard sandwich variance estimator in Equation (5) is conservative. In such cases, a bootstrap estimator of the asymptotic variance is recommended.

AIPWCC estimator versus pseudo-observation-based DR estimator

The pseudo-observation-based DR estimator proposed by Wang (2018) can only be applied when loss to follow-up is completely at random since the pseudo outcomes are estimated using the KM estimator where covariate adjustment is not allowed
Even though the pseudo-observation-based DR estimator can be extended to handle covariate-dependent censoring (Binder, Gerds, and Andersen 2014), it models the censoring mechanism and outcome within the same model
In contrast, the AIPWCC estimator allows to model the censoring and outcome separately
In addition, the AIPWCC estimator includes an augmentation term for censoring to make it robust to censoring model misspecification and gain efficiency when the censoring model is correct

Implementation of the AIPWCC estimator

We now demonstrate how to implement the AIPWCC estimator under the informative censoring scenario.

Data Generation

We first simulate the data according to the variable relationships depicted in the directed acyclic graph (DAG) below.

dag <- dagitty("dag {
  Z1 -> Treatment
  Z2 -> Treatment
  Z3 -> Treatment
  Z4 -> Treatment
  Z1 -> Death
  Z2 -> Death
  Z6 -> Death
  Z7 -> Death
  Treatment -> Death
}")

coordinates(dag) <- list(
  x = c(Z3 = 0, Z4 = 0, Z1 = 1, Z2 = 1, Treatment = 1.2, Death = 3.7, Z6 = 4, Z7 = 4),
  y = c(Z3 = 0.8, Z4 = 0.4, Z1 = 0.7, Z2 = 0.3, Treatment = 0.5, Death = 0.5, Z6 = 0.7, Z7 = 0.3)
)

plot(dag)

alpha_0 <- -0.5
alpha_Z1 <- 1.2
alpha_Z2 <- 1.1
alpha_Z3 <- 0.2
alpha_Z4 <- 0.2 

beta_0 <- -2.55      
beta_D <- -0.8 
beta_Z1 <- 1.25 
beta_Z2 <- 1.8
beta_Z6 <- 1.4 
beta_Z7 <- 1

gamma_0 <- -2.1
gamma_D <- 0.25
gamma_Z1 <- 0.2
gamma_Z2 <- 1.8
gamma_Z8 <- 0.7
gamma_Z9 <- 0.17

generate_data_inf <- function(n = 500) {
  Z1 <- rtruncnorm(n, mean = 0, sd = 1, a = -0.5, b = 0.5)
  Z2 <- rnorm(n, 0, 1)
  Z3 <- rnorm(n, 0, 1)
  Z4 <- rnorm(n, 0, 1)
  Z5 <- rnorm(n, 0, 1)
  Z6 <- rnorm(n, 0, 1)
  Z7 <- rnorm(n, 0, 1)
  Z8 <- rnorm(n, 0, 1)
  Z9 <- rnorm(n, 0, 1)
  
  logit_p <- alpha_0 + alpha_Z1 * Z1 + alpha_Z2 * Z2 + alpha_Z3 * Z3 + alpha_Z4 * Z4
  p <- plogis(logit_p) # exp(logit_p) / (1 + exp(logit_p))
  D <- rbinom(n, 1, p) # Treatment indicator
  
  hazard <- exp(beta_0 + beta_D * D + beta_Z1 * Z1 + beta_Z2 * Z2 + beta_Z6 * Z6 + beta_Z7 * Z7)
  T_event <- rexp(n, rate = hazard)
  
  hazard_cen <- exp(gamma_0 + gamma_D * D + gamma_Z1 * Z1 + gamma_Z2 * Z2 + gamma_Z8 * Z8 + gamma_Z9 * Z9)
  T_censor <- rexp(n, rate = hazard_cen)
  T_admin <- rep(5, n)
  
  Time <- pmin(T_event, T_censor, T_admin)
  Event <- as.numeric(T_event <= T_censor & T_event <= T_admin)
  Censored <- 1 - Event
  
  return(data.frame(Z1 = Z1, Z2 = Z2, Z3 = Z3, Z4 = Z4, Z5 = Z5, Z6 = Z6, Z7 = Z7, Z8 = Z8, Z9 = Z9, D = D, Time = Time, Event = Event, Censored = Censored))
}

data <- generate_data_inf()

attr(data$Z1, "label") <- "Baseline value of Z1"
attr(data$Z2, "label") <- "Baseline value of Z2"
attr(data$Z3, "label") <- "Baseline value of Z3"
attr(data$Z4, "label") <- "Baseline value of Z4"
attr(data$Z5, "label") <- "Baseline value of Z5"
attr(data$Z6, "label") <- "Baseline value of Z6"
attr(data$Z7, "label") <- "Baseline value of Z7"
attr(data$Z8, "label") <- "Baseline value of Z8"
attr(data$Z9, "label") <- "Baseline value of Z9"
attr(data$D, "label") <- "Treatment indicator (1 = treated, 0 = control)"
attr(data$Time, "label") <- "Time index for longitudinal records"
attr(data$Event, "label") <- "Event indicator (1 = event occurred, 0 = otherwise)"
attr(data$Censored, "label") <- "Censoring indicator (1 = censored, 0 = otherwise)"

get_label <- function(x) {
  lbl <- attr(x, "label", exact = TRUE)
  if (is.null(lbl)) "" else as.character(lbl)
}

dict <- data.frame(
  Variable = names(data),
  Meaning  = vapply(data, get_label, character(1)),
  check.names = FALSE
)

knitr::kable(dict, caption = "Data Dictionary", row.names = FALSE)

Data Dictionary
Variable	Meaning
Z1	Baseline value of Z1
Z2	Baseline value of Z2
Z3	Baseline value of Z3
Z4	Baseline value of Z4
Z5	Baseline value of Z5
Z6	Baseline value of Z6
Z7	Baseline value of Z7
Z8	Baseline value of Z8
Z9	Baseline value of Z9
D	Treatment indicator (1 = treated, 0 = control)
Time	Time index for longitudinal records
Event	Event indicator (1 = event occurred, 0 = otherwise)
Censored	Censoring indicator (1 = censored, 0 = otherwise)

summary(data$Time)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> 0.001226 0.444392 1.844927 2.504917 5.000000 5.000000

table(data$Event)
#> 
#>   0   1 
#> 364 136

table(data$D)
#> 
#>   0   1 
#> 307 193

censored_subjects <- data[data$Censored == 1, ]
(early_dropout <- sum(censored_subjects$Time < 5))
#> [1] 198
(admin_censored <- sum(censored_subjects$Time >= 5))
#> [1] 166

Among the 364 subjects with no events, 198 are early dropouts and 166 are administratively censored. Administrative censoring is random, so no additional modeling is needed. Early dropout could be covariate-dependent or occur completely at random, and we assume it is informative in this data generation.

With the simulation above, we can calculate the true survival probability.

get_true_survival <- function(data, t, D) {
  hazard <- exp(beta_0 + beta_D * D + beta_Z1 * data$Z1 + beta_Z2 * data$Z2 + beta_Z6 * data$Z6 + beta_Z7 * data$Z7)
  survival_probs <- exp(-hazard * t) # Baseline hazard is constant under an exponential survival model
  return(mean(survival_probs)) 
}

We generate multiple datasets to compute the true ATEs at prespecified time points and corresponding 95% confidence intervals (CIs) using percentiles.

n_datasets <- 100
datasets <- replicate(n_datasets, generate_data_inf(), simplify = FALSE)

times <- c(1, 2, 3)
true_ate_samples <- matrix(NA, nrow = n_datasets, ncol = length(times))

for (i in 1:n_datasets) {
  # Compute potential survival probabilities by setting treatment to 1 and 0 for all individuals
  true_surv_1 <- sapply(times, function(t) get_true_survival(datasets[[i]], t, D = 1))
  true_surv_0 <- sapply(times, function(t) get_true_survival(datasets[[i]], t, D = 0))
  
  # Calculate ATE as difference in survival probabilities
  true_ate_samples[i, ] <- true_surv_1 - true_surv_0
}

true_ate_mean <- colMeans(true_ate_samples)
true_ate_ci_lower <- apply(true_ate_samples, 2, quantile, probs = 0.025)
true_ate_ci_upper <- apply(true_ate_samples, 2, quantile, probs = 0.975)

results_true <- list(ATE = true_ate_mean, ci_lower = true_ate_ci_lower, ci_upper = true_ate_ci_upper)

Using R package

First, we implement the AIPWCC estimator using the adjustedsurv function in the adjustedCurves R package.

aipwcc_dr_ate_pkg <- function(data, times, ps_correct = TRUE, or_correct = TRUE, 
                              cen_correct = TRUE, variance = "asymptotic") {
  
  data$D <- factor(data$D, levels = c(0, 1))
  
  # Propensity score model
  if (ps_correct) {
    ps_model <- glm(D ~ Z1 + Z2 + Z3 + Z4, data = data, family = "binomial")
  } else {
    ps_model <- glm(D ~ Z5, data = data, family = "binomial")
  }

  # Censoring model 
  # Note: Event == 0 (Censored == 1) includes both loss to follow-up and administrative censoring
  if (cen_correct) {
    cen_model <- coxph(Surv(Time, Event == 0) ~ D + Z1 + Z2 + Z8 + Z9, data = data, x = TRUE)
  } else {
    cen_model <- coxph(Surv(Time, Event == 0) ~ Z6 + Z7, data = data, x = TRUE)
  }
  
  # Outcome regression (Cox proportional hazard model)
  if (or_correct) {
    or_model <- coxph(Surv(Time, Event) ~ D + Z1 + Z2 + Z6 + Z7, data = data, x = TRUE)
  } else {
    or_model <- coxph(Surv(Time, Event) ~ D + Z5, data = data, x = TRUE)
  }
  
  # AIPWCC estimator for survival probabilities 
  if (variance == "asymptotic") { # Use asymptotic variance
    adj_surv <- adjustedsurv(
      data = data,
      variable = "D",
      ev_time = "Time",
      event = "Event",
      method = "aiptw",
      treatment_model = ps_model,
      censoring_model = cen_model,
      outcome_model = or_model,
      times = times, 
      conf_int = TRUE,
      conf_level = 0.95,
      bootstrap = FALSE) 
  } else if (variance == "bootstrap"){ # Use bootstrap variance
    adj_surv <- adjustedsurv(
      data = data,
      variable = "D",
      ev_time = "Time",
      event = "Event",
      method = "aiptw",
      treatment_model = ps_model,
      censoring_model = cen_model,
      outcome_model = or_model,
      times = times,
      conf_int = TRUE,
      conf_level = 0.95,
      bootstrap = TRUE,
      n_boot = 100) 
  }
    
  # Treatment effect is defined as the difference in survival probabilities.
  ate_results <- adjusted_curve_diff(
    adj = adj_surv,
    group_1 = "1",  # Treated group
    group_2 = "0",  # Control group
    conf_int = TRUE,
    conf_level = 0.95)
  
  # Extract survival probabilities at times of interest
  ate_at_times <- ate_results[ate_results$time %in% times, ]
  
  return(list(ATE = ate_at_times$diff, ci_lower = ate_at_times$ci_lower, ci_upper = ate_at_times$ci_upper))
}

We estimate the ATEs at prespecified time points and corresponding 95% confidence intervals (CIs).

results_asymp <- aipwcc_dr_ate_pkg(data, times, variance = "asymptotic")
results_boot <- aipwcc_dr_ate_pkg(data, times, variance = "bootstrap")

results_package_tbl <- data.frame(
  Time = paste0("T = ", times),
  True = sprintf("%.3f (%.3f, %.3f)", results_true$ATE, results_true$ci_lower, results_true$ci_upper),
  Asymptotic = sprintf("%.3f (%.3f, %.3f)", results_asymp$ATE, results_asymp$ci_lower, results_asymp$ci_upper),
  Bootstrap = sprintf("%.3f (%.3f, %.3f)", results_boot$ATE, results_boot$ci_lower, results_boot$ci_upper)
)

kable(results_package_tbl,
      caption = "ATE (95% CI)",
      col.names = c("Time", "True", "Asymptotic", "Bootstrap"),
      align = c('c', 'c', 'c', 'c')) %>%
  kable_styling(full_width = FALSE, font_size = 12)

ATE (95% CI)
Time	True	Asymptotic	Bootstrap
T = 1	0.077 (0.069, 0.084)	0.103 (0.036, 0.169)	0.103 (0.036, 0.169)
T = 2	0.092 (0.085, 0.100)	0.102 (0.025, 0.178)	0.102 (0.025, 0.178)
T = 3	0.100 (0.093, 0.109)	0.086 (0.000, 0.171)	0.086 (0.000, 0.171)

Crude estimation

Fit a crude Cox model with no covariate or propensity score adjustment to show the amount of confounding bias in our data

crude_cox_estimator <- function(data, times) {
  crude_cox <- coxph(Surv(Time, Event) ~ D, data = data)
  treatment_effect <- coef(crude_cox)["D"]
  
  base_surv <- survfit(crude_cox)
  s0 <- summary(base_surv, times = times, extend = TRUE)$surv
  
  # Survival for treated group: S_baseline^exp(beta_D)
  surv_1 <- s0^exp(treatment_effect)
  
  # Survival for control group: S_baseline^exp(0) 
  surv_0 <- s0

  ATE <- surv_1 - surv_0
  return(list(ATE = ATE))
}

results_crude <- crude_cox_estimator(data, times)

Evaluate the performance of AIPWCC estimator

We use the datasets generated previously to calculate RMSE, bias, and variance.

plan(multisession, workers = n_cores)

calculate_metrics <- function(datasets, true_ate_samples, times, method = "package", variance = "asymptotic", 
                              ps_correct = TRUE, or_correct = TRUE, cen_correct = TRUE) {
  
  n_datasets <- length(datasets)
  
  ate <- future_lapply(1:n_datasets, function(i) {
    library(survival)
    library(adjustedCurves)
    
    est_results <- switch(
      method,
      "package" = aipwcc_dr_ate_pkg(datasets[[i]], times, variance = variance, 
                                    ps_correct = ps_correct, or_correct = or_correct, cen_correct = cen_correct),
      "crude" = crude_cox_estimator(datasets[[i]], times))
    
    est_results$ATE}, future.seed = TRUE)
  
  ate_all <- do.call(rbind, ate)
  
  bias <- colMeans(ate_all - true_ate_samples)
  variance <- apply(ate_all, 2, var)
  rmse <- sqrt(colMeans((ate_all - true_ate_samples)^2))
  
  return(list(bias = bias, variance = variance, rmse = rmse))
}

Model misspecification scenarios

To illustrate the robustness of the AIPWCC estimator to model misspecification, we compare the following scenarios:

Scenario 1: The outcome model, the propensity score model, and the censoring model are all correctly specified
Scenario 2: The outcome model and propensity score model are correctly specified, while the censoring model is misspecified
Scenario 3: The outcome model and censoring model are correctly specified, while the propensity score model is misspecified
Scenario 4: The outcome model is correctly specified, while the propensity score model and censoring model are both misspecified
Scenario 5: The propensity score model and censoring model are correctly specified, but the outcome model is misspecified
Scenario 6: The propensity score model is correctly specified, but the outcome model and censoring model are both misspecified
Scenario 7: The censoring model is correctly specified, but the outcome model and propensity score model are both misspecified
Scenario 8: All three models are misspecified

# Scenario 1: All models correct
metrics_1_asymp <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "asymptotic",
                                     ps_correct = TRUE, cen_correct = TRUE, or_correct = TRUE)
metrics_1_boot <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "bootstrap",
                                    ps_correct = TRUE, cen_correct = TRUE, or_correct = TRUE)

# Scenario 2: PS correct, Cen incorrect, OR correct
metrics_2_asymp <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "asymptotic",
                                     ps_correct = TRUE, cen_correct = FALSE, or_correct = TRUE)
metrics_2_boot <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "bootstrap",
                                    ps_correct = TRUE, cen_correct = FALSE, or_correct = TRUE)

# Scenario 3: PS incorrect, Cen correct, OR correct
metrics_3_asymp <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "asymptotic",
                                     ps_correct = FALSE, cen_correct = TRUE, or_correct = TRUE)
metrics_3_boot <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "bootstrap",
                                    ps_correct = FALSE, cen_correct = TRUE, or_correct = TRUE)

# Scenario 4: PS incorrect, Cen incorrect, OR correct
metrics_4_asymp <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "asymptotic",
                                     ps_correct = FALSE, cen_correct = FALSE, or_correct = TRUE)
metrics_4_boot <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "bootstrap",
                                    ps_correct = FALSE, cen_correct = FALSE, or_correct = TRUE)

# Scenario 5: PS correct, Cen correct, OR incorrect
metrics_5_asymp <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "asymptotic",
                                     ps_correct = TRUE, cen_correct = TRUE, or_correct = FALSE)
metrics_5_boot <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "bootstrap",
                                    ps_correct = TRUE, cen_correct = TRUE, or_correct = FALSE)

# Scenario 6: PS correct, Cen incorrect, OR incorrect
metrics_6_asymp <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "asymptotic",
                                     ps_correct = TRUE, cen_correct = FALSE, or_correct = FALSE)
metrics_6_boot <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "bootstrap",
                                    ps_correct = TRUE, cen_correct = FALSE, or_correct = FALSE)

# Scenario 7: PS incorrect, Cen correct, OR incorrect
metrics_7_asymp <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "asymptotic",
                                     ps_correct = FALSE, cen_correct = TRUE, or_correct = FALSE)
metrics_7_boot <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "bootstrap",
                                    ps_correct = FALSE, cen_correct = TRUE, or_correct = FALSE)

# Scenario 8: All models incorrect
metrics_8_asymp <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "asymptotic",
                                     ps_correct = FALSE, cen_correct = FALSE, or_correct = FALSE)
metrics_8_boot <- calculate_metrics(datasets, true_ate_samples, times, method = "package", variance = "bootstrap",
                                    ps_correct = FALSE, cen_correct = FALSE, or_correct = FALSE)

# Scenario 0: Crude Cox
metrics_crude <- calculate_metrics(datasets, true_ate_samples, times, method = "crude")

scenarios <- c("1: All Correct",
  "2: PS Correct, Cen Incorrect, OR Correct",
  "3: PS Incorrect, Cen Correct, OR Correct",
  "4: PS Incorrect, Cen Incorrect, OR Correct",
  "5: PS Correct, Cen Correct, OR Incorrect",
  "6: PS Correct, Cen Incorrect, OR Incorrect",
  "7: PS Incorrect, Cen Correct, OR Incorrect",
  "8: All Incorrect")

# Combine values across all time points into a single string
format_values <- function(vec) {
  paste(sprintf("%.4f", vec), collapse = ", ")
}

# RMSE table
rmse_asymp <- c(format_values(metrics_1_asymp$rmse), format_values(metrics_2_asymp$rmse), 
                format_values(metrics_3_asymp$rmse), format_values(metrics_4_asymp$rmse),
                format_values(metrics_5_asymp$rmse), format_values(metrics_6_asymp$rmse), 
                format_values(metrics_7_asymp$rmse), format_values(metrics_8_asymp$rmse))

rmse_boot <- c(format_values(metrics_1_boot$rmse), format_values(metrics_2_boot$rmse),
               format_values(metrics_3_boot$rmse), format_values(metrics_4_boot$rmse),
               format_values(metrics_5_boot$rmse), format_values(metrics_6_boot$rmse),
               format_values(metrics_7_boot$rmse), format_values(metrics_8_boot$rmse))

rmse_table <- data.frame(Scenario = scenarios, Asymptotic = rmse_asymp, Bootstrap = rmse_boot)

kable(rmse_table,
      col.names = c("Scenario", "T = 1, 2, 3", "T = 1, 2, 3"),
      caption = "Performance of the estimator (RMSE)",
      align = c('c', 'c', 'c')) %>%
  kable_styling(full_width = FALSE, font_size = 12) %>%
  add_header_above(c(" " = 1, "Asymptotic" = 1, "Bootstrap" = 1)) %>%
  footnote(general = paste0("Crude Cox: ", format_values(metrics_crude$rmse)), general_title = "")

Performance of the estimator (RMSE)
	Asymptotic	Bootstrap
Scenario	T = 1, 2, 3	T = 1, 2, 3
1: All Correct	0.0292, 0.0330, 0.0345	0.0292, 0.0330, 0.0345
2: PS Correct, Cen Incorrect, OR Correct	0.0285, 0.0320, 0.0344	0.0285, 0.0320, 0.0344
3: PS Incorrect, Cen Correct, OR Correct	0.0297, 0.0316, 0.0326	0.0297, 0.0316, 0.0326
4: PS Incorrect, Cen Incorrect, OR Correct	0.0282, 0.0299, 0.0313	0.0282, 0.0299, 0.0313
5: PS Correct, Cen Correct, OR Incorrect	0.0391, 0.0414, 0.0419	0.0391, 0.0414, 0.0419
6: PS Correct, Cen Incorrect, OR Incorrect	0.0367, 0.0384, 0.0399	0.0367, 0.0384, 0.0399
7: PS Incorrect, Cen Correct, OR Incorrect	0.1409, 0.1614, 0.1738	0.1409, 0.1614, 0.1738
8: All Incorrect	0.1347, 0.1560, 0.1696	0.1347, 0.1560, 0.1696
Crude Cox: 0.1232, 0.1524, 0.1688

# Bias table
bias_asymp <- c(format_values(metrics_1_asymp$bias), format_values(metrics_2_asymp$bias), 
                format_values(metrics_3_asymp$bias), format_values(metrics_4_asymp$bias),
                format_values(metrics_5_asymp$bias), format_values(metrics_6_asymp$bias), 
                format_values(metrics_7_asymp$bias), format_values(metrics_8_asymp$bias))

bias_boot <- c(format_values(metrics_1_boot$bias), format_values(metrics_2_boot$bias),
               format_values(metrics_3_boot$bias), format_values(metrics_4_boot$bias),
               format_values(metrics_5_boot$bias), format_values(metrics_6_boot$bias),
               format_values(metrics_7_boot$bias), format_values(metrics_8_boot$bias))

bias_table <- data.frame(Scenario = scenarios, Asymptotic = bias_asymp, Bootstrap = bias_boot)

kable(bias_table,
      col.names = c("Scenario", "T = 1, 2, 3", "T = 1, 2, 3"),
      caption = "Performance of the estimator (Bias)",
      align = c('c', 'c', 'c')) %>%
  kable_styling(full_width = FALSE, font_size = 12) %>%
  add_header_above(c(" " = 1, "Asymptotic" = 1, "Bootstrap" = 1)) %>%
  footnote(general = paste0("Crude Cox: ", format_values(metrics_crude$bias)), general_title = "")

Performance of the estimator (Bias)
	Asymptotic	Bootstrap
Scenario	T = 1, 2, 3	T = 1, 2, 3
1: All Correct	-0.0004, 0.0008, -0.0021	-0.0004, 0.0008, -0.0021
2: PS Correct, Cen Incorrect, OR Correct	-0.0002, 0.0010, -0.0020	-0.0002, 0.0010, -0.0020
3: PS Incorrect, Cen Correct, OR Correct	0.0002, 0.0006, -0.0019	0.0002, 0.0006, -0.0019
4: PS Incorrect, Cen Incorrect, OR Correct	0.0000, -0.0004, -0.0028	0.0000, -0.0004, -0.0028
5: PS Correct, Cen Correct, OR Incorrect	-0.0006, 0.0036, 0.0021	-0.0006, 0.0036, 0.0021
6: PS Correct, Cen Incorrect, OR Incorrect	-0.0010, 0.0038, 0.0035	-0.0010, 0.0038, 0.0035
7: PS Incorrect, Cen Correct, OR Incorrect	-0.1352, -0.1550, -0.1668	-0.1352, -0.1550, -0.1668
8: All Incorrect	-0.1298, -0.1509, -0.1641	-0.1298, -0.1509, -0.1641
Crude Cox: -0.1194, -0.1475, -0.1632

# Variance table
var_asymp <- c(format_values(metrics_1_asymp$variance), format_values(metrics_2_asymp$variance), 
               format_values(metrics_3_asymp$variance), format_values(metrics_4_asymp$variance),
               format_values(metrics_5_asymp$variance), format_values(metrics_6_asymp$variance), 
               format_values(metrics_7_asymp$variance), format_values(metrics_8_asymp$variance))

var_boot <- c(format_values(metrics_1_boot$variance), format_values(metrics_2_boot$variance),
              format_values(metrics_3_boot$variance), format_values(metrics_4_boot$variance),
              format_values(metrics_5_boot$variance), format_values(metrics_6_boot$variance),
              format_values(metrics_7_boot$variance), format_values(metrics_8_boot$variance))

var_table <- data.frame(Scenario = scenarios, Asymptotic = var_asymp, Bootstrap = var_boot)

kable(var_table,
      col.names = c("Scenario", "T = 1, 2, 3", "T = 1, 2, 3"),
      caption = "Performance of the estimator (Variance)",
      align = c('c', 'c', 'c')) %>%
  kable_styling(full_width = FALSE, font_size = 12) %>%
  add_header_above(c(" " = 1, "Asymptotic" = 1, "Bootstrap" = 1)) %>%
  footnote(general = paste0("Crude Cox: ", format_values(metrics_crude$variance)), general_title = "")

Performance of the estimator (Variance)
	Asymptotic	Bootstrap
Scenario	T = 1, 2, 3	T = 1, 2, 3
1: All Correct	0.0009, 0.0011, 0.0012	0.0009, 0.0011, 0.0012
2: PS Correct, Cen Incorrect, OR Correct	0.0008, 0.0010, 0.0012	0.0008, 0.0010, 0.0012
3: PS Incorrect, Cen Correct, OR Correct	0.0009, 0.0010, 0.0011	0.0009, 0.0010, 0.0011
4: PS Incorrect, Cen Incorrect, OR Correct	0.0008, 0.0009, 0.0010	0.0008, 0.0009, 0.0010
5: PS Correct, Cen Correct, OR Incorrect	0.0016, 0.0017, 0.0018	0.0016, 0.0017, 0.0018
6: PS Correct, Cen Incorrect, OR Incorrect	0.0014, 0.0015, 0.0016	0.0014, 0.0015, 0.0016
7: PS Incorrect, Cen Correct, OR Incorrect	0.0015, 0.0020, 0.0024	0.0015, 0.0020, 0.0024
8: All Incorrect	0.0013, 0.0016, 0.0019	0.0013, 0.0016, 0.0019
Crude Cox: 0.0009, 0.0015, 0.0019

Funding

This work was supported by Utah Clinical & Translational Science Institute (CTSI) Translational Innovation Pilot (TIP) Program Award (NCATS UM1TR004409).

References

Bai, Xiaofei, Anastasios A. Tsiatis, and Sean M. O’Brien. 2013. “Doubly-Robust Estimators of Treatment-Specific Survival Distributions in Observational Studies with Stratified Sampling.” Biometrics 69 (December): 830–39. https://doi.org/10.1111/biom.12076.

Binder, Nadine, Thomas A. Gerds, and Per Kragh Andersen. 2014. “Pseudo-Observations for Competing Risks with Covariate Dependent Censoring.” Lifetime Data Analysis 20 (April): 303–15. https://doi.org/10.1007/s10985-013-9247-7.

Tsiatis, Anastasios A. 2006. Semiparametric Theory and Missing Data. Springer Series in Statistics. New York, NY: Springer. https://doi.org/10.1007/0-387-37345-4.

Wang, Jixian. 2018. “A Simple, Doubly Robust, Efficient Estimator for Survival Functions Using Pseudo Observations.” Pharmaceutical Statistics 17 (1): 38–48.

Corresponding author: yizhe.xu@hsc.utah.edu ↩︎

Augmented Inverse Probability Weighted Complete Case Estimator for Survival Analysis

Yizhe Xu¹, Peirong Hao, Ravinder Singh, Kevin Ying, Adam Bress, Tom Greene

2026-01-29

Background

Augmented Inverse Probability Weighted Complete Case Estimator

AIPWCC estimator versus pseudo-observation-based DR estimator

Implementation of the AIPWCC estimator

Data Generation

Using R package

Crude estimation

Evaluate the performance of AIPWCC estimator

Model misspecification scenarios

Funding

References

Augmented Inverse Probability Weighted Complete Case Estimator for Survival Analysis

Yizhe Xu1, Peirong Hao, Ravinder Singh, Kevin Ying, Adam Bress, Tom Greene

2026-01-29

Background

Augmented Inverse Probability Weighted Complete Case Estimator

AIPWCC estimator versus pseudo-observation-based DR estimator

Implementation of the AIPWCC estimator

Data Generation

Using R package

Crude estimation

Evaluate the performance of AIPWCC estimator

Model misspecification scenarios

Funding

References

Yizhe Xu¹, Peirong Hao, Ravinder Singh, Kevin Ying, Adam Bress, Tom Greene