Correlation Matrix

dummy_rwe_training <- observational_rwe_training %>%
  dplyr::mutate(
    race = (race == "White"),
    app = (app == "used app"),
    infection = (infection == "1. infected")
  )

knitr::kable(round(cor(dummy_rwe_training), 4))

corrplot::corrplot.mixed(cor(dummy_rwe_training))

Even in this simple correlation matrix, there are already signs of the challenges ahead. Can you spot them?

Infections by App Usage (Marginal Model)

Let’s examine our main relationship of interest.

observational_rwe_training %>%
  ggplot2::ggplot(ggplot2::aes(x = app, fill = infection)) +
  ggplot2::theme_classic() +
  ggplot2::geom_bar(position = "dodge") +
  ggplot2::ggtitle("Infections by App Usage")

df_rwe_training <- with(
  observational_rwe_training,
  cbind(
    table(app, infection),
    prop.table(table(app, infection), margin = 1) # row proportions
  )
)

knitr::kable(df_rwe_training) # row proportions

There was a lower infection rate (i.e., empirical or observed infection risk) among app users: Only 9.1% of users had infections, compared to 27.7% of non-users. The empirical RD is 0.091 \(-\) 0.277 \(=\) -0.186—a marginal quantity in statistical terms because it does not account for (i.e., “is marginalized over”) any other variables. A model that does account for other variables (here, race) in addition to the would-be intervention (here, app usage) is called a conditional model.

out_fisher_rwe_training <- with(
  observational_rwe_training,
  fisher.test(app, infection)
)

out_fisher_rwe_training

## 
##  Fisher's Exact Test for Count Data
## 
## data:  app and infection
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.2493633 0.2749802
## sample estimates:
## odds ratio 
##  0.2618449

In addition, there is strong statistical evidence (i.e., statistical significance) that infections varied by app usage (\(p_{\text{EDA}} <<\) 0.001). Here, the estimated odds of infection for app users were only 0.262 (i.e., roughly a quarter) that of non-users, with an EDA 95% confidence interval (CI) of (0.249, 0.275).

• Statistical significance does not imply importance: This strong statistical evidence is not surprising. Bigger samples typical of RWE data provide better estimates of progressively smaller differences in outcome means (i.e., for any given set of outcome distributions with finite means and variances). This says nothing about the domain (e.g., scientific, clinical, research, business) usefulness of those differences.

• Correlation does not imply causation: In no way does this finding imply that using the app reduces infection risk. Statistical findings like this cannot imply causation. Statistical findings are implications of causal relationships, not vice versa. But it should prompt us to investigate further, and should be assessed alongside subsequent findings and hypotheses to help explain how this association arose.

Infections by Race

Race seems to have been associated with infection.

observational_rwe_training %>%
  ggplot2::ggplot(ggplot2::aes(x = race, fill = infection)) +
  ggplot2::theme_classic() +
  ggplot2::geom_bar(position = "dodge") +
  ggplot2::ggtitle("Infections by Race")

df_rwe_training_race_infection <- with(
  observational_rwe_training,
  cbind(
    table(race, infection),
    prop.table(table(race, infection), margin = 1) # row proportions
  )
)

knitr::kable(df_rwe_training_race_infection) # row proportions

out_fisher_rwe_training_race_infection <- with(
  observational_rwe_training,
  fisher.test(race, infection)
)

out_fisher_rwe_training_race_infection

## 
##  Fisher's Exact Test for Count Data
## 
## data:  race and infection
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.1147125 0.1272759
## sample estimates:
## odds ratio 
##  0.1208232

Blacks were more likely to have been infected than Whites (\(p_{\text{EDA}} <<\) 0.001): 42.1% of Blacks were infected, compared to 8.1% of Whites. Put differently, the estimated odds of infection for Whites were 0.121 (95% \(\text{CI}_{\text{EDA}}\): 0.115, 0.127) times that of Blacks. This finding is plausible due to structural inequality: African Americans tend to have worse access to care, lower socioeconomic status (SES), and more underlying health conditions.

App Usage by Race

Importantly, it looks like race was strongly associated with app usage.

observational_rwe_training %>%
  ggplot2::ggplot(ggplot2::aes(x = race, fill = app)) +
  ggplot2::theme_classic() +
  ggplot2::geom_bar(position = "dodge") +
  ggplot2::ggtitle("App Usage by Race")

