Chapter 13 Causality-roadblocks on the way to it

This chapter contains a discussion of the differences between causality and association. It also covers specific issues that may be encountered in this area when investigating the effects of environmental hazards on health. From this chapter, the reader will have gained an understanding of the following topics:

  • Issues with causality in observational studies.
  • The Bradford–Hill criteria which are a group of minimal conditions necessary to provide adequate evidence of a causal relationship.
  • Ecological bias which may occur when inferences about the nature of individuals are made using aggregated data.
  • The role of exposure variability in determining the extent of ecological bias.
  • Approaches to acknowledging ecological bias in ecological studies.
  • Concentration and exposure response functions.
  • Models for estimating personal exposures including micro-environments.

Solutions to Selected Exercises

Question 13.16: Load the Italy and China data, described in Example 13.3, and included in the supplementary Bookdown material into R.

The data for this exercise follows that of Exercise 13.3, we can load it as follows.

library(tidyverse)
library(janitor)

vonK_Data <- read_csv("data/vonK_China_Italy.csv")
# Aggregate these datasets to have just two age ranges---$<70$ and $\geq 70$.
vonK_Data_70 <-
  mutate(
    vonK_Data,
    .keep = "unused",
    Age_Range = fct_collapse(
      `Age group`,
      Over_70 = c("70-79", "80+"),
      other_level = "Under_70"
    )
  )

# Sum up all death data grouped by age range and country
vonK_Data_70 <- vonK_Data_70 |> group_by(Age_Range, Country) |>
  summarise(
    Fatalities = sum(`Confirmed fatalities`),
    Total_Cases = sum(`Confirmed cases`),
    .groups = "drop_last"
  ) |>
  mutate(Survived = Total_Cases - Fatalities, .before = Total_Cases) |>  ## Convert from fatalities to survivals
  arrange(Country) |> ungroup()

i): Display a 2 by 2 contingency table of deaths/survivors vs. China/Italy aggregated over all age groups.

## Sum up all deaths by country
vonK_Aggregated_Table <- vonK_Data_70 |>
  dplyr::select(-Age_Range) |>
  group_by(Country) |>
  summarise_all(.funs = c(sum)) |>
  adorn_totals("row")
Table 13.1: Aggregated Contingency Table
Country Fatalities Survived Total_Cases
China 1023 43649 44672
Italy 357 7669 8026
Total 1380 51318 52698

ii): Display 2 by 2 contingency tables stratified by the age groups \(<70\) and \(\geq 70\).

## Filter by age range and then sum up deaths
vonK_Under70_Table <- vonK_Data_70 |>
  filter(Age_Range == "Under_70") |>
  dplyr::select(-Age_Range) |>
  adorn_totals("row")
vonK_Over70_Table <- vonK_Data_70 |>
  filter(Age_Range == "Over_70") |>
  dplyr::select(-Age_Range) |>
  adorn_totals("row")
Table 13.2: Under 70 Contingency Table
Country Fatalities Survived Total_Cases
China 503 38843 39346
Italy 41 4668 4709
Total 544 43511 44055
Table 13.3: Over 70 Contingency Table
Country Fatalities Survived Total_Cases
China 520 4806 5326
Italy 316 3001 3317
Total 836 7807 8643

iii): Does a Simpson’s disaggregation (SDis) empirically occur in these two tables?

Following the notation of the textbook (Example 13.3), let \(S\) denote alive and \(\bar{S}\) denote death. Let \(T\) and \(\bar{T}\) denote China and Italy. Finally, let \(C\) and \(\bar{C}\) denote the binarized age groups under 70 and over 70. Recall that \[ P_1 = Pr(S|T), \; P_2 = Pr(S|\tilde{T}), \; R_1 = Pr(S|\bar{C}T), \\ R_2 = Pr(S|\bar{C}\bar{T}), \; Q_1 = Pr(S|CT), \; Q_2 = Pr(S|C\bar{T}). \] A Simpson’s disaggregation occurs when \(P_1 > P_2\), while \(R_1 < R_2\) and \(Q_1 < Q_2\).

We compute the empirical values denoted \(\hat{P}_1, \hat{P}_2, \hat{Q}_1, \hat{Q}_2, \hat{R}_1, \hat{R}_2\):

P_hats = head(vonK_Aggregated_Table$Survived / vonK_Aggregated_Table$Total_Cases,
              -1)
Q_hats = head(vonK_Under70_Table$Survived / vonK_Under70_Table$Total_Cases,
              -1)
R_hats = head(vonK_Over70_Table$Survived / vonK_Over70_Table$Total_Cases,
              -1)

Solution <-
  tibble(
    index = c(1, 2),
    P_hat = P_hats,
    Q_hat = Q_hats,
    R_hat = R_hats
  )

knitr::kable(Solution)
index P_hat Q_hat R_hat
1 0.9770997 0.9872160 0.9023658
2 0.9555196 0.9912933 0.9047332

