给女朋友写的生统资料_Part16

2019-06-14 本文已影响0人城管大队哈队长

多元线性回归

多元线性回归的方程写为：
$Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdot\cdot\cdot + \beta_p X_p + \epsilon$
其中 $X_j$ 代表第 $j$ 个预测变量， $\beta_j$ 是对应的模型参数。 $\beta_j$ 可以解释为在所有其他预测变量保持不变的情况下， $X_j$ 增加一个单位对 $Y$ 产生的平均影响。（这段话其实很重要，我觉得有可能要考的）

和前面的一元线性回归类似，我们需要用 $\hat{\beta}_0,\hat{\beta}_1,……,\hat{\beta}_p$ 来估计 $\beta_0,\beta_1,……,\beta_p$ 。对于给定的 $\hat{\beta}_0,\hat{\beta}_1,……,\hat{\beta}_p$ ，我们的回归模型就变成了
$\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \cdot\cdot\cdot + \hat{\beta}_p x_p$
与一元线性回归类似，多元线性回归的参数也是用最小二乘法进行估计的，也是使得残差平方和（RSS）最小。计算方法的话，似乎是矩阵代数最为方便。一来我也不会╮（╯＿╰）╭，二来老师PPT也没有放，所以这里就不算了。

然后我们还是按照一元线性回归的步骤，首先去检验整个方程的显著性。

模型的显著性检验（模型的意义）

假设检验第一步首先是建立零假设和备则假设（这里跟一元就不一样了）。
$H_{0}: \beta_1 = \beta_2 =... = \beta_p = 0$

$H_{1}: 至少一个\beta_j不等于0$

统计检验的F统计量为
$F = \frac{(TSS-RSS)/p}{RSS/(n-p-1)} \sim F(p,n-p-1)$
然后计算p-value就可以了。

模型参数的显著性

还记得我们在简单线性回归中提到的，模型参数的准确性么。我们利用了 t 统计量，即可以用来做模型参数的置信区间，也可以用来做模型参数的统计检验。但在一元线性回归中，t 检验跟 F 检验几乎可以说是等价的，因为只有一个参数，你检验了那一个参数的显著性，就是等于检验了整个方程的显著性。但在这里，t 检验就非常有用了。其可以检验每个模型参数的统计显著性，而F检验做的，是整个模型的统计检验性。

假设检验：
$H0:\beta_j=0,其他模型参数不为0$

$H1:\beta_j \ne0，其他模型参数不为0$

然后统计量为：
$t = \frac{\hat{\beta_1} - 0}{SE(\hat{\beta_1})} \sim t(n-p-1)$
然后计算 p-value 就行了。

模型的比较

我们先前讲过了 F 统计量是检验整体的统计量，t 统计量可以检验单参数的统计量，其实我们还可以检验部分参数的统计量（这个考试的时候，应该是用anova函数来比较含有不同参数的模型）。

在前面 F 统计量的假设检验是假设所有系数都是零。但这里我们想要检验的是规模为q的特定子集为0。
$H0:\beta_{p-q+1}=\beta_{p-q+2}……=\beta_{p}=0，其他剩余的p-q个不等于0$

$H1:规模为q的特定子集不都为0，其他剩余的p-q个不等于0$

为方便起见，我们用除了q个变量之外的所有变量建立第二个模型。假设该模型的残差平方和为 $RSS_0$ ，相应的F统计量为
$F=\frac{(RSS_0-RSS)/q}{RSS/(n-p-1)} \sim F(q,n-p-1)$
这部分我觉得不太需要掌握。唯一需要掌握的是利用anova函数比较去掉变量前后的模型，这个后面会给出例子。

评估模型的准确性

这部分还是用残差标准误和 $R^2$ 统计量暂时应该就可以了。残差标准误稍微注意下
$RSE = \sqrt{\frac{RSS}{n-p-1}}$
参考文章

《统计学习导论》

Ｒ语言实现

R 语言实现的函数就是一个 lm 函数，用法见《R语言实战》第二版162页。我这里放个截图

15_1.png

15_2.png

我这里用的是我们生统第六次作业的第二题，我将这个数据集命名为test3（数据集已放入part0部分）。

某农场通过试验取得早稻收获量与春季降雨量和春季温度的数据如下：

# 准备数据
> test3 <- read.table("rawdata/homework-6.2.txt",header = T)
> colnames(test3) <- c("weight","rainfall","tempature")
> test3
  weight rainfall tempature
1   2250       25         6
2   3450       33         8
3   4500       45        10
4   6750      105        13
5   7200      110        14
6   7500      115        16
7   8250      120        17

建立早稻收获量对春季降雨量和春季温度的二元线性回归方程。

用的是lm函数，用summary函数展现结果

> test3_lm <- lm(weight ~ ., data = test3)
> summary(test3_lm)

Call:
lm(formula = weight ~ ., data = test3)

Residuals:
       1        2        3        4        5        6        7 
-275.101   90.464  216.483  140.280  150.676 -316.599   -6.203 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   -0.591    505.004  -0.001   0.9991  
rainfall      22.387      9.601   2.332   0.0801 .
tempature    327.672     98.798   3.317   0.0295 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 261.4 on 4 degrees of freedom
Multiple R-squared:  0.9913,  Adjusted R-squared:  0.987 
F-statistic: 228.4 on 2 and 4 DF,  p-value: 7.532e-05

写出方程
$y= -0.591+22.387x1+327.672x2$

检验模型的显著性（模型的意义）**

第一步 提出假设，设置显著性水平：
$H_{0}: \beta_1 = \beta_2 = 0$

$H1:不都为0$

显著性水平 $\alpha$ 设置为0.05

第二步 计算检验统计量 F

