2  线性回归

2.1 简单线性回归

\[Y = \beta_0 + \beta_1 X + \epsilon\]

其中，\(\epsilon\) 是随机误差项，通常假设误差项是独立于\(X\)的，\(\beta_0\) 和 \(\beta_1\) 是系数或参数。

我们希望找到\(\beta_0\) 和 \(\beta_1\) 的最优估计，使得预测的直线尽可能接近数据点。测量接近程度的方法有很多，最长用的是残差平方和最小化准则。而最小二乘估计就是满足这个准则的估计方法。

残差平方和（Residual Sum of Squares, RSS）：

\[RSS = \sum_{i=1}^n (y_i - \hat{y}_i)^2\]

其中，\(\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i\) 是预测值，\(y_i\) 是实际值。

通过计算\(\beta_0\) 和 \(\beta_1\)的置信区间（通常用标准误差计算）以及对两者进行假设检验来提供估计的准确性。

一个无偏估计不会系统的高估或低估真实参数

评价模型的准确性：残差标准误（RSE）,\(R^2\)统计量

RSE被认为是对模型失拟的度量，RSE越小，模型拟合越好。

\(R^2\)统计量，也称为决定系数，衡量的是\(Y\)的变异中能被\(X\)解释的部分所占的比例。\(R^2\)越接近于1，模型拟合越好。

2.2 多元线性回归

\[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p + \epsilon\]

响应变量与预测变量之间是否有关系？也就是说，所有回归系数是否都为0？即\(\beta_0 = beta_1 = \cdots = beta_p = 0\)。

\[H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0\]

\(H_1: 至少有一个\beta_j 不全为0\)。

如果\(H_0\)为真，那么\(Y\)与\(X_i\)之间没有线性关系。这里可以对系数进行F检验。

如果\(p=100\)，那么与每个变量对应的\(p\)值，有约5%将碰巧小于0.05，而F检验则不会有这个问题

变量选择：向前选择、向后选择、混合选择

加入新的变量必然会使\(R^2\)增大，但可能不会使\(RSE\)减小，这种情况容易出现过拟合。

线性模型拓展

标准的线性回归模型假设预测变量和响应变量的关系是可加和线性的，如果预测变量和响应变量之间的关系不是可加的，我们就需要放宽这种假设。

去除可加性假设 \(Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_1 X_2 + \epsilon\), 其中，\(X_1 X_2\) 是交互项。在\(\beta_3\)的\(p值\)很小时，也就是\(X_1\)和\(X_2\)之间的交互作用重要时，\(\beta_1\)和\(\beta_2\)的\(p\)值即使很大，也不应该去除这两个变量（实验分层原则）。

多项式回归：\(Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \cdots + \beta_p X^p + \epsilon\)。

2.3 线性模型的潜在问题

1. 数据非线性：如果数据的分布不是线性关系，可以考虑使用变量转换，如\(logX, \sqrt{X}, X^2\)等，或者使用非线性模型。可以使用残差图来考察数据分布是否线性，理想情况下，残差图中显示不出明显的规律。残差图，也就是残差与\(x_i\)的散点图，残差计算为\(e_i=y_i - \hat{y}_i\)。在多元线性回归中，残差图是残差（\(e_i\)）与\(\hat{y_i}\)的散点图。
2. 误差项自相关：通常出现在时间序列的数据中。如果误差项是自相关的，那么残差图会显示出某种模式。
3. 误差项方差非恒定：方差非恒定通常出现在数据存在异方差性的情况下。异方差性是指误差项的方差随着预测变量\(X\)的值变化。如果存在异方差性，那么残差图会呈现漏斗形。可以使用凹函数对\(Y\)进行转换，如\(\sqrt{Y}和log Y\)。
4. 离群点：给定\(x_i\)的值，\(y_i\)异常的点，离群点会显著影响模型的拟合效果。可以用学生化残差（studentized residuals）来检测离群点，|学生化残差| > 3 则可能是离群点。
5. 高杠杆点：\(x_i\)的值异常。可以用杠杆统计量来检测高杠杆点， 杠杆统计量大于\(3(p+1)/n\)则可能是高杠杆点。
6. 共线性：当两个或多个预测变量高度相关时，称为共线性。如果存在共线性，那么系数估计值会不稳定，置信区间会很大，F检验会不显著，可能导致检验效能降低，无法拒绝\(H_0\)。可以使用方差膨胀因子（VIF）来检测共线性。VIF > 5 则可能存在共线性。解决办法可以是删除其中一个预测变量，或者合并两个共线的预测变量。

