Correlation Matrix

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

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

	race	app	infection
race	1.0000	0.0054	-0.3156
app	0.0054	1.0000	-0.0497
infection	-0.3156	-0.0497	1.0000

corrplot::corrplot.mixed(cor(dummy_rct))

Unlike the RWE correlation matrix in Part 1, race is no longer correlated with app. This is because we randomized app usage.

Infections by App Usage (Marginal Model)

Let’s first examine our main relationship of interest, as we did with the training data.

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

df_rct <- with(
  experimental_rct,
  cbind(
    table(app, infection),
    prop.table(table(app, infection), margin = 1) # row proportions
  )
)
rct_rd <- df_rct[2,4] - df_rct[1,4]

rct_rd # empirical RD

## [1] -0.03392808

knitr::kable(df_rct) # row proportions

	0. uninfected	1. infected	0. uninfected	1. infected
didn’t use app	4056	726	0.8481807	0.1518193
used app	4250	568	0.8821088	0.1178912

As in Part 1, there was a lower infection rate among app users: Only 11.8% of users had infections, compared to 15.2% of non-users. However, the empirical RD is -0.034, much closer to the true ATE of -0.031. This is markedly smaller than the “false ATE estimate” of -0.169 (i.e., the RWE holdout empirical RD) that would have misled public health authorities on the true effectiveness of our app at reducing coronavirus infections.

out_fisher_rct <- with(
  experimental_rct,
  fisher.test(app, infection)
)

out_fisher_rct

## 
##  Fisher's Exact Test for Count Data
## 
## data:  app and infection
## p-value = 1.254e-06
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.6623291 0.8415144
## sample estimates:
## odds ratio 
##  0.7466837

In addition, there is strong statistical evidence (i.e., statistical significance) that infections varied by app usage (p << 0.001). (Note that this corresponds to a two-sided hypothesis test; our RCT hypothesis is one-sided.) Here, the estimated odds of infection for app users were 0.747 (i.e., roughly three quarters) that of non-users, with a 95% confidence interval (CI) of (0.662, 0.842).

out_epi.2by2 <- epiR::epi.2by2(
  with(
    experimental_rct,
    table(app == "didn't use app", infection == "0. uninfected")
  )
)

out_epi.2by2$res$ARisk.crude.wald / 100

##           est       lower       upper
## 1 -0.03392808 -0.04757939 -0.02027677

The corresponding estimated RD with 95% CI is -0.034 (95% CI: -0.048, -0.02). (Note that this corresponds to a two-sided hypothesis test; our RCT hypothesis is one-sided.) That is, we’re 95% confident that using the app lowered the risk of infection by somewhere between 0.02 and 0.048. Note that this is now a statement about a causal effect, not a statistical association.

We could stop here because we designed our RCT in order to properly estimate the ATE. To review from Part 1, the ATE is 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.

But in the Explanatory Model section, we’ll go ahead and fit the Part 1 conditional model that included race as a confounder. In the next section (Causal Inference: It’s the Law), we’ll see how to “roll up” the predictions from this model to equal our estimated RD. We’ll then learn how to use this same procedure to estimate the ATE in a statistically consistent way using our original RWE data.

Infections by Race

As in Part 1, race seems to have been associated with infection.

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

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

knitr::kable(df_rct_race_infection) # row proportions

	0. uninfected	1. infected	0. uninfected	1. infected
Black	794	537	0.5965440	0.4034560
White	7512	757	0.9084533	0.0915467

out_fisher_rct_race_infection <- with(
  experimental_rct,
  fisher.test(race, infection)
)

out_fisher_rct_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.1302266 0.1705730
## sample estimates:
## odds ratio 
##  0.1490586

As before, Blacks were more likely to have been infected than Whites (p << 0.001): 40.3% of Blacks were infected, compared to 9.2% of Whites. Put differently, the estimated odds of infection for Whites were 0.149 (95% CI: 0.13, 0.171) times that of Blacks.

