Add HERS dataset linear regression example with interactions, interspersed with birthweight analysis (Vittinghoff Ch. 4)

DESCRIPTION

@@ -61,6 +61,8 @@ Suggests:
     sjPlot,
     equatiomatic,
     broom (>= 1.0.8),
+    GGally,
+    regress3d,
     lmtest,
     gh,
     lintr,

_sec_hers_data.qmd

@@ -0,0 +1,139 @@
+### Motivating example: `hers` data {.smaller}
+
+:::{.callout-note}
+This section is based on @vittinghoff2e, Chapter 4.
+:::
+
+::: notes
+
+{{< include _subfiles/shared/_sec_hers_intro.qmd >}}
+
+:::
+
+```{r}
+#| eval: false
+#| code-fold: show
+library(haven)
+hers <- haven::read_dta(
+  paste0(
+    "https://regression.ucsf.edu/sites/g/files",
+    "/tkssra6706/f/wysiwyg/home/data/hersdata.dta"
+  )
+)
+```
+
+```{r}
+#| include: false
+library(haven)
+hers <-
+  fs::path_package("rme", "extdata/hersdata.dta") |>
+  read_dta() |>
+  dplyr::mutate(
+    HT = as_factor(HT) |>
+      relevel(ref = "placebo"),
+    statins = as_factor(statins) |>
+      relevel(ref = "no")
+  )
+```
+
+::::: {.panel-tabset}
+
+#### Data as table
+
+```{r}
+#| label: tbl-hers-ch4
+#| tbl-cap: "`hers` data"
+hers |> head()
+```
+
+#### Data as graph
+
+```{r}
+#| label: fig-hers-scatter
+#| fig-cap: >
+#|   `hers` data (@vittinghoff2e):
+#|   age (years) and BMI (kg/m²) vs. baseline LDL (mg/dL),
+#|   colored by statin use.
+library(plotly)
+hers_scatter_data <- hers |>
+  dplyr::filter(!is.na(age), !is.na(BMI), !is.na(LDL))
+plotly::plot_ly(
+  x = hers_scatter_data[["age"]],
+  y = hers_scatter_data[["BMI"]],
+  z = hers_scatter_data[["LDL"]],
+  color = as.character(hers_scatter_data[["statins"]]),
+  colors = c("no" = "steelblue", "yes" = "darkorange"),
+  type = "scatter3d",
+  mode = "markers",
+  marker = list(size = 3, opacity = 0.5)
+) |>
+  plotly::layout(
+    scene = list(
+      xaxis = list(title = "Age (yr)"),
+      yaxis = list(title = "BMI (kg/m²)"),
+      zaxis = list(title = "LDL (mg/dL)")
+    ),
+    legend = list(title = list(text = "Statins"))
+  )
+```
+
+#### Key variables
+
+```{r}
+#| label: fig-hers-key-vars
+#| fig-cap: >
+#|   Key variables in hers data: outcome (LDL),
+#|   treatment (HT), and covariates (BMI, statins, age)
+#| fig-height: 7
+#| fig-width: 8
+library(GGally)
+hers |>
+  dplyr::select(LDL, HT, BMI, statins, age) |>
+  ggpairs(
+    mapping = aes(col = statins),
+    lower = list(continuous = GGally::wrap("points", alpha = 0.3)),
+    columnLabels = c(
+      "LDL (mg/dL)",
+      "HT",
+      "BMI (kg/m²)",
+      "Statins",
+      "Age (yr)"
+    )
+  ) +
+  theme_bw() +
+  theme(legend.position = "bottom")
+```
+
+:::::
+
+---
+
+#### Data notation {.smaller}
+
+::: notes
+Let's define some notation to represent this data:
+:::
+
+-   $Y$: LDL cholesterol (mg/dL)
+-   $A$: age (years)
+-   $B$: BMI (kg/m²)
+-   $T$: hormone therapy treatment assignment
+    ("placebo" or "hormone therapy")
+-   $H$: indicator variable for $T$ = "hormone therapy"
+    -   $H = 0$ if $T$ = "placebo"
+    -   $H = 1$ if $T$ = "hormone therapy"
+-   $U$: statin use ("no" or "yes")
+-   $V$: indicator variable for $U$ = "yes"
+    -   $V = 0$ if $U$ = "no"
+    -   $V = 1$ if $U$ = "yes"
+
+::: notes
+"Placebo" is the **reference level** for the categorical variable $T$,
+and "no" is the **reference level** for statin use $U$.
+The choice of reference level is arbitrary;
+it only affects the interpretation of the intercept and corresponding indicator coefficients.
+
+Since LDL is measured at **baseline** (before the hormone therapy was administered),
+$H$ is not included as a predictor in our regression models for LDL.
+We instead focus on statin use $U$ (and its indicator $V$) as the key grouping variable.
+:::

