User-defined target population

Introduction

See other vignettes for details.

library(boot)      # non-parametric bootstrap in MAIC and ML G-computation
library(copula)    # simulating BC covariates from Gaussian copula
library(rstanarm)  # fit outcome regression, draw outcomes in Bayesian G-computation
library(tidyr)
library(dplyr)
library(MASS)
library(outstandR)
library(simcovariates)

Data

N <- 200
allocation <- 2/3      # active treatment vs. placebo allocation ratio (2:1)
b_trt <- log(0.17)     # conditional effect of active treatment vs. common comparator
b_X <- -log(0.5)       # conditional effect of each prognostic variable
b_EM <- -log(0.67)     # conditional interaction effect of each effect modifier
meanX_AC <- c(0.45, 0.45)       # mean of normally-distributed covariate in AC trial
meanX_BC <- c(0.6, 0.6)         # mean of each normally-distributed covariate in BC
meanX_EM_AC <- c(0.45, 0.45)    # mean of normally-distributed EM covariate in AC trial
meanX_EM_BC <- c(0.6, 0.6)      # mean of each normally-distributed EM covariate in BC
sdX <- c(0.4, 0.4)     # standard deviation of each covariate (same for AC and BC)
sdX_EM <- c(0.4, 0.4)  # standard deviation of each EM covariate
corX <- 0.2            # covariate correlation coefficient  
b_0 <- -0.6            # baseline intercept coefficient  ##TODO: fixed value

ipd_trial <- gen_data(N, b_trt, b_X, b_EM, b_0,
                      meanX_AC, sdX, 
                      meanX_EM_AC, sdX_EM, 
                      corX, allocation,
                      family = binomial("logit"))

ipd_trial$trt <- factor(ipd_trial$trt, labels = c("C", "A"))

BC.IPD <- gen_data(N, b_trt, b_X, b_EM, b_0,
                   meanX_BC, sdX, 
                   meanX_EM_BC, sdX_EM, 
                   corX, allocation,
                   family = binomial("logit"))

BC.IPD$trt <- factor(BC.IPD$trt, labels = c("C", "B"))

# covariate summary statistics
# assume same between treatments
cov.X <- 
  BC.IPD %>%
  as.data.frame() |> 
  dplyr::select(X1, X2, X3, X4, trt) %>%
  pivot_longer(cols = starts_with("X"), names_to = "variable", values_to = "value") %>%
  group_by(variable) %>%
  summarise(
    mean = mean(value),
    sd = sd(value)
  ) %>%
  pivot_longer(cols = c("mean", "sd"), names_to = "statistic", values_to = "value") %>%
  ungroup() |> 
  mutate(trt = NA)

# outcome
summary.y <- 
  BC.IPD |> 
  as.data.frame() |> 
  dplyr::select(y, trt) %>%
  pivot_longer(cols = "y", names_to = "variable", values_to = "value") %>%
  group_by(variable, trt) %>%
  summarise(
    mean = mean(value),
    sd = sd(value),
    sum = sum(value)
  ) %>%
  pivot_longer(cols = c("mean", "sd", "sum"),
               names_to = "statistic", values_to = "value") %>%
  ungroup()

# sample sizes
summary.N <- 
  BC.IPD |> 
  group_by(trt) |> 
  count(name = "N") |> 
  pivot_longer(cols = "N", names_to = "statistic", values_to = "value") |> 
  mutate(variable = NA_character_) |> 
  dplyr::select(variable, statistic, value, trt)
  
ald_trial <- rbind.data.frame(cov.X, summary.y, summary.N)

lin_form <- as.formula("y ~ X3 + X4 + trt + trt:X1 + trt:X2")

Parametric G-computation with maximum-likelihood estimation

outstandR_gcomp_ml <-
  outstandR(ipd_trial, ald_trial,
            strategy = strategy_gcomp_ml(
              formula = lin_form,
              family = binomial(link = "logit")))

outstandR_gcomp_ml

Once again, let us change the outcome scale to LRR

outstandR_gcomp_ml_lrr <-
  outstandR(ipd_trial, ald_trial,
            strategy = strategy_gcomp_ml(
              formula = lin_form,
              family = binomial(link = "logit")),
            scale = "log_relative_risk")

outstandR_gcomp_ml_lrr

Bayesian G-computation with MCMC

outstandR_gcomp_stan <-
  outstandR(ipd_trial, ald_trial,
            strategy = strategy_gcomp_stan(
              formula = lin_form,
              family = binomial(link = "logit")))

outstandR_gcomp_stan

As before, we can change the outcome scale to LRR.

outstandR_gcomp_stan_lrr <-
  outstandR(ipd_trial, ald_trial,
            strategy = strategy_gcomp_stan(
              formula = lin_form,
              family = binomial(link = "logit")),
            scale = "log_relative_risk")

outstandR_gcomp_stan_lrr

Multiple imputation marginalisation

outstandR_mim <-
  outstandR(ipd_trial, ald_trial,
            strategy = strategy_mim(
              formula = lin_form,
              family = binomial(link = "logit")))

outstandR_mim

Change the outcome scale again for LRR,

outstandR_mim_lrr <-
  outstandR(ipd_trial, ald_trial,
            strategy = strategy_mim(
              formula = lin_form,
              family = binomial(link = "logit")),
            scale = "log_relative_risk")

outstandR_mim_lrr