多项式数据与R的glm()

Question

如何让R的glm()匹配多项式数据？我已经尝试了几次'family=AAA(link="BBB")‘的迭代，但我似乎无法得到与之匹配的微不足道的预测。

例如，请使用R的glm帮助匹配多项式数据

x=seq(-6,6,2)
y=x*x
parabola=data.frame(x,y)
plot(parabola)
model=glm(y~x,dat=parabola)
test=data.frame(x=seq(-5,5,2))
test$y=predict(model,test)
plot(test)

plot(抛物线)看起来和预期的一样，但我可以找到glm()的咒语，它将使plot(test)看起来像抛物线。

Answer 1

我认为您需要后退一步，开始考虑模型以及如何在R中表示这一点。在您的示例中，y是x的二次函数，因此您需要在模型公式中包括x和x^2，即作为预测器，您需要在给定手头数据的情况下估计x和x^2对响应的影响。

如果y是高斯的，取决于模型，那么您可以使用lm()实现这一点，或者

y ~ x + I(x^2)

或

y ~ poly(x, 2)

首先，我们将二次项包装在I()中，因为^运算符在R模型公式中具有特殊含义(而不是其数学含义)。后一版本给出了正交多项式，因此x和x^2项不会相互关联，这有助于拟合，但在某些情况下，使用poly()解释系数会更棘手。

把所有这些放在一起，我们有(请注意，我在y中添加了一些随机误差，以便不能完美地预测它，因为我使用的示例在现实中更常见)：

x <- seq(-6 ,6 ,2)
y <- x^2 + rnorm(length(x), sd = 2)
parabola <- data.frame(x = x, y = y)

mod <- lm(y ~ poly(x, 2), data = parabola)

plot(parabola)
lines(fitted(mod) ~ x, data = parabola, col = "red")

生成的图是：

​

​

另一个问题是y是否是高斯的？如果y不能为负(即计数)，并且/或者是离散的，那么使用lm()建模将是错误的。这就是glm()的用武之地，通过它，您可以在不需要x^2的情况下拟合曲线(尽管如果数据真的是抛物线，那么x本身不会拟合响应)，因为存在从线性预测器到响应尺度的数据的显式转换。

最好考虑数据的属性和您想要拟合的模型类型，然后在该建模框架内建立多项式的次数，而不是跳入各种咒语中简单地对数据进行曲线拟合。

Answer 2

匹配现在是完美的。一个稍微有趣的抛物线：

x=seq(-16,16,2)
y= 4*x*x + 10*x + 6
parabola=data.frame(x,y)
plot(parabola)
model=lm(y~poly(x,2),dat=parabola)
summary(model)
test=data.frame(x=seq(-15,15,2))
test$y=predict(model,test)
points(test,pch=3)

一个业余爱好者(像我一样)可能希望模型的系数是(4,10,6)匹配4*x*x + 10*x +6

Call:
lm(formula = y ~ poly(x, 2), data = parabola)

Residuals:
       Min         1Q     Median         3Q        Max 
-3.646e-13 -8.748e-14 -3.691e-14  4.929e-14  6.387e-13 

Coefficients:
             Estimate Std. Error   t value Pr(>|t|)    
(Intercept) 3.900e+02  5.192e-14 7.511e+15   <2e-16 ***
poly(x, 2)1 4.040e+02  2.141e-13 1.887e+15   <2e-16 ***
poly(x, 2)2 1.409e+03  2.141e-13 6.581e+15   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.141e-13 on 14 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 2.343e+31 on 2 and 14 DF,  p-value: < 2.2e-16

为什么系数会是(390,404,1409)？