App Usage by Race

Unlike our RWE findings in Part 1, race is no longer correlated with app. Again, this is because we randomized app usage.

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

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

knitr::kable(df_rct_race_app) # row proportions

	didn’t use app	used app	didn’t use app	used app
Black	672	659	0.5048835	0.4951165
White	4110	4159	0.4970371	0.5029629

out_fisher_rct_race_app <- with(
  experimental_rct,
  fisher.test(race, app)
)

out_fisher_rct_race_app

## 
##  Fisher's Exact Test for Count Data
## 
## data:  race and app
## p-value = 0.6156
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.917487 1.160591
## sample estimates:
## odds ratio 
##    1.03189

Causal Model

From Part 1, recall that an explanatory model consists of a causal model and a statistical model. The causal model is often specified as a directed acyclic graph (DAG) (Pearl, 2009). We had assumed that the true DAG was:

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.

Check Sample Size

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

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

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

knitr::kable(df_rct_race_app) # row proportions

	didn’t use app	used app	didn’t use app	used app
Black	672	659	0.5048835	0.4951165
White	4110	4159	0.4970371	0.5029629

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

Fit Statistical Model

We fit our RCT logistic model as follows.

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

knitr::kable(summary(glm_rct)$coefficients)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-0.2383012	0.0637708	-3.73684	0.0001863
appused app	-0.3133894	0.0634329	-4.94049	0.0000008
raceWhite	-1.9089811	0.0678516	-28.13467	0.0000000

After controlling for race, there is strong statistical evidence (p << 0.001) for the estimated effect of app usage on infection risk. 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.909\() =\) 0.148 (95% CI: 0.13, 0.169) times that of Blacks (regardless of app usage).

The estimated odds of infection for app users were \(\exp(\)-0.313\() =\) 0.731 (95% CI: 0.645, 0.828) times that of non-users (regardless of race). The corresponding estimated infection risks for app usage by race are:

risk_didnt_use_app_black_rct <- plogis(coef(glm_rct) %*% c(1, 0, 0))
risk_used_app_black_rct <- plogis(coef(glm_rct) %*% c(1, 1, 0))
risk_didnt_use_app_white_rct <- plogis(coef(glm_rct) %*% c(1, 0, 1))
risk_used_app_white_rct <- plogis(coef(glm_rct) %*% c(1, 1, 1))

• 0.441 for Blacks not using the app
• 0.365 for Blacks using the app
• 0.105 for Whites not using the app
• 0.079 for Whites using the app

The estimated RDs by race are:

rct_rd_black <- risk_used_app_black_rct - risk_didnt_use_app_black_rct
rct_rd_white <- risk_used_app_white_rct - risk_didnt_use_app_white_rct
confint_glm_rct <- confint(glm_rct) # 95% CIs: odds ratios of infection
rct_rd_ci_black <- c(
  plogis(confint_glm_rct[, 1] %*% c(1, 1, 0)) - plogis(confint_glm_rct[, 1] %*% c(1, 0, 0)),
  plogis(confint_glm_rct[, 2] %*% c(1, 1, 0)) - plogis(confint_glm_rct[, 2] %*% c(1, 0, 0))
)
rct_rd_ci_white <- c(
  plogis(confint_glm_rct[, 1] %*% c(1, 1, 1)) - plogis(confint_glm_rct[, 1] %*% c(1, 0, 1)),
  plogis(confint_glm_rct[, 2] %*% c(1, 1, 1)) - plogis(confint_glm_rct[, 2] %*% c(1, 0, 1))
)

• -0.075 (95% CI: -0.1, -0.047) for Blacks
• -0.026 (95% CI: -0.028, -0.02) for Whites

Randomized Controlled Trial

We can use our explanatory model’s estimated coefficients to predict the infection risk for each person.

