Simulation Workflow for Primary Analysis

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Lei Shi, Matthew Secrest

2026-08-20

Primary analysis

This vignette demonstrates Monte Carlo simulation for evaluating the EC-IPW and EC-AIPW weighting estimators proposed in Zhou et al. (2024) for the primary (placebo-controlled) phase.

1 Simulate a list of datasets for primary analysis

1.1 Create basic setup

# Initialize an empty data list
set.seed(2023)

data_matrix_list_null <- list()
data_matrix_list_alt <- list()
# a small ntrial keeps this vignette fast to build; use hundreds to
# thousands of trials for stable operating characteristic estimates
ntrial <- 20

# Specify the significance level alpha
alpha <- 0.05

# Specify the true effect size at the end of the study
true_effect <- 0
alt_effect <- 2.0 # tune this

# Specify column names
covariates_col_name <- c("x1", "x2", "x3", "x4", "x5")
outcome_col_name <- c("y1", "y2")
treatment_col_name <- "A"
trial_status_col_name <- "S"

1.2 Generate simulation data list for null hypothesis

# Sequentially adding in datasets
for (trial_iter in 1:ntrial) {
  # simulate 300 sample
  normal <- copula::normalCopula(param = c(0.8), dim = 4, dispstr = "ar1")

  # ========== generate internal covariates =============
  X_int <- simulate_X_copula(
    n = 200,
    p = 4,
    cp = normal, # copula
    margins = c("binom", "binom", "binom", "exp"), # specify marginal distributions
    paramMargins = list(
      list(size = 1, prob = 0.7), # specify parameters for marginals
      list(size = 1, prob = 0.9),
      list(size = 1, prob = 0.3),
      list(rate = 1 / 10)
    )
  )

  X_int$x4 <- round(X_int$x4) + 1
  X_int$x5 <- 30 + 10 * X_int$x1 + (7) * X_int$x2 + (-6) * X_int$x3 + (-0.5) * X_int$x4 + rnorm(200, mean = 0, sd = 10)

  varnames <- c("1", paste0("x", 1:5))

  # ============ generate external covariates ==============
  X_ext <- simulate_X_copula(
    n = 100,
    p = 4,
    cp = normal, # copula
    margins = c("binom", "binom", "binom", "exp"), # specify marginal distributions
    paramMargins = list(
      list(size = 1, prob = 0.7), # specify parameters for marginals
      list(size = 1, prob = 0.9),
      list(size = 1, prob = 0.3),
      list(rate = 1 / 10)
    )
  )

  X_ext$x4 <- round(X_ext$x4) + 1
  X_ext$x5 <- 50 + 10 * X_ext$x1 + (2) * X_ext$x2 + (-1) * X_ext$x3 + (-0.3) * X_ext$x4 + rnorm(100, mean = 0, sd = 10)

  varnames <- c("1", paste0("x", 1:5))

  # ============ Specify outcome models ==============
  model_form_x_t1 <- setNames(c(10.0, 0.05, -1.5, -1.0, -0.2, -0.1), varnames) # 1.5*A, sigma = 4.0
  model_form_x_t2 <- setNames(c(6.0, 0.5, -0.5, -1.0, -0.3, -0.06), varnames) # 1.8*A, sigma = 4.0
  # model_form_x_t3 = setNames(c(5.0, 1.9, 1.4, -1.3, -0.4, -0.15), varnames) # 1.6*A, sigma = 4.0
  # model_form_x_t4 = setNames(c(1.2, 1.0, 2.0, -0.5, -0.4, -0.10), varnames) # 2.5*A, sigma = 5.0

  outcome_model_specs <- list(
    list(
      effect = 0, model_form_x = model_form_x_t1, # from data: true_effect = 1.5
      noise_mean = 0, noise_sd = 4
    ), # model form for the first time point, given by model_form1
    list(
      effect = true_effect, model_form_x = model_form_x_t2, # from data: true_effect = 1.8
      noise_mean = 0, noise_sd = 4
    ) # model form for the second time point, given by model_form2
  )

  # =========== generate trial data ============
  Data <- simulate_trial(X_int,
    X_ext,
    num_treated = 100,
    OLE_flag = FALSE,
    T_cross = 2,
    outcome_model_specs
  )

  data_matrix_list_null[[trial_iter]] <- Data
}


head(data_matrix_list_null[[1]])
#>   x1 x2 x3 x4       x5 A S         y1         y2
#> 1  0  0  0  2 37.56031 1 1  3.0752036  7.4562185
#> 2  1  1  0 15 32.11467 0 1  0.3814777  1.4237345
#> 3  1  1  0  5 29.65879 1 1 -2.1309058  4.8196265
#> 4  1  1  0 14 48.56377 1 1 -0.7990459 -3.6642685
#> 5  1  1  1 21 15.37401 1 1  1.7543228 -0.5165716
#> 6  1  1  1 13 23.47198 1 1 -2.1335935 -4.0533044

1.3 Generate simulation data list for alternative hypothesis

