R语言poly拟合立方模型系数异常但预测准确问题咨询

问题描述 #

在R中拟合立方模型并预测时，出现模型预测结果准确，但直接用summary()输出的系数手动计算预测值与predict()结果完全不符的情况：

1. 拟合模型的代码及summary输出：

summary(lm(Blur ~ poly(logMAR, 3), data = df))

输出结果：

Call:
lm(formula = Blur ~ poly(logMAR, 3), data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.7838  -2.0661  -0.0868   2.3009  23.8918 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       35.0000     0.7775  45.015   <2e-16 ***
poly(logMAR, 3)1 230.6529     6.7782  34.029   <2e-16 ***
poly(logMAR, 3)2  94.3878     6.7782  13.925   <2e-16 ***
poly(logMAR, 3)3  19.6777     6.7782   2.903   0.0049 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.778 on 72 degrees of freedom
Multiple R-squared:  0.9497,    Adjusted R-squared:  0.9476 
F-statistic: 453.4 on 3 and 72 DF,  p-value: < 2.2e-16

2. 提取系数结果：

coef(model)

输出：

(Intercept) poly(logMAR, 3)1 poly(logMAR, 3)2 poly(logMAR, 3)3 
        35.00000        230.65294         94.38777         19.67767 

3. predict()输出的合理预测值：

predict(model, newdata = data.frame(logMAR = c(0.37, 0.74, 1.15)))

输出：

1         2         3 
 3.492585 11.469869 31.099493 

4. 手动用系数计算的错误结果：

# 手动计算示例：35 + 230.65*0.37 +94.39*0.37² +19.68*0.37³
[1] 134.2600 265.3438 455.0060

对比图表可知，predict()的结果更合理，且与其他软件拟合的原始多项式系数（β0=1.3415、β1=0.9073、β2=9.2517、β3=10.8354）手动计算结果一致。

问题原因 #

这是因为R中poly()函数默认使用正交多项式（orthogonal polynomials），而非原始的单项式多项式（raw polynomials）：

• 正交多项式的系数是对应正交基的权重，不是直接对应logMAR、logMAR²、logMAR³的系数，因此不能直接用原始x值去乘这些系数计算预测值。
• predict()函数会自动使用模型中生成的正交基计算预测值，因此结果是正确的。

解决方法 #

如果需要得到原始多项式的系数（即对应x、x²、x³的系数），可以在poly()中添加raw=TRUE参数：

model_raw <- lm(Blur ~ poly(logMAR, 3, raw=TRUE), data = df)
summary(model_raw)

此时coef(model_raw)输出的系数可以直接用原始x值手动计算，结果会和predict(model_raw, ...)一致。

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