experimental_rct_preds <- experimental_rct %>%
  dplyr::mutate(predicted_risk = predict(object = glm_rct, type = "response"))
tbl_estimated_risks_rct <- experimental_rct_preds %>%
  dplyr::select(race, app, predicted_risk) %>%
  dplyr::distinct() %>%
  dplyr::arrange(race, app)

The predicted risks for each race and app combination are:

knitr::kable(tbl_estimated_risks_rct)

race	app	predicted_risk
Black	didn’t use app	0.4407050
Black	used app	0.3654723
White	didn’t use app	0.1045855
White	used app	0.0786616

These predicted risks are just the estimated infection risks for app usage by race we’d calculated earlier. (Note that while predicted probabilities such as these are used in classification, they aren’t “predictions” in the classification sense. The latter can only take on the values of “used app” or “didn’t use app”, not continuous probabilities. The predicted risks are “predictions” by the standard statistics—not machine learning—definition.)

We “roll up” these estimated risks to the level of app usage by taking their averages over both racial categories for app users and non-users.

tbl_aac_rct <- experimental_rct_preds %>%
  dplyr::group_by(app) %>%
  dplyr::summarize(mean_preds = mean(predicted_risk))
estimated_AAC_rct <- tbl_aac_rct$mean_preds[2] - tbl_aac_rct$mean_preds[1]

knitr::kable(tbl_aac_rct)

app	mean_preds
didn’t use app	0.1518193
used app	0.1178912

But the estimated RD we calculate from these average estimated risks is exactly the same as the empirical RD we’d calculated earlier!

rct_rd # RCT: empirical RD

## [1] -0.03392808

estimated_AAC_rct # RCT: RD from average estimated risks

## [1] -0.03392808

Real-world Evidence

So can’t we just do this with our RWE estimated risks to come up with a similar estimate?

# RWE: empirical RD
df_rwe_holdout <- with(
  observational_rwe_holdout, # see Part 1 to generate this:
    # https://towardsdatascience.com/coronavirus-telemedicine-and-race-part-1-simulated-real-world-evidence-9971f553194d?source=email-8430d9f1992d-1586971057342-layerCake.autoLayerCakeWriterNotification-------------------------90bb4612_8a36_4b93_81dd_1656d841e715&sk=32bfb03e1ca157e423ecf8cb69835ef3
  cbind(
    table(app, infection),
    prop.table(table(app, infection), margin = 1) # row proportions
  )
)
rwe_holdout_rd <- df_rwe_holdout[2,4] - df_rwe_holdout[1,4]

# RWE: RD from average estimated risks
observational_rwe_preds <- observational_rwe_holdout %>%
  dplyr::mutate(predicted_risk = predict(object = glm_rwe_holdout, type = "response"))
tbl_estimated_risks_rwe <- observational_rwe_preds %>%
  dplyr::select(race, app, predicted_risk) %>%
  dplyr::distinct() %>%
  dplyr::arrange(race, app)
tbl_aac_rwe <- observational_rwe_preds %>%
  dplyr::group_by(app) %>%
  dplyr::summarize(mean_preds = mean(predicted_risk))
estimated_AAC_rwe <- tbl_aac_rwe$mean_preds[2] - tbl_aac_rwe$mean_preds[1]

knitr::kable(tbl_estimated_risks_rwe)

race	app	predicted_risk
Black	didn’t use app	0.4201799
Black	used app	0.3650084
White	didn’t use app	0.0947902
White	used app	0.0766925

knitr::kable(tbl_aac_rwe)

app	mean_preds
didn’t use app	0.2598114
used app	0.0904586

rwe_holdout_rd # RWE: empirical RD

## [1] -0.1693528

estimated_AAC_rwe # RWE: RD from average estimated risks

## [1] -0.1693528

Apparently not! What happened?

knitr::kable(tbl_lte_example)

Overall Average

