热线电话:13121318867

登录
首页精彩阅读R语言多项式回归
R语言多项式回归
2016-10-17
收藏

R语言多项式回归

含有x和y这两个变量的线性回归是所有回归分析中最常见的一种;而且,在描述它们关系的时候,也是最有效、最容易假设的一种模型。然而,有些时候,它的实际情况下某些潜在的关系是非常复杂的,不是二元分析所能解决的,而这时,我们需要多项式回归分析来找到这种隐藏的关系。

让我们看一下经济学里的一个例子:假设你要买一个具体的产品,而你要买的个数是q。如果产品的单价是p,然后,你要给y元。其实,这就是一个很典型的线性关系。而总价和产品数量呈正比例关系。下面,根据这个实例,我们敲击行代码来作它们的线性关系图:

  1. p <- 0.5
  2. q <- seq(0,100,1)
  3. y <- p*q
  4. plot(q,y,type='l',col='red',main='Linear relationship')

下面是它的线性关系图:

现在,我们看到这确实是一个不错的估计,这个图很好的模拟成q和y的线性关系。然而,当我们在做买卖要考虑别的因素的时候,诸如这种商品要买多少,很有可能,我们可以通过询问和讨价赚得折扣,或者,当我们越来越多的买一种具体的商品的时候,我们也可能让这种商品升价了。
这样,我们根据上面的条件,我们在写脚本的时候,我们要注意,总价与产品的数量不再具有线性关系了:

  1. y <- 450 + p*(q-10)^3
  2. plot(q,y,type='l',col='navy',main='Nonlinear relationship',lwd=3)

利用多项式回归,我们可以拟合n>1张订单所产生的数据的模型,并且能试着建一个非线性模型。

怎样拟合一个多项式回归

首先,当我们要创建一串虚拟随机数的时候,我们必须总要记得写set.seed(n)。这样做,随机数生成器总能产生同等数目的数据。

  1. set.seed(20)

预测变量q:使用seq来快速产生等间距的序列:

  1. q <- seq(from=0, to=20, by=0.1)

预测y值:

  1. y <- 500 + 0.4 * (q-10)^3

我们现在产生一些噪音并把它添加到模型中:

  1. noise <- rnorm(length(q), mean=10, sd=80)
  2. noisy.y <- y + noise

对噪声数据进行画图:

  1. plot(q,noisy.y,col='deepskyblue4',xlab='q',main='Observed data')
  2. lines(q,y,col='firebrick1',lwd=3)

下面的这个图根据观测数据进行模拟。其中,模拟的图的散点是蓝色的,而红色线则是信号(信号是一种术语,它通常用于表示我们感兴趣的东西的通常变化趋势)。

我们得出的模型应当是 y = aq + bq2 + c*q3 + cost。

现在,我们用R对此进行模拟。要拟合一个多项式模型,你也可以这样用:

  1. model <- lm(noisy.y ~ poly(q,3))

或者:

  1. model <- lm(noisy.y ~ x + I(X^2) + I(X^3))

然而,我们要知道q,I(q^2),I(q^3)存在相关的关系,而这些相关变量很有可能引起某些问题的产生。这时,使用poly()可以避免这个问题,因为它是创建一个垂直的多项式。因此,我喜欢第一种方法:

  1. summary(model)
  2. Call:
  3. lm(formula = noisy.y ~ poly(q, 3))
  4. Residuals:
  5.      Min       1Q   Median       3Q      Max
  6. -212.326  -51.186    4.276   61.485  165.960
  7. Coefficients:
  8.             Estimate Std. Error t value Pr(>|t|)
  9. (Intercept)  513.615      5.602   91.69   <2e-16 ***
  10. poly(q, 3)1 2075.899     79.422   26.14   <2e-16 ***
  11. poly(q, 3)2 -108.004     79.422   -1.36    0.175
  12. poly(q, 3)3  864.025     79.422   10.88   <2e-16 ***
  13. ---
  14. Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  15. Residual standard error: 79.42 on 197 degrees of freedom
  16. Multiple R-squared:  0.8031,    Adjusted R-squared:  0.8001
  17. F-statistic: 267.8 on 3 and 197 DF,  p-value:  0

我们可以使用confint()来获得一个模型的参数的置信区间。

以下是模型参数的置信区间:

  1. confint(model, level=0.95)
  2.                 2.5 %     97.5 %
  3. (Intercept)  502.5676  524.66261
  4. poly(q, 3)1 1919.2739 2232.52494
  5. poly(q, 3)2 -264.6292   48.62188
  6. poly(q, 3)3  707.3999 1020.65097

现在,我们要作一个拟合VS残差图。如果这是一个拟合效果比较不错的模型,我们应该看不到任何一种模型的模式特征

  1. plot(fitted(model),residuals(model))

整体来说,这个模型的拟合效果还是不错的,毕竟残差为0.8。第一和第三个订单序列的系数,在统计学当中,是相当这样的,这样在我们的意料之中。现在,我们可以使用predict()函数来获得拟合数据以及置信区间,这样,我们可以不按照数据来作图。 下面是预测值和预测置信区间:

  1. predicted.intervals <- predict(model,data.frame(x=q),interval='confidence',level=0.99)

在已有的图像中添加拟合线:

  1. lines(q,predicted.intervals[,1],col='green',lwd=3)
  2. lines(q,predicted.intervals[,2],col='black',lwd=1)
  3. lines(q,predicted.intervals[,3],col='black',lwd=1)

添加图例:

  1. legend("bottomright",c("Observ.","Signal","Predicted"),
  2. col=c("deepskyblue4","red","green"), lwd=3)

下面是它的拟合图像:

我们可以看到我们的模型在数据的拟合方面做的不错,我们也因此感到非常满意。

注意:多项式回归是一种更能强大的工具。可是,我们也可能得到事与愿违的结果:在这个例子中,我们知道我们的信号是使用三次多项式而产生的,然而,当我们在分析实际数据的时候,我们通常对此不知情,因此,正因为多项式次数n大于4的时候会产生过度拟合的情况,我们要在这里注意一下。但你的模型取了噪音而不是信号的时候会产生过拟合的情况;甚至,当你在现有的数据进行模型优化的时候,当你要尝试预测新的数据的时候就不好了,它会导致缺失值的产生。

数据分析咨询请扫描二维码

最新资讯
更多
客服在线
立即咨询