可视化含两个连续预测变量且交互作用显著的线性混合效应模型

Great to hear you’ve uncovered a significant interaction in your linear mixed-effects model! Visualizing interactions between two continuous predictors is key to understanding exactly how they combine to influence your outcome variable. Let’s walk through a few straightforward approaches using the packages you already have loaded, plus some extra tips to make your plots clear and informative.

1. Visualize the Overall Interaction with the effects Package #

Since you’re already using effects, this is the easiest starting point—it extracts the marginal effects of your interaction term and handles the heavy lifting of calculating predicted values and standard errors.

First, extract the interaction effect and convert it to a data frame for ggplot2:

# Extract the interaction effect between the two continuous predictors
int_effect <- Effect(c("AIP_s_parent.z", "Q_mean.z"), mod_father_son)

# Convert to a data frame compatible with ggplot2
int_effect_df <- as.data.frame(int_effect)

Now, plot the relationship between parent AIP and child AIP at three meaningful levels of Q_mean (since your variables are z-scored, we’ll use -1 SD, mean, and +1 SD to represent low, average, and high Q_mean):

ggplot(int_effect_df, aes(x = AIP_s_parent.z, y = fit)) +
  # Add lines for each Q_mean level
  geom_line(aes(color = factor(Q_mean.z)), linewidth = 1) +
  # Add shaded confidence bands (±1 SE)
  geom_ribbon(aes(ymin = fit - se, ymax = fit + se, fill = factor(Q_mean.z)), 
              alpha = 0.2, color = NA) +
  # Customize labels and legend
  scale_color_discrete(
    name = "Q_mean (z-score)",
    labels = c("Low (-1 SD)", "Mean (0)", "High (+1 SD)")
  ) +
  scale_fill_discrete(
    name = "Q_mean (z-score)",
    labels = c("Low (-1 SD)", "Mean (0)", "High (+1 SD)")
  ) +
  labs(
    x = "Parent AIP (z-score)",
    y = "Child AIP (z-score)",
    title = "Interaction Between Parent AIP and Q_mean on Child AIP"
  ) +
  # Use clean theme and custom font (since you loaded extrafont)
  theme_minimal() +
  theme(text = element_text(family = "Arial"))

This plot will show you how the slope of parent AIP’s effect on child AIP changes depending on Q_mean—for example, maybe the effect is stronger when Q_mean is high, or reverses direction at low Q_mean.

2. Explore Simple Slopes with emmeans #

To quantify and visualize how the slope of parent AIP changes across Q_mean levels, you can use the emmeans package (a common companion for LME models):

library(emmeans)

# Calculate the slope of parent AIP at key Q_mean levels (-1, 0, +1 SD)
simple_slopes <- emtrends(
  mod_father_son, 
  ~ Q_mean.z, 
  var = "AIP_s_parent.z", 
  at = list(Q_mean.z = c(-1, 0, 1))
)

# Convert to data frame for plotting
slopes_df <- as.data.frame(simple_slopes)

# Plot the slopes with confidence intervals
ggplot(slopes_df, aes(x = Q_mean.z, y = AIP_s_parent.z.trend)) +
  geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL), size = 1) +
  # Add a dashed line at 0 to highlight non-significant slopes
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    x = "Q_mean (z-score)",
    y = "Slope of Parent AIP on Child AIP",
    title = "Simple Slopes of Parent AIP at Different Q_mean Levels"
  ) +
  theme_minimal()

This plot tells you whether the effect of parent AIP is significant at each Q_mean level (if the confidence interval doesn’t cross the red dashed line, the slope is statistically significant).

3. Optional: Visualize Family-Specific Trends #

Your model includes family-level random effects for the interaction, so you might want to see how the interaction varies across families. Note that this works best if you don’t have hundreds of families (otherwise the plot will be too cluttered):

# Generate predicted values for each observation, including random effects
data_father_son$predicted_child_AIP <- predict(mod_father_son)

# Plot family-specific trends overlaid with the overall average trend
ggplot(data_father_son, aes(x = AIP_s_parent.z, y = predicted_child_AIP)) +
  # Add faint lines for each family
  geom_line(aes(color = factor(Family_number)), alpha = 0.3) +
  # Overlay the overall average trend from the effects package
  geom_line(data = int_effect_df, aes(y = fit), color = "black", linewidth = 1.5) +
  labs(
    x = "Parent AIP (z-score)",
    y = "Predicted Child AIP (z-score)",
    title = "Family-Specific vs. Overall Interaction Patterns"
  ) +
  theme_minimal() +
  theme(legend.position = "none") # Hide family legend to avoid clutter

This helps you see how much individual families deviate from the overall interaction effect.

内容的提问来源于stack exchange，提问作者Julia M