The overall average height is just 5.364. Let \(y_i\) represent the height for individual \(i = 1, ..., n\), where \(n = 8\). This overall average height is explicitly calculated as \(\frac{1}{n} \sum_i^n y_i\), or:

(woman1 + woman2 + man1 + man2 + man3 + man4 + man5 + man6) / 8

## [1] 5.364126

Overall Average as Weighted Average of Components

But the overall average height is also:

((woman1 + woman2) / 2) * (2/8) + ((man1 + man2 + man3 + man4 + man5 + man6) / 6) * (6/8)

## [1] 5.364126

Put differently, but using code similar to what we’d used earlier:

tbl_lte_example_components <- tbl_lte_example %>%
  dplyr::group_by(gender) %>%
  dplyr::summarize(mean_height = mean(height))

knitr::kable(tbl_lte_example_components)

gender	mean_height
0	5.459415
1	5.078258

tbl_lte_example_components[2, 2] * prop.table(table(tbl_lte_example$gender))[2] +
  tbl_lte_example_components[1, 2] * prop.table(table(tbl_lte_example$gender))[1]

##   mean_height
## 1    5.364126

Relationship of Overall Average to Component Weighted Average

Let \(z_i = 1\) for women, and \(z_i = 0\) for men. The average height is:

• \(\left( \sum_i^n y_i z_i \right) / \left( \sum_i^n z_i \right) =\) 5.078 for women
• \(\left\{ \sum_i^n y_i (1 - z_i) \right\} / \left\{ \sum_i^n (1 - z_i) \right\} =\) 5.459 for men

The proportion of each gender is:

• \(\frac{1}{n} \sum_i^n z_i =\) 0.25 for women
• \(\frac{1}{n} \sum_i^n (1 - z_i) =\) 0.75 for men

The calculation above adds these two component averages, but it weights each component average by its respective component proportions. Hence, this weighted average is equal to the overall average height:

\[ \left\{ \frac{\sum_i^n y_i z_i}{\sum_i^n z_i} \times \frac{\sum_i^n z_i}{n} \right\} + \left\{ \frac{\sum_i^n y_i (1 - z_i)}{\sum_i^n (1 - z_i)} \times \frac{\sum_i^n (1 - z_i)}{n} \right\} \] \[ = \frac{\sum_i^n y_i z_i}{n} + \frac{\sum_i^n y_i (1 - z_i)}{n} = \frac{1}{n} \left( \sum_i^n y_i z_i + \sum_i^n y_i - y_i z_i \right) \] \[ = \frac{1}{n} \sum_i^n \left( y_i z_i + y_i - y_i z_i \right) = \frac{1}{n} \sum_i^n y_i \]

## [1] "(5.078 * 0.25) + (5.459 * 0.75) = 5.364"

## [1] "((woman1 + woman2) / 2) * (2/8) + ((man1 + man2 + man3 + man4 + man5 + man6) / 6) * (6/8) = 5.364"

This statement of equivalence is an expression of the Law of Total Expectation (LTE): The overall average is a sum of the component averages weighted by the proportion of each component.

See the Appendix for the definition of the LTE.

Real-world Evidence: App Users

Let’s now apply the LTE to our RWE holdout data. For simplicity, we’ll just look at app users; similar reasoning holds for app non-users.

• We have \(n =\) 9600 individuals.
• Let \(Y_i = 1\) if individual \(i\) was infected, and \(Y_i = 0\) if uninfected.
• Let \(X_i = 1\) if individual \(i\) used the app, and \(X_i = 0\) if they didn’t.
• Let \(Z_i = 1\) if individual \(i\) is White, and \(Z_i = 0\) if they are Black.
• Let \(\hat{R}_i = \widehat{\text{Pr}} ( Y_i=1 | X_i, Z_i )\) represent the estimated infection risk given \(X_i\) and \(Z_i\).

