By

文章目录
  1. 回归
    1. 1、回归的多面性
    2. 2、最小二乘回归
      1. 2.1 使用lm()拟合回归模型
      2. 2.2 简单线性回归
      3. 2.3 多项式回归
      4. 2.4 多元线性回归
      5. 2.5 有交互项的多元线性回归
    3. 3、回归诊断
      1. 3.1 标准方法
      2. 3.2 改进的方法
      3. 3.3 线性模型假设的综合验证
      4. 3.4 多重共线性

回归

1、回归的多面性

  • 简单线性：用一个量化的解释变量预测一个量化的响应变量
  • 多项式：用一个量化的解释变量预测一个量化的响应变量，模型的关系是n阶多项式
  • 多元线性：用两个或多个量化的解释变量预测一个量化的响应变量
  • 多变量：用一个或多个解释变量预测多个响应变量
  • Logistic：用一个或多个解释变量预测一个类别型响应变量
  • 泊松：用一个或多个解释变量预测一个代表频数的响应变量
  • Cox比例风险：用一个或多个解释变量预测一个事件（死亡、失败或旧病复发）发生的时间
  • 时间序列：对误差项相关的时间序列数据建模
  • 非线性：用一个或多个量化的解释变量预测一个量化的响应变量，不过模型是非线性的
  • 非参数：用一个或多个量化的解释变量预测一个量化的响应变量，模型的形式源自数据形式，不事先设定
  • 稳健：用一个或多个量化的解释变量预测一个量化的响应变量，能抵御强影响点的干扰

    2、最小二乘回归


    为了能够恰当地解释OLS模型的系数，数据必须满足以下统计假设。
  • 正态性: 对于固定的自变量值，因变量值成正态分布。
  • 独立性: Yi值之间相互独立。
  • 线性: 因变量与自变量之间为线性相关。
  • 同方差性: 因变量的方差不随自变量的水平不同而变化。也可称作不变方差，但是说同方差性感觉上更犀利。

2.1 使用lm()拟合回归模型

1
myfit <- lm(formula, data)

常用的参数和解释

对拟合线性模型非常有用的其他函数:

2.2 简单线性回归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
> #OLS回归
> 
> #简单线性回归 女性身高体重数据
> 
> fit<-lm(weight~height, data=women)
> 
> summary(fit)
Call:
lm(formula = weight ~ height, data = women)
Residuals:
    Min      1Q  Median      3Q     Max 
-1.7333 -1.1333 -0.3833  0.7417  3.1167 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
height        3.45000    0.09114   37.85 1.09e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared:  0.991,	Adjusted R-squared:  0.9903 
F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14
> 
> women$weight
 [1] 115 117 120 123 126 129 132 135 139 142 146 150 154 159 164
> 
> fitted(fit)
       1        2        3        4        5        6        7        8        9       10       11 
112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 133.2833 136.7333 140.1833 143.6333 147.0833 
      12       13       14       15 
150.5333 153.9833 157.4333 160.8833 
> 
> residuals(fit)
          1           2           3           4           5           6           7           8 
 2.41666667  0.96666667  0.51666667  0.06666667 -0.38333333 -0.83333333 -1.28333333 -1.73333333 
          9          10          11          12          13          14          15 
-1.18333333 -1.63333333 -1.08333333 -0.53333333  0.01666667  1.56666667  3.11666667 
> 
> plot(women)
> abline(fit)
# 此模型的得到的方程式是：weight= -87.51667+3.45*height

2.3 多项式回归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
> #增加一个二次项提高精度
> 
> fit2<-lm(weight~height+I(height^2),data=women)
> 
> #此处I()的含义是height^2表示为一个常规的R表达式，因为^在在表达式中有特殊含义
> 
> summary(fit2)
Call:
lm(formula = weight ~ height + I(height^2), data = women)
Residuals:
     Min       1Q   Median       3Q      Max 
-0.50941 -0.29611 -0.00941  0.28615  0.59706 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 261.87818   25.19677  10.393 2.36e-07 ***
height       -7.34832    0.77769  -9.449 6.58e-07 ***
I(height^2)   0.08306    0.00598  13.891 9.32e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3841 on 12 degrees of freedom
Multiple R-squared:  0.9995,	Adjusted R-squared:  0.9994 
F-statistic: 1.139e+04 on 2 and 12 DF,  p-value: < 2.2e-16
> 
> women$weight
 [1] 115 117 120 123 126 129 132 135 139 142 146 150 154 159 164
