Exposure and outcome

Let’s assume that we are interested in studying the comparative effectiveness of two antineoplastic systemic drug treatment regimens (exposure (0/1)) on the time to death due to any reason as the outcome (overall survival). The anticipated time of follow-up is truncated to 5 years.

The desired strength of effectiveness of the exposure of interest is defined with a hazard ratio (HR) of 1.0, i.e. there is no difference in overall survival among patients who are treated with the exposure of interest as compared to the comparator regimen. The proportional hazards assumption is fulfilled for this dataset.

Confounders

We further assume that there are a some confounders which we need to specify to estimate our outcome model. Most of the covariates are associated with both the probability of treatment initiation and the outcome but there are also some that are not predictive of the exposure and just the outcome or not associated with any of the exposure or outcome, whatsoever.

Despite the low dimensionality, the dataset is simulated as realistically as possible with varying strengths of associations between covariates and treatment initiation and the outcome.

Missingness

Most importantly for testing the functionality of this package, some of the above mentioned confounders are just partially observed according to the missingness mechanisms and proportions specified below in the table below.

Overview covariates/confounder structure

To get an overview of the dataset, this table provides a summary of the different covariate-exposure-outcome-missingness correlations.

Simulate time-to-event

Next, we simulate a time-to-event outcome for overall survival. For this, the simsurv package is used with the following assumptions:

• Parametric event times, following an exponential distribution
• Max event times: 5 years of follow-up
• Event times depend on some (time-fixed/baseline) covariate effects as listed in the table above
• beta coefficients for the outcome model as defined in the following

betas_os <- c(
  exposure = log(1),
  age_num = log(1.05),
  female_cat = log(0.94),
  ecog_cat = log(1.25),
  smoking_cat = log(1.3),
  physical_cat = log(0.79),
  egfr_cat = log(0.5),
  alk_cat = log(0.91),
  pdl1_num = log(0.98),
  histology_cat = log(1.15)
  )

betas_os %>% 
  as.data.frame() %>% 
  transmute(logHR = round(`.`, 2)) %>% 
  rownames_to_column(var = "Covariate") %>% 
  mutate(HR = round(exp(logHR), 2)) %>% 
  gt()

set.seed(seed_value)

sim_df <- sim_covar %>% bind_cols(
  simsurv(
    dist = "exponential",
    lambdas = 0.05,
    betas = betas_os,
    x = sim_covar,
    maxt = 5 
    )
  ) %>% 
  select(-id)

km_overall <- survfit(Surv(eventtime, status) ~ 1, data = sim_df)
km_overall
#> Call: survfit(formula = Surv(eventtime, status) ~ 1, data = sim_df)
#> 
#>         n events median 0.95LCL 0.95UCL
#> [1,] 2500   2017   1.57    1.46    1.68

km_exposure <- survfit(Surv(eventtime, status) ~ exposure, data = sim_df)
km_exposure
#> Call: survfit(formula = Surv(eventtime, status) ~ exposure, data = sim_df)
#> 
#>               n events median 0.95LCL 0.95UCL
#> exposure=0 1502   1277   1.28    1.17    1.41
#> exposure=1  998    740   2.18    1.91    2.43

km_exposure <- survfit(Surv(eventtime, status) ~ exposure, data = sim_df)

ggsurvplot(
  km_exposure, 
  data = sim_df,
  conf.int = TRUE,
  surv.median.line = "hv",
  palette = "jco",
  xlab = "Time [Years]",
  legend.labs = c("Comparator", "Exposure of interest")
  )
#> Warning in geom_segment(aes(x = 0, y = max(y2), xend = max(x1), yend = max(y2)), : All aesthetics have length 1, but the data has 2 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.
#> All aesthetics have length 1, but the data has 2 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.

cox_lhs <- "survival::Surv(eventtime, status)"
cox_rhs <- paste(colnames(sim_covar), collapse = " + ")
cox_form = as.formula(paste(cox_lhs, "~ exposure +", cox_rhs))
  
cox_fit <- coxph(cox_form, data = sim_df)

