---
title: "Visualizing Linear Models and Mean Comparisons"
author: "Michael Hallquist, PSYC 859"
date: last-modified
week: 12
download: "materials/downloads/w12_models/visualizing_lm_and_means.qmd"
format:
  html:
    code-fold: show
execute:
  warning: false
  message: false
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.width = 10, fig.height = 7, echo = TRUE)
```

# Overview

Statistical models are tools for summarizing data: they use parameters to capture the processes and individual differences that give rise to observed values. But a table of coefficients and *p*-values is only part of the story. Graphical methods let us *see* what a model predicts, compare those predictions to the observed data, and spot places where the model breaks down.

A useful rule of thumb from the class slides: **If you don't know how to visualize the fit and misfit of a statistical model, you shouldn't move ahead with interpretation.**

This walkthrough demonstrates how to visualize predictions, coefficients, and diagnostics for the models that appear most frequently in psychology: linear regression, *t*-tests, and ANOVA. We use three complementary R packages that each address a different part of the visualization problem:

| Package | What it does | Key functions |
|---|---|---|
| **marginaleffects** | Model predictions, comparisons, and slopes on the response scale | `predictions()`, `plot_predictions()`, `comparisons()`, `plot_slopes()` |
| **sjPlot** | Publication-ready coefficient plots and model tables | `plot_model()`, `tab_model()` |
| **performance** | Assumption diagnostics in a single dashboard | `check_model()`, `check_normality()`, `check_heteroscedasticity()` |

## Learning goals

After working through this document, you should be able to:

1. Explain why visualizing model predictions alongside raw data is essential for interpretation.
2. Use `marginaleffects` to generate and plot predicted values, contrasts, and marginal effects.
3. Use `sjPlot` to create forest plots of coefficients and publication-quality model tables.
4. Use `performance::check_model()` to produce a multi-panel diagnostic dashboard for any `lm` object.
5. Recognize common violations (nonlinearity, heteroscedasticity, influential observations) in diagnostic plots.

## Packages

```{r}
# Install any missing packages (uncomment if needed):
# install.packages(c(
#   "ggplot2", "dplyr", "tidyr", "psych", "MASS",
#   "marginaleffects", "sjPlot", "performance", "see",
#   "car", "patchwork"
# ))

library(ggplot2)
library(dplyr)
library(tidyr)
library(psych)
library(MASS)
library(marginaleffects)
library(sjPlot)
library(performance)
library(see) # required for performance plot methods
library(car)
library(patchwork)
```

# Data preparation: The Big Five Inventory

The `bfi` dataset from the `psych` package contains 2,800 responses to 25 personality items measuring the Big Five traits: Agreeableness (A1–A5), Conscientiousness (C1–C5), Extraversion (E1–E5), Neuroticism (N1–N5), and Openness (O1–O5). It also includes `gender` (1 = Male, 2 = Female), `education` (1–5), and `age`.

We will score the Agreeableness and Neuroticism scales and use them as outcomes.

```{r}
data(bfi, package = "psych")

# Score agreeableness and neuroticism.
# Some items are reverse-keyed; the psych package documents which ones.
# A1 is reverse-keyed; so are E1, E2, O2, O5, C4, C5.
bfi_scored <- bfi %>%
  mutate(
    # Reverse-key relevant items (on 1-6 scale)
    A1r = 7 - A1,
    # Score agreeableness as the mean of the 5 items
    agreeableness = rowMeans(pick(A1r, A2, A3, A4, A5), na.rm = TRUE),
    neuroticism = rowMeans(pick(N1, N2, N3, N4, N5), na.rm = TRUE),
    gender_f = factor(gender, levels = 1:2, labels = c("Male", "Female")),
    education_f = factor(education, levels = 1:5, labels = c(
      "Some HS", "HS grad", "Some college", "College grad", "Graduate"
    ))
  ) %>%
  filter(!is.na(agreeableness), !is.na(age), !is.na(gender_f))