> 
> fitted(fit2)
       1        2        3        4        5        6        7        8        9       10       11       12       13 
115.1029 117.4731 120.0094 122.7118 125.5804 128.6151 131.8159 135.1828 138.7159 142.4151 146.2804 150.3118 154.5094 
      14       15 
158.8731 163.4029 
> 
> residuals(fit2)
           1            2            3            4            5            6            7            8            9 
-0.102941176 -0.473109244 -0.009405301  0.288170653  0.419618617  0.384938591  0.184130575 -0.182805430  0.284130575 
          10           11           12           13           14           15 
-0.415061409 -0.280381383 -0.311829347 -0.509405301  0.126890756  0.597058824 
> 
> plot(women)
> lines(women$height,fitted(fit2))

1
2
3
4
5
6
7
8
9
10
11
#car包中的scatterplot方便的绘制二元关系图
library(car)
scatterplot(weight~height,
            data=women,
            spread=FALSE,#spread=FALSE选项删除了残差正负均方根在平滑
            #曲线上的展开和非对称信息
            pch=19,
            main="Women Age 30-39",
            xlab="Height (inches)",
            ylab="Weight (lbs.)"
            )

2.4 多元线性回归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
> #多项式回归
> #观察谋杀率 人口 文盲率 平均收入 结霜天数关系
> states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])
> #观察二元关系
> cor(states)
               Murder Population Illiteracy     Income      Frost
Murder      1.0000000  0.3436428  0.7029752 -0.2300776 -0.5388834
Population  0.3436428  1.0000000  0.1076224  0.2082276 -0.3321525
Illiteracy  0.7029752  0.1076224  1.0000000 -0.4370752 -0.6719470
Income     -0.2300776  0.2082276 -0.4370752  1.0000000  0.2262822
Frost      -0.5388834 -0.3321525 -0.6719470  0.2262822  1.0000000
> 
> library(car)
> 
> #大致观察数据图示关系(见下图)
> scatterplotMatrix(states, spread = FALSE, lty = 2, main="Scatter Plot Matrix")
> 
> #回归
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
> 
> summary(fit)
Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
    data = states)
Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
Population  2.237e-04  9.052e-05   2.471   0.0173 *  
Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
Income      6.442e-05  6.837e-04   0.094   0.9253    
Frost       5.813e-04  1.005e-02   0.058   0.9541    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared:  0.567,	Adjusted R-squared:  0.5285 
F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08
> 
> #当预测变量不止一个时，回归系数的含义为，一个预测变量增加一个单位，其他预测变量保
> #持不变时，因变量将要增加的数量。
> #例如本例中，文盲率的回归系数为4.14，表示控制人口、收入和温度不变时，文盲率上升1%，
> #谋杀率将会上升4.14%，它的系数在p<0.001的水平下显著不为0。

2.5 有交互项的多元线性回归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
> #有交互项的回归  y=x-y+xy
> fit <- lm(mpg ~ hp + wt +hp:wt, data=mtcars)
> summary(fit)
Call:
lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars)
Residuals:
    Min      1Q  Median      3Q     Max 
-3.0632 -1.6491 -0.7362  1.4211  4.5513 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 49.80842    3.60516  13.816 5.01e-14 ***
hp          -0.12010    0.02470  -4.863 4.04e-05 ***
wt          -8.21662    1.26971  -6.471 5.20e-07 ***
hp:wt        0.02785    0.00742   3.753 0.000811 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.153 on 28 degrees of freedom
Multiple R-squared:  0.8848,	Adjusted R-squared:  0.8724 
F-statistic: 71.66 on 3 and 28 DF,  p-value: 2.981e-13
> 
> #若两个预测
> #变量的交互项显著，说明响应变量与其中一个预测变量的关系依赖于另外一个预测变量的水平。
> #因此此例说明，每加仑汽油行驶英里数与汽车马力的关系依车重不同而不同。
> 
> #通过effects包中的effect()函数，你可以用图形展示交互项的结果。
> library(effects)
> plot(effect("hp:wt", fit, xlevels=list(wt=c(2.2,3.2,4.2))),multiline=TRUE)


