双/三因素ANOVA模型绘图、高阶交互拟合及不平衡数据处理问询

我正在用R复现Montgomery《Design and Analysis of Experiments》（2013年第8版）中的例5.5（双因素案例），以此学习三因素析因实验的模型拟合方法。目前已成功生成ANOVA表并提取模型系数（推测为β₀、β₁、β₂、β₁₂），但遇到以下问题：

• 不知如何绘制该模型对应的曲面；
• 能否在R中加入β₁₁₂、β₁₂₂、β₁₁₂₂等高阶交互项以获得更优拟合效果；
• 注意到aov()提示“Estimated effects may be unbalanced”，已知aov()更适用于水平和重复数相同的因素，想了解针对此类不平衡数据是否有更合适的ANOVA函数。

我熟悉Python但未用matplotlib绘制过曲面，刚接触R且无绘图经验，未找到相关解决资料。

我的R代码及输出 #

代码 #

lewensduur_data <- data.frame(A = rep(c(15, 20, 25), each = 2),
                         B = rep(c(125, 150, 175), each = 6),
                         lewe = c(-2, -1, 0, 2, -1, 0,
                                  -3, 0, 1, 3, 5, 6,
                                  2, 3, 4, 6, 0, -1))

lewensduur_anova <- aov(lewe ~ A * B, data = lewensduur_data)
lewensduur_anova

ANOVA表输出 #

Call:
   aov(formula = lewe ~ A * B, data = lewensduur_data)

Terms:
                       A        B      A:B Residuals
Sum of Squares   8.33333 21.33333  8.00000  86.33333
Deg. of Freedom        1        1        1        14

Residual standard error: 2.483277
Estimated effects may be unbalanced

模型系数输出 #

通过coefficients(lewensduur_anova)提取的系数：

(Intercept)  -34
A            1.36666666666667
B            0.213333333333333
A:B         -0.008

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