glimpse(bfi_scored)
```

# Simple linear regression

We start with the simplest possible model: does agreeableness change with age?

## Fit the model

```{r}
m_simple <- lm(agreeableness ~ age, data = bfi_scored)
summary(m_simple)
```

## Using ggplot to plot the linear relationship 

Most people start here: `geom_smooth(method = "lm")` draws a regression line directly from the data, while `geom_point()` shows the raw data for comparison.

```{r}
ggplot(bfi_scored, aes(x = age, y = agreeableness)) +
  geom_point(alpha = 0.15, size = 1.5) +
  geom_smooth(method = "lm", color = "#2c7fb8", fill = "#7fcdbb") +
  labs(
    title = "Agreeableness by age",
    subtitle = "Simple OLS with 95% confidence band",
    x = "Age (years)",
    y = "Agreeableness (mean of 5 items)"
  ) +
  theme_minimal(base_size = 14)
```

This is quick and useful, but it has limits. The band shows uncertainty about the *regression line* (mean prediction), not about individual data points. And as the model gets more complex — multiple predictors, interactions — `geom_smooth` can't keep up.

## Using `marginaleffects::plot_predictions()`

The `marginaleffects` package provides a general-purpose approach that works for *any* model:

```{r}
plot_predictions(m_simple, condition = "age") +
  labs(
    title = "Model-predicted agreeableness by age",
    subtitle = "marginaleffects::plot_predictions()",
    y = "Predicted agreeableness"
  ) +
  theme_minimal(base_size = 14)
```

The strength of `plot_predictions()` is that it works identically whether the model is a simple `lm`, a generalized linear model, or a multilevel model. The same code, the same logic.

## Residual diagnostics with `performance::check_model()`

The `check_model()` function produces a multi-panel dashboard of the most important diagnostic plots for any regression model. It replaces the need to remember which base R `plot()` panel is which.

```{r}
#| fig-width: 12
#| fig-height: 10
check_model(m_simple)
```

What to look for in each panel:

Even for ordinary least squares models, `performance` includes a posterior predictive-style panel by simulating replicated outcome data from the fitted model and comparing those simulated outcomes to the observed outcome distribution.

1. **Posterior predictive check**: Do simulated datasets from the model look like the observed data?
2. **Linearity**: Are residuals scattered randomly around zero, or is there a pattern (e.g., curve)?
3. **Homogeneity of variance (scale-location)**: Does the spread of residuals change as fitted values increase?
4. **Influential observations**: In this plot, the x-axis shows leverage (how unusual an observation's predictor values are) and the y-axis shows standardized residuals (how far the observed outcome is from the model's prediction, in standard deviation units). Observations with both high leverage (far right) and large standardized residuals (far from zero) are potentially influential because they can disproportionately affect the fitted model. The dashed lines indicate Cook's distance thresholds; observations beyond these lines may be unduly influential on the model estimates.
5. **Normality of residuals (flattened QQ plot)**: This plot checks the normality of the residuals. The points should tightly hug the horizontal line. The x-axis shows the theoretical quantiles of a standard normal distribution, and the y-axis shows the deviations of the sample residual quantiles from those theoretical values. Departures at the tails are common, but systematic deviations from the line (e.g., curvature) may indicate that the model does not conform to the normal residuals assumption.

Look carefully at the **Linearity** panel. There seems to be a slight curve in the reference line (rather than a flat horizontal line at zero), which suggests the relationship between age and agreeableness may not be purely linear. Let's investigate.

## Diagnosing and addressing curvilinearity

### Zooming into the residuals-vs-fitted plot

The `check_model()` dashboard packs a lot into a small space. Let's isolate the residuals-vs-fitted panel so we can examine the pattern more clearly:

```{r}
#| fig-width: 10
#| fig-height: 6

# Extract residuals and fitted values from the linear model
diag_simple <- data.frame(
  fitted = fitted(m_simple),
  residuals = residuals(m_simple)
)