F统计量为
$F = \frac{(TSS-RSS)/p}{RSS/(n-p-1)} \sim F(p,n-p-1)$
第三步 做出决策。按 R 输出结果为 228.4，显著性水平 P=7.532e-05 < 0.05。所以拒绝 $H_0$ ，即收获量y与降雨量和温度之间的线性关系显著。
判断每个自变量对因变量的影响是否都是显著的

需要对各回归系数 $\beta_j$ 分别进行 t 检验，假设检验和t统计量见前面的模型参数的显著性。根据前面线性回归输出的结果，可以看到降雨量和温度的回归系数对应的显著性水平分别为0.0801和0.0295。只有温度对应的显著性水平小于0.05通过检验，这表明影响收获量的自变量中，只有温度对收获量的影响显著，而降雨量对收获量的影响不显著。

各回归系数的置信区间（根据前面截图的函数）

> confint(test3_lm)
                   2.5 %     97.5 %
(Intercept) -1402.707516 1401.52552
rainfall       -4.268921   49.04184
tempature      53.364699  601.97873

检验两个模型的区别

因为先前我们发现降雨量对收获量的影响不显著，所以我们剔除到降雨量，再重新构建一个模型。然后拿新模型和旧模型比较，看去掉的那个变量是否重要

# 重新构建模型
> test3_lm_new <- lm(weight ~ tempature, data = test3)
> summary(test3_lm_new)

Call:
lm(formula = weight ~ tempature, data = test3)

Residuals:
   1    2    3    4    5    6    7 
-150  -50 -100  500  400 -400 -200 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -900.00     447.82   -2.01    0.101    
tempature     550.00      35.56   15.47 2.05e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 359.2 on 5 degrees of freedom
Multiple R-squared:  0.9795,  Adjusted R-squared:  0.9754 
F-statistic: 239.2 on 1 and 5 DF,  p-value: 2.052e-05

# 利用anova函数比较
> anova(test3_lm,test3_lm_new)
Analysis of Variance Table

Model 1: weight ~ rainfall + tempature
Model 2: weight ~ tempature
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1      4 273385                              
2      5 645000 -1   -371615 5.4372 0.08009 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

稍微要注意的一点就是，在后面逻辑斯蒂回归的（广义线性模型）比较中，anova要设置test="Chisq"。

其余一些结果
- 残差标准误（似乎也有叫剩余标准差的）：261.4
- $R^2$ 统计量：0.9913。
  
  这里你会发现还有一个 $adjust\ R^2$ 。因为实际上只要你的变量加的足够多，你的 $R^2$ 就肯定会增加，哪怕你的变量跟你的Y值是无关的。所以就会引入adjust，来适当调整 $R^2$ 。我觉得考试时候两个值可以都写。

置信区间与预测区间

这部分我觉得大家不用看，我只是写一下，考试肯定也不会考。

先提下可约和不可约误差。

我们利用统计模型，对Y做了预测，得到预测值 $\hat Y$ 。预测值的精确度取决于两个方面，一方面是可约误差，另一方面是不可约误差。其中可约误差是来自于我们所构建的模型。如果我们构建的模型不是那么好，可约误差就会大一点。但我们可以通过优化我们的模型，比如增加一些变量或者把线性变成非线性等等，来降低这种误差。然而，尽管我们可以得到一个精确的模型，我们还是不能减少不可约误差。因为不可约误差是跟误差项，即我们在最一开始的时候写过
$Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdot\cdot\cdot + \beta_p X_p + \epsilon$
里面的 $\epsilon$ 有关的。按照定义， $\epsilon$ 是不能用X去预测的。所以我们是无法去减少这部分误差的。那么其实最后的误差来源就表示为
$E(Y-\hat{Y})^2=E[f(X)+\epsilon-\hat f(X)]^2$

$=[f(X)-\hat f(X)]^2+Var(\epsilon)$

我为啥要提上面这个呢。因为其实在针对线性回归有一个函数叫predict，里面有个参数叫interval。大家可以看下。

## S3 method for class 'lm'
predict(object, newdata, se.fit = FALSE, scale = NULL, df = Inf,
        interval = c("none", "confidence", "prediction"),
        level = 0.95, type = c("response", "terms"),
        terms = NULL, na.action = na.pass,
        pred.var = res.var/weights, weights = 1, ...)

书里面一般写的是prediction。prediction其实就是在构建区间的时候考虑了误差项的影响，而confidence则没有考虑误差项的影响。这就会导致预测区间比置信区间要宽很多。

顺便提一下， $y_0=E(y_0)+\epsilon$ 。 $E(y_0)$ 就是我们通常用线性回归预测出来的值。而confidence构建的就是 $E(y_0)$ 的区间，而prediction构建的就是 $y_0$ 的区间。

# 利用先前的结果，看下区别。可以看到预测区间更宽一点。
> predict(test3_lm,test3,interval = "confidence")
       fit      lwr      upr
1 2525.101 1992.745 3057.457
2 3359.536 2928.396 3790.676
3 4283.517 3795.736 4771.298
4 6609.720 6096.023 7123.417
5 7049.324 6620.305 7478.342
6 7816.599 7407.001 8226.198
7 8256.203 7748.618 8763.789
> predict(test3_lm,test3,interval = "prediction")
       fit      lwr      upr
1 2525.101 1624.957 3425.245
2 3359.536 2515.298 4203.774
3 4283.517 3408.996 5158.038
4 6609.720 5720.483 7498.956
5 7049.324 6206.167 7892.481
6 7816.599 6983.156 8650.043
7 8256.203 7370.483 9141.923

参考文献

《统计学习导论-基于R应用》
Y叔的统计笔记