R中含交互项的线性回归：拟合数据与原始数据绘图问题

1. 为什么geom_smooth()未采用自定义模型？ #

geom_smooth()的默认逻辑是：根据当前aes()中指定的x、y变量，自动拟合一个新模型（小样本默认用loess）。你当前代码里仅指定了y = .fitted，但没有明确告诉ggplot要复用预先拟合的lm模型，因此它会基于x=Wind和y=.fitted重新拟合loess模型，而非使用你包含交互项的线性回归模型。

2. 基于自定义模型绘制拟合线与原始数据的方法 #

以下两种方法可实现需求：

方法一：手动生成预测数据（灵活可控） #

通过predict()生成模型预测值，可灵活控制其他变量的取值（如均值、分位数），再结合原始数据绘图：

library(tidyverse)
library(broom)

# 拟合含交互项的线性模型
model <- lm(Ozone ~ Solar.R + Wind + Solar.R:Temp, data = airquality)

# 生成预测数据集：固定Solar.R和Temp为均值，Wind取序列值
pred_data <- expand.grid(
  Wind = seq(min(airquality$Wind, na.rm = TRUE), max(airquality$Wind, na.rm = TRUE), length.out = 100),
  Solar.R = mean(airquality$Solar.R, na.rm = TRUE),
  Temp = mean(airquality$Temp, na.rm = TRUE)
) %>%
  mutate(.fitted = predict(model, newdata = .))

# 绘制原始数据点 + 模型拟合线
ggplot() +
  geom_point(data = airquality, aes(x = Wind, y = Ozone, size = Solar.R), na.rm = TRUE) +
  geom_line(data = pred_data, aes(x = Wind, y = .fitted), size = 1, color = "purple") +
  labs(x = "Wind", y = "Ozone", title = "Ozone vs. Wind (基于自定义线性模型)")

若需展示不同Solar.R水平下的拟合线，可调整预测数据的生成逻辑：

# 生成多水平Solar.R的预测数据
pred_data_multi <- expand.grid(
  Wind = seq(min(airquality$Wind, na.rm = TRUE), max(airquality$Wind, na.rm = TRUE), length.out = 100),
  Solar.R = quantile(airquality$Solar.R, na.rm = TRUE, c(0.25, 0.5, 0.75)),
  Temp = mean(airquality$Temp, na.rm = TRUE)
) %>%
  mutate(.fitted = predict(model, newdata = .),
         Solar.R_level = factor(Solar.R, labels = c("低Solar.R", "中Solar.R", "高Solar.R")))

# 绘图
ggplot() +
  geom_point(data = airquality, aes(x = Wind, y = Ozone, size = Solar.R), na.rm = TRUE) +
  geom_line(data = pred_data_multi, aes(x = Wind, y = .fitted, color = Solar.R_level), size = 1) +
  labs(x = "Wind", y = "Ozone", title = "不同Solar.R水平下Ozone与Wind的拟合关系", color = "Solar.R水平") +
  scale_color_viridis_d()

方法二：使用ggpredict包简化操作 #

ggpredict可直接提取模型预测值，快速生成可视化所需数据：

library(tidyverse)
library(broom)
library(ggpredict)

model <- lm(Ozone ~ Solar.R + Wind + Solar.R:Temp, data = airquality)

# 提取模型预测：控制Temp为均值，按Solar.R四分位数分组
pred <- ggpredict(model, terms = c("Wind", "Solar.R [quartile]"))

# 绘图
ggplot(pred, aes(x = x, y = predicted)) +
  geom_point(data = airquality, aes(x = Wind, y = Ozone, size = Solar.R), na.rm = TRUE) +
  geom_line(aes(color = group), size = 1) +
  labs(x = "Wind", y = "Ozone", title = "Ozone vs. Wind (基于自定义模型)", color = "Solar.R水平")

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