ggplot(diag_simple, aes(x = fitted, y = residuals)) +
  geom_point(alpha = 0.12, size = 1.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_smooth(method = "loess", se = TRUE, color = "#e41a1c", linewidth = 1) +
  labs(
    title = "Residuals vs. fitted values: Linear model",
    subtitle = "Red LOESS curve should be flat if the linearity assumption holds",
    x = "Fitted values",
    y = "Residuals"
  ) +
  theme_minimal(base_size = 14)
```

If the red LOESS curve deviates from the dashed zero line — especially showing a U-shape or inverted-U — then a linear specification is missing curvature in the data.

### Centering age and fitting a polynomial model

The standard remedy for mild curvilinearity is to add a quadratic (squared) term. We **center** age before squaring it to reduce nonessential collinearity between the linear and quadratic terms. Centering means subtracting the mean, so that $\text{age\_c} = \text{age} - \overline{\text{age}}$. Without centering, `age` and `age^2` are highly correlated simply because squaring a variable that is far from zero inflates shared variance.

```{r}
# Center age
bfi_scored <- bfi_scored %>%
  mutate(
    age_c = age - mean(age, na.rm = TRUE),
    age_c2 = age_c^2
  )

# Fit the polynomial model
m_quad <- lm(agreeableness ~ age_c + age_c2, data = bfi_scored)
summary(m_quad)
```

The coefficient for `age_c2` tells us whether there is significant curvature. A positive coefficient means the relationship is U-shaped (agreeableness accelerates at older ages); a negative one means the curve flattens or turns over.

### Why centering matters

To see the effect of centering on collinearity, compare the correlation between the linear and quadratic terms with and without centering:

```{r}
cat("Correlation WITHOUT centering (age, age²):",
    round(cor(bfi_scored$age, bfi_scored$age^2), 3), "\n")
cat("Correlation WITH centering (age_c, age_c²):",
    round(cor(bfi_scored$age_c, bfi_scored$age_c2), 3), "\n")
```

Centering dramatically reduces the correlation, which makes the coefficient estimates for the linear and quadratic terms more stable and interpretable.

### Formal model comparison: F-test

Is the quadratic model significantly better than the linear one? The `anova()` function compares nested models using an incremental F-test:

```{r}
# Refit the linear model using centered age for a fair comparison
m_linear_c <- lm(agreeableness ~ age_c, data = bfi_scored)

anova(m_linear_c, m_quad)
```

This tests the null hypothesis that the quadratic term adds no explanatory power. A significant *p*-value (typically *p* < .05) indicates that the curved model fits the data better than the straight line.

### Visualizing the quadratic prediction

Let's overlay both the linear and quadratic predictions on the raw data to see the practical difference:

```{r}
# Generate predictions from both models
pred_grid <- data.frame(
  age = seq(min(bfi_scored$age), max(bfi_scored$age), length.out = 200)
) %>%
  mutate(
    age_c = age - mean(bfi_scored$age),
    age_c2 = age_c^2
  )

pred_grid$pred_linear <- predict(m_linear_c, newdata = pred_grid)
pred_grid$pred_quad <- predict(m_quad, newdata = pred_grid)

pred_long <- pred_grid %>%
  pivot_longer(
    cols = c(pred_linear, pred_quad),
    names_to = "model",
    values_to = "predicted"
  ) %>%
  mutate(model = ifelse(model == "pred_linear", "Linear", "Quadratic"))

ggplot() +
  geom_point(
    data = bfi_scored, aes(x = age, y = agreeableness),
    alpha = 0.12, size = 1.5
  ) +
  geom_line(
    data = pred_long, aes(x = age, y = predicted, color = model),
    linewidth = 1.3
  ) +
  scale_color_manual(
    values = c("Linear" = "#2c7fb8", "Quadratic" = "#e41a1c"),
    name = "Model"
  ) +
  labs(
    title = "Linear vs. quadratic fit: Agreeableness by age",
    subtitle = "Quadratic model uses centered age to reduce collinearity",
    x = "Age (years)",
    y = "Agreeableness (mean of 5 items)"
  ) +
  theme_minimal(base_size = 14)
```

### Comparing residual diagnostics

The real test of improvement is whether the residual plots look better. Here we compare the residuals-vs-fitted plots side by side:

```{r}
#| fig-width: 12
#| fig-height: 5

diag_quad <- data.frame(
  fitted = fitted(m_quad),
  residuals = residuals(m_quad)
)

p_resid_linear <- ggplot(diag_simple, aes(x = fitted, y = residuals)) +
  geom_point(alpha = 0.12, size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_smooth(method = "loess", se = TRUE, color = "#2c7fb8", linewidth = 1) +
  labs(
    title = "Linear model",
    x = "Fitted values",
    y = "Residuals"
  ) +
  theme_minimal(base_size = 13)

p_resid_quad <- ggplot(diag_quad, aes(x = fitted, y = residuals)) +
  geom_point(alpha = 0.12, size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_smooth(method = "loess", se = TRUE, color = "#e41a1c", linewidth = 1) +
  labs(
    title = "Quadratic model",
    x = "Fitted values",
    y = "Residuals"
  ) +
  theme_minimal(base_size = 13)

p_resid_linear + p_resid_quad +
  plot_annotation(
    title = "Residuals vs. fitted: Linear vs. quadratic model",
    subtitle = "A flatter LOESS curve indicates the model captures the trend better"
  )
```

The quadratic model seems to have addressed the curvilinearity: its LOESS curve should be closer to the dashed zero line than the linear model's curve. That said, notice how few observations there are of lower Agreeableness individuals. Likewise, the curvilinearity seems like it may rely on a few older people -- the dataset does not have many people above age 60. Thus, in the full diagnostic dashboard, we should ensure that there are no unduly influential observations (points outside the Cook's distance reference lines).

### Full diagnostic dashboard for the quadratic model

```{r}
#| fig-width: 12
#| fig-height: 10
check_model(m_quad)
```

### Model comparison summary

```{r}
tab_model(m_linear_c, m_quad,
  dv.labels = c("Linear (age only)", "Quadratic (age + age²)"),
  show.ci = 0.95, show.aic = TRUE
)
```

The take-away: when diagnostic plots suggest nonlinearity, adding a polynomial term is a simple and interpretable remedy. Centering the predictor before squaring ensures that the linear and quadratic terms are not artificially correlated, producing more stable estimates and cleaner interpretation.

# Multiple regression

Because the earlier diagnostics suggested that age is better modeled with a quadratic term, we will **carry forward** the centered linear and quadratic age terms and now add gender and education as predictors.

## Fit the model

```{r}
age_m <- mean(bfi_scored$age, na.rm = TRUE)

m_multi <- lm(
  agreeableness ~ I(age - age_m) + I((age - age_m)^2) +
    gender_f + education_f,
  data = bfi_scored
)
summary(m_multi)
```

## Visualizing predictions across predictors

With multiple predictors, we need to think about which variable(s) to display on the x-axis and which to condition on. `plot_predictions()` makes this straightforward.

### Predicted agreeableness by age, separated by gender

```{r}
plot_predictions(m_multi, condition = c("age", "gender_f")) +
  labs(
    title = "Predicted agreeableness by age and gender",
    subtitle = "Holding education at a typical value",
    y = "Predicted agreeableness",
    color = "Gender",
    fill = "Gender"
  ) +
  theme_minimal(base_size = 14)
```

> **Notice the matching curvature:** Because we fit an additive model with centered linear and quadratic age terms plus `gender_f`, the model assumes the age trajectory has the same shape for males and females. The curves can be shifted up or down, but they cannot bend differently by gender. If we wanted to allow the age pattern to differ by gender, we would need to add interaction terms between gender and the age terms. Visualizing predictions makes the model's structural assumptions instantly visible.

### Predicted agreeableness by education level

```{r}
plot_predictions(m_multi, condition = "education_f") +
  labs(
    title = "Predicted agreeableness by education level",
    subtitle = "Holding age and gender at typical values",
    y = "Predicted agreeableness",
    x = "Education"
  ) +
  theme_minimal(base_size = 14)
```

## Coefficient forest plot with `sjPlot`

A forest plot shows each coefficient as a point estimate with a confidence interval. This is a much more informative summary than a table of numbers, because it makes the relative magnitude and precision of effects immediately visible.

```{r}
plot_model(
  m_multi,
  show.values = TRUE,
  value.offset = 0.3,
  axis.labels = c(
    "Centered age",
    "Centered age^2",
    "Female (vs Male)",
    "HS grad (vs Some HS)",
    "Some college (vs Some HS)",
    "College grad (vs Some HS)",
    "Graduate (vs Some HS)"
  )
) +
  labs(title = "Coefficient plot: Agreeableness model") +
  theme_minimal(base_size = 14)
```

The coefficients are on the unstandardized scale by default. To see standardized coefficients (useful for comparing the relative importance of predictors), use `type = "std"`:

```{r}
plot_model(
  m_multi,
  type = "std",
  show.values = TRUE,
  value.offset = 0.3,
  axis.labels = c(
    "Centered age",
    "Centered age^2",
    "Female (vs Male)",
    "HS grad (vs Some HS)",
    "Some college (vs Some HS)",
    "College grad (vs Some HS)",
    "Graduate (vs Some HS)"
  )
) +
  labs(title = "Standardized coefficients: Agreeableness model") +
  theme_minimal(base_size = 14)
```

## Publication-quality table with `sjPlot`

```{r}
tab_model(
  m_multi,
  show.intercept = FALSE,
  show.std = TRUE,
  show.ci = 0.95,
  pred.labels = c(
    "Centered age",
    "Centered age^2",
    "Female (vs Male)",
    "HS grad (vs Some HS)",
    "Some college (vs Some HS)",
    "College grad (vs Some HS)",
    "Graduate (vs Some HS)"
  )
)
```

## Diagnostics

```{r}
#| fig-width: 12
#| fig-height: 10
check_model(m_multi)
```

# Mean comparisons: t-tests and ANOVA

## Framing: ANOVA as `lm()`

A one-way ANOVA is just an `lm()` with a categorical predictor. By fitting the model with `lm()` rather than `aov()`, we stay within the same framework and can use all the same visualization tools.

## Data: The anorexia dataset

The `anorexia` dataset from the `MASS` package records pre- and post-treatment weights (in pounds) for 72 young women receiving one of three treatments for anorexia:

- **CBT**: Cognitive-behavioral therapy
- **Cont**: Control (no treatment)
- **FT**: Family therapy

```{r}
data(anorexia, package = "MASS")

anorexia <- anorexia %>%
  mutate(
    weight_change = Postwt - Prewt,
    Treat = factor(Treat, levels = c("Cont", "CBT", "FT"))
  )

# Quick look at the data
anorexia %>%
  group_by(Treat) %>%
  summarize(
    n = n(),
    mean_change = mean(weight_change),
    sd_change = sd(weight_change),
    .groups = "drop"
  )
```

## Visualizing the raw data first

Before fitting any model, always look at the data. Here we layer individual observations, a boxplot, and group means.

```{r}
ggplot(anorexia, aes(x = Treat, y = weight_change, fill = Treat)) +
  geom_boxplot(alpha = 0.4, outlier.shape = NA, width = 0.5) +
  geom_jitter(width = 0.15, alpha = 0.6, size = 2) +
  stat_summary(
    fun = mean, geom = "point",
    shape = 18, size = 5, color = "black"
  ) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "Weight change by treatment group",
    subtitle = "Diamonds show group means; dashed line = no change",
    x = "Treatment",
    y = "Weight change (lbs)"
  ) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none")
```

## Fit the ANOVA as `lm()`

```{r}
m_anova <- lm(weight_change ~ Treat, data = anorexia)
summary(m_anova)
```

The intercept is now the mean of the reference group (Cont). The coefficients for CBT and FT represent *differences from Control*.

## Estimated marginal means with `marginaleffects`

`predictions()` gives us the estimated marginal means — the model's expected value for each group — complete with confidence intervals.

To do this, we use the `newdata = datagrid(Treat = unique)` argument. The `datagrid()` function is a powerful tool: it creates a synthetic, balanced dataset containing exactly one row for each unique level of `Treat`. If we had other predictors in the model (like age or baseline weight), `datagrid()` would automatically hold them at their mean or typical values, ensuring we make a clean, "apples-to-apples" comparison between treatment groups.

```{r}
# Estimated marginal means
emm <- predictions(m_anova, newdata = datagrid(Treat = unique))
emm
```

```{r}
plot_predictions(m_anova, condition = "Treat") +
  labs(
    title = "Estimated marginal means by treatment group",
    subtitle = "From lm(weight_change ~ Treat)",
    x = "Treatment",
    y = "Predicted weight change (lbs)"
  ) +
  theme_minimal(base_size = 14)
```

## Pairwise contrasts

Which groups differ from each other? `comparisons()` computes all pairwise differences between group means.

```{r}
# Pairwise contrasts between treatment groups
avg_comparisons(m_anova, variables = "Treat")
```

## Diagnostics for the ANOVA model

A one-way ANOVA fit as `lm(weight_change ~ Treat)` has only a categorical predictor, so the fitted values take on just a few group means. In this situation, the residual-vs-fitted and scale-location smoothers inside the full `check_model()` dashboard can behave erratically and are not very informative. We therefore keep the most useful panels and then inspect the key assumptions with dedicated plots below.

```{r}
#| fig-width: 12
#| fig-height: 6
check_model(m_anova, check = c("qq", "pp_check", "outliers"))
check_model(m_anova)
```

### Normality of residuals

```{r}
#| fig-width: 7
#| fig-height: 5
plot(check_normality(m_anova))
```

### Testing homogeneity of variance

A key assumption of the standard ANOVA is equal variances across groups. The `performance` package offers both a dedicated test and a groupwise visualization:

```{r}
check_homogeneity(m_anova)
```

```{r}
#| fig-width: 8
#| fig-height: 6
plot(check_homogeneity(m_anova))
```

### Influential observations

```{r}
check_outliers(m_anova)
```

# Interaction effects

## Does the relationship between pre-treatment weight and weight change differ by treatment group?

This is a classic ANCOVA design: a continuous covariate (Prewt) interacting with a categorical factor (Treat).

```{r}
m_interact <- lm(weight_change ~ Prewt * Treat, data = anorexia)
summary(m_interact)
```

## Visualizing the interaction

`plot_predictions()` can display the predicted weight change as a function of pre-treatment weight, separately for each treatment group:

```{r}
plot_predictions(m_interact, condition = c("Prewt", "Treat")) +
  labs(
    title = "Interaction: Weight change by pre-treatment weight and treatment",
    subtitle = "Separate regression lines for each treatment group",
    x = "Pre-treatment weight (lbs)",
    y = "Predicted weight change (lbs)",
    color = "Treatment"
  ) +
  theme_minimal(base_size = 14)
```

## Marginal effects (slopes) by group

How does the effect of pre-treatment weight vary by treatment group? `plot_slopes()` answers this directly:

```{r}
plot_slopes(m_interact, variables = "Prewt", condition = "Treat") +
  labs(
    title = "Slope of pre-treatment weight by treatment group",
    subtitle = "How does the Prewt → weight change relationship differ?",
    x = "Treatment",
    y = "Marginal effect of Prewt on weight change"
  ) +
  theme_minimal(base_size = 14)
```

## Added-variable plots (partial regression plots)

Fox (2018) describes added-variable (AV) plots as a way to visualize the unique contribution of a predictor after controlling for other variables. The `car` package provides `avPlots()`:

```{r}
#| fig-width: 10
#| fig-height: 8
avPlots(m_interact, id = FALSE)
```

The AV plot for a predictor $x_j$ is constructed by:

1. Regressing the outcome $y$ on all predictors *except* $x_j$ and taking the residuals.
2. Regressing $x_j$ on all other predictors and taking the residuals.
3. Plotting those two sets of residuals against each other.

The slope of the line in the AV plot is exactly the coefficient of $x_j$ in the full model. This makes AV plots useful for judging whether a term is making an **incremental contribution** beyond the rest of the model.

How to read them:

1. If the fitted line has a clear positive or negative slope, that term is contributing unique information after controlling for the other predictors.
2. If the fitted line is nearly flat and the points look like an undirected cloud, that suggests little incremental contribution from that term.
3. Points that sit far from the rest of the cloud can signal observations that are disproportionately affecting the apparent contribution of that term.

For this interaction model, note that the treatment factor is represented by two dummy-coded contrasts (`TreatCont` and `TreatFT`), and the interaction is represented by two product terms. So these AV plots are showing the incremental contribution of each **coefficient-coded term**, not a single omnibus test of the entire factor or the entire interaction.

## Coefficient plot for the interaction model

For coefficient displays, the centered interaction model is easier to interpret because the treatment coefficients are defined at the **average pre-treatment weight** rather than at `Prewt = 0`.

```{r}
anorexia_c <- anorexia %>%
  mutate(Prewt_c = Prewt - mean(Prewt, na.rm = TRUE))

m_interact_c <- lm(weight_change ~ Prewt_c * Treat, data = anorexia_c)

plot_model(m_interact_c, show.values = TRUE, value.offset = 0.3) +
  labs(title = "Coefficient plot: Centered interaction model") +
  theme_minimal(base_size = 14)
```

## Diagnostics for the interaction model

```{r}
#| fig-width: 12
#| fig-height: 10
check_model(m_interact)
```

## Follow-up: Why are the Treat and Treat x Prewt VIFs so high?

One thing worth noticing in the `check_model()` output is that the treatment main effect and the `Treat x Prewt` interaction can show extremely high VIFs. At first glance, this looks alarming. But in interaction models, a high VIF is often caused by **nonessential collinearity**: overlap created by the coding of the predictors, not by a substantive flaw in the data.

Let's inspect the collinearity directly:

```{r}
check_collinearity(m_interact)
```

The problem is structural. `Prewt` is coded in its raw units, so the main effect of `Treat` is defined at `Prewt = 0`, which is far outside the observed range of pre-treatment weights. That coding choice makes the treatment dummy columns and their interaction columns overlap much more than necessary.

This is the same issue we saw earlier with quadratic terms and the same issue that often appears in multilevel interaction models: when a continuous predictor is left uncentered, the parameterization can induce avoidable collinearity. The remedy is to **center** the continuous predictor at a meaningful value, often its sample mean.

```{r}
anorexia_c <- anorexia %>%
  mutate(Prewt_c = Prewt - mean(Prewt, na.rm = TRUE))

mean(anorexia_c$Prewt_c)
```

Now refit the same interaction model using centered pre-treatment weight:

```{r}
m_interact_c <- lm(weight_change ~ Prewt_c * Treat, data = anorexia_c)

check_collinearity(m_interact_c)
```

Centering dramatically reduces the VIFs for the treatment main effect and the interaction terms. Importantly, this does **not** change the fitted values or the substantive interaction. It changes the parameterization so that the coefficients are easier to interpret and the model matrix is less redundant.

```{r}
rbind(
  uncentered = coef(m_interact),
  centered = coef(m_interact_c)
)
```

Notice what stays the same and what changes:

1. The slope for pre-treatment weight and the interaction terms are unchanged apart from relabeling (`Prewt` becomes `Prewt_c`).
2. The coefficients for `Treat` change because they now refer to group differences at the **average pre-treatment weight**, rather than at `Prewt = 0`, which is outside the observed data.
3. The intercept also changes for the same reason: it now refers to the predicted weight change for a person at the sample mean pre-treatment weight in the reference treatment group.

We can verify that this is only a reparameterization by comparing fitted values:

```{r}
cat("Maximum absolute difference in fitted values:",
    max(abs(fitted(m_interact) - fitted(m_interact_c))), "\n")
```

For completeness, here is the diagnostic dashboard for the centered interaction model:

```{r}
#| fig-width: 12
#| fig-height: 10
check_model(m_interact_c)
```

The broader lesson is:

1. When a model includes a continuous predictor and its interaction with a factor, check whether a high VIF reflects **nonessential collinearity**.
2. Centering the continuous predictor often reduces this problem substantially.
3. Centering does not change the fitted lines or the substantive interaction; it changes the reference point for the main effects.
4. Main effects in interaction models are easiest to interpret when the continuous predictor is centered at a meaningful value.

# Take-home summary

## The three-step workflow

1. **Visualize predictions**: Use `marginaleffects::plot_predictions()` to see what the model expects at different values of the predictors. This is the most direct way to understand what a model *says*.
2. **Summarize coefficients**: Use `sjPlot::plot_model()` for forest plots and `sjPlot::tab_model()` for publication tables. Coefficients answer "how much does the outcome change per unit of the predictor?"
3. **Diagnose misfit**: Use `performance::check_model()` to see whether the model's assumptions are met. If the diagnostic plots show problems, the predictions and coefficients may not be trustworthy.

## Tool reference

| Task | Function | Package |
|---|---|---|
| Predicted values with CIs | `predictions()`, `plot_predictions()` | marginaleffects |
| Group contrasts (pairwise comparisons) | `comparisons()`, `avg_comparisons()` | marginaleffects |
| Marginal effects / slopes | `slopes()`, `plot_slopes()` | marginaleffects |
| Coefficient forest plot | `plot_model()` | sjPlot |
| Model summary table | `tab_model()` | sjPlot |
| Multi-panel diagnostics | `check_model()` | performance |
| Collinearity diagnostics | `check_collinearity()` | performance |
| Homogeneity of variance test | `check_homogeneity()` | performance |
| Added-variable plots | `avPlots()` | car |

## Further reading

- Arel-Bundock, V. *Model to Meaning*. <https://marginaleffects.com>
- Fox, J. (2018). Visualizing fit and misfit in regression models. *Handbook of Data Visualization*.
- Lüdecke, D. et al. (2021). performance: An R Package for Assessment, Comparison and Testing of Statistical Models. *JOSS*, 6(60), 3139.
