一元线性回归

1 介绍

  本小节主要介绍一元线性回归，并使用R语言去分析某联合循环电厂中发电量与锅炉温度关系。介绍一元线性回归时主要涉及到以下几个部分：

• 一元线性回归的简介及要点
• 一元线性回归模型的参数估计
• 一元线性回归方程的诊断

了解完我们将要讲解的主要内容之后，接下来正式开始本小节的内容。

2 一元线性回归简介及要点

  在讲一元线性回归之前，先说一下什么是回归分析。回归分析就是研究\(xy\)相关性的分析，比如说研究税收和财政收入的关系、研究身高和体重的关系、研究房价和面积、城区的关系等等。回归分析包括一元线性回归、多元线性回归、多项式回归、逻辑回归、定序回归等等，这些回归都有相同点和不同点。

  我们要讲的一元线性回归模型可以说是最简单的一种回归模型，它就是通过回归分析研究两变量之间的依存关系，将变量区分出来自变量和因变量，并研究确定自变量和因变量之间的具体关系的方程式。

一元线性回归有两个要点必须得注意：

• 第一，因变量\(y\)必须是连续型变量，或者是近似连续的变量。当然，也不用严格上的连续，至少从理论上讲是可正可负的。

• 第二，它研究的是\(xy\)线性相关关系的模型。非线性相关关系不属于线性回归的研究范畴。

3 回归模型各参数简介

  通过阅读一元线性回归的简介和要点，我们知道了一元线性回归就涉及到了两个变量：因变量\(y\)、自变量\(x\)。那么下面就给大家展示一元线性回归方程表达式： \[y=\beta_0+\beta_1x+\mu\] 我们可以用下图去形象的理解上述表达式：

  大家看到表达式之后会发现里面还多出了\(\beta_0\)、\(\beta_1\)和\(\mu\)。下面就给大家去解释表达式中各个系数的含义：

  # \(x\)：自变量
  # \(y\)：因变变
  # \(\beta_0\)：俗称截距项，是当\(x=0\)时，\(y\)的值
  # \(\beta_1\)：俗称斜率，是整个回归的核心
  # \(\mu\)：干扰项，是其他因素对\(y\)的影响的大小，是统计学的大智慧

  从上面的表达式中不难看出，要想得到具体的这个回归模型，我们得估计出\(\beta_0\)和\(\beta_1\)。常见的系数估计方法有：极大似然估计、最小二乘法、梯度下降法、最小一法、M-估计等等，那么就有同学问了：哪个是最好的估计方法呢？答：这些估计方法各有千秋，没有任何一种是最优的。下面我们使用的是最常用的最小二乘法。

4 参数估计之最小二乘法

  使用最小二乘法拟合回归系数，并且得到一个最佳的一元回归线的前提就是：估计出的\(\beta_0\)和\(\beta_1\)得使上图中的所有\(\mu\)(残差平方和)的平方和最小。在这里可能有同学会问：为什么要用残差平方和，而不用残差和呢？答案很简单：因为残差有正有负，如果单纯的对所有残差求和，那么正负抵消的情况会导致估计系数出现较大的偏差。下面给大家演示一下最小二乘法的推导过程：

1、拟合方程\((y\sim{x})\)

我们将要拟合的方程一元线性方程是：\(y=\beta_0+\beta_1x\)

2、样本点\((x,y)\)

我们用来拟合方程的样本点为：\((x_1, y_1),(x_2, y_2)...(x_n, y_n)\)

3、残差(Residual)

残差就是拟合值与真实值的距离，我们在此假设残差为\(\mu_i\)，则\(\mu_i=y_i-(\beta_0+\beta_1x_i)\)

4、残差平方和(SSR)

\(D=\sum\limits_{i=1}^{n}\mu_i^2=\sum\limits_{i=1}^{n}(y_i-\beta_0-\beta_1x_i)^2\)

