如何在ggplot2中正确拟合Random Intercept-Fixed Slope模型？

我用以下代码生成模拟数据：

# Create an empty data frame to store the simulated data
data <- data.frame(
  Lot = rep(1:num_lots, each = 9),
  Time = rep(3 * 0:8, times = num_lots),
  Purity = numeric(num_lots * 9)
)

# Simulate purity data for each lot and time point
for (lot in 1:num_lots) {
  # Generate random intercept and slope for each lot
  intercept <- rnorm(1, mean = 95, sd = 2)
  slope <- runif(1, min = -.7, max = 0)
  
  for (month in 0:8) {
    # Simulate purity data with noise
    data[data$Lot == lot &
           data$Time == month * 3, "Purity"] <-
      intercept + slope * month * 3 + rnorm(1, mean = 0, sd = .35)
  }
}

我想拟合随机截距-固定斜率的混合效应模型，用了下面的ggplot2代码却没得到预期结果：

ggplot(data, aes(x = Time, y = Purity)) +
  geom_point(aes(color = as.factor(Lot), shape = as.factor(Lot))) +
  geom_smooth(
    method = "lm",
    formula = y ~ x  + (1 | Lot), data = data,
    se = FALSE,
    size = 0.5,
    show.legend = FALSE
  ) +
  labs(
    title = "Random Intercept and Fixed Slope ",
    x = "month",
    y = "Purity",
    color = "Lot",
    shape = "Lot"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = c(0, 3, 6, 9, 12, 15, 18, 21, 24))

我期望每个Lot对应一条斜率相同、截距不同的回归线，请问代码哪里出错了？

你的核心问题是geom_smooth的method = "lm"不支持混合效应模型的公式语法(1 | Lot)——这个随机效应项会被直接忽略，最终只会画出一条全局的普通线性回归线，而非每个Lot单独的回归线。

下面提供两种符合需求的解决方案：

方案1：基于混合效应模型预测绘图（推荐） #

先用lme4包拟合标准的随机截距-固定斜率模型，再生成每个Lot的预测值，最后用ggplot绘制回归线：

library(lme4)
library(ggplot2)

# 拟合混合效应模型
model <- lmer(Purity ~ Time + (1 | Lot), data = data)

# 构建包含所有Time和Lot组合的预测数据集
pred_data <- expand.grid(
  Time = unique(data$Time),
  Lot = unique(data$Lot)
)
# 生成预测值
pred_data$Purity_pred <- predict(model, newdata = pred_data)

# 绘图
ggplot(data, aes(x = Time, y = Purity)) +
  geom_point(aes(color = as.factor(Lot), shape = as.factor(Lot))) +
  geom_line(data = pred_data, aes(x = Time, y = Purity_pred, color = as.factor(Lot)), size = 0.5) +
  labs(
    title = "Random Intercept and Fixed Slope",
    x = "Month",
    y = "Purity",
    color = "Lot",
    shape = "Lot"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = c(0, 3, 6, 9, 12, 15, 18, 21, 24))

方案2：强制全局斜率的分组拟合 #

如果不想依赖lme4，可以先计算全局线性模型的斜率，再让每个Lot仅拟合截距，从而实现固定斜率、不同截距的效果：

library(ggplot2)

# 计算全局斜率
global_lm <- lm(Purity ~ Time, data = data)
global_slope <- coef(global_lm)[["Time"]]

ggplot(data, aes(x = Time, y = Purity)) +
  geom_point(aes(color = as.factor(Lot), shape = as.factor(Lot))) +
  geom_smooth(
    aes(group = Lot, color = as.factor(Lot)),
    method = "lm",
    formula = y ~ I(x * global_slope) + 1, # 固定斜率，仅拟合截距
    se = FALSE,
    size = 0.5
  ) +
  labs(
    title = "Random Intercept and Fixed Slope",
    x = "Month",
    y = "Purity",
    color = "Lot",
    shape = "Lot"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = c(0, 3, 6, 9, 12, 15, 18, 21, 24))

说明 #

• 方案1严格遵循混合效应模型的统计逻辑，生成的回归线是模型的预测结果，更严谨；
• 方案2是一种近似实现，适合快速验证可视化效果，结果和混合模型的预测线会略有差异。

内容的提问来源于stack exchange，提问作者Joe the Second