cox_fit %>% 
  tbl_regression(exponentiate = T)

Export smdi_data_complete

To provide flexibility to play with the complete data before introducing missingness, we offer two datasets:

• smdi_data_complete: complete data
• smdi_data: data with certain covariates missing (see next chapter)

smdi_data_complete <- sim_df
use_data(smdi_data_complete, overwrite = TRUE)

Missing complete at random

ecog_cat¹ will be set to missing according to the following specification:

• Missingness pattern: Only ecog_cat will be set to missing
• Weight vector: There are no weights, since the missingness is not dependent on any other observed or unobserved covariates
• Type: Not applicable

# specify missingness pattern
# (0 = set to missing, 1 = remains complete)
mcar_col <- which(colnames(tmp)=="ecog_cat")
miss_pattern_mcar <- replace(miss_pattern, mcar_col, 0)

miss_prop_mcar <- .35

set.seed(42)
smdi_data_mcar <- ampute(
  data = tmp,
  prop = miss_prop_mcar,
  mech = "MCAR",
  patterns = miss_pattern_mcar,
  bycases = TRUE
  )$amp %>% 
  select(ecog_cat)

Missing at random

egfr_cat will be set to missing according to the following specification:

• Missingness pattern: Only egfr_cat will be set to missing
• Weight vector: all covariates (except egfr_cat itself and ses_cat) contribute with equal weights as linear predictors for the probability of observations becoming missing
• Type: higher weighted sum scores (i.e. higher odds of having an egfr_cat mutation) will have a larger probability of becoming incomplete
• Define an overall missingness proportion of ~40%

# specify missingness pattern
# (0 = set to missing, 1 = remains complete)
mar_col <- which(colnames(tmp)=="egfr_cat")
miss_pattern_mar <- replace(miss_pattern, mar_col, 0)

# weights to compute missingness probability 
# by assigning a non-zero value
miss_weights_mar <- rep(1, ncol(tmp))
miss_weights_mar <- replace(miss_weights_mar, mar_col, 0)

miss_prop_mar <- .4

set.seed(42)
smdi_data_mar <- ampute(
  data = tmp,
  prop = miss_prop_mar,
  mech = "MAR",
  patterns = miss_pattern_mar,
  weights = miss_weights_mar,
  bycases = TRUE,
  type = "RIGHT"
  )$amp

Missing not at random - value

• Missingness pattern: Only pdl1_num will be set to missing
• Weight vector: Only pdl1_num itself (by a non-zero value) is the linear predictor for the probability of observations becoming missing
• Type: Specify a type of logistic probability distribution for the missingness weights, so that cases with low weighted sum scores (i.e. lower pdl1_num expression) will have a larger probability of becoming incomplete
• Define an overall missingness proportion of ~20%

# determine missingness pattern
mnar_v_col <- which(colnames(tmp)=="pdl1_num")
miss_pattern <- rep(1, ncol(tmp))
miss_pattern_mnar_v <- replace(miss_pattern, mnar_v_col, 0)

# weights to compute missingness probability 
# by assigning a non-zero value
# MNAR_v: covariate itself is only predictor
miss_weights_mnar_v <- rep(0, ncol(tmp))
miss_weights_mnar_v <- replace(miss_weights_mnar_v, mnar_v_col, 1)

miss_prop_mnar_v <- .2

set.seed(42)
smdi_data_mnar_v <- ampute(
  data = tmp,
  prop = miss_prop_mnar_v,
  mech = "MNAR",
  patterns = miss_pattern_mnar_v,
  weights = miss_weights_mnar_v,
  bycases = TRUE,
  type = "LEFT"
  )$amp

Assemble final dataset

smdi_data <- smdi_data_complete %>% 
  select(-c(ecog_cat, egfr_cat, pdl1_num)) %>% 
  bind_cols(ecog_cat = smdi_data_mcar$ecog_cat, egfr_cat = smdi_data_mar$egfr_cat, pdl1_num = smdi_data_mnar_v$pdl1_num) %>% 
  mutate(across(ends_with("cat"), as.factor))