• Posts tagged "线性回归"

Blog Archives

Posted:
Jul 27, 2016
Tags:
Comments:
15 Comments

R语言解读多元线性回归模型

R的极客理想系列文章，涵盖了R的思想，使用，工具，创新等的一系列要点，以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言，一直在小众领域闪耀着光芒。直到大数据的爆发，R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入，R语言的社区在迅速扩大成长。现在已不仅仅是统计领域，教育，银行，电商，互联网….都在使用R语言。

要成为有理想的极客，我们不能停留在语法上，要掌握牢固的数学，概率，统计知识，同时还要有创新精神，把R语言发挥到各个领域。让我们一起动起来吧，开始R的极客理想。

关于作者：

  • 张丹(Conan), 程序员R,Nodejs,Java
  • weibo：@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处：
http://blog.fens.me/r-multi-linear-regression

reg-multi-liner

前言

本文接上一篇R语言解读一元线性回归模型。在许多生活和工作的实际问题中，影响因变量的因素可能不止一个，比如对于知识水平越高的人，收入水平也越高，这样的一个结论。这其中可能包括了因为更好的家庭条件，所以有了更好的教育；因为在一线城市发展，所以有了更好的工作机会；所处的行业赶上了大的经济上行周期等。要想解读这些规律，是复杂的、多维度的，多元回归分析方法更适合解读生活的规律。

由于本文为非统计的专业文章，所以当出现与教课书不符的描述，请以教课书为准。本文力求用简化的语言，来介绍多元线性回归的知识，同时配合R语言的实现。

目录

  1. 多元线性回归介绍
  2. 元线性回归建模
  3. 模型优化
  4. 案例：黑色系期货日K线数据验证

1. 多元线性回归介绍

对比一元线性回归，多元线性回归是用来确定2个或2个以上变量间关系的统计分析方法。多元线性回归的基本的分析方法与一元线性回归方法是类似的，我们首先需要对选取多元数据集并定义数学模型，然后进行参数估计，对估计出来的参数进行显著性检验，残差分析，异常点检测，最后确定回归方程进行模型预测。

由于多元回归方程有多个自变量，区别于一元回归方程，有一项很重要的操作就是自变量的优化，挑选出相关性最显著的自变量，同时去除不显著的自变量。在R语言中，有很方便地用于优化函数，可以很好的帮助我们来改进回归模型。

下面就开始多元线性回归的建模过程。

2. 多元线性回归建模

做过商品期货研究的人，都知道黑色系品种是具有产业链上下游的关系。铁矿石是炼钢的原材料，焦煤和焦炭是炼钢的能源资源，热卷即热轧卷板是以板坯为原料经加热后制成的钢板，螺纹钢是表面带肋的钢筋。

由于有产业链的关系，假设我们想要预测螺纹钢的价格，那么影响螺纹钢价格的因素可以会涉及到原材料，能源资源和同类材料等。比如，铁矿石价格如果上涨，螺纹钢就应该要涨价了。

2.1 数据集和数学模型

先从数据开始介绍，这次的数据集，我选择的期货黑色系的品种的商品期货，包括了大连期货交易所的 焦煤(JM)，焦炭(J)，铁矿石(I)，上海期货交易所的 螺纹钢(RU) 和 热卷(HC)。

数据集为2016年3月15日，当日白天开盘的交易数据，为黑色系的5个期货合约的分钟线的价格数据。


# 数据集已存在df变量中
> head(df,20)
                       x1    x2    x3   x4    y
2016-03-15 09:01:00 754.5 616.5 426.5 2215 2055
2016-03-15 09:02:00 752.5 614.5 423.5 2206 2048
2016-03-15 09:03:00 753.0 614.0 423.0 2199 2044
2016-03-15 09:04:00 752.5 613.0 422.5 2197 2040
2016-03-15 09:05:00 753.0 615.5 424.0 2198 2043
2016-03-15 09:06:00 752.5 614.5 422.0 2195 2040
2016-03-15 09:07:00 752.0 614.0 421.5 2193 2036
2016-03-15 09:08:00 753.0 615.0 422.5 2197 2043
2016-03-15 09:09:00 754.0 615.5 422.5 2197 2041
2016-03-15 09:10:00 754.5 615.5 423.0 2200 2044
2016-03-15 09:11:00 757.0 616.5 423.0 2201 2045
2016-03-15 09:12:00 756.0 615.5 423.0 2200 2044
2016-03-15 09:13:00 755.5 615.0 423.0 2197 2042
2016-03-15 09:14:00 755.5 615.0 423.0 2196 2042
2016-03-15 09:15:00 756.0 616.0 423.5 2200 2045
2016-03-15 09:16:00 757.5 616.0 424.0 2205 2052
2016-03-15 09:17:00 758.5 618.0 424.0 2204 2051
2016-03-15 09:18:00 759.5 618.5 424.0 2205 2053
2016-03-15 09:19:00 759.5 617.5 424.5 2206 2053
2016-03-15 09:20:00 758.5 617.5 423.5 2201 2050