2.4 KNN回归

当X和Y之间存在非线性关系时，可以使用KNN回归。KNN回归的基本思想是，给定一个新的观测值，找到离它最近的K个观测值，然后使用这K个观测值的Y值来预测新观测值的Y值。

当X和Y的真是关系为线性时，KNN回归的效果略逊于线性回归，而为非线性时，KNN回归的效果会大大优于线性回归。

当X较多时（高维），KNN往往不如线性回归。X仅有少量观测时，参数方法更优秀。

2.5 R语言实现

导入需要用到的包

library(MASS)
library(ISLR)
library(car)

Loading required package: carData

查看数据集列名

names(Boston)

 [1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"    
 [8] "dis"     "rad"     "tax"     "ptratio" "black"   "lstat"   "medv"   

2.5.1 简单线性回归

lm.fit = lm(medv~lstat, data=Boston)

# to see the detailed information about the model.
summary(lm.fit)

Call:
lm(formula = medv ~ lstat, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.168  -3.990  -1.318   2.034  24.500 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 34.55384    0.56263   61.41   <2e-16 ***
lstat       -0.95005    0.03873  -24.53   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.216 on 504 degrees of freedom
Multiple R-squared:  0.5441,    Adjusted R-squared:  0.5432 
F-statistic: 601.6 on 1 and 504 DF,  p-value: < 2.2e-16

系数估计的置信区间

# confidence interval for the coefficient estimates
confint(lm.fit)

                2.5 %     97.5 %
(Intercept) 33.448457 35.6592247
lstat       -1.026148 -0.8739505

系数估计值

# extractor functions to acess the information stored in lm.fit
coef(lm.fit)

(Intercept)       lstat 
 34.5538409  -0.9500494 

给定lstat值，预测medv时，predict()函数计算置信区间和预测区间

# prediction interval
predict(lm.fit, data.frame(lstat=c(5, 10, 15)), interval="prediction")

       fit       lwr      upr
1 29.80359 17.565675 42.04151
2 25.05335 12.827626 37.27907
3 20.30310  8.077742 32.52846

#| label: confidence interval

# confidence interval
predict(lm.fit, data.frame(lstat=c(5, 10, 15)), interval="confidence")

       fit      lwr      upr
1 29.80359 29.00741 30.59978
2 25.05335 24.47413 25.63256
3 20.30310 19.73159 20.87461

medv和lstat的散点图以及最小二乘回归直线

attach(Boston)
plot(lstat,medv)
detach()
abline(lm.fit,lwd=3,col="red")

诊断图

par(mfrow=c(2,2))
plot(lm.fit)

残差对拟合值的散点图

# compute the residuals from a linear regression fit
plot(predict(lm.fit),residuals(lm.fit))

学生化残差对拟合值的散点图

# compute the studentized residuals from a linear regression fit
plot(predict(lm.fit),rstudent(lm.fit))

杠杆统计量

# compute the leverage statistics from a linear regression fit
plot(hatvalues(lm.fit))

# which observation has the largest leverage statistic
which.max(hatvalues(lm.fit))

375 
375 

2.5.2 多元线性回归

lm.fit_lstat_age <- lm(medv ~ lstat+age,data=Boston)
summary(lm.fit_lstat_age)

Call:
lm(formula = medv ~ lstat + age, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.981  -3.978  -1.283   1.968  23.158 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 33.22276    0.73085  45.458  < 2e-16 ***
lstat       -1.03207    0.04819 -21.416  < 2e-16 ***
age          0.03454    0.01223   2.826  0.00491 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.173 on 503 degrees of freedom
Multiple R-squared:  0.5513,    Adjusted R-squared:  0.5495 
F-statistic:   309 on 2 and 503 DF,  p-value: < 2.2e-16

对所有预测变量进行线性回归

# use all the variables to perform a regresssion
lm.fit_all <- lm(medv~.,data=Boston)
summary(lm.fit_all)

Call:
lm(formula = medv ~ ., data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.595  -2.730  -0.518   1.777  26.199 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
zn           4.642e-02  1.373e-02   3.382 0.000778 ***
indus        2.056e-02  6.150e-02   0.334 0.738288    
chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
age          6.922e-04  1.321e-02   0.052 0.958229    
dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
black        9.312e-03  2.686e-03   3.467 0.000573 ***
lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406,    Adjusted R-squared:  0.7338 
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

vif()函数计算方差膨胀因子

# compute variance inflation factors
vif(lm.fit_all)

    crim       zn    indus     chas      nox       rm      age      dis 
1.792192 2.298758 3.991596 1.073995 4.393720 1.933744 3.100826 3.955945 
     rad      tax  ptratio    black    lstat 
7.484496 9.008554 1.799084 1.348521 2.941491 

age变量具有很高的p值，这里使用所有变量，除了age，进行线性回归

# use all the variables to perform a regresssion except age
lm.fit_all_except_age <- lm(medv~.-age,data=Boston)
summary(lm.fit_all_except_age)

Call:
lm(formula = medv ~ . - age, data = Boston)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.6054  -2.7313  -0.5188   1.7601  26.2243 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  36.436927   5.080119   7.172 2.72e-12 ***
crim         -0.108006   0.032832  -3.290 0.001075 ** 
zn            0.046334   0.013613   3.404 0.000719 ***
indus         0.020562   0.061433   0.335 0.737989    
chas          2.689026   0.859598   3.128 0.001863 ** 
nox         -17.713540   3.679308  -4.814 1.97e-06 ***
rm            3.814394   0.408480   9.338  < 2e-16 ***
dis          -1.478612   0.190611  -7.757 5.03e-14 ***
rad           0.305786   0.066089   4.627 4.75e-06 ***
tax          -0.012329   0.003755  -3.283 0.001099 ** 
ptratio      -0.952211   0.130294  -7.308 1.10e-12 ***
black         0.009321   0.002678   3.481 0.000544 ***
lstat        -0.523852   0.047625 -10.999  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.74 on 493 degrees of freedom
Multiple R-squared:  0.7406,    Adjusted R-squared:  0.7343 
F-statistic: 117.3 on 12 and 493 DF,  p-value: < 2.2e-16

2.5.3 交互项

lstat, age 和交互项 lstat × age 同时作为预测变量

# The syntax lstat:black tells R to include an interaction term between
# lstat and black. The syntax lstat*age simultaneously includes lstat,age,
# and the interaction term lstat × age as predictors
lm.fit_ia_lstat_age <- lm(medv~lstat*age,data=Boston)
summary(lm.fit_ia_lstat_age)

Call:
lm(formula = medv ~ lstat * age, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.806  -4.045  -1.333   2.085  27.552 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 36.0885359  1.4698355  24.553  < 2e-16 ***
lstat       -1.3921168  0.1674555  -8.313 8.78e-16 ***
age         -0.0007209  0.0198792  -0.036   0.9711    
lstat:age    0.0041560  0.0018518   2.244   0.0252 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.149 on 502 degrees of freedom
Multiple R-squared:  0.5557,    Adjusted R-squared:  0.5531 
F-statistic: 209.3 on 3 and 502 DF,  p-value: < 2.2e-16

2.5.4 预测变量的非线性变换

建立medv对lstat和lstat平方的线性回归模型

# we can create a predictor X 2 using I(X^2)
lm.fit2 <- lm(medv~lstat+I(lstat^2),data=Boston)
summary(lm.fit2)

Call:
lm(formula = medv ~ lstat + I(lstat^2), data = Boston)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.2834  -3.8313  -0.5295   2.3095  25.4148 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 42.862007   0.872084   49.15   <2e-16 ***
lstat       -2.332821   0.123803  -18.84   <2e-16 ***
I(lstat^2)   0.043547   0.003745   11.63   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.524 on 503 degrees of freedom
Multiple R-squared:  0.6407,    Adjusted R-squared:  0.6393 
F-statistic: 448.5 on 2 and 503 DF,  p-value: < 2.2e-16

使用anova()函数进一步量化二次拟合优于线性拟合的程度

# further quantify the extent to which the quadratic fit is superior to the linear fit.
anova(lm.fit,lm.fit2)

Analysis of Variance Table

Model 1: medv ~ lstat
Model 2: medv ~ lstat + I(lstat^2)
  Res.Df   RSS Df Sum of Sq     F    Pr(>F)    
1    504 19472                                 
2    503 15347  1    4125.1 135.2 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

这里的模型1代表只包含lstat的线性模型，模型2代表包含lstat和lstat平方的二次模型。从F统计量来看，二次模型显著优于线性模型。

par(mfrow=c(2,2))
plot(lm.fit2)

当模型中包含交互项时，残差中可识别的规律很少

建立5阶多项式拟合

# create the polynomial within lm()
lm.fit5 <- lm(medv~poly(lstat,5),data=Boston)
summary(lm.fit5)

Call:
lm(formula = medv ~ poly(lstat, 5), data = Boston)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.5433  -3.1039  -0.7052   2.0844  27.1153 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       22.5328     0.2318  97.197  < 2e-16 ***
poly(lstat, 5)1 -152.4595     5.2148 -29.236  < 2e-16 ***
poly(lstat, 5)2   64.2272     5.2148  12.316  < 2e-16 ***
poly(lstat, 5)3  -27.0511     5.2148  -5.187 3.10e-07 ***
poly(lstat, 5)4   25.4517     5.2148   4.881 1.42e-06 ***
poly(lstat, 5)5  -19.2524     5.2148  -3.692 0.000247 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.215 on 500 degrees of freedom
Multiple R-squared:  0.6817,    Adjusted R-squared:  0.6785 
F-statistic: 214.2 on 5 and 500 DF,  p-value: < 2.2e-16

对数变换

# create the log transformation within lm()
lm.fit_log <- lm(medv~log(rm),data=Boston)
summary(lm.fit_log)

Call:
lm(formula = medv ~ log(rm), data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.487  -2.875  -0.104   2.837  39.816 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -76.488      5.028  -15.21   <2e-16 ***
log(rm)       54.055      2.739   19.73   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.915 on 504 degrees of freedom
Multiple R-squared:  0.4358,    Adjusted R-squared:  0.4347 
F-statistic: 389.3 on 1 and 504 DF,  p-value: < 2.2e-16

2.5.5 定性预测变量

lm.fit_qp = lm ( Sales~.+Income:Advertising+Price:Age,data = Carseats )
summary(lm.fit_qp)

Call:
lm(formula = Sales ~ . + Income:Advertising + Price:Age, data = Carseats)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9208 -0.7503  0.0177  0.6754  3.3413 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         6.5755654  1.0087470   6.519 2.22e-10 ***
CompPrice           0.0929371  0.0041183  22.567  < 2e-16 ***
Income              0.0108940  0.0026044   4.183 3.57e-05 ***
Advertising         0.0702462  0.0226091   3.107 0.002030 ** 
Population          0.0001592  0.0003679   0.433 0.665330    
Price              -0.1008064  0.0074399 -13.549  < 2e-16 ***
ShelveLocGood       4.8486762  0.1528378  31.724  < 2e-16 ***
ShelveLocMedium     1.9532620  0.1257682  15.531  < 2e-16 ***
Age                -0.0579466  0.0159506  -3.633 0.000318 ***
Education          -0.0208525  0.0196131  -1.063 0.288361    
UrbanYes            0.1401597  0.1124019   1.247 0.213171    
USYes              -0.1575571  0.1489234  -1.058 0.290729    
Income:Advertising  0.0007510  0.0002784   2.698 0.007290 ** 
Price:Age           0.0001068  0.0001333   0.801 0.423812    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.011 on 386 degrees of freedom
Multiple R-squared:  0.8761,    Adjusted R-squared:  0.8719 
F-statistic:   210 on 13 and 386 DF,  p-value: < 2.2e-16

返回哑变量值

# returns the coding that R uses for the dummy variables
attach(Carseats)
contrasts(ShelveLoc)

       Good Medium
Bad       0      0
Good      1      0
Medium    0      1