Indeed, we see that \(\hat{P}_1 = 0.977 > \hat{P}_2 = 0.956\), while \(\hat{Q}_1 = 0.987 < \hat{Q}_2 = 0.991\) and \(\hat{R}_1 = 0.902 < \hat{R}_2 = 0.905\), hence a Simpson’s disaggregation has occured.

Question 13.17: Using the Italy and China data stratified by the age groups \(<70\) and \(\geq 70\), obtain empirical estimates of \(P_1\) and \(P_2\), as defined in an SDis. In the following questions, assume that \(P_1\) and \(P_2\) are the population values.

i): What are the possible ranges of \(\gamma_1\) and \(\gamma_2\) for a controlled SDis to exist (can you visualize it in the rectangle \([0,1]^2\))? Keep all intermediate computations to three decimal precision.

We compute \(R\) as defined in Theorem 2, with \(P_i = \hat{P}_i\).

\[ R = \frac{P_2/(1-P_2)}{P_1/(1-P_1)} = \frac{0.956/0.044}{0.977/0.023} = 0.511. \] The allowable region defined by Theorem 2 is hence:

\[ \Gamma < 0.511 \implies \frac{\gamma_2/(1-\gamma_2)}{\gamma_1/(1-\gamma_1)} < 0.511, \; \gamma_2 \leq \gamma_1. \] One can plug the above inequality into freeware such as (Desmos)[https://www.desmos.com/] or (Wolfram-Alpha)[https://www.wolframalpha.com/] to visualize the region.

ii): Suppose a centered interval for the \(\gamma\)’s are desired, i.e., \(\gamma_1 = 1 - \gamma_2\). What are the allowable values of \(\gamma_1\) and \(\gamma_2\)? Furthermore, what choices (considering three decimal points) yield the tightest interval for \(\beta_i\)?

We can re-write the above by substituting \(\gamma_1 = 1 - \gamma_2\) as

\[ \frac{\gamma_2/(1-\gamma_2)}{(1-\gamma_2)(1-(1-\gamma_2))} = \left(\frac{\gamma_2}{1-\gamma_2}\right)^2 < 0.511, \] which, taking only the positive solution of the square root, yields the allowable interval for \(\gamma_2\) \[ \gamma_2 < 0.417, \] implying that \(\gamma_1 > 0.583\). Being conservative up to three decimal points, the smallest interval is hence \([0.416, 0.584]\).

iii): Compute empirical estimates of \(\beta_1\) and \(\beta_2\), and interpret these values.

Recall that \(\beta_1 = Pr(T|\bar{C}), \beta_2 = Pr(T|C)\). In our context, an empirical estimate of \(\beta_1\) is hence the proportion of Chinese patients over 70, and \(\beta_2\) the proportion of Chinese patients under 70.

beta_hat_1 <-
  vonK_Over70_Table$Total_Cases[1] / vonK_Over70_Table$Total_Cases[3]
beta_hat_2 <-
  vonK_Under70_Table$Total_Cases[2] / vonK_Over70_Table$Total_Cases[3]
beta_hats <- c(beta_hat_1, beta_hat_2)
Solution <- Solution |> add_column(beta_hat = beta_hats)
knitr::kable(Solution)
index P_hat Q_hat R_hat beta_hat
1 0.9770997 0.9872160 0.9023658 0.6162212
2 0.9555196 0.9912933 0.9047332 0.5448340

We obtain \(\hat{\beta}_1 = 0.616\) and \(\hat{\beta}_2 = 0.545\).

iv) (Bonus): What is the tightest centered interval at \(\frac{1}{2}(\hat{\beta}_1 - \hat{\beta}_2)\) for \(\gamma_1\) and \(\gamma_2\) such that a controlled SDis can theoretically exist (assuming population parameters)? How well does it cover the estimated \(\hat{\beta}_i\) (if at all)? Was the empirically observed SDis surprising given these results?

We compute the average:

mean(c(beta_hat_1, beta_hat_2))
## [1] 0.5805276

We wish to have the tightest interval centered at \(0.581\) such that \(\gamma_1, \gamma_2\) satisfy \[ \frac{\gamma_2/(1-\gamma_2)}{\gamma_1/(1-\gamma_1)} < 0.511, \; \gamma_2 \leq \gamma_1. \] We can re-write this as follows. The desired interval is \([0.581-\gamma, 0.581+\gamma]\), setting \(\gamma_2 = 0.581-\gamma, \gamma_1 = 0.581+\gamma\). Using numerical methods accessed via (Wolfram-Alpha)[https://www.wolframalpha.com/], we obtain the inequality \[ \gamma > 0.081 \] This indicates that \(\gamma = 0.080\) is a boundary value (at 3 decimal points), and hence the tightest interval is \([0.501, 0.661]\), which well-covers the estimated \(\hat{\beta}\). This indicates that perhaps the empirical existence of an SDis is not so surprising.

Question 13.18: Stratify the Italy and China data by the age groups \(<50\) and \(\geq 50\). Display the corresponding contingency tables.

Re-loading the data and re-stratifying:

vonK_Data_50 <-
  mutate(
    vonK_Data,
    .keep = "unused",
    Age_Range = fct_collapse(
      `Age group`,
      Over_50 = c("50-59", "60-69", "70-79", "80+"),
      other_level = "Under_50"
    )
  )

i): Does an SDis still occur? Why or why not?