数据集包括有6列:

  • 索引, 为时间
  • x1, 为焦炭(j1605)合约的1分钟线的报价数据
  • x2, 为焦煤(jm1605)合约的1分钟线的报价数据
  • x3, 为铁矿石(i1605)合约的1分钟线的报价数据
  • x4, 为热卷(hc1605)合约的1分钟线的报价数据
  • y, 为螺纹钢(rb1605)合约的1分钟线的报价数据

假设螺纹钢的价格与其他4个商品的价格有线性关系，那么我们建立以螺纹钢为因变量，以焦煤、焦炭、铁矿石和热卷的为自变量的多元线性回归模型。用公式表示为：

y = a + b * x1 + c * x2 + d * x3 + e * x4 + ε
  • y，为因变量，螺纹钢
  • x1，为自变量，焦煤
  • x2，为自变量，焦炭
  • x3，为自变量，铁矿石
  • x4，为自变量，热卷
  • a，为截距
  • b,c,d,e，为自变量系数
  • ε, 为残差，是其他一切不确定因素影响的总和，其值不可观测。假定ε服从正态分布N(0,σ^2)。

通过对多元线性回归模型的数学定义，接下来让我们利用数据集做多元回归模型的参数估计。

2.2. 回归参数估计

上面公式中，回归参数 a, b, c, d,e都是我们不知道的，参数估计就是通过数据来估计出这些参数，从而确定自变量和因变量之前的关系。我们的目标是要计算出一条直线，使直线上每个点的Y值和实际数据的Y值之差的平方和最小，即（Y1实际-Y1预测）^2+（Y2实际-Y2预测）^2+ …… +（Yn实际-Yn预测）^2 的值最小。参数估计时，我们只考虑Y随X自变量的线性变化的部分，而残差ε是不可观测的，参数估计法并不需要考虑残差。

类似于一元线性回归，我们用R语言来实现对数据的回归模型的参数估计，用lm()函数来实现多元线性回归的建模过程。


# 建立多元线性回归模型
> lm1<-lm(y~x1+x2+x3+x4,data=df)

# 打印参数估计的结果
> lm1

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4  
   212.8780       0.8542       0.6672      -0.6674       0.4821

这样我们就得到了y和x关系的方程。

y = 212.8780 + 0.8542 * x1 + 0.6672 * x2 - 0.6674 * x3 + 0.4821 * x4

2.3. 回归方程的显著性检验

参考一元线性回归的显著性检验，多元线性回归的显著性检验，同样是需要经过 T检验，F检验，和R^2(R平方)相关系统检验。在R语言中这三种检验的方法都已被实现，我们只需要把结果解读，我们可以summary()函数来提取模型的计算结果。


> summary(lm1)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.9648 -1.3241 -0.0319  1.2403  5.4194 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 212.87796   58.26788   3.653 0.000323 ***
x1            0.85423    0.10958   7.795 2.50e-13 ***
x2            0.66724    0.12938   5.157 5.57e-07 ***
x3           -0.66741    0.15421  -4.328 2.28e-05 ***
x4            0.48214    0.01959  24.609  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.028 on 221 degrees of freedom
Multiple R-squared:  0.9725,	Adjusted R-squared:  0.972 
F-statistic:  1956 on 4 and 221 DF,  p-value: < 2.2e-16
  • T检验：所自变量都是非常显著***
  • F检验：同样是非常显著，p-value < 2.2e-16
  • 调整后的R^2：相关性非常强为0.972

最后，我们通过的回归参数的检验与回归方程的检验，得到最后多元线性回归方程为：


y = 212.87796 + 0.85423 * x1 + 0.66724 * x2 - 0.66741 * x3 + 0.48214 * x4

即

螺纹钢 = 212.87796 + 0.85423 * 焦炭 + 0.66724 * 焦煤 - 0.66741 * 铁矿石 + 0.48214 * 热卷

2.4 残差分析和异常点检测

在得到的回归模型进行显著性检验后，还要在做残差分析(预测值和实际值之间的差)，检验模型的正确性，残差必须服从正态分布N(0,σ^2)。直接用plot()函数生成4种用于模型诊断的图形，进行直观地分析。


> par(mfrow=c(2,2))
> plot(lm1)