# Sequentially adding in datasets
for (trial_iter in 1:ntrial) {
  # simulate 300 sample
  normal <- copula::normalCopula(param = c(0.8), dim = 4, dispstr = "ar1")

  # ========== generate internal covariates =============
  X_int <- simulate_X_copula(
    n = 200,
    p = 4,
    cp = normal, # copula
    margins = c("binom", "binom", "binom", "exp"), # specify marginal distributions
    paramMargins = list(
      list(size = 1, prob = 0.7), # specify parameters for marginals
      list(size = 1, prob = 0.9),
      list(size = 1, prob = 0.3),
      list(rate = 1 / 10)
    )
  )

  X_int$x4 <- round(X_int$x4) + 1
  X_int$x5 <- 30 + 10 * X_int$x1 + (7) * X_int$x2 + (-6) * X_int$x3 + (-0.5) * X_int$x4 + rnorm(200, mean = 0, sd = 10)

  varnames <- c("1", paste0("x", 1:5))

  # ============ generate external covariates ==============
  X_ext <- simulate_X_copula(
    n = 100,
    p = 4,
    cp = normal, # copula
    margins = c("binom", "binom", "binom", "exp"), # specify marginal distributions
    paramMargins = list(
      list(size = 1, prob = 0.7), # specify parameters for marginals
      list(size = 1, prob = 0.9),
      list(size = 1, prob = 0.3),
      list(rate = 1 / 10)
    )
  )

  X_ext$x4 <- round(X_ext$x4) + 1
  X_ext$x5 <- 50 + 10 * X_ext$x1 + (2) * X_ext$x2 + (-1) * X_ext$x3 + (-0.3) * X_ext$x4 + rnorm(100, mean = 0, sd = 10)

  varnames <- c("1", paste0("x", 1:5))

  # ============ Specify outcome models ==============
  model_form_x_t1 <- setNames(c(10.0, 0.05, -1.5, -1.0, -0.2, -0.1), varnames) # 1.5*A, sigma = 4.0
  model_form_x_t2 <- setNames(c(6.0, 0.5, -0.5, -1.0, -0.3, -0.06), varnames) # 1.8*A, sigma = 4.0
  # model_form_x_t3 = setNames(c(5.0, 1.9, 1.4, -1.3, -0.4, -0.15), varnames) # 1.6*A, sigma = 4.0
  # model_form_x_t4 = setNames(c(1.2, 1.0, 2.0, -0.5, -0.4, -0.10), varnames) # 2.5*A, sigma = 5.0

  outcome_model_specs <- list(
    list(
      effect = 0, model_form_x = model_form_x_t1, # from data: true_effect = 1.5
      noise_mean = 0, noise_sd = 4
    ), # model form for the first time point, given by model_form1
    list(
      effect = alt_effect, model_form_x = model_form_x_t2, # from data: true_effect = 1.8
      noise_mean = 0, noise_sd = 4
    ) # model form for the second time point, given by model_form2
  )

  # =========== generate trial data ============
  Data <- simulate_trial(X_int,
    X_ext,
    num_treated = 100,
    OLE_flag = FALSE,
    T_cross = 2,
    outcome_model_specs
  )

  data_matrix_list_alt[[trial_iter]] <- Data
}

head(data_matrix_list_alt[[1]])
#>   x1 x2 x3 x4       x5 A S         y1        y2
#> 1  0  1  0  2 36.90508 1 1  3.6260347  6.265349
#> 2  1  1  0 14 25.44926 1 1  5.6809759  2.566003
#> 3  1  1  1 28 27.19396 1 1 -0.5518921 -2.074243
#> 4  0  1  0 10 32.65669 0 1 -0.1818891  5.195533
#> 5  0  1  0  6 31.22264 1 1  4.8591983  3.665490
#> 6  1  1  1 35 22.77644 0 1  0.4364034 -9.368366

2 Parametric inference

method_IPW_optimal_weight <- ec_ipw(
  ps_formula = "S ~ x1 + x2 + x3 + x4 + x5"
)

method_AIPW_optimal_weight <- ec_aipw(
  ps_formula = "S ~ x1 + x2 + x3 + x4 + x5",
  outcome_formula = c(
    "y1 ~ x1 + x2 + x3 + x4 + x5",
    "y2 ~ x1 + x2 + x3 + x4 + x5"
  )
)

method_IPW_zero_weight <- ec_ipw(
  ps_formula = "S ~ x1 + x2 + x3 + x4 + x5",
  weight = 0
)

method_AIPW_zero_weight <- ec_aipw(
  ps_formula = "S ~ x1 + x2 + x3 + x4 + x5",
  outcome_formula = c(
    "y1 ~ x1 + x2 + x3 + x4 + x5",
    "y2 ~ x1 + x2 + x3 + x4 + x5"
  ),
  weight = 0
)

method_obj_list <- list(
  method_IPW_optimal_weight,
  method_AIPW_optimal_weight,
  method_IPW_zero_weight,
  method_AIPW_zero_weight
)