5、估计系数\(\beta_0\)和\(\beta_1\)

 (1)对\(\beta_0\)求偏导

  令 \(G=y_i-\beta_0-\beta_1x_i\)，则 \(D=\sum\limits_{i=1}^nG^2\)

\[\begin{aligned} \frac{\partial D}{\partial \beta_0} &=2GG' \\ &=\sum\limits_{i=1}^{n}2(y_i-\beta_0-\beta_1x_i)(-1)\\ &=-2\sum\limits_{i=1}^{n}(y_i-\beta_0-\beta_1x_i)\\ &=-2(\sum\limits_{i=1}^{n}y_i-\sum\limits_{i=1}^{n}a-b\sum\limits_{i=1}^{n}x_i)\\ &=-2(n\bar{y}-n\beta_0-n\beta_1\bar{x}) \end{aligned}\]

 (2)对\(\beta_1\)求偏导 \[\begin{aligned} \frac{\partial D}{\partial \beta_1} &=2GG'\\ &=\sum\limits_{i=1}^{n}2(y_i-\beta_0-\beta_1x_i)(-x_i)\\ &=-2\sum\limits_{i=1}^{n}(x_iy_i-\beta_0x_i-\beta_1x_i^2)\\ &=-2(\sum\limits_{i=1}^{n}x_iy_i-\sum\limits_{i=1}^{n}\beta_0x_i-\sum\limits_{i=1}^{n}\beta_1x_i^2)\\ &=-2(\sum\limits_{i=1}^{n}x_iy_i-n\beta_0\bar{x}-\beta_1\sum\limits_{i=1}^{n}x_i^2) \end{aligned}\]

 (3)令 \(\frac{\partial D}{\partial \beta_0}=0\)

\[\begin{alignat}{2} &-2(n\bar{y}-n\beta_0-n\beta_1\bar{x})=0\\ &\Rightarrow\color{red}{\beta_0=\bar{y}-\beta_1\bar{x}} \end{alignat}\]

 (4)令 \(\frac{\partial D}{\partial \beta_1}=0\) \[\begin{alignat}{2} &-2(\sum\limits_{i=1}^{n}x_iy_i-n\color{red}{\beta_0}\bar{x}-\beta_1\sum\limits_{i=1}^{n}x_i^2)=0\\ &\Rightarrow\sum\limits_{i=1}^nx_1y_i-n\color{red}{\beta_0}\bar{x}-\beta_1\sum\limits_{i=1}^nx_i^2=0\\ &\Rightarrow\sum\limits_{i=1}^nx_iy_i-n\bar{x}(\color{red}{\bar{y}-\beta_1\bar{x}})-\beta_1\sum\limits_{i=1}^nx_i^2=0\\ &\Rightarrow\sum\limits_{i=1}^{n}x_iy_i-n\bar{x}\bar{y}=\beta_1(\sum\limits_{i=1}^{n}x_i^2-n\bar{x}^2)\\ &\Rightarrow\beta_1=\frac{\sum\limits_{i=1}^{n}x_iy_i-n\bar{x}\bar{y}}{\sum\limits_{i=1}^{n}x_i^2-n\bar{x}^2} \end{alignat}\]

 又因为 \[\begin{alignat}{2} \sum\limits_{i=1}^n(x_i-\bar{x})(y_i-\bar{y}) &= \sum\limits_{i=1}^nx_iy_i-n\bar{x}\bar{y}\\ \sum\limits_{i=1}^n(x_i-\bar{x}) &= \sum\limits_{i=1}^nx_i^2-n\bar{x}^2 \end{alignat}\]  所以\(\color{red}{\beta_1=\frac{\sum\limits_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum\limits_{i=1}^n(x_i-\bar{x})^2}}\)

  以上就是推导的全过程，看起来一大堆，但实际上依旧是浅显易懂的。最小二乘法的部分到此结束，开始下面的模型拟合。

5 案例实战

  本部分主要利用某联合循环发电厂的数据去拟合一条回归线。该数据一共有两列，分别是自变量锅炉温度（AT）和因变量每小时发电量（EP）。温度单位是摄氏度，发电量单位兆瓦。案例主要涉及到以下5小部分：

• 导入外部数据
• 数据简单探索
• 拟合模型
• 模型解释
• 模型诊断

5.1 导入数据

  案例使用的数据的格式是csv，所以下面我们使用R中read.csv()函数读取此数据。该函数里面有很多的参数，一般情况下只需要一个主要的参数就可以方便的导入数据，那就是指定数据所在位置的那个参数file。

## 导入数据
mydata <- read.csv(file = "E:/BeiDa_task/CCPP.csv")

5.2 简单探索数据

数据的简单探索主要包括：

• 查看数据的前5行和后5行，一般使用head()和tail()函数，为了输出的美观可以嵌套kable()函数
• 数据集的基本结构，一般使用str()函数
• 数据集的主要描述性统计量，一般使用summary()函数

5.2.1 mydata

library(knitr)          # 加载knitr包，为了使用kable函数
kable(head(mydata, 5))  # 查看数据前5行

AT	PE
8.34	480.48
23.64	445.75
29.74	438.76
19.07	453.09
11.80	464.43

kable(tail(mydata, 5))  # 查看数据后5行

	AT	PE
181	32.17	430.64
182	26.23	447.16
183	11.69	471.72
184	14.73	470.41
185	12.36	468.37

5.2.2 str()

## 探索数据集的基本结构
str(mydata)

## 'data.frame':    185 obs. of  2 variables:
##  $ AT: num  8.34 23.64 29.74 19.07 11.8 ...
##  $ PE: num  480 446 439 453 464 ...

简单的一条命令，返回了数据的大量信息。

• 语句data.frame告诉我们数据集是数据框类型
• 语句185 obs. of 2 variabels告诉我们数据一共包含185个观测、2个变量
• 左下角的AT、PE和num告诉我们每个变量的变量名和所对应的类型（num表示数值型）
• 右下角的那些数值是每个变量所对应的部分数值

5.2.3 summary()

  此函数可以返回各个变量（数值型、整型）的主要描述性统计量：最小值、第四分之一位数、中位数、均值、第四分之三位数和最大值。这几个主要描述性统计量可以帮助我们很快的了解各变量的分布情况。

注：该函数还可以返回各个变量包含的缺失值个数。

## 查看变量的主要描述性统计量
summary(mydata)

##        AT              PE       
##  Min.   : 5.12   Min.   :426.2  
##  1st Qu.:13.97   1st Qu.:438.5  
##  Median :20.01   Median :453.1  
##  Mean   :19.78   Mean   :454.1  
##  3rd Qu.:26.23   3rd Qu.:468.4  
##  Max.   :34.20   Max.   :489.8

5.3 拟合一元回归模型

  在R中常用的回归函数是lm()，下面我们也使用这个函数数据进行构建一元线性回归模型。该函数在构建一元线性回归模型时，主要使用以下两个参数：

• formula：一个回归方程的表达式
• data：表达式中的变量所属的数据集（即用来拟合模型的数据）

## 拟合回归模型
fit <- lm(formula = PE ~ AT, data = mydata)

5.4 模型解释

  我们在上个操作针对数据建立了一个一元回归模型，然后可以使用summary()函数去查看模型的详细信息。

## 查看模型的详细信息
summary(fit)

## 
## Call:
## lm(formula = PE ~ AT, data = mydata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.0969  -3.5455   0.5229   3.4218  12.5321 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 495.88937    1.02383  484.35   <2e-16 ***
## AT           -2.11226    0.04828  -43.75   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.012 on 183 degrees of freedom
## Multiple R-squared:  0.9127, Adjusted R-squared:  0.9123 
## F-statistic:  1914 on 1 and 183 DF,  p-value: < 2.2e-16

上面返回的模型的细节，大致可以分为三个部分来看：

• Residuals（残差）：该部分提供了预测误差的主要描述性统计量。

  从返回的残差值（真实值减去预测值）可知，拟合误差最大是12.5321，说明该模型至少对一个观测的发电量少预测了12.53兆瓦左右，相对于实际值的区间[426.2, 489.8]而言，误差还是比较小的。另一方面，误差值的50%落在了1Q（第四分之一位数）和3Q（第四分之三位数）内，所以大部分的预测值在超过真实值3.55与低于真实值3.42之间。根据残差来看，整体的上的拟合还是很不错的。

• Confficients：该部分提供了线性模型的各参数的系数的详细信息。主要关注以下两点：

  Estimate（估计值）：各个变量所对应的估计值，比如AT的系数是-2.11226。
  Pr（>|t|）（P值）：通常与预设的0.05做对比来判定自变量的显著性（检验的原假设是，该洗漱是否显著为0，若P值小于0.05，则拒绝原假设，即对应的变量显著不为0）。后面的星号与之对应，Signif.codes部分提供了一种度量方式，即真实系数有多大可能是0。比如3颗星表示显著性水平是0，意味着该变量极不可能与因变量无关。

• Multiple R-squared（多元R方值）和Adjusted R-squared（调整后的R方值）：

  提供一种度量模型性能的方式，即从整体上，模型能多大程度解释因变量的值。它类似于相关系数，越接近于1，模型的解释数据的性能越好。我们这个模型的调整后的R方值超过了0.9，可以说是已经是相当好了。

讲到这里，大家还记得我在之前说过\(\beta_1\)是整个回归的核心吗？接下来就给大家讲一下为什么它是核心：

  就拿咱们拟合出的模型而言，通过返回的模型信息可以得到此一元回归模型是： \[y=-2.11226AT+495.88937\]   这里面的\(\beta_1\)是-2.11226，说明锅炉温度与发电量负相关，锅炉温度每增长\(1^{\circ}\)C，发电量会降低2.11226兆瓦。我们根据这个结论就可以采取某些措施去改变锅炉的温度，从而降低其对发电量的消极影响，提高发电量。

5.5 模型诊断

  一元线性回归的模型诊断主要包括残差是否符合正态分布、数据异常值检测和是否存在异常差性，一般情况下使用plot()函数就可以返回上述我们的所需信息。

library(lmtest)        # 加载程序包，为了使用bptest函数

par(mfrow = c(2, 2))   # 把绘图区域分割为“2*2”四份
plot(fit, whic = 1:4)  # 绘制残差图、QQ图、位置尺度图和Cook距离图

bptest(fit)            # 判断异方差是否存在

## 
##  studentized Breusch-Pagan test
## 
## data:  fit
## BP = 1.3294, df = 1, p-value = 0.2489

• 左上角的残差图可以看出，所有的残差点并不是均匀分布在中间的虚线两侧，也并不存在特别明显的趋势，不能精准的判断是否存在异方差。但是经过Breush Pagan检验，发现P值是0.2489>0.05，所以可以说明检验拒绝原假设（残差的方差是恒定的），即不存在异常差。
• 右上角的QQ图得出，大部分的点均在对角线上，说明标准化的残差符合正态分布。
• 右下角的Cook距离相对较均衡，表明数据中不存在异常值。

综述所述，无论是数据表现还是模型表现都不错，整体上可以说是一个不错的回归模型。

6 总结

  一元线性回归的介绍也就到此结束了，以上所有内容看起来还是不少，但是并不是一元线性回归的全部，仅仅是一部分而已。

简单总结一下一元线性回归的几个注意点：

• 确定变量之间是线性关系
• 确定因变量的连续性
• 构建模型之后要进行模型诊断
• 根据诊断结果采取相应措施