df_rwe_training_race_app <- with(
  observational_rwe_training,
  cbind(
    table(race, app),
    prop.table(table(race, app), margin = 1) # row proportions
  )
)

knitr::kable(df_rwe_training_race_app) # row proportions

out_fisher_rwe_training_race_app <- with(
  observational_rwe_training,
  fisher.test(race, app)
)

out_fisher_rwe_training_race_app

## 
##  Fisher's Exact Test for Count Data
## 
## data:  race and app
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  18.95585 21.10802
## sample estimates:
## odds ratio 
##   20.00128

Blacks were less likely to have used the app than Whites (\(p_{\text{EDA}} <<\) 0.001): Only 27.2% of Blacks used the app, compared to 88.2% of Whites. Put differently, the estimated odds of using the app for Whites were 20.001 (95% \(\text{CI}_{\text{EDA}}\): 18.956, 21.108) times that of Blacks.

This finding is plausible also due to structural inequality. For example, suppose our app was only available on iPhone. Blacks are less likely than Whites to own iPhones (e.g., page 7 of Smith, 2013), which would have limited their access to the app. African Americans may also have worse access to telemedicine options in general.

Implications for Explanatory Modeling

What do these findings imply about how we should approach estimating the ATE? Could race affect both app usage and infection risk? If so, race would be a confounder—meaning our marginal RD estimate earlier will not suffice as an estimate of the ATE. (We’ll see how later on, and why this is so in Part 2.)

Prediction Modeling

Suppose we decide to fit a prediction model. We think this may help us decide if we should include race in our final explanatory model (in addition to app usage). We’ll fit a 10-fold cross-validated ridge-regularized generalized linear model (GLM) with logit link (i.e., corresponding to a logistic regression).

set.seed(2004110254)
glmnet_rwe_training <- glmnet::cv.glmnet(
  x = model.matrix(
    object = infection ~ app + race,
    data = observational_rwe_training
  ),
  y = observational_rwe_training$infection,
  family = "binomial",
  type.measure = "auc",
  alpha = 0 # ridge: https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html
)
glmnet_rwe_training$glmnet_preds <- as.numeric(
  predict(
    object = glmnet_rwe_training,
    newx = model.matrix(
      object = infection ~ app + race,
      data = observational_rwe_training
    ),
    s = glmnet_rwe_training$lambda.min,
    type = "response"
  )
)
glmnet_ROC_rwe_training <- pROC::roc(
  response = observational_rwe_training$infection,
  predictor = glmnet_rwe_training$glmnet_preds
)
coef_glmnet_rwe_training <- coef(glmnet_rwe_training)

plot(glmnet_ROC_rwe_training)

glmnet_ROC_rwe_training$auc # training set cross-validated AUC

## Area under the curve: 0.6954

coef_glmnet_rwe_training

## 4 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept) -1.900132655
## (Intercept)  .          
## appused app -0.001874018
## raceWhite   -0.003428336

We’ll select app regardless of its coefficient size because it is our potenial intervention. Still, we note that an automated search for important predictors ranked by coefficient size would have ranked race higher than app. We wonder how such a procedure could potentially have excluded true causes with smaller coefficients from other analyses we’ve done—causes that were important for recommending courses of action.

Importantly, note that the area under the receiver operating characteristic (ROC) curve (AUC) is not particularly large even though there exists a true ATE of app usage on infection risk by design. We use the simulation parameters and models (see Appendix) to calculate this ATE as -0.031. (We’ll derive this calculation in Part 2.) This ATE is only modestly sized; hence the modest AUC.

“To explain or to predict?” These findings underscore the important points in (Shmueli, 2010) about the different objectives, procedures, and properties of prediction and explanation.

Implications for Explanatory Modeling

We decide to keep race in the final explanatory model due to its relatively large prediction coefficient, even after controlling for app.

Causal Model

From our findings, we propose the following causal model, formalized as a directed acyclic graph (DAG) (Pearl, 2009):

1. App Usage \(\rightarrow\) Infection
2. Race \(\rightarrow\) Infection
3. Race \(\rightarrow\) App Usage