# create a simulation object for primary analysis
simulation_primary_obj <- setup_simulation_primary(
  data_matrix_list_null = data_matrix_list_null, # two scenarios
  data_matrix_list_alt = data_matrix_list_alt,
  trial_status_col_name = trial_status_col_name,
  treatment_col_name = treatment_col_name,
  outcome_col_name = outcome_col_name,
  covariates_col_name = covariates_col_name,
  method_obj_list = method_obj_list,
  true_effect = true_effect,
  alt_effect = alt_effect,
  alpha = alpha,
  method_description = c(
    "IPW, optimal weight",
    "AIPW, optimal weight",
    "IPW, zero weight",
    "AIPW, zero weight"
  )
)

simulation_report <- run_simulation(simulation_primary_obj, quiet = TRUE)

simulation_report # Type I error and Power
#>     method_description         bias  variance       mse coverage type_I_error
#> 1  IPW, optimal weight -0.003333374 0.5948377 0.5948488     0.90         0.10
#> 2 AIPW, optimal weight -0.105637948 0.2263140 0.2374734     0.95         0.05
#> 3     IPW, zero weight  0.009894312 0.6239639 0.6240618     0.95         0.05
#> 4    AIPW, zero weight -0.147679223 0.2401897 0.2619988     0.95         0.05
#>   power
#> 1  0.75
#> 2  1.00
#> 3  0.75
#> 4  1.00

3 Bootstrap inference

method_IPW_optimal_weight <- ec_ipw(
  ps_formula = "S ~ x1 + x2 + x3 + x4 + x5",
  bootstrap = 50,
  bootstrap_ci_type = "perc"
)

method_AIPW_optimal_weight <- ec_aipw(
  ps_formula = "S ~ x1 + x2 + x3 + x4 + x5",
  outcome_formula = c(
    "y1 ~ x1 + x2 + x3 + x4 + x5",
    "y2 ~ x1 + x2 + x3 + x4 + x5 + y1"
  ),
  bootstrap = 50,
  bootstrap_ci_type = "perc"
)

method_IPW_zero_weight <- ec_ipw(
  ps_formula = "S ~ x1 + x2 + x3 + x4 + x5",
  weight = 0,
  bootstrap = 50,
  bootstrap_ci_type = "perc"
)

method_AIPW_zero_weight <- ec_aipw(
  ps_formula = "S ~ x1 + x2 + x3 + x4 + x5",
  outcome_formula = c(
    "y1 ~ x1 + x2 + x3 + x4 + x5",
    "y2 ~ x1 + x2 + x3 + x4 + x5 + y1"
  ),
  weight = 0,
  bootstrap = 50,
  bootstrap_ci_type = "perc"
)

method_obj_list <- list(
  method_IPW_optimal_weight,
  method_AIPW_optimal_weight,
  method_IPW_zero_weight,
  method_AIPW_zero_weight
)

# create a simulation object for primary analysis
simulation_primary_obj <- setup_simulation_primary(
  data_matrix_list_null = data_matrix_list_null,
  data_matrix_list_alt = data_matrix_list_alt,
  trial_status_col_name = trial_status_col_name,
  treatment_col_name = treatment_col_name,
  outcome_col_name = outcome_col_name,
  covariates_col_name = covariates_col_name,
  method_obj_list = method_obj_list,
  true_effect = true_effect,
  alt_effect = alt_effect,
  alpha = alpha,
  method_description = c(
    "IPW, optimal weight, bootstrap",
    "AIPW, optimal weight, bootstrap",
    "IPW, zero weight, bootstrap",
    "AIPW, zero weight, bootstrap"
  )
)

simulation_report_bs <- run_simulation(simulation_primary_obj, quiet = TRUE)

simulation_report_bs
#>                method_description         bias  variance       mse coverage
#> 1  IPW, optimal weight, bootstrap -0.003333374 0.5948377 0.5948488     0.95
#> 2 AIPW, optimal weight, bootstrap -0.108635705 0.2250218 0.2368236     0.95
#> 3     IPW, zero weight, bootstrap  0.009894312 0.6239639 0.6240618     0.90
#> 4    AIPW, zero weight, bootstrap -0.150771625 0.2366983 0.2594304     0.95
#>   type_I_error power
#> 1         0.05   0.7
#> 2         0.05   0.9
#> 3         0.10   0.7
#> 4         0.05   1.0

References

Zhou X, Zhu J, Drake C, Pang H (2024). “Causal estimators for incorporating external controls in randomized trials with longitudinal outcomes.” Journal of the Royal Statistical Society Series A: Statistics in Society. doi: 10.1093/jrsssa/qnae075.
Shi L, Pang H, Chen C, Zhu J (2025). “rdborrow: an R package for causal inference incorporating external controls in randomized controlled trials with longitudinal outcomes.” Journal of Biopharmaceutical Statistics, 35(6), 1043-1066. doi: 10.1080/10543406.2025.2489283.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.