从图中可以很清晰地看出，随着车重的增加，马力与每加仑汽油行驶英里数的关系减弱了。
当wt = 4.2时，直线几乎是水平的，表明随着hp的增加， mpg不会发生改变。

3、回归诊断

  • 对回归结果的信念，都只建立在你的数据满足统计假设的前提之上。

    3.1 标准方法

    1
    2
    3
    4
    5
    #回归诊断
    fit<- lm(weight~height,data=women)
    #将plot()函数绘制的四幅图形组合在一个大的2×2的图中。
    par(mfrow=c(2,2))
    plot(fit)

  • 正态性：当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正
    态分布。正态Q-Q图（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概率
    图。若满足正态假设，那么图上的点应该落在呈45度角的直线上；若不是如此，那么就
    违反了正态性的假设。

  • 独立性：由样本的数据收集去确保

  • 线性：若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。
    换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“残差图与拟合图”
    （Residuals vs Fitted，左上）中可以清楚的看到一个曲线关系，这暗示着你可能需要对回
    归模型加上一个二次项。

  • 同方差性 ：若满足不变方差假设，那么在位置尺度图（Scale-Location Graph，左下）中，
    水平线周围的点应该随机分布。该图似乎满足此假设。

  • 最后一幅“残差与杠杆图”（Residuals vs Leverage，右下）提供了你可能关注的单个观测点的信息。从图形可以鉴别出离群点、高杠杆值点和强影响点。

看下二次拟合的图：

1
2
3
4
5
6
7
fit<- lm(weight~height+I(height^2),data=women)
#将plot()函数绘制的四幅图形组合在一个大的2×2的图中。
par(mfrow=c(2,2))
plot(fit)
#这第二组图表明多项式回归拟合效果比较理想，基本符合了线性假设、残差正态性（除了
#观测点13）和同方差性（残差方差不变）。观测点15看起来像是强影响点（根据是它有较大的
#Cook距离值），删除它将会影响参数的估计。事实上，删除观测点13和15，模型会拟合得会更好。

3.2 改进的方法

car包提供多种函数用以评价回归模型：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
qqPlot() 分位数比较图
durbinWatsonTest() 对误差自相关性做Durbin-Watson检验
crPlots() 成分与残差图
ncvTest() 对非恒定的误差方差做得分检验
spreadLevelPlot() 分散水平检验
outlierTest() Bonferroni离群点检验
avPlots() 添加的变量图形
inluencePlot() 回归影响图
scatterplot() 增强的散点图
scatterplotMatrix() 增强的散点图矩阵
vif() 方差膨胀因子

  • 正态性：qqPlot()
    1
    2
    3
    4
    5
    6
    7
    #它画出了在n-p-1个自由度的t分布下的学生化残差（studentized residual，也称学生化删除残差或折叠
    #化残差）图形，其中n是样本大小， p是回归参数的数目（包括截距项）。
    library(car)
    fit<-lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
    par(mfrow=c(1,1))
    qqPlot(fit,labels=row.names(states),id.method = "identify",simulate=TRUE, main="Q-Q Plot")
    #当simulate=TRUE时， 95%的置信区间将会用参数自助法

另一种和可视化误差的方法