_sec_hers_lm_diagnostics.qmd

@@ -0,0 +1,68 @@
+### Diagnostics for `hers` parallel-planes model
+
+#### Residuals vs fitted for `hers_lm1`
+
+```{r}
+#| label: fig-hers-resid-fitted
+#| fig-cap: "Residuals vs fitted values for `hers_lm1` (parallel planes model)"
+library(ggplot2)
+hers_diag <- hers |>
+  dplyr::mutate(
+    .fitted = fitted(hers_lm1),
+    .resid = residuals(hers_lm1)
+  )
+
+ggplot(hers_diag, aes(x = .fitted, y = .resid)) +
+  geom_point(alpha = 0.3) +
+  geom_hline(yintercept = 0, linetype = "dashed") +
+  facet_wrap(~statins, labeller = label_both) +
+  xlab("Fitted values") +
+  ylab("Residuals") +
+  theme_bw()
+```
+
+#### QQ plot for `hers_lm1`
+
+```{r}
+#| label: fig-hers-qq
+#| fig-cap: "QQ plot of residuals for `hers_lm1` (parallel planes model)"
+ggplot(hers_diag, aes(sample = .resid)) +
+  stat_qq() +
+  stat_qq_line() +
+  facet_wrap(~statins, labeller = label_both) +
+  theme_bw()
+```
+
+### Diagnostics for `hers` interaction model
+
+#### Residuals vs fitted for `hers_lm2`
+
+```{r}
+#| label: fig-hers-resid-fitted-lm2
+#| fig-cap: "Residuals vs fitted values for `hers_lm2` (interaction model)"
+hers_diag2 <- hers |>
+  dplyr::mutate(
+    .fitted = fitted(hers_lm2),
+    .resid = residuals(hers_lm2)
+  )
+
+ggplot(hers_diag2, aes(x = .fitted, y = .resid)) +
+  geom_point(alpha = 0.3) +
+  geom_hline(yintercept = 0, linetype = "dashed") +
+  facet_wrap(~statins, labeller = label_both) +
+  xlab("Fitted values") +
+  ylab("Residuals") +
+  theme_bw()
+```
+
+#### QQ plot for `hers_lm2`
+
+```{r}
+#| label: fig-hers-qq-lm2
+#| fig-cap: "QQ plot of residuals for `hers_lm2` (interaction model)"
+ggplot(hers_diag2, aes(sample = .resid)) +
+  stat_qq() +
+  stat_qq_line() +
+  facet_wrap(~statins, labeller = label_both) +
+  theme_bw()
+```

_sec_hers_lm_gof.qmd

@@ -0,0 +1,15 @@
+### Goodness of fit for `hers` models
+
+#### AIC and BIC for `hers` models
+
+```{r}
+AIC(hers_lm1, hers_lm2)
+BIC(hers_lm1, hers_lm2)
+```
+
+#### Deviance for `hers` models
+
+```{r}
+deviance(hers_lm1)
+deviance(hers_lm2)
+```

_sec_hers_lm_interact.qmd