m01

  • 残差和拟合值(左上)，残差和拟合值之间数据点均匀分布在y=0两侧，呈现出随机的分布，红色线呈现出一条平稳的曲线并没有明显的形状特征。
  • 残差QQ图(右上)，数据点按对角直线排列，趋于一条直线，并被对角直接穿过，直观上符合正态分布。
  • 标准化残差平方根和拟合值(左下)，数据点均匀分布在y=0两侧，呈现出随机的分布，红色线呈现出一条平稳的曲线并没有明显的形状特征。
  • 标准化残差和杠杆值(右下)，没有出现红色的等高线，则说明数据中没有特别影响回归结果的异常点。

结论，没有明显的异常点，残差符合假设条件。

2.5. 模型预测

我们得到了多元线性回归方程的公式，就可以对数据进行预测了。我们可以用R语言的predict()函数来计算预测值y0和相应的预测区间，并把实际值和预测值一起可视化化展示。


> par(mfrow=c(1,1))  #设置画面布局

# 预测计算
> dfp<-predict(lm1,interval="prediction")

# 打印预测时
> head(dfp,10)
                fit      lwr      upr
2014-03-21 3160.526 3046.425 3274.626
2014-03-24 3193.253 3078.868 3307.637
2014-03-25 3240.389 3126.171 3354.607
2014-03-26 3228.565 3114.420 3342.710
2014-03-27 3222.528 3108.342 3336.713
2014-03-28 3262.399 3148.132 3376.666
2014-03-31 3291.996 3177.648 3406.344
2014-04-01 3305.870 3191.447 3420.294
2014-04-02 3275.370 3161.018 3389.723
2014-04-03 3297.358 3182.960 3411.755

# 合并数据
> mdf<-merge(df$y,dfp)	 

# 画图
> draw(mdf)

m02

图例说明

  • y, 实际价格，红色线
  • fit, 预测价格，绿色线
  • lwr，预测最低价，蓝色线
  • upr，预测最高价，紫色线

从图中看出，实际价格y和预测价格fit，在大多数的时候都是很贴近的。我们的一个模型就训练好了！

3. 模型优化

上文中，我们已经很顺利的发现了一个非常不错的模型。如果要进行模型优化，可以用R语言中update()函数进行模型的调整。我们首先检查一下每个自变量x1,x2,x3,x4和因变量y之间的关系。

pairs(as.data.frame(df))

m03

从图中，我们可以发现x2与Y的关系，可能是最偏离线性的。那么，我们尝试对多元线性回归模型进行调整，从原模型中去掉x2变量。



# 模型调整
> lm2<-update(lm1, .~. -x2)

> summary(lm2)