1
2
3
4
5
6
7
8
9
10
residplot<-function(fit,nbreaks=10){
  z<-rstudent(fit)
  hist(z,breaks=nbreaks,freq = FALSE, xlab="Studentized Residual",main="Distribution of Errors")
  rug(jitter(z), col = "brown")
  curve(dnorm(x,mean=mean(z),sd=sd(z)),
        add=TRUE, col="blue",lwd=2)
  lines(density(z)$x,density(z)$y,col="red",lwd=2,lty=2)
  legend("topright",legend = c("Normal Curve","Kernel Density Curve"),lty=1:2,col=c("blue","red"),cex=.7)
}
residplot(fit)

  • 独立性

    1
    2
    3
    4
    5
    6
    7
    8
    9
    > #误差的独立性
    > durbinWatsonTest(fit)
     lag Autocorrelation D-W Statistic p-value
       1      -0.2006929      2.317691   0.258
     Alternative hypothesis: rho != 0
    > #p值不显著说明无相关性，误差项之间独立，滞后项（lag=1）表明数据集中
    > #每个数据都是与其后一个数据进行比较的。该检验适用于时间独立的数据，对于非聚集型的数据
    > #并不适用。注意， durbinWatsonTest()函数使用自助法（参见第12章）来导出p值，如果添加
    > #了选项simulate=FALSE，则每次运行测试时获得的结果都将略有不同。

  • 线性

    1
    2
    3
    4
    crPlots(fit)
    #若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分，
    #那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代
    #替X），或用其他回归变体形式而不是线性回归。

  • 同方差性
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    > #线性
    > crPlots(fit)
    > #若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分，
    > #那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代
    > #替X），或用其他回归变体形式而不是线性回归。
    > 
    > #同方差性
    > #ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。
    > #若检验显著，则说明存在异方差性（误差方差不恒定）
    > ncvTest(fit)
    Non-constant Variance Score Test 
    Variance formula: ~ fitted.values 
    Chisquare = 1.746514    Df = 1     p = 0.1863156 
    > p=0.18不显著接受零假设
    > #spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。
    > spreadLevelPlot(fit)
    Suggested power transformation:  1.209626 
    #你将会看到一个非水平的曲线。代码结果建议幂次变换（suggested power transformation）的含义
    #是，经过p次幂（Y^p）变换，非恒定的误差方差将会平稳。
    #建议幂次转换为0.5，在回归等式中用 Y 代替Y，可能会使模型满足同方差性。若建议幂次为0，
    #则使用对数变换。对于当前例子，异方差性很不明显，因此建议幂次接近1（不需要进行变换）。

3.3 线性模型假设的综合验证

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
> library(gvlma)
> gvmodel<-gvlma(fit)
> summary(gvmodel)
Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
    data = states)
Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
Population  2.237e-04  9.052e-05   2.471   0.0173 *  
Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
Income      6.442e-05  6.837e-04   0.094   0.9253    
Frost       5.813e-04  1.005e-02   0.058   0.9541    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared:  0.567,	Adjusted R-squared:  0.5285 
F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 
Call:
 gvlma(x = fit) 
                    Value p-value                Decision
Global Stat        2.7728  0.5965 Assumptions acceptable.
Skewness           1.5374  0.2150 Assumptions acceptable.
Kurtosis           0.6376  0.4246 Assumptions acceptable.
Link Function      0.1154  0.7341 Assumptions acceptable.
Heteroscedasticity 0.4824  0.4873 Assumptions acceptable.
#从输出项（Global Stat中的文字栏）我们可以看到数据满足OLS回归模型所有的统计假设
#（p=0.5965597） 。若Decision下的文字表明违反了假设条件（比如p<0.05） ，你可以使用前几节讨论
#的方法来判断哪些假设没有被满足。

3.4 多重共线性

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
> #多重共线性
> #例如 握力~年龄+DOB(Date Of Birth)
> #F检验显著， p<0.001。但是当你观察DOB和年龄的回归系数时，却
> #发现它们都不显著（也就是说无法证明它们与握力相关）
> #原因是DOB与年龄在四舍五入后相关性极大。回归系数测量的是当其他预测变量不变时，某
> #个预测变量对响应变量的影响。那么此处就相当于假定年龄不变，然后测量握力与年龄的关系，
> #这种问题就称作多重共线性（multicollinearity）。它会导致模型参数的置信区间过大，使单个系数
> #解释起来很困难。
> #使用VIF度量（Variance Inflation Factor，方差膨胀因子），
> #VIF的平方根表示变量回归参数的置信区间能膨胀为与模型无关的预测变量的程度（因此而得名）
> library(car)
> vif(fit)
Population Illiteracy     Income      Frost 
  1.245282   2.165848   1.345822   2.082547 
> sqrt(vif(fit))>2
Population Illiteracy     Income      Frost 
     FALSE      FALSE      FALSE      FALSE
文章目录
  1. 回归
    1. 1、回归的多面性
    2. 2、最小二乘回归
      1. 2.1 使用lm()拟合回归模型
      2. 2.2 简单线性回归
      3. 2.3 多项式回归
      4. 2.4 多元线性回归
      5. 2.5 有交互项的多元线性回归
    3. 3、回归诊断
      1. 3.1 标准方法
      2. 3.2 改进的方法
      3. 3.3 线性模型假设的综合验证
      4. 3.4 多重共线性