DiagrammeR::grViz("
digraph causal {

  # Nodes
  node [shape = plaintext]
  Z [label = 'Race']
  X [label = 'App \n Usage']
  Y [label = 'Infection']

  # Edges
  edge [color = black,
        arrowhead = vee]
  rankdir = LR
  X -> Y
  Z -> X
  Z -> Y

  # Graph
  graph [overlap = true]
}")

Here, race confounds the effect of app usage on infection.

Statistical Model and Sample Size

Based on our DAG, it makes sense to model the infection risk as a function of race and app usage. We propose the following logistic model.

glm_rwe_training <- glm(
  data = observational_rwe_training,
  formula = as.factor(infection) ~ app + race,
  family = "binomial"
)

summary(glm_rwe_training)$coefficients

##               Estimate Std. Error   z value     Pr(>|z|)
## (Intercept) -0.2280908 0.02278461 -10.01074 1.367266e-23
## appused app -0.3432987 0.03249127 -10.56588 4.289439e-26
## raceWhite   -1.9073638 0.03260492 -58.49926 0.000000e+00

The estimated odds of infection for app users were \(\exp(\)-0.343\() =\) 0.709 (95% \(\text{CI}_{\text{EDA}}\): 0.666, 0.756) times that of non-users (regardless of race). The corresponding estimated infection risks for app usage by race are:

risk_didnt_use_app_black_training <- plogis(coef(glm_rwe_training) %*% c(1, 0, 0))
risk_used_app_black_training <- plogis(coef(glm_rwe_training) %*% c(1, 1, 0))
risk_didnt_use_app_white_training <- plogis(coef(glm_rwe_training) %*% c(1, 0, 1))
risk_used_app_white_training <- plogis(coef(glm_rwe_training) %*% c(1, 1, 1))

• 0.443 for Blacks not using the app
• 0.361 for Blacks using the app
• 0.106 for Whites not using the app
• 0.077 for Whites using the app

In order to estimate the RD for Blacks with 80% statistical power at an \(\alpha = 0.05\) overall level of statistical evidence, the required sample sizes are:

# https://cran.r-project.org/web/packages/pwr/vignettes/pwr-vignette.html
sample_size_black <- ceiling(
  pwr::pwr.2p.test(
    h = pwr::ES.h(p1 = risk_used_app_black_training, p2 = risk_didnt_use_app_black_training),
    sig.level = 0.05,
    power = .80,
    alternative = "less"
  )$n
)
sample_size_white <- ceiling(
  pwr::pwr.2p.test(
    h = pwr::ES.h(p1 = risk_used_app_white_training, p2 = risk_didnt_use_app_white_training),
    sig.level = 0.05,
    power = .80,
    alternative = "less"
  )$n
)

• 556 Black individuals in each app usage group (i.e., users or non-users)
• 1617 White individuals in each app usage group

Note that these figures are for one-sided hypothesis tests that the infection risk is lower for app users versus non-users within each racial group. They adjust for multiple testing (i.e., multiple comparisons, multiplicity) by requiring each of the two individual statistical hypothesis tests to meet an \(\alpha = 0.025\) level of statistical evidence.

Check Sample Size

Do we have a large enough sample for each app usage and racial group combination to meet our statistical power and evidence requirements?

observational_rwe_holdout %>%
  ggplot2::ggplot(ggplot2::aes(x = race, fill = app)) +
  ggplot2::theme_classic() +
  ggplot2::geom_bar(position = "dodge") +
  ggplot2::ggtitle("App Usage by Race")

df_rwe_holdout_race_app <- with(
  observational_rwe_holdout,
  cbind(
    table(race, app),
    prop.table(table(race, app), margin = 1) # row proportions
  )
)

knitr::kable(df_rwe_holdout_race_app) # row proportions

Yes: We have at least 556 Black individuals in each app usage group, and at least 1617 White individuals likewise.

Fit Statistical Model

Based on our proposed causal and statistical models, we fit the following logistic model.

glm_rwe_holdout <- glm(
  data = observational_rwe_holdout,
  formula = as.factor(infection) ~ app + race,
  family = "binomial"
)

summary(glm_rwe_holdout)$coefficients

##               Estimate Std. Error    z value      Pr(>|z|)
## (Intercept) -0.3220348 0.04608849  -6.987316  2.801938e-12
## appused app -0.2316567 0.06680031  -3.467899  5.245447e-04
## raceWhite   -1.9344662 0.06675370 -28.979159 1.204740e-184

After controlling for race, there is good statistical evidence (p < 0.001) for the estimated association between app usage and infection. This confirms our training-data findings. There is also very strong statistical evidence (p <<< 0.001) for race’s association with infection risk. Specifically, the estimated odds of infection for Whites were \(\exp(\)-1.934\() =\) 0.145 (95% CI: 0.127, 0.165) times that of Blacks (regardless of app usage).

The estimated odds of infection for app users were \(\exp(\)-0.232\() =\) 0.793 (95% CI: 0.696, 0.905) times that of non-users (regardless of race). The corresponding estimated infection risks for app usage by race are:

risk_didnt_use_app_black_holdout <- plogis(coef(glm_rwe_holdout) %*% c(1, 0, 0))
risk_used_app_black_holdout <- plogis(coef(glm_rwe_holdout) %*% c(1, 1, 0))
risk_didnt_use_app_white_holdout <- plogis(coef(glm_rwe_holdout) %*% c(1, 0, 1))
risk_used_app_white_holdout <- plogis(coef(glm_rwe_holdout) %*% c(1, 1, 1))

• 0.42 for Blacks not using the app
• 0.365 for Blacks using the app
• 0.095 for Whites not using the app
• 0.077 for Whites using the app

The estimated RDs by race are:

rwe_holdout_rd_black <- risk_used_app_black_holdout - risk_didnt_use_app_black_holdout
rwe_holdout_rd_white <- risk_used_app_white_holdout - risk_didnt_use_app_white_holdout
confint_glm_rwe_holdout <- confint(glm_rwe_holdout) # 95% CIs: odds ratios of infection
rwe_holdout_rd_ci_black <- c(
  plogis(confint_glm_rwe_holdout[, 1] %*% c(1, 1, 0)) - plogis(confint_glm_rwe_holdout[, 1] %*% c(1, 0, 0)),
  plogis(confint_glm_rwe_holdout[, 2] %*% c(1, 1, 0)) - plogis(confint_glm_rwe_holdout[, 2] %*% c(1, 0, 0))
)
rwe_holdout_rd_ci_white <- c(
  plogis(confint_glm_rwe_holdout[, 1] %*% c(1, 1, 1)) - plogis(confint_glm_rwe_holdout[, 1] %*% c(1, 0, 1)),
  plogis(confint_glm_rwe_holdout[, 2] %*% c(1, 1, 1)) - plogis(confint_glm_rwe_holdout[, 2] %*% c(1, 0, 1))
)

• -0.055 (95% CI: -0.083, -0.025) for Blacks
• -0.018 (95% CI: -0.022, -0.01) for Whites

Feature distribution

Race: Based on the U.S. Census, the probability of being White was set as \(\Pr(Z=1) = \pi_Z =\) 0.8560091.

Outcome Model

DiagrammeR::grViz("
digraph causal {

  # Nodes
  node [shape = plaintext]
  X [label = 'App \n Usage \n (X)']
  Z [label = 'Race \n (Z)']
  Y [label = 'Infection \n (Y)']

  # Edges
  edge [color = black,
        arrowhead = vee]
  rankdir = LR
  X -> Y
  Z -> Y

  # Graph
  graph [overlap = true]
}")

a <- round((1 - prInfectionGivenBlack) * prBlack, 3)
b <- round(prBlackGivenInfection * prInfection, 3)
c <- round((1 - prInfectionGivenWhite) * prWhite, 3)
d <- round(prWhiteGivenInfection * prInfection, 3)
knitr::kable(
  dplyr::tribble(
    ~Race, ~Uninfected, ~Infected, ~`Row Total`,
    "Black", a, b, a + b,
    "White", c, d, c + d,
    "Column Total", a + c, b + d, a + b + c + d
  )
)

Propensity Model: App Usage

DiagrammeR::grViz("
digraph causal {

  # Nodes
  node [shape = plaintext]
  Z [label = 'Race \n (Z)']
  X [label = 'App \n Usage \n (X)']

  # Edges
  edge [color = black,
        arrowhead = vee]
  rankdir = LR
  Z -> X

  # Graph
  graph [overlap = true]
}")

Causal Model

The full DAG is:

1. App Usage (\(X\)) \(\rightarrow\) Infection (\(Y\))
2. Race (\(Z\)) \(\rightarrow\) Infection (\(Y\))
3. Race (\(Z\)) \(\rightarrow\) App Usage (\(X\))

DiagrammeR::grViz("
digraph causal {

  # Nodes
  node [shape = plaintext]
  X [label = 'App \n Usage \n (X)']
  Z [label = 'Race \n (Z)']
  Y [label = 'Infection \n (Y)']

  # Edges
  edge [color = black,
        arrowhead = vee]
  rankdir = LR
  X -> Y
  Z -> Y
  Z -> X

  # Graph
  graph [overlap = true]
}")