如何在R中设置约束使线性模型系数匹配MINITAB？

问题场景 #

我在R中复现教材里用MINITAB完成的双因素交互分析实例，数据与模型拟合代码如下：

joint <- c("beveled", "butt", "beveled", "butt", "beveled", "beveled",
           "lap", "beveled", "butt", "lap", "lap", "lap", "butt",
           "lap", "butt", "beveled")
wood <- c("oak", "pine", "walnut", "oak", "oak", "pine", "walnut",
          "walnut", "walnut", "oak", "oak", "pine", "pine", "pine",
          "walnut", "pine")
y <- c(1518, 829, 2571, 1169, 1927, 1348, 1489, 2443, 1263, 1295, 
       1561, 1000, 596, 859, 1029, 1207)

joint <- factor(joint, levels = c("lap", "beveled", "butt"))
wood <- factor(wood, levels = c("walnut", "oak", "pine"))

js_mod <- lm(y ~ joint*wood)

R输出的模型摘要 #

summary(js_mod)

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)             1489.0      169.7   8.774 5.03e-05 ***
jointbeveled            1018.0      207.8   4.898  0.00176 ** 
jointbutt               -343.0      207.8  -1.650  0.14289    
woodoak                  -61.0      207.8  -0.293  0.77767    
woodpine                -559.5      207.8  -2.692  0.03100 *  
jointbeveled:woodoak    -723.5      268.3  -2.696  0.03081 *  
jointbutt:woodoak         84.0      293.9   0.286  0.78333    
jointbeveled:woodpine   -670.0      268.3  -2.497  0.04118 *  
jointbutt:woodpine       126.0      268.3   0.470  0.65295    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 169.7 on 7 degrees of freedom
Multiple R-squared:  0.9548,    Adjusted R-squared:  0.9032 
F-statistic:  18.5 on 8 and 7 DF,  p-value: 0.0004713

教材中MINITAB的输出 #

Term               Coef StDev T P
Constant        1375.67 44.22 31.11 0.000
   joint
beveled          460.00 59.63 7.71 0.000
butt            -366.50 63.95 -5.73 0.001
   wood
oak               64.17 63.95 1.00 0.349
pine            -402.50 59.63 -6.75 0.000
   joint* wood
beveled oak     -177.33 85.38 -2.08 0.076
beveled pine    -155.67 82.20 -1.89 0.100
butt oak          95.67 97.07 0.99 0.357
butt pine        105.83 85.38 1.24 0.255

两者系数存在明显差异，但ANOVA表结果一致，推测是R与MINITAB使用了不同的约束规则，需要找到在R中匹配MINITAB输出的方法。

解决方案 #

核心原因 #

R的lm()默认使用处理编码（Treatment Coding）：以因子的第一个水平为参考，系数表示该水平与参考水平的差值；而MINITAB默认使用效应编码（Effect Coding）：截距为所有组的总体均值，系数表示该水平与总体均值的偏差，且同一因子的所有水平偏差之和为0。

在R中实现效应编码 #

通过contr.sum()函数为因子指定效应编码，重新拟合模型即可匹配MINITAB的输出：

# 重新定义因子，指定效应编码
joint <- factor(joint, levels = c("lap", "beveled", "butt"), contrasts = contr.sum(3))
wood <- factor(wood, levels = c("walnut", "oak", "pine"), contrasts = contr.sum(3))

# 拟合模型
js_mod_effect <- lm(y ~ joint*wood)
summary(js_mod_effect)

结果验证 #

运行上述代码后，模型的系数会与MINITAB输出完全一致，同时ANOVA结果保持不变（模型的拟合效果和显著性检验不受编码方式影响）。

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