We reproduce the contingency tables of Exercise 13.16.

vonK_Data_50 <- vonK_Data_50 |> group_by(Age_Range, Country) |>
  summarise(
    Fatalities = sum(`Confirmed fatalities`),
    Total_Cases = sum(`Confirmed cases`),
    .groups = "drop_last"
  ) |>
  mutate(Survived = Total_Cases - Fatalities, .before = Total_Cases) |>
  arrange(Country) |> ungroup()

vonK_Aggregated_Table <- vonK_Data_50 |>
  dplyr::select(-Age_Range) |>
  group_by(Country) |>
  summarise_all(.funs = c(sum)) |>
  adorn_totals("row")
vonK_Under50_Table <- vonK_Data_50 |>
  filter(Age_Range == "Under_50") |>
  dplyr::select(-Age_Range) |> adorn_totals("row")
vonK_Over50_Table <- vonK_Data_50 |>
  filter(Age_Range == "Over_50") |>
  dplyr::select(-Age_Range) |>
  adorn_totals("row")

Displaying the tables:

Table 13.4: Aggregated Contingency Table
Country Fatalities Survived Total_Cases
China 1023 43649 44672
Italy 357 7669 8026
Total 1380 51318 52698
Table 13.5: Under 50 Contingency Table
Country Fatalities Survived Total_Cases
China 64 20691 20755
Italy 1 1784 1785
Total 65 22475 22540
Table 13.6: Over 50 Contingency Table
Country Fatalities Survived Total_Cases
China 959 22958 23917
Italy 356 5885 6241
Total 1315 28843 30158

Now, checking if an SDis still occurs.

P_hats <-
  head(vonK_Aggregated_Table$Survived / vonK_Aggregated_Table$Total_Cases,
       -1)
Q_hats <-
  head(vonK_Under50_Table$Survived / vonK_Under50_Table$Total_Cases,
       -1)
R_hats <-
  head(vonK_Over50_Table$Survived / vonK_Over50_Table$Total_Cases,
       -1)

Solution <-
  tibble(
    index = c(1, 2),
    P_hat = P_hats,
    Q_hat = Q_hats,
    R_hat = R_hats
  )

knitr::kable(Solution)
index P_hat Q_hat R_hat
1 0.9770997 0.9969164 0.9599030
2 0.9555196 0.9994398 0.9429579

An SDis no longer occurs. \(\hat{P}_1 = 0.977 > \hat{P}_2 = 0.956\), but \(\hat{R}_1 = 0.960 > \hat{R}_2 = 0.943\) also. Probabilistically: different definitions of the random variables involved change the probabilities (empirical or population), and an SDis is not at all guaranteed to appear, even under the same dataset. Conceptually: see Figure 1 in Von Kügelen et al. Although each 10-year range results in a higher CFR for China, the aggregated data shows a higher CFR for Italy. This is exactly Simpson’s paradox. The same aggregation phenomenon happens in our binarization, resulting in a smaller version of Simpson’s paradox. By aggregating the bins 50-80+, we have produced an estimate where the Italian CFR is higher than the Chinese, resulting in a similar bin to the “Total” one given in Figure 1.

ii): Compute \(\hat{\beta}_1\) and \(\hat{\beta}_2\), and repeat the bonus question of Exercise 13.17. What has changed? Do the changes match your observations?

Note that because the aggregated data does not change from Exercise 13.17, our estimates for \(P_1\) and \(P_2\) also remain the same. Therefore, the allowable region for \(\gamma_i\) does not change: \[ \frac{\gamma_2/(1-\gamma_2)}{\gamma_1/(1-\gamma_1)} < 0.511, \; \gamma_2 \leq \gamma_1. \]

We first compute \(\hat{\beta_i}\) again.

beta_hat_1 <-
  vonK_Over50_Table$Total_Cases[1] / vonK_Over50_Table$Total_Cases[3]
beta_hat_2 <-
  vonK_Under50_Table$Total_Cases[1] / vonK_Over50_Table$Total_Cases[3]
beta_hats <- c(beta_hat_1, beta_hat_2)

print(beta_hats)
## [1] 0.7930566 0.6882088

We obtain \(\hat{\beta}_1 = 0.793\), while \(\hat{\beta}_2 = 0.688\). We compute their mid-point to be \(0.741\). Hence, we look for an interval of the form \([0.741 - \gamma, 0.741 + \gamma]\). Repeating the exercise, we obtain \[ \gamma > 0.063. \] Taking \(\gamma = 0.064\), this results in the interval \([0.677, 0.805]\). our estimated \(\hat{\beta}_i\) does still fall in this interval, but it is near the boundary values. This perhaps suggests that an SDis was less likely to occur in this stratification (and indeed, it does not occur empirically).