R中模型与原始数据同图可视化、验证及二阶响应面拟合咨询

DevOps 小满

2026-8-12

析因实验模型拟合相关问题解答

背景说明

我是R语言新手，仅掌握基础统计学知识，正在学习析因实验及模型拟合，参考《Design and Analysis of Experiments》（Montgomery 2013，ISBN：9781118097939）中的例5.3完成了如下操作：

构建数据集：

bottel_vul_data <- data.frame(A = rep(c(10, 12, 14), each = 2),
                              B = rep(c(25, 30), each = 12),
                              C = rep(c(200, 250), each = 6),        
                              vul = c(-3, -1, 0, 1, 5, 4,
                                      -1,  0, 2, 1, 7, 6,   
                                      -1,  0, 2, 3, 7, 9,
                                       1,  1, 6, 5, 10, 11))

执行ANOVA分析：

bottel_anova <- aov(vul ~ factor(A) * B * C, data = bottel_vul_data)
summary(bottel_anova)

结果与教材一致，随后对比了4个线性回归模型，通过模型间ANOVA确定**Model4（包含A、B、C及A:B交互项）**为合适模型。现针对提出的技术问题逐一解答：

1. 如何在R中将拟合模型与原始数据绘制在同一张图上进行可视化？

可以用ggplot2或基础绘图系统实现，先拟合Model4并生成预测值，再将原始数据与预测结果合并展示：

步骤1：拟合模型并生成预测值

model4 <- lm(vul ~ A + B + C + A:B, data = bottel_vul_data)
bottel_vul_data$pred_vul <- predict(model4)

步骤2：用ggplot2绘图（推荐）

以A为横轴，按B、C的水平区分数据，同时展示原始点与拟合线：

library(ggplot2)

ggplot(bottel_vul_data, aes(x = A, y = vul, color = factor(B))) +
  geom_point(size = 3, alpha = 0.7) +  # 原始数据点
  geom_line(aes(y = pred_vul), linewidth = 1) +  # 模型拟合线
  facet_wrap(~factor(C)) +  # 按C的水平分面展示
  labs(title = "原始数据与Model4拟合结果对比",
       x = "因子A", y = "响应vul", color = "因子B") +
  theme_minimal()

步骤3：用基础绘图系统实现

par(mfrow = c(1,2))  # 分两列展示C的两个水平
for (c_val in unique(bottel_vul_data$C)) {
  sub_data <- subset(bottel_vul_data, C == c_val)
  plot(sub_data$A, sub_data$vul, col = sub_data$B, pch = 16,
       main = paste("C =", c_val), xlab = "A", ylab = "vul")
  # 按B的水平绘制拟合线
  for (b_val in unique(sub_data$B)) {
    fit_line <- subset(sub_data, B == b_val)
    lines(fit_line$A, fit_line$pred_vul, col = b_val, lwd = 2)
  }
  legend("topleft", legend = unique(sub_data$B), col = unique(sub_data$B), lwd = 2, pch = 16)
}
par(mfrow = c(1,1))

2. 如何验证该模型是否能准确代表数据？

从残差分析、拟合优度、显著性检验、预测精度四个维度验证：

(1) 残差分析

通过残差诊断图判断模型假设（线性、方差齐性、正态性）是否成立：

bottel_vul_data$residuals <- resid(model4)

# 生成标准残差诊断图
par(mfrow = c(2,2))
plot(model4)
par(mfrow = c(1,1))

残差vs拟合值：点随机分布在0线附近说明方差齐性，无遗漏非线性关系；
正态Q-Q图：点接近直线说明残差符合正态分布；
尺度-位置图：点分布均匀说明方差稳定；
残差vs杠杆值：识别Cook距离大的异常点。

(2) 拟合优度

查看$R^{2$和调整$R}2$，值越接近1拟合效果越好：

summary(model4)$r.squared
summary(model4)$adj.r.squared

(3) 显著性检验

验证整体模型及各变量的显著性：

anova(model4)  # 整体模型及变量的F检验
summary(model4)  # 各系数的t检验

(4) 预测精度验证

用留一交叉验证评估模型预测能力：

library(boot)
cv.err <- cv.glm(bottel_vul_data, model4)$delta[1]
cat("留一交叉验证均方误差：", cv.err, "\n")

均方误差越小，模型预测精度越高。

3. 如何在R中拟合二阶模型以生成响应面，是否有类似lm()的可用函数？

可以直接用lm()拟合二阶模型，也可以用专门的rsm包简化操作：

(1) 用lm()拟合二阶模型

二阶模型包含线性项、平方项和所有二阶交互项，公式中用I()表示数学运算：

rs_model <- lm(vul ~ A + B + C + I(A^2) + I(B^2) + I(C^2) + A:B + A:C + B:C, data = bottel_vul_data)
summary(rs_model)

(2) 用rsm包拟合并生成响应面

rsm包专为响应面分析设计，语法更简洁，绘图更便捷：

install.packages("rsm")
library(rsm)

# 拟合二阶模型，SO()表示包含所有二阶项
rs_model_rsm <- rsm(vul ~ SO(A, B, C), data = bottel_vul_data)
summary(rs_model_rsm)

# 绘制等高线图和3D响应面图（固定C为均值）
contour(rs_model_rsm, ~ A + B, at = list(C = mean(bottel_vul_data$C)))
persp(rs_model_rsm, ~ A + B, at = list(C = mean(bottel_vul_data$C)), theta = 30, phi = 30)

内容的提问来源于stack exchange，提问作者Philip de Bruin