Call:
lm(formula = y ~ x1 + x3 + x4, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.0039 -1.3842  0.0177  1.3513  4.8028 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 462.47104   34.26636   13.50  < 2e-16 ***
x1            1.08728    0.10543   10.31  < 2e-16 ***
x3           -0.40788    0.15394   -2.65  0.00864 ** 
x4            0.42582    0.01718   24.79  < 2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.142 on 222 degrees of freedom
Multiple R-squared:  0.9692,	Adjusted R-squared:  0.9688 
F-statistic:  2330 on 3 and 222 DF,  p-value: < 2.2e-16

当把自变量x2去掉后，自变量x3的T检验反而变大了，同时Adjusted R-squared变小了，所以我们这次调整是有问题的。

如果通过生产和原材料的内在逻辑分析，焦煤与焦炭属于上下游关系。焦煤是生产焦炭的一种原材料，焦炭是焦煤与其他炼焦煤经过配煤焦化形成的产品，一般生产 1 吨焦炭需要1.33 吨炼焦煤，其中焦煤至少占 30% 。

我们把焦煤 和 焦炭的关系改变一下，增加x1*x2的关系匹配到模型，看看效果。


# 模型调整
> lm3<-update(lm1, .~. + x1*x2)
> summary(lm3)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x1:x2, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8110 -1.3501 -0.0595  1.2019  5.3884 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7160.32231 7814.50048   0.916    0.361    
x1            -8.45530   10.47167  -0.807    0.420    
x2           -10.58406   12.65579  -0.836    0.404    
x3            -0.64344    0.15662  -4.108 5.63e-05 ***
x4             0.48363    0.01967  24.584  < 2e-16 ***
x1:x2          0.01505    0.01693   0.889    0.375    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.029 on 220 degrees of freedom
Multiple R-squared:  0.9726,	Adjusted R-squared:  0.972 
F-statistic:  1563 on 5 and 220 DF,  p-value: < 2.2e-16

从结果中发现，增加了x1*x2列后，原来的x1,x2和Intercept的T检验都不显著。继续调整模型，从模型中去掉x1,x2两个自变量。


# 模型调整
> lm4<-update(lm3, .~. -x1-x2)
> summary(lm4)

Call:
lm(formula = y ~ x3 + x4 + x1:x2, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.9027 -1.2516 -0.0167  1.2748  5.8683 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.950e+02  1.609e+01  43.183  < 2e-16 ***
x3          -6.284e-01  1.530e-01  -4.108 5.61e-05 ***
x4           4.959e-01  1.785e-02  27.783  < 2e-16 ***
x1:x2        1.133e-03  9.524e-05  11.897  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.035 on 222 degrees of freedom
Multiple R-squared:  0.9722,	Adjusted R-squared:  0.9718 
F-statistic:  2588 on 3 and 222 DF,  p-value: < 2.2e-16

从调整后的结果来看，效果还不错。不过，也并没有比最初的模型有所提高。

对于模型调整的过程，如果我们手动调整测试时，一般都会基于业务知识来操作。如果是按照数据指标来计算，我们可以用R语言中提供的逐步回归的优化方法，通过AIC指标来判断是否需要参数优化。


#对lm1模型做逐步回归
> step(lm1)
Start:  AIC=324.51
y ~ x1 + x2 + x3 + x4

       Df Sum of Sq    RSS    AIC
               908.8 324.51
- x3    1     77.03  985.9 340.90
- x2    1    109.37 1018.2 348.19
- x1    1    249.90 1158.8 377.41
- x4    1   2490.56 3399.4 620.65

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4  
   212.8780       0.8542       0.6672      -0.6674       0.4821

通过计算AIC指标，lm1的模型AIC最小时为324.51，每次去掉一个自变量都会让AIC的值变大，所以我们还是不调整比较好。

对刚才的lm3模型做逐步回归的模型调整。


#对lm3模型做逐步回归
> step(lm3)
Start:  AIC=325.7               #当前AIC
y ~ x1 + x2 + x3 + x4 + x1:x2

        Df Sum of Sq    RSS    AIC
- x1:x2  1      3.25  908.8 324.51
                905.6 325.70
- x3     1     69.47  975.1 340.41
- x4     1   2487.86 3393.5 622.25

Step:  AIC=324.51               #去掉x1*x2项的AIC
y ~ x1 + x2 + x3 + x4

       Df Sum of Sq    RSS    AIC
               908.8 324.51
- x3    1     77.03  985.9 340.90
- x2    1    109.37 1018.2 348.19
- x1    1    249.90 1158.8 377.41
- x4    1   2490.56 3399.4 620.65

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4  
   212.8780       0.8542       0.6672      -0.6674       0.4821

通过AIC的判断，去掉X1*X2项后AIC最小，最后的检验结果告诉我们，还是原初的模型是最好的。

4. 案例：黑色系期货日K线数据验证

最后，我们用上面5个期货合约的日K线数据测试一下，找到多元回归关系。


> lm9<-lm(y~x1+x2+x3+x4,data=df)  # 日K线数据
> summary(lm9)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-173.338  -37.470    3.465   32.158  178.982 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 386.33482   31.07729  12.431  < 2e-16 ***
x1            0.75871    0.07554  10.045  < 2e-16 ***
x2           -0.62907    0.14715  -4.275 2.24e-05 ***
x3            1.16070    0.05224  22.219  < 2e-16 ***
x4            0.46461    0.02168  21.427  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 57.78 on 565 degrees of freedom
Multiple R-squared:  0.9844,	Adjusted R-squared:  0.9843 
F-statistic:  8906 on 4 and 565 DF,  p-value: < 2.2e-16

数据集的基本统计信息。


> summary(df)
     Index                           x1               x2       
 Min.   :2014-03-21 00:00:00   Min.   : 606.5   Min.   :494.0  
 1st Qu.:2014-10-21 06:00:00   1st Qu.: 803.5   1st Qu.:613.1  
 Median :2015-05-20 12:00:00   Median : 939.0   Median :705.8  
 Mean   :2015-05-21 08:02:31   Mean   : 936.1   Mean   :695.3  
 3rd Qu.:2015-12-16 18:00:00   3rd Qu.:1075.0   3rd Qu.:773.0  
 Max.   :2016-07-25 00:00:00   Max.   :1280.0   Max.   :898.0  

       x3              x4             y       
 Min.   :284.0   Min.   :1691   Min.   :1626  
 1st Qu.:374.1   1st Qu.:2084   1st Qu.:2012  
 Median :434.0   Median :2503   Median :2378  
 Mean   :476.5   Mean   :2545   Mean   :2395  
 3rd Qu.:545.8   3rd Qu.:2916   3rd Qu.:2592  
 Max.   :825.0   Max.   :3480   Max.   :3414

m04

对于日K线数据，黑色系的5个品种，同样具有非常强的相关关系，那么我们就可以把这个结论应用到实际的交易中了。

本文通过多元回归的统计分析方法，介绍多元回归在金融市场的基本应用。我们通过建立因变量和多个自变量的模型，从而发现生活中更复杂的规律，并建立有效的验证指标。让我们我们的技术优势，去金融市场抢钱吧。

转载请注明出处：
http://blog.fens.me/r-multi-linear-regression

打赏作者

Posted:
Jul 25, 2016
Tags:
Comments:
6 Comments

R语言解读一元线性回归模型

R的极客理想系列文章，涵盖了R的思想，使用，工具，创新等的一系列要点，以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言，一直在小众领域闪耀着光芒。直到大数据的爆发，R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入，R语言的社区在迅速扩大成长。现在已不仅仅是统计领域，教育，银行，电商，互联网….都在使用R语言。

要成为有理想的极客，我们不能停留在语法上，要掌握牢固的数学，概率，统计知识，同时还要有创新精神，把R语言发挥到各个领域。让我们一起动起来吧，开始R的极客理想。

关于作者：

  • 张丹(Conan), 程序员R,Nodejs,Java
  • weibo：@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处：
http://blog.fens.me/r-linear-regression/

reg-liner

前言

在我们的日常生活中，存在大量的具有相关性的事件，比如大气压和海拔高度，海拔越高大气压强越小；人的身高和体重，普遍来看越高的人体重也越重。还有一些可能存在相关性的事件，比如知识水平越高的人，收入水平越高；市场化的国家经济越好，则货币越强势，反而全球经济危机，黄金等避险资产越走强。

如果我们要研究这些事件，找到不同变量之间的关系，我们就会用到回归分析。一元线性回归分析是处理两个变量之间关系的最简单模型，是两个变量之间的线性相关关系。让我们一起发现生活中的规律吧。

由于本文为非统计的专业文章，所以当出现与教课书不符的描述，请以教课书为准。本文力求用简化的语言，来介绍一元线性回归的知识，同时配合R语言的实现。

目录

  1. 一元线性回归介绍
  2. 数据集和数学模型
  3. 回归参数估计
  4. 回归方程的显著性检验
  5. 残差分析和异常点检测
  6. 模型预测

1. 一元线性回归介绍

回归分析（Regression Analysis)是用来确定2个或2个以上变量间关系的一种统计分析方法。如果回归分析中，只包括一个自变量X和一个因变量Y时，且它们的关系是线性的，那么这种回归分析称为一元线性回归分析。

回归分析属于统计学的基本模型，涉及统计学基础，就会有一大堆的名词和知识点需要介绍。

在回归分析中，变量有2类：因变量 和 自变量。因变量通常是指实际问题中所关心的指标，用Y表示。而自变量是影响因变量取值的一个变量，用X表示，如果有多个自变量则表示为X1, X2, …, Xn。

回归分析研究的主要步骤：

  1. 确定因变量Y 与 自变量X1, X2, …, Xn 之间的定量关系表达式，即回归方程。
  2. 对回归方程的置信度检查。
  3. 判断自变量Xn(n=1,2,…,m)对因变量的影响。
  4. 利用回归方程进行预测。

本文会根据回归分析的的主要步骤，进行结构梳理，介绍一元线性回归模型的使用方法。

reg

2. 数据集和数学模型

先让我们通过一个例子开始吧，用一组简单的数据来说明一元线性回归分析的数学模型的原理和公式。找出下面数据集中Y与X的定量关系。

数据集为2016年3月1日，白天开盘的交易数据，为锌的2个期货合约的分钟线的价格数据。数据集包括有3列，索引列为时间，zn1.Close为ZN1604合约的1分钟线的报价数据，zn2.Close为ZN1605合约的1分钟线的报价数据。

数据集如下：


                    zn1.Close zn2.Close
2016-03-01 09:01:00     14075     14145
2016-03-01 09:02:00     14095     14160
2016-03-01 09:03:00     14095     14160
2016-03-01 09:04:00     14095     14165
2016-03-01 09:05:00     14120     14190
2016-03-01 09:06:00     14115     14180
2016-03-01 09:07:00     14110     14170
2016-03-01 09:08:00     14110     14175
2016-03-01 09:09:00     14105     14170
2016-03-01 09:10:00     14105     14170
2016-03-01 09:11:00     14120     14180
2016-03-01 09:12:00     14105     14170
2016-03-01 09:13:00     14105     14170
2016-03-01 09:14:00     14110     14175
2016-03-01 09:15:00     14105     14175
2016-03-01 09:16:00     14120     14185
2016-03-01 09:17:00     14125     14190
2016-03-01 09:18:00     14115     14185
2016-03-01 09:19:00     14135     14195
2016-03-01 09:20:00     14125     14190
2016-03-01 09:21:00     14135     14205
2016-03-01 09:22:00     14140     14210
2016-03-01 09:23:00     14140     14200
2016-03-01 09:24:00     14135     14205
2016-03-01 09:25:00     14140     14205
2016-03-01 09:26:00     14135     14205
2016-03-01 09:27:00     14130     14205

我们以zn1.Close列的价格为X，zn2.Close列的价格为Y，那么试试找到自变量X和因变量Y的关系的表达式。

为了直观起见，我们可以先画出一张散点图，以X为横坐标，Y为纵坐标，每个点对应一个X和一个Y。


# 数据集已存在df变量中
> head(df)
                    zn1.Close zn2.Close
2016-03-01 09:01:00     14075     14145
2016-03-01 09:02:00     14095     14160
2016-03-01 09:03:00     14095     14160
2016-03-01 09:04:00     14095     14165
2016-03-01 09:05:00     14120     14190
2016-03-01 09:06:00     14115     14180

# 分别给x,y赋值
> x<-as.numeric(df[,1])
> y<-as.numeric(df[,2])

# 画图
> plot(y~x+1)

01

从散点图上发现 X和Y 的排列基本是在一条直线附近，那么我们可以假设X和Y的关系是线性，可以用公式表式为。

Y = a + b * X + c
  • Y，为因变量
  • X，为自变量
  • a，为截距
  • b，为自变量系数
  • a+b*X, 表示Y随X的变化而线性变化的部分
  • c, 为残差或随机误差，是其他一切不确定因素影响的总和，其值不可观测。假定c是符合均值为0方差为σ^2的正态分布 ，记作c~N(0,σ^2)

对于上面的公式，称函数f(X) = a + b * X 为一元线性回归函数，a为回归常数，b为回归系数，统称回归参数。X 为回归自变量或回归因子，Y 为回归因变量或响应变量。如果(X1,Y1)，(X2,Y2)…(Xn,Yn)是(X,Y)的一组观测值，则一元线性回归模型可表示为


Yi = a + b * X + ci，     i= 1,2,...n

其中E(ci)=0, var(ci)=σ^2, i=1,2,...n

通过对一元线性回归模型的数学定义，接下来让我们利用数据集做回归模型的参数估计。

3. 回归参数估计

对于上面的公式，回归参数a,b是我们不知道的，我们需要用参数估计的方法来计算出a,b的值，而从得到数据集的X和Y的定量关系。我们的目标是要计算出一条直线，使直接线上每个点的Y值和实际数据的Y值之差的平方和最小，即（Y1实际-Y1预测）^2+（Y2实际-Y2预测）^2+ …… +（Yn实际-Yn预测）^2 的值最小。参数估计时，我们只考虑Y随X的线性变化的部分，而残差c是不可观测的，参数估计法并不需要考虑残差，对于残差的分析在后文中介绍。

令公式变形为a和b的函数Q(a,b), 即 (Y实际-Y测试)的平方和，变成到(Y实际 – (a+b*X))的平方和。

reg2

公式一 回归参数变形公式

通过最小二乘估计推导出a和b的求解公式，详细的推导过程请参考文章一元线性回归的细节

reg3

公式二 回归参数计算公式

其中 x和y的均值，计算方法如下

reg4

公式三 均值计算公式

有了这个公式，我们就可以求出a和b两个的回归参数的解了。

接下来，我们用R语言来实现对上面数据的回归模型的参数估计，R语言中可以用lm()函数来实现一元线性回归的建模过程。


# 建立线性回归模型
> lm.ab<-lm(y ~ 1+x)

# 打印参数估计的结果
> lm.ab

Call:
lm(formula = y ~ 1 + x)

Coefficients:
(Intercept)            x  
   -349.493        1.029

如果你想动手来计算也可以自己实现公式。


# x均值
> Xm<-mean(x);Xm 
[1] 14034.82

# y均值
> Ym<-mean(y);Ym
[1] 14096.76

# 计算回归系数
> b <- sum((x-Xm)*(y-Ym)) / sum((x-Xm)^2) ;b
[1] 1.029315

# 计算回归常数
> a <- Ym - b * Xm;a
[1] -349.493

回归参数a和b的计算结果，与lm()函数的计算结果是相同的。有了a和b的值，我们就可以画出这条近似的直接线。

计算公式为：

Y= a + b * X = -349.493 + 1.029315 * X

画出回归线。


> plot(y~x+1)
> abline(lm.ab)

02

这条直线是我们用数据拟合出来的，是一个近似的值。我们看到有些点在线上，有些点不在线上。那么要评价这条回归线拟合的好坏，我们就需要对回归模型进行显著性检验。

4. 回归方程的显著性检验

从回归参数的公式二可知，在计算过程中并不一定要知道Y和X是否有线性相关的关系。如果不存相关关系，那么回归方程就没有任何意义了，如果Y和X是有相关关系的，即Y会随着X的变化而线性变化，这个时候一元线性回归方程才有意义。所以，我们需要用假设检验的方法，来验证相关性的有效性。

通常会采用三种显著性检验的方法。

  • T检验法：T检验是检验模型某个自变量Xi对于Y的显著性，通常用P-value判断显著性，小于0.01更小时说明这个自变量Xi与Y相关关系显著。
  • F检验法：F检验用于对所有的自变量X在整体上看对于Y的线性显著性，也是用P-value判断显著性，小于0.01更小时说明整体上自变量与Y相关关系显著。
  • R^2(R平方)相关系统检验法：用来判断回归方程的拟合程度，R^2的取值在0，1之间，越接近1说明拟合程度越好。

在R语言中，上面列出的三种检验的方法都已被实现，我们只需要把结果解读。上文中，我们已经通过lm()函数构建一元线性回归模型，然后可以summary()函数来提取模型的计算结果。


> summary(lm.ab)      # 计算结果

Call:
lm(formula = y ~ 1 + x)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.9385  -2.2317  -0.1797   3.3546  10.2766 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.495e+02  7.173e+01  -4.872 2.09e-06 ***
x            1.029e+00  5.111e-03 201.390  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.232 on 223 degrees of freedom
Multiple R-squared:  0.9945,	Adjusted R-squared:  0.9945 
F-statistic: 4.056e+04 on 1 and 223 DF,  p-value: < 2.2e-16

模型解读：

  • Call，列出了回归模型的公式。
  • Residuals，列出了残差的最小值点，1/4分位点，中位数点，3/4分位点，最大值点。
  • Coefficients，表示参数估计的计算结果。
  • Estimate，为参数估计列。Intercept行表示常数参数a的估计值 ，x行表示自变量x的参数b的估计值。
  • Std. Error，为参数的标准差，sd(a), sd(b)
  • t value，为t值，为T检验的值
  • Pr(>|t|) ，表示P-value值，用于T检验判定，匹配显著性标记
  • 显著性标记，***为非常显著，**为高度显著, **为显著，·为不太显著，没有记号为不显著。
  • Residual standard error，表示残差的标准差，自由度为n-2。
  • Multiple R-squared，为相关系数R^2的检验，越接近1则越显著。
  • Adjusted R-squared，为相关系数的修正系数，解决多元回归自变量越多，判定系数R^2越大的问题。
  • F-statistic，表示F统计量，自由度为(1,n-2)，p-value:用于F检验判定，匹配显著性标记。

通过查看模型的结果数据，我们可以发现通过T检验的截距和自变量x都是非常显著，通过F检验判断出整个模型的自变量是非常显著，同时R^2的相关系数检验可以判断自变量和因变量是高度相关的。

最后，我们通过的回归参数的检验与回归方程的检验，得到最后一元线性回归方程为：

Y = -349.493 + 1.029315 * X

5. 残差分析和异常点检测

在得到的回归模型进行显著性检验后，还要在做残差分析(预测值和实际值之间的差)，检验模型的正确性，残差必须服从正态分布N(0,σ^2)。

我们可以自己计算数据残差，并进行正态分布检验。


# 残差
> y.res<-residuals(lm.ab)

# 打印前6条数据
> head(y.res)
        1         2         3         4         5         6 
6.8888680 1.3025744 1.3025744 6.3025744 5.5697074 0.7162808 

# 正态分布检验
> shapiro.test(y.res)

	Shapiro-Wilk normality test

data:  y.res
W = 0.98987, p-value = 0.1164

# 画出残差散点图
> plot(y.res)

09

对残差进行Shapiro-Wilk正态分布检验，W接近1，p-value>0.05，证明数据集符合正态分布！关于正态分布的介绍，请参考文章常用连续型分布介绍及R语言实现。

同时，我们也可以用R语言的工具生成4种用于模型诊断的图形，简化自己写代码计算的操作。


# 画图，回车展示下一张
> plot(lm.ab)    
Hit  to see next plot:   # 残差拟合图
Hit  to see next plot:   # 残差QQ图
Hit  to see next plot:   # 标准化的残差对拟合值 
Hit  to see next plot:   # 标准化残差对杠杆值

04

图1，残差和拟合值对比图

对残差和拟合值作图，横坐标是拟合值，纵坐标是残差。残差和拟合值之间，数据点均匀分布在y=0两侧，呈现出随机的分布，红色线呈现出一条平稳的曲线并没有明显的形状特征，说明残差数据表现非常好。

05

图2，残差QQ图

残差QQ图，用来描述残差是否符合正态分布。图中的数据点按对角直线排列，趋于一条直线，并被对角直接穿过，直观上符合正态分布。对于近似服从正态分布的标准化残差，应该有 95% 的样本点落在 [-2,2] 区间内。

06

图3，标准化残差平方根和拟合值对比图

对标准化残差平方根和拟合值作图，横坐标是拟合值，纵坐标是标准化后的残差平方根。与残差和拟合值对比图(图1)的判断方法类似，数据随机分布，红色线呈现出一条平稳的曲线，无明显的形状特征。

07

图4，标准残差和杠杆值对比图

对标准化残差和杠杆值作图，虚线表示的cooks距离等高线，通常用Cook距离度量的回归影响点。本图中没有出现红色的等高线，则说明数据中没有特别影响回归结果的异常点。

如果想把把4张图画在一起进行展示，可以改变画布布局。


> par(mfrow=c(2,2))
> plot(lm.ab)

08

看到上面4幅中，每幅图上都有一些点被特别的标记出来了，这些点是可能存在的异常值点，如果要对模型进行优化，我们可以从这些来入手。但终于本次残差分析的结果已经很好了，所以对于异常点的优化，可能并不能明显的提升模型的效果。

从图中发现，索引编号为27和192的2个点在多幅图中出现。我们假设这2个点为异常点，从数据中去掉这2个点，再进行显著性检验和残差分析。


# 查看27和192
> df[c(27,192),]
                    zn1.Close zn2.Close
2016-03-01 09:27:00     14130     14205
2016-03-01 14:27:00     14035     14085

# 新建数据集，去掉27和192
> df2<-df[-c(27,192),]

回归建模和显著性检验。


> x2<-as.numeric(df2[,1])
> y2<-as.numeric(df2[,2])
> lm.ab2<-lm(y2 ~ 1+x2)
> summary(lm.ab2)

Call:
lm(formula = y2 ~ 1 + x2)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.0356 -2.1542 -0.2727  3.3336  9.5879 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.293e+02  7.024e+01  -4.688 4.83e-06 ***
x2           1.028e+00  5.004e-03 205.391  < 2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.117 on 221 degrees of freedom
Multiple R-squared:  0.9948,	Adjusted R-squared:  0.9948 
F-statistic: 4.219e+04 on 1 and 221 DF,  p-value: < 2.2e-16

对比这次的显著性检验结果和之前结果，T检验，F检验 和 R^2检验，并没有明显的效果提升，结果和我预想的是一样的。所以，通过残差分析和异常点分析，我认为模型是有效的。

6. 模型预测

最后，我们获得了一元线性回归方程的公式，就可以对数据进行预测了。比如，对给定X=x0时，计算出y0=a+b*x0的值，并计算出置信度为1-α的预测区间。

当X=x0,Y=y0时，置信度为1-α的预测区间为

reg5

即

reg6

我们可以用R语言的predict()函数来计算预测值y0，和相应的预测区间。程序算法如下。


> new<-data.frame(x=14040)
> lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95)

# 预测结果
> lm.pred
       fit      lwr      upr
1 14102.09 14093.73 14110.44

当x0=14040时，在预测区间为0.95的概率时，y0的值为 14102，预测区间为[14093.73,14110.44]。

我们通过图形来表示。


> plot(y~x+1)
> abline(lm.ab,col='red')
> points(rep(newX$x,3),y=lm.pred,pch=19,col=c('red','blue','green'))

03

其中，红色点为y0的值，蓝色点为预测区间最小值，绿色点为预测区间最大值。

对于统计模型中最核心部分就在结果解读，本文介绍了一元回归模型的基本的建模过程和模型的详细解读方法。在我们掌握了这种方法以后，就可以更容易地理解和学习 多元回归，非线性回归 等更多的模型，并把这些模型应用到实际的工作中了。下一篇文章R语言解读多元线性回归模型，请继续阅读。

转载请注明出处：
http://blog.fens.me/r-linear-regression/

打赏作者