Among app users, the overall estimated risk was \(\left( \sum_i^n Y_i X_i \right) / \left( \sum_i^n X_i \right) =\) 0.09. The average estimated risk was:

• \(\left( \sum_i^n \hat{R}_i X_i Z_i \right) / \left( \sum_i^n X_i Z_i \right) = \widehat{\text{Pr}} ( Y=1 | X=1, Z=1 ) =\) 0.077 for Whites
• \(\left\{ \sum_i^n \hat{R}_i X_i (1 - Z_i) \right\} / \left\{ \sum_i^n X_i (1 - Z_i) \right\} = \widehat{\text{Pr}} ( Y=1 | X=1, Z=0 ) =\) 0.365 for Blacks

The proportion of each race was:

• \(\left( \sum_i^n X_i Z_i \right) / \left( \sum_i^n X_i \right) =\) 0.952 White
• \(\left\{ \sum_i^n X_i (1 - Z_i) \right\} / \left\{ \sum_i^n X_i \right\} =\) 0.048 Black

The weighted average of average estimated risks taken over both race categories is equal to the overall estimated risk:

\[ \left\{ \frac{\sum_i^n \hat{R}_i X_i Z_i}{\sum_i^n X_i Z_i} \times \frac{\sum_i^n X_i Z_i}{\sum_i^n X_i} \right\} + \left\{ \frac{\sum_i^n \hat{R}_i X_i (1 - Z_i)}{\sum_i^n X_i (1 - Z_i)} \times \frac{\sum_i^n X_i (1 - Z_i)}{\sum_i^n X_i} \right\} \] \[ = \frac{\sum_i^n \hat{R}_i X_i}{\sum_i^n X_i} = \frac{\sum_i^n Y_i X_i}{\sum_i^n X_i} \]

## [1] "(0.077 * 0.952) + (0.365 * 0.048) = 0.09 = 0.09"

The last equivalence follows from the definition of a logistic regression. Here’s some code that does this explicitly:

observational_rwe_preds %>%
  dplyr::filter(app == "used app") %>%
  dplyr::mutate(
    weight_Black = (race == "Black") * mean(observational_rwe_holdout$race[observational_rwe_holdout$app == "used app"] == "Black"),
    weight_White = (race == "White") * mean(observational_rwe_holdout$race[observational_rwe_holdout$app == "used app"] == "White"),
    predicted_risk_weighted = predicted_risk * (race == "Black") * weight_Black +
      predicted_risk * (race == "White") * weight_White
  ) %>%
  dplyr::group_by(race) %>%
  dplyr::summarize(mean_preds_weighted_rwe = mean(predicted_risk_weighted)) %>%
  dplyr::summarize(sum_mean_preds_weighted_rwe = sum(mean_preds_weighted_rwe))

## # A tibble: 1 x 1
##   sum_mean_preds_weighted_rwe
##                         <dbl>
## 1                      0.0905

Compare this to our earlier code we used to calculate the empirical RD, which did not explicitly average of race-specific average estimated risks:

observational_rwe_preds %>%
  dplyr::filter(app == "used app") %>%
  dplyr::summarize(mean_preds = mean(predicted_risk))

## # A tibble: 1 x 1
##   mean_preds
##        <dbl>
## 1     0.0905

Randomized Controlled Trial: App Users

Now let’s apply the LTE to our RCT data. Again, we’ll just look at app users. We’ll use the exact same formulas.

Among app users, the overall estimated risk was 0.118. The average estimated risk was:

• 0.079 for Whites
• 0.365 for Blacks

The proportion of each race was:

• 0.863 White
• 0.137 Black

The weighted average of average estimated risks taken over both race categories is equal to the overall estimated risk:

## [1] "(0.079 * 0.863) + (0.365 * 0.137) = 0.118 = 0.118"

Here’s some code that does this explicitly:

