Merge pull request #201 from r-causal/revert-189-update-dag-ch

DESCRIPTION

@@ -27,7 +27,6 @@ Imports:
     kableExtra,
     lubridate,
     MatchIt,
-    patchwork,
     propensity (>= 0.0.0.9000),
     PSW,
     ragg,

R/setup.R

@@ -31,7 +31,7 @@ theme_dag <- function() {
 }
 
 geom_dag_label_repel <- function(..., seed = 10) {
-  ggdag_geom_dag_label_repel(
+  ggdag::geom_dag_label_repel(
     aes(x, y, label = label),
     box.padding = 3.5,
     inherit.aes = FALSE,

_freeze/chapters/chapter-13/execute-results/html.json

@@ -1,10 +1,8 @@
 {
-  "hash": "f63c78940c149a1b3d9159a74137a66a",
+  "hash": "1d5c11e4eb1e56d478902dc918373b89",
   "result": {
-    "markdown": "# Continuous exposures\n\n\n\n\n\n## Calculating propensity scores for continuous exposures\n\nPropensity scores generalize to many other types of exposures, including continuous exposures.\nAt its heart, the workflow is the same: we fit a model where the exposure is the outcome and then use that model to weight a second outcome model.\nFor continuous exposures, linear regression is the simplest way to create propensities.\nInstead of probabilities, we use the cumulative density function.\nThen, we use this density to weight the outcome model.\n\nLet's take a look at an example.\nIn the `touringplans` data set, we have information about the posted waiting times for rides.\nWe also have a limited amount of data on the observed actual times.\nThe question we will consider is this: Do posted wait times for the Seven Dwarfs Mine Train at 8 am affect actual wait times at 9 am?\nHere's our DAG:\n\n\n::: {.cell}\n::: {.cell-output-display}\n![Proposed DAG for the relationship between posted wait in the morning at a particular park and the average wait time between 5 pm and 6 pm.](chapter-13_files/figure-html/fig-dag-avg-wait-1.png){#fig-dag-avg-wait width=672}\n:::\n:::\n\n\nIn @fig-dag-avg-wait, we're assuming that our primary confounders are when the park closes, historic high temperatures, whether or not the ride has extra magic morning hours, and the ticket season.\nThis is the only minimal adjustment set in the DAG, as well.\nThe confounders precede the exposure and outcome, and (by definition) the exposure precedes the outcome.\nThe average posted wait time is, in theory, a manipulable exposure because the park could post a time different from what they expect.\nThe adjustment set\n\nThe model is similar to the binary exposure case, but we need to use linear regression, as the posted time is a continuous variable.\nSince we're not using probabilities, we'll calculate denominators for our weights from a normal density.\nWe then calculate the denominator using the `dnorm()` function, which calculates the normal density for the `exposure`, using `.fitted` as the mean and `mean(.sigma)` as the SD.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlm(\n  exposure ~ confounder_1 + confounder_2,\n  data = df\n) |>\n  augment(data = df) |>\n  mutate(\n    denominator = dnorm(exposure, .fitted, mean(.sigma, na.rm = TRUE))\n  )\n```\n:::\n\n\n## Diagnostics and stabilization\n\nContinuous exposure weights, however, are very sensitive to modeling choices.\nOne problem, in particular, is the existence of extreme weights, an issue that can also affect other types of exposures.\nWhen some observations have extreme weights, the propensities are *destabilized,* which results in wider confidence intervals.\nWe can stabilize them using the marginal distribution of the exposure.\nA common way to calculate the marginal distribution for propensity scores is to use a regression model with no predictors.\n\n::: callout-caution\nExtreme weights destabilize estimates, resulting in wider confidence intervals.\nExtreme weights can be an issue for any time of weight (including those for binary and other types of exposures) that is not bounded.\nBounded weights like the ATO (which are bounded to 0 and 1) do not have this problem, however---one of their many benefits.\n:::\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# for continuous exposures\nlm(\n  exposure ~ 1,\n  data = df\n) |>\n  augment(data = df) |>\n  transmute(\n    numerator = dnorm(exposure, .fitted, mean(.sigma, na.rm = TRUE))\n  )\n\n# for binary exposures\nglm(\n  exposure ~ 1,\n  data = df,\n  family = binomial()\n) |>\n  augment(type.predict = \"response\", data = df) |>\n  select(numerator = .fitted)\n```\n:::\n\n\nThen, rather than inverting them, we calculate the weights as `numerator / denominator`.\nLet's try it out on our posted wait times example.\nFirst, let's wrangle our data to address our question: do posted wait times at 8 affect actual weight times at 9?\nWe'll join the baseline data (all covariates and posted wait time at 8) with the outcome (average actual time).\nWe also have a lot of missingness for `wait_minutes_actual_avg`, so we'll drop unobserved values for now.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(touringplans)\neight <- seven_dwarfs_train_2018 |>\n  filter(wait_hour == 8) |>\n  select(-wait_minutes_actual_avg)\n\nnine <- seven_dwarfs_train_2018 |>\n  filter(wait_hour == 9) |>\n  select(park_date, wait_minutes_actual_avg)\n\nwait_times <- eight |>\n  left_join(nine, by = \"park_date\") |>\n  drop_na(wait_minutes_actual_avg)\n```\n:::\n\n\nFirst, let's calculate our denominator model.\nWe'll fit a model using `lm()` for `wait_minutes_posted_avg` with our covariates, then use the fitted predictions of `wait_minutes_posted_avg` (`.fitted`) to calculate the density using `dnorm()`.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(broom)\ndenominator_model <- lm(\n  wait_minutes_posted_avg ~\n    park_close + park_extra_magic_morning + park_temperature_high + park_ticket_season,\n  data = wait_times\n)\n\ndenominators <- denominator_model |>\n  augment(data = wait_times) |>\n  mutate(\n    denominator = dnorm(\n      wait_minutes_posted_avg,\n      .fitted,\n      mean(.sigma, na.rm = TRUE)\n    )\n  ) |>\n  select(park_date, denominator, .fitted)\n```\n:::\n\n\nWhen we only use the inverted values of `denominator`, we end up with several extreme weights:\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndenominators |>\n  mutate(wts = 1 / denominator) |>\n  ggplot(aes(wts)) +\n  geom_histogram(fill = \"#E69F00\", color = \"white\", bins = 50) +\n  scale_x_log10(name = \"weights\")\n```\n\n::: {.cell-output-display}\n![A histogram of the inverse probability weights for posted waiting time. Weights for continuous exposures are prone to extreme values, which can unstabilize estimates and variance.](chapter-13_files/figure-html/fig-hist-sd-unstable-1.png){#fig-hist-sd-unstable width=672}\n:::\n:::\n\n\nIn @fig-hist-sd-unstable, we see several weights over 100 and one over 10,000; these extreme weights will put undue stress on specific points, complicating the results we will estimate.\n\nLet's now fit the marginal density to use for stabilized weights:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nnumerator_model <- lm(\n  wait_minutes_posted_avg ~ 1,\n  data = wait_times\n)\n\nnumerators <- numerator_model |>\n  augment(data = wait_times) |>\n  mutate(\n    numerator = dnorm(\n      wait_minutes_posted_avg,\n      .fitted,\n      mean(.sigma, na.rm = TRUE)\n    )\n  ) |>\n  select(park_date, numerator)\n```\n:::\n\n\nWe also need to join the fitted values back to our original data set by date, then calculate the stabilized weights (`swts`) using `numerator / denominator`.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nwait_times_wts <- wait_times |>\n  left_join(numerators, by = \"park_date\") |>\n  left_join(denominators, by = \"park_date\") |>\n  mutate(swts = numerator / denominator)\n```\n:::\n\n\nThe stabilized weights are much less extreme.\nStabilized weights should have a mean close to 1 (in this example, it is `round(mean(wait_times_wts$swts), digits = 2)`); when that is the case, then the pseudo-population (that is, the equivalent number of observations after weighting) is equal to the original sample size.\nIf the mean is far from 1, we may have issues with model misspecification or positivity violations [@hernan2021].\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(wait_times_wts, aes(swts)) +\n  geom_histogram(fill = \"#E69F00\", color = \"white\", bins = 50) +\n  scale_x_log10(name = \"weights\")\n```\n\n::: {.cell-output-display}\n![A histogram of the stabilized inverse probability weights for posted waiting time. These weights are much more reasonable and will allow the outcome model to behave better.](chapter-13_files/figure-html/fig-hist-sd-stable-1.png){#fig-hist-sd-stable width=672}\n:::\n:::\n\n\nWhen we compare the exposure---average posted wait times---to the standardized weights, we still have one exceptionally high weight.\nIs this a problem, or is this a valid data point?\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(wait_times_wts, aes(wait_minutes_posted_avg, swts)) +\n  geom_point(size = 3, color = \"grey80\", alpha = 0.7) +\n  geom_point(\n    data = function(x) filter(x, swts > 10),\n    color = \"firebrick\",\n    size = 3\n  ) +\n  geom_text(\n    data = function(x) filter(x, swts > 10),\n    aes(label = park_date),\n    size = 5,\n    hjust = 0,\n    nudge_x = -15.5,\n    color = \"firebrick\"\n  ) +\n  scale_y_log10() +\n  labs(x = \"Average Posted Wait\", y = \"Stabilized Weights\")\n```\n\n::: {.cell-output-display}\n![A scatter of the stabilized inverse probability weights for posted waiting time vs. posted waiting times. Days with more values of `wait_minutes_posted_avg` farther from the mean appear to be downweighted, with a few exceptions. The most unusual weight is for June 23, 2018.](chapter-13_files/figure-html/fig-stabilized-wts-scatter-1.png){#fig-stabilized-wts-scatter width=672}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nwait_times_wts |>\n  filter(swts > 10) |>\n  select(\n    park_date, wait_minutes_posted_avg, .fitted,\n    park_close, park_extra_magic_morning,\n    park_temperature_high, park_ticket_season\n  ) |>\n  knitr::kable()\n```\n\n::: {.cell-output-display}\n\n\n|park_date  | wait_minutes_posted_avg| .fitted|park_close | park_extra_magic_morning| park_temperature_high|park_ticket_season |\n|:----------|-----------------------:|-------:|:----------|------------------------:|---------------------:|:------------------|\n|2018-06-23 |                      81|    28.1|24:00:00   |                        0|                 91.36|regular            |\n\n\n:::\n:::\n\n\nOur model predicted a much lower posted wait time than observed, so this date was upweighted.\nWe don't know why the posted time was so high (the actual time was much lower), but we did find an artist rendering from that date of [Pluto digging for Seven Dwarfs Mine Train treasure](https://disneyparks.disney.go.com/blog/2018/06/disney-doodle-pluto-sniffs-out-fun-at-seven-dwarfs-mine-train/).\n\n## Fitting the outcome model for continuous exposures\n",
-    "supporting": [
-      "chapter-13_files/figure-html"
-    ],
+    "markdown": "# Categorical exposures\n\n\n\n\n\n## Calculating propensity scores for categorical exposures\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] -0.6711  0.9282  0.9044 -0.1768  1.4206\n```\n\n\n:::\n:::\n\n\n## Diagnostics with many categories\n\n## Fitting the outcome model again\n",
+    "supporting": [],
     "filters": [
       "rmarkdown/pagebreak.lua"
     ],

_freeze/chapters/chapter-14/execute-results/html.json

@@ -1,8 +1,10 @@
 {
-  "hash": "1d5c11e4eb1e56d478902dc918373b89",
+  "hash": "22032a0e6e14411b69033ee514dc2692",
   "result": {
-    "markdown": "# Categorical exposures\n\n\n\n\n\n## Calculating propensity scores for categorical exposures\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1]  0.55045 -0.09999  1.27700 -0.90324 -0.65693\n```\n\n\n:::\n:::\n\n\n## Diagnostics with many categories\n\n## Fitting the outcome model again\n",
-    "supporting": [],
+    "markdown": "# Missingness\n\n\n\n\n\n## Missingness as structural bias\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] -0.5199 -1.1160  1.4866  1.0149 -1.4796\n```\n\n\n:::\n:::\n\n\n## Multiple Imputation\n\n## Combining MICE and IPW\n",
+    "supporting": [
+      "chapter-14_files"
+    ],
     "filters": [
       "rmarkdown/pagebreak.lua"
     ],

_freeze/chapters/chapter-15/execute-results/html.json

@@ -1,8 +1,10 @@
 {
   "hash": "d24262a4d5a1dbb0a621f98b43effcb4",
   "result": {
-    "markdown": "# G-Computation {#sec-g-comp}\n\n\n\n\n\n## The Parametric G-Formula\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] -1.4226 -1.1423  1.3277  2.0461  0.3256\n```\n\n\n:::\n:::\n\n\n## Calculating estimates with G-Computation\n\n## The Natural Course\n",
-    "supporting": [],
+    "markdown": "# G-Computation {#sec-g-comp}\n\n\n\n\n\n## The Parametric G-Formula\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1]  0.5788 -1.2622 -0.9627 -0.9578  0.9935\n```\n\n\n:::\n:::\n\n\n## Calculating estimates with G-Computation\n\n## The Natural Course\n",
+    "supporting": [
+      "chapter-15_files"
+    ],
     "filters": [
       "rmarkdown/pagebreak.lua"
     ],

_freeze/chapters/chapter-16/execute-results/html.json

@@ -1,8 +1,10 @@
 {
-  "hash": "05f0f21481f04bb0be23917f8f78e338",
+  "hash": "0c4033870d4c2666b8e9a897aaa34986",
   "result": {
-    "markdown": "# Interaction {#sec-interaction}\n\n\n\n\n\n## Functional form, hetereogenous effects, and joint causes\n\n## Fitting interaction terms in causal models\n\n<!-- Maybe also: -->\n<!-- ## type ii selection bias, induced by EM [@lu2022] -->",
-    "supporting": [],
+    "markdown": "# Causal survival analysis\n\n\n\n\n\n## Preparing data for survival analysis\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] -0.1047 -0.7075 -0.4694  1.2174 -1.6054\n```\n\n\n:::\n:::\n\n\n## Pooled logistic regression\n\n## Confidence intervals for causal survival models\n",
+    "supporting": [
+      "chapter-16_files"
+    ],
     "filters": [
       "rmarkdown/pagebreak.lua"
     ],

_freeze/chapters/chapter-17/execute-results/html.json

@@ -1,8 +1,10 @@
 {
-  "hash": "ddc7311c8eb85246136331989d1ca532",
+  "hash": "4295ad1482e0a56a886c07d3b85f2c3b",
   "result": {
-    "markdown": "# Missingness {#sec-missingness}\n\n\n\n\n\n## Missingness as structural bias\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1]  2.6365  1.8457 -0.4574  0.1535  1.1835\n```\n\n\n:::\n:::\n\n\n## Multiple Imputation\n\n## Combining MICE and IPW\n",
-    "supporting": [],
+    "markdown": "# Causal mediation analysis\n\n\n\n\n\n## Estimating direct effects\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 1.64179 1.65583 0.22422 0.02398 0.86649\n```\n\n\n:::\n:::\n\n\n## Estimating indirect effects\n\n## Allowing for interactions\n",
+    "supporting": [
+      "chapter-17_files"
+    ],
     "filters": [
       "rmarkdown/pagebreak.lua"
     ],

_freeze/chapters/chapter-18/execute-results/html.json

@@ -1,8 +1,10 @@
 {
-  "hash": "b9dc2626f7c608a76f7c7caba0addd09",
+  "hash": "f1d3d81eeb5540167668bf71aebc6749",
   "result": {
-    "markdown": "# Causal inference across time {#sec-longitudinal}\n\n\n\n\n\n## Loss to follow-up\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] -0.3767 -0.6036 -0.6856  0.2084 -0.2844\n```\n\n\n:::\n:::\n\n\n## Time-varying confounding and exposures\n\n",
-    "supporting": [],
+    "markdown": "# Sensitivity analysis\n\n## Quantitative bias analyses\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] -0.6136 -0.3190  0.7672  0.3058  0.4258\n```\n\n\n:::\n:::\n\n\n## Tipping point analyses\n",
+    "supporting": [
+      "chapter-18_files"
+    ],
     "filters": [
       "rmarkdown/pagebreak.lua"
     ],

_freeze/chapters/chapter-19/execute-results/html.json

@@ -1,8 +1,10 @@
 {
-  "hash": "b054d57b42451b94833be81a7144a2c7",
+  "hash": "e048fceeeb7a0c5e7ada95e14bb9e9e4",
   "result": {
-    "markdown": "# Causal survival analysis {#sec-survival}\n\n\n\n\n\n## Preparing data for survival analysis\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1]  0.5082 -0.3546 -0.2759  0.2368  0.1820\n```\n\n\n:::\n:::\n\n\n## Pooled logistic regression\n\n## Confidence intervals for causal survival models\n",
-    "supporting": [],
+    "markdown": "# Other approaches to causal inference\n\n\n\n\n\n## Targeted Learning\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1]  1.0362  2.1558 -0.2929  0.8677  0.2885\n```\n\n\n:::\n:::\n\n\n## Instrumental variable analysis\n\n## Regression discontinuity\n\n## Difference-in-Difference\n",
+    "supporting": [
+      "chapter-19_files"
+    ],
     "filters": [
       "rmarkdown/pagebreak.lua"
     ],

_freeze/chapters/chapter-20/execute-results/html.json

@@ -1,7 +1,7 @@
 {
-  "hash": "35e28ba55bdcd7229ce521d2f04664ff",
+  "hash": "e048fceeeb7a0c5e7ada95e14bb9e9e4",
   "result": {
-    "markdown": "# Causal mediation analysis {#sec-mediation}\n\n\n\n\n\n## Estimating direct effects\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] -0.65327  1.14621 -1.26733 -0.25579  0.04294\n```\n\n\n:::\n:::\n\n\n## Estimating indirect effects\n\n## Allowing for interactions\n",
+    "markdown": "# Other approaches to causal inference\n\n\n\n\n\n## Targeted Learning\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n[1] 0.9632 1.2420 0.6826 1.3024 0.6905\n```\n:::\n:::\n\n\n## Instrumental variable analysis\n\n## Regression discontinuity\n\n## Difference-in-Difference\n",
     "supporting": [],
     "filters": [
       "rmarkdown/pagebreak.lua"