@@ -0,0 +1,123 @@
+### Interactions in `hers` data {.smaller}
+
+::: notes
+What if the slope of LDL with respect to BMI
+differs depending on age?
+Then we need an "interaction" between age $A$ and BMI $B$:
+:::
+
+$$
+\ba
+Y|A,B &\sciid N(\mu(A,B), \sigma^2)\\
+\mu(a,b) &= \beta_0 + \beta_A a + \beta_B b + \beta_{AB}(a \cdot b)
+\ea
+$$ {#eq-hers-interact}
+
+::: notes
+Now the slope of mean LDL with respect to BMI $B$
+depends on age $A$:
+
+$$
+\ba
+\deriv{b}\mu(A=\red{0}, B=b) &= \beta_B + \beta_{AB} \cdot \red{0} = \beta_B \\
+\deriv{b}\mu(A=\red{a}, B=b) &= \beta_B + \beta_{AB} \red{a}
+\ea
+$$
+
+So the slope of LDL with respect to BMI
+changes by $\beta_{AB}$ for each one-year increase in age.
+:::
+
+```{r}
+#| label: tbl-hers-lm2
+#| tbl-cap: "HERS interaction model"
+hers_lm2 <- lm(
+  LDL ~ age + BMI + age:BMI,
+  data = hers,
+  na.action = na.exclude
+)
+hers_plot_data <- hers |>
+  dplyr::filter(!is.na(age), !is.na(BMI), !is.na(LDL))
+
+hers_lm2 |>
+  parameters::parameters() |>
+  parameters::print_md(
+    select = "{estimate}"
+  )
+```
+
+::::: {.panel-tabset}
+
+#### Statins: No
+
+:::{#fig-hers-interact-fit-no}
+
+```{r}
+#| code-fold: true
+hers_no <- hers_plot_data |> dplyr::filter(statins == "no")
+
+plotly::plot_ly(
+  x = hers_no[["age"]],
+  y = hers_no[["BMI"]],
+  z = hers_no[["LDL"]],
+  type = "scatter3d",
+  mode = "markers",
+  name = "No statins",
+  marker = list(size = 3, opacity = 0.3, color = "steelblue")
+) |>
+  regress3d::add_3d_surface(
+    model = hers_lm2,
+    data = hers_plot_data,
+    showlegend = TRUE
+  ) |>
+  plotly::layout(
+    scene = list(
+      xaxis = list(title = "Age (yr)"),
+      yaxis = list(title = "BMI (kg/m²)"),
+      zaxis = list(title = "LDL (mg/dL)")
+    )
+  )
+```
+
+Interaction model regression surface for `hers` data
+(patients not taking statins)
+
+:::
+
+#### Statins: Yes
+
+:::{#fig-hers-interact-fit-yes}
+
+```{r}
+#| code-fold: true
+hers_yes <- hers_plot_data |> dplyr::filter(statins == "yes")
+
+plotly::plot_ly(
+  x = hers_yes[["age"]],
+  y = hers_yes[["BMI"]],
+  z = hers_yes[["LDL"]],
+  type = "scatter3d",
+  mode = "markers",
+  name = "Yes statins",
+  marker = list(size = 3, opacity = 0.3, color = "darkorange")
+) |>
+  regress3d::add_3d_surface(
+    model = hers_lm2,
+    data = hers_plot_data,
+    showlegend = TRUE
+  ) |>
+  plotly::layout(
+    scene = list(
+      xaxis = list(title = "Age (yr)"),
+      yaxis = list(title = "BMI (kg/m²)"),
+      zaxis = list(title = "LDL (mg/dL)")
+    )
+  )
+```
+
+Interaction model regression surface for `hers` data
+(patients taking statins)
+
+:::
+
+:::::

_sec_hers_lm_model_selection.qmd

@@ -0,0 +1,8 @@
+### Model selection for `hers` data
+
+#### Comparing HERS models using LRT
+
+```{r}
+library(lmtest)
+lrtest(hers_lm1, hers_lm2)
+```