experimental_rct_preds %>%
  dplyr::filter(app == "used app") %>%
  dplyr::mutate(
    weight_Black = (race == "Black") * mean(experimental_rct$race[experimental_rct$app == "used app"] == "Black"),
    weight_White = (race == "White") * mean(experimental_rct$race[experimental_rct$app == "used app"] == "White"),
    predicted_risk_weighted = predicted_risk * (race == "Black") * weight_Black +
      predicted_risk * (race == "White") * weight_White
  ) %>%
  dplyr::group_by(race) %>%
  dplyr::summarize(mean_preds_weighted_rct = mean(predicted_risk_weighted)) %>%
  dplyr::summarize(sum_mean_preds_weighted_rct = sum(mean_preds_weighted_rct))

## # A tibble: 1 x 1
##   sum_mean_preds_weighted_rct
##                         <dbl>
## 1                       0.118

Can you spot the main difference between the four corresponding RWE and RCT components? This difference unlocks the key insight of the g-formula.

RCT Substitution

In the RWE data, race influenced the propensity of using the app. As a result, app users were more likely to be White than non-users:

knitr::kable(df_rwe_holdout_race_app_colprop) # column proportions

	didn’t use app	used app	didn’t use app	used app
Black	1667	607	0.5071494	0.0477464
White	1620	12106	0.4928506	0.9522536

However, the proportion of each race among app users and non-users were very similar in the RCT data:

knitr::kable(df_rct_race_app_colprop) # column proportions

	didn’t use app	used app	didn’t use app	used app
Black	672	659	0.140527	0.1367787
White	4110	4159	0.859473	0.8632213

This is because randomizing app usage made it statistically independent of race. Hence, the app-usage-specific proportions approximated the overall race proportions, which were:

• \(\frac{1}{n} \sum_i^n (1 - Z_i) =\) 0.139 Black
• \(\frac{1}{n} \sum_i^n Z_i =\) 0.861 White

That is, app users were:

• \(\left( \sum_i^n X_i Z_i \right) / \left( \sum_i^n X_i \right) =\) 0.863 \(\approx\) 0.861 \(= \frac{1}{n} \sum_i^n Z_i\) White
• \(\left\{ \sum_i^n X_i (1 - Z_i) \right\} / \left\{ \sum_i^n X_i \right\} =\) 0.137 \(\approx\) 0.139 \(= \frac{1}{n} \sum_i^n (1 - Z_i)\) Black

So we can set \(X_i=1\) for all \(i\), such that the race proportions we use are:

• \(\left( \sum_i^n 1 Z_i \right) / \left( \sum_i^n 1 \right) = \frac{1}{n} \sum_i^n Z_i =\) 0.861 White
• \(\left\{ \sum_i^n 1 (1 - Z_i) \right\} / \left\{ \sum_i^n 1 \right\} = \frac{1}{n} \sum_i^n (1 - Z_i) =\) 0.139 Black

This is the same as removing [experimental_rct$app == "used app"] from our earlier code like so:

(
  experimental_rct_preds %>%
    dplyr::filter(app == "used app") %>%
    dplyr::mutate(
      weight_Black = (race == "Black") * mean(experimental_rct$race == "Black"),
      weight_White = (race == "White") * mean(experimental_rct$race == "White"),
      predicted_risk_weighted = predicted_risk * (race == "Black") * weight_Black +
        predicted_risk * (race == "White") * weight_White
    ) %>%
    dplyr::group_by(race) %>%
    dplyr::summarize(mean_preds_weighted_rct = mean(predicted_risk_weighted)) %>%
    dplyr::summarize(sum_mean_preds_weighted_rct = sum(mean_preds_weighted_rct))
)$sum_mean_preds_weighted_rct

## [1] 0.1184267

Statistically, this substitution is acceptable: We are simply replacing one statistically consistent estimate of the true race proportions with another. But what happens if we do this with our RWE data?

RWE Re-weighting

