---
title: "Data generation"
subtitle: "A simulated dataset for package functionality demonstration"
author: "Janick Weberpals"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Data generation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dpi = 150,
  fig.width = 6,
  fig.height = 4.5
  )
```

```{r setup, message=FALSE}
library(smdi)
library(dplyr)
library(tibble)
library(gtsummary)
library(gt)
library(survival)
library(simsurv)
library(survminer)
library(usethis)
library(mice)
library(cardx)

# some global simulation parameters
seed_value <- 42
n <- 2500
```

# `smdi` dataset background

To get acquainted with the functionality and usage of the `smdi` package, the package includes a simulated example dataset. The dataset is an exemplary low-dimensional electronic health records (EHR) dataset depicting a cohort of `r formatC(n, format = "fg", big.mark = ",")` lung cancer patients. The dataset follows the general *one-row-per-patient* structure, in which one row stands for an individual patient and the columns represent the variables.

## 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 {#overview}

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

| Variable      | Description                                      | Associated with exposure/outcome | Missingness [%]      |
|------------------|-------------------|------------------|------------------|
| age_num       | Age at baseline (continuous)                     | Yes/Yes                          |                      |
| female_cat    | Female gender (binary)                           | Yes/Yes                          |                      |
| ecog_cat      | ECOG performance score 0/1 or \>1 (binary)       | Yes/Yes                          | MCAR [35%]           |
| smoking_cat   | Smoker vs. non-smoker at baseline (binary)       | Yes/Yes                          |                      |
| physical_cat  | Physically active vs not active (binary)         | Yes/Yes                          |                      |
| egfr_cat      | EGFR alteration (binary)                         | Yes/Yes                          | MAR [40%]            |
| alk_cat       | ALK translocation (binary)                       | No/Yes                           |                      |
| pdl1_num      | PD-L1 expression in % (continuous)               | Yes/Yes                          | MNAR(value) [20%]    |
| histology_cat | Tumor histology squamous vs nonsquamous (binary) | No/Yes                           |                      |
| ses_cat       | Socio-economic status (multi-categorical)        | No/No                            |                      |
| copd_cat      | History of COPD (binary)                         | No/No                            | Auxiliary to smoking |

# Simulation of covariates and exposure

```{r, covar_generation}
set.seed(seed_value)

# start with basic dataframe, covariates and their association with exposure
sim_covar <- tibble(
  exposure = rbinom(n = n, size = 1, prob = 0.4),
  age_num = rnorm(n, mean = 64 - 7.5*exposure, sd = 13.7),
  female_cat = rbinom(n, size = 1, prob = 0.39 - 0.05*exposure),
  ecog_cat = rbinom(n, size = 1, prob = 0.63 - 0.04*exposure),
  smoking_cat = rbinom(n, size = 1, prob = 0.45 + 0.1*exposure), 
  physical_cat = rbinom(n, size = 1, prob = 0.35 + 0.02*exposure),
  egfr_cat = rbinom(n, size = 1, prob = 0.20 + 0.07*exposure), 
  alk_cat = rbinom(n, size = 1, prob = 0.03),
  pdl1_num = rnorm(n, mean = 40 + 10*exposure, sd = 10.5),
  histology_cat = rbinom(n, size = 1, prob = 0.2),
  ses_cat = sample(x = c("1_low", "2_middle", "3_high"), size = n, replace = TRUE, prob = c(0.2 , 0.4, 0.4)),
  copd_cat = rbinom(n, size = 1, prob = 0.3 + 0.5*smoking_cat)
  ) %>%  
  # bring data in right format
  mutate(across(ends_with("num"), as.numeric)) %>% 
  mutate(across(ends_with("num"), function(x) round(x, digits = 2)))
```

In the first step, we create a dataset with `r formatC(nrow(sim_covar), format = "fg", big.mark = ",")` patients and `r ncol(sim_covar)` variables.

```{r, distributions_covars, include = FALSE, eval = FALSE}
sim_covar %>% 
  tbl_summary(by = "exposure") %>% 
  add_difference()
```

The following table illustrates the odds of exposure assignment.

```{r, pr(exposure)_tbl}
exposure_form <- as.formula(paste("exposure ~ ", paste(colnames(sim_covar %>% select(-exposure)), collapse = " + ")))

exposure_fit <- glm(
  exposure_form,
  data = sim_covar,
  family = "binomial"
  )

exposure_fit %>% 
  tbl_regression(exponentiate = T)
