在R中绘制多层模型的二次效应：倒U型关系可视化方法咨询

Hey Gina, nice work getting that mixed-effects model set up with the quadratic term! Visualizing that inverted U-shape is key to verifying your relationship, and I’ll walk you through two robust, easy-to-implement methods using R’s most popular plotting tools.

方法1：使用predict() + ggplot2（最灵活） #

This approach lets you generate population-level (or individual-level) predictions from your model, then plot them alongside your raw data to highlight the quadratic trend.

First, load the necessary packages and create a dataset for predictions:

# Load required packages
library(lme4)
library(ggplot2)
library(dplyr)

# Create a dataset covering the full range of your MTsqZ variable
# We'll focus on population-level fit first (ignoring random intercepts for the trend)
pred_data <- tibble(
  MTsqZ = seq(min(MT$MTsqZ, na.rm = TRUE), max(MT$MTsqZ, na.rm = TRUE), length.out = 100),
  VP04_01 = factor(rep(unique(MT$VP04_01)[1], 100)) # Dummy value for random effect
)

# Generate predictions with 95% confidence intervals
# Use re.form = ~0 to get population-level (fixed effects only) predictions
preds <- predict(Pfad, newdata = pred_data, re.form = ~0, interval = "confidence")

# Merge predictions with our dataset
pred_data <- pred_data %>% bind_cols(as_tibble(preds))

Now plot the raw data and fitted curve:

ggplot(MT, aes(x = MTsqZ, y = FlowZ)) +
  # Raw data points (transparent to avoid clutter)
  geom_point(alpha = 0.3, color = "gray50") +
  # Fitted quadratic curve
  geom_line(data = pred_data, aes(y = fit), color = "#2c3e50", linewidth = 1) +
  # 95% confidence interval ribbon
  geom_ribbon(data = pred_data, aes(ymin = lwr, ymax = upr), fill = "#2c3e50", alpha = 0.2) +
  # Clean labels and theme
  labs(
    x = "Standardized Quadratic MT (MTsqZ)",
    y = "Standardized Flow (FlowZ)",
    title = "Population-Level Inverted U-Shape Relationship",
    subtitle = "Adjusted for random participant intercepts"
  ) +
  theme_minimal()

方法2：使用emmeans包（简化边际均值计算） #

The emmeans package is designed for extracting marginal means from mixed models, making it even easier to plot your quadratic trend without manually creating prediction datasets.

# Load the package
library(emmeans)

# Extract marginal means for MTsqZ across its full range
emm_trend <- emmeans(Pfad, ~ MTsqZ, at = list(MTsqZ = seq(min(MT$MTsqZ), max(MT$MTsqZ), length.out = 100)))

# Convert to a data frame for plotting
emm_df <- as.data.frame(emm_trend)

# Plot with ggplot2
ggplot(MT, aes(x = MTsqZ, y = FlowZ)) +
  geom_point(alpha = 0.3, color = "gray50") +
  geom_line(data = emm_df, aes(y = emmean), color = "#e74c3c", linewidth = 1) +
  geom_ribbon(data = emm_df, aes(ymin = lower.CL, ymax = upper.CL), fill = "#e74c3c", alpha = 0.2) +
  labs(
    x = "Standardized Quadratic MT (MTsqZ)",
    y = "Standardized Flow (FlowZ)",
    title = "Inverted U-Shape Relationship (Marginal Means)",
    subtitle = "With 95% confidence intervals"
  ) +
  theme_minimal()

可选：展示个体水平的拟合曲线 #

If you want to see how the quadratic trend varies across individual participants (VP04_01), adjust the prediction dataset to include all participant levels:

# Create dataset with all participants and MTsqZ range
ind_pred_data <- expand.grid(
  MTsqZ = seq(min(MT$MTsqZ), max(MT$MTsqZ), length.out = 100),
  VP04_01 = unique(MT$VP04_01)
)

# Generate individual-level predictions
ind_preds <- predict(Pfad, newdata = ind_pred_data, interval = "confidence")
ind_pred_data <- cbind(ind_pred_data, ind_preds)

# Plot individual curves
ggplot(MT, aes(x = MTsqZ, y = FlowZ)) +
  geom_point(alpha = 0.2, color = "gray50") +
  geom_line(data = ind_pred_data, aes(y = fit, color = VP04_01), linewidth = 0.8, alpha = 0.7) +
  labs(
    x = "Standardized Quadratic MT (MTsqZ)",
    y = "Standardized Flow (FlowZ)",
    title = "Individual-Level Inverted U-Shape Trends",
    color = "Participant ID"
  ) +
  theme_minimal()

关键注意点 #

• Always use model-generated predictions (not raw data regression lines) for mixed models—raw lines won’t account for the random intercept structure of your model.
• The re.form = ~0 argument in predict() ensures you’re plotting the population-level trend (average across all participants). Omit this if you want to include individual random effects.

内容的提问来源于stack exchange，提问作者Gina W.