Using our RWE dataset, let’s re-weight the average estimated risks for app users using the overall race proportions instead of the app-usage-specific proportions by removing [observational_rwe_holdout$app == "used app"] from our earlier code like so:

(
  observational_rwe_preds %>%
    dplyr::filter(app == "used app") %>%
    dplyr::mutate(
      weight_Black = (race == "Black") * mean(observational_rwe_holdout$race == "Black"),
      weight_White = (race == "White") * mean(observational_rwe_holdout$race == "White"),
      predicted_risk_weighted = predicted_risk * (race == "Black") * weight_Black +
        predicted_risk * (race == "White") * weight_White
    ) %>%
    dplyr::group_by(race) %>%
    dplyr::summarize(mean_preds_weighted_rwe = mean(predicted_risk_weighted)) %>%
    dplyr::summarize(sum_mean_preds_weighted_rwe = sum(mean_preds_weighted_rwe))
)$sum_mean_preds_weighted_rwe

## [1] 0.1176694

This quantity is much closer to our RCT estimate. We apply this same re-weighting to the app non-users using the following code:

tbl_ate_gf <- observational_rwe_preds %>%
  dplyr::mutate(
    predicted_risk = predict(object = glm_rwe_holdout, type = "response"),
    predicted_risk_X1 = predict(
      object = glm_rwe_holdout,
      newdata = observational_rwe_holdout %>%
        dplyr::mutate(app = "used app"),
      type = "response"
    ),
    predicted_risk_X0 = predict(
      object = glm_rwe_holdout,
      newdata = observational_rwe_holdout %>%
        dplyr::mutate(app = "didn't use app"),
      type = "response"
    ),
    gformula_weight = (race == "Black") * mean(observational_rwe_holdout$race == "Black") +
      (race == "White") * mean(observational_rwe_holdout$race == "White"),
    predicted_risk_weighted = predicted_risk * gformula_weight
  ) %>%
  dplyr::group_by(app, race) %>%
  dplyr::summarize(mean_preds_weighted_rwe = mean(predicted_risk_weighted))
estimated_ATE_gf <- (tbl_ate_gf$mean_preds_weighted_rwe[3] + tbl_ate_gf$mean_preds_weighted_rwe[4]) -
  (tbl_ate_gf$mean_preds_weighted_rwe[1] + tbl_ate_gf$mean_preds_weighted_rwe[2])

We calculate the estimated RD as -0.023. Compare this to our RCT estimated RD of -0.034. Both estimates are much closer to the true ATE of -0.031 than our original RWE estimated RD of -0.169.

The re-weighting procedure can be more generally stated as:

1. Calculate the proportion of each combination of confounder values (i.e., the empirical joint distribution of confounders).
2. Weight the estimated average outcomes (conditional on all confounders) by their corresponding confounder proportions from Step 1.
3. For each would-be intervention group, sum these weighted estimates.

This general procedure is known as the g-formula. The contrast (e.g., difference, ratio) of these standardized sums is a statistically consistent estimate of the ATE.

Causal and Statistical Models

We had assumed that race affects app usage. The relevant part of our earlier DAG is:

1. Race (\(Z\)) \(\rightarrow\) App Usage (\(X\))

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]
}")

Race is said to influence the propensity to use the app. Hence, the statistical model of this relationship is often called a propensity model in the causal inference literature. That is, it models how the would-be intervention is statistically related to its inputs.

The true model for app usage as a function of race is the logistic model listed in the Appendix of Part 1. The propensity to use the app is measured as the probability of using the app. In the context of this causal model, this probability is therefore called a propensity score (Rosenbaum and Rubin, 1983).

In reality, we wouldn’t know this true model. But suppose we’d correctly modeled it and estimated its parameters as follows:

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

knitr::kable(summary(glm_rwe_holdout_ps)$coefficients)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-1.010252	0.0474059	-21.31067	0
raceWhite	3.021527	0.0542881	55.65728	0

Now what?