```

Fitting a generalized linear model and assessing the probability of treatment assignment, the above constellation of odds results in the following simulated distributions depicting propensities of treatment initiation (aka propensity scores).

```{r, pr_treatment_assignment, fig.cap="Treatment assignment probabilities."}
# compute propensity score
exposure_plot <- sim_covar %>% 
  mutate(ps = fitted(exposure_fit))

# plot density
exposure_plot %>% 
  ggplot(aes(x = ps, fill = factor(exposure))) +
  geom_density(alpha = .5) +
  theme_bw() +
  labs(
    x = "Pr(exposure)",
    y = "Density",
    fill = "Exposed"
  )
```

## 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

```{r, betas_outcome_generation}
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()
```

```{r, outcome_generation, message=FALSE}
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)
```

### Kaplan-Meier estimates

The simulation resulted in the following crude Kaplan-Meier estimates for overall survival.

```{r, km_estimates, include = FALSE, eval = FALSE}
km_overall <- survfit(Surv(eventtime, status) ~ 1, data = sim_df)
km_exposure <- survfit(Surv(eventtime, status) ~ exposure, data = sim_df)

tbl_survfit(
  list(km_overall, km_exposure),
  times = c(1, 5),
  label_header = "**{time} Years**"
  )
```

```{r, km_estimates2}
km_overall <- survfit(Surv(eventtime, status) ~ 1, data = sim_df)
km_overall

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

Given, that the true exposure effect is null, the crude model is severely biased as we can see even more clearly in the crude Kaplan-Meier curve.

```{r}
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")
  )
```

### Cox proportional hazards

After adjusting, the simulated data results in the following hazard ratio (HR) estimates.

```{r, Cox_estimates}
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)

```{r, export_complete_data, message=FALSE}
smdi_data_complete <- sim_df
use_data(smdi_data_complete, overwrite = TRUE)
```

# Introduce missingness

In many different quantitative disciplines from classic epidemiology to machine and deep learning there is an increasing interest in utilizing electronic health records (EHR) to develop prognostic/predictive models or study the comparative effectiveness and safety of medical interventions. Especially information on variables which are not readily available in other datasets (e.g. administrative claims) are of high interest, including vital signs, biomarkers and lab data. However, these covariates are often just partially observed for various reasons:

-   Physician did simply not perform/order a certain test
-   Certain measurements are just collected for particularly sick patients
-   Information is 'hiding' in unstructured records, e.g. clinical notes

To illustrate `smdi's` main functions using this dataset, we introduce some missingness to relevant covariates, which are critical confounders of the causal exposure-outcome relationship in the `smdi_data` dataset.

In order to introduce missingness following different missingness mechanisms, we use the `ampute` function of the [mice](https://amices.org/mice/index.html) package. An excellent tutorial on this very flexible and elegant function can be found [here](https://rianneschouten.github.io/mice_ampute/vignette/ampute.html).

In brief, we introduce missingness by...

-   Determining a [missingess pattern](https://rianneschouten.github.io/mice_ampute/vignette/ampute.html#Patterns)
-   Creating a [weight vector](https://rianneschouten.github.io/mice_ampute/vignette/ampute.html#Weights), using fitted linear models for the probability of observations becoming missing
-   Specifying a type of logistic probability distribution for the missingness weights
-   Defining an overall missingness proportion

```{r}
# prepare a placeholder df for missing simulation
# we do not consider ses_cat
tmp <- smdi_data_complete %>% 
  select(-c(ses_cat))

# determine missingness pattern template
miss_pattern <- rep(1, ncol(tmp))
```

## Missing complete at random

`ecog_cat`[^1] will be set to missing according to the following specification:

[^1]: The Eastern Cooperative Oncology Group (ECOG) performance score is a clinical measure for how the cancer disease affects the daily living abilities of the patient and is often used as a patient inclusion criterion for clinical trials. The scale ranges from 0 to 5 and typically only patients with 0 and 1 are eligible. Let's assume for our example, we already have a subset of such a pre-selected clinical trial-like cohort, but we still want to adjust for a baseline ECOG of 0 and 1 (definitions are taken from <https://ecog-acrin.org/resources/ecog-performance-status/>

-   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

```{r, mcar}
# 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%

```{r, mar}
# 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%

```{r, create_mnar_v}
# 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

```{r}
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))
```

# Export `smdi_data`

Exporting the data to `data/`, so it can be used for educative purposes and for playing around with the package.

```{r, export_missing_data, message=FALSE}
use_data(smdi_data, overwrite = TRUE)
```