一元线性回归模型及诊断（原理+实例+代码）

🍉CSDN小墨&晓末:https://blog.csdn.net/jd1813346972

   个人介绍: 研一｜统计学｜干货分享
         擅长Python、Matlab、R等主流编程软件
         累计十余项国家级比赛奖项，参与研究经费10w、40w级横向

1 目的

  利用最小二乘法构建一保险公司营业部每周加班时间和签发的新保单数目关系的一元线性回归模型，通过数理统计方法对回归方程及回归系数进行显著性检验和经济 学意义解释，并通过拟合结果对未知变量进行预测及置信度为 95%的区间估计

2 数据背景

  一家保险公司十分担心其总公司营业部加班的程度，决定认真调查一 下现状。经过 10 周时间，收集了每周加班时间的数据和签发的新保单数目，x 为 每周签发的新保单数目，y 为每周加班时间（小时），数据见表 1。

表1 保险公司营业部加班及业务情况

3 数据基本情况

  根据实际情况，因变量每周加班时间 y ,程序中对应变量为 overtime。 自变量为每周签发的新保单数目 x ,程序中对应变量为 np_number。，首先绘制该保险公司每周签发的新保单数和每周加班时间（小时）的散点图，见图1。
  运行程序：

1. np_number<-c(825,215,1070,550,480,920,1350,325,670,1215) #读取新保单数目 
2. overtime<-c(3.5,1.0,4.0,2.0,1.0,3.0,4.5,1.5,3.0,5.0) #读取加班时间 
3. Com<-data.frame('新保单数目'=np_number,'加班时间'=overtime) #合并数据 
4. Com #数据展示 
5. plot(Com$新保单数目,Com$加班时间,xlab="新保单数目",ylab="加班时间"，type="p") 
6. #散点图

  运行结果：

图1 原始数据散点图

  然后绘制根据数据拟合该保险公司每周签发的新保单数和每周加班时间（小时）之间的关系，判断 二者之间是否大致呈线性关系，拟合直线图见图2所示。
  运行程序：

1. abline(lm(Com$加班时间~Com$新保单数目)) #拟合直线 

  运行结果：

图2 拟合直线图

  通过图 2 可以看出每周加班时间和新保单数目大致呈线性关系。

4 建模分析

4.1 最小二乘法回归

  设回归方程为：
$$y={\hat{\beta }}_0+{\hat{\beta }}_{1}x$$
求解方法：

1.lsfit 函数

  运行程序：

1. lsfit(Com$新保单数目,Com$加班时间) #最小二乘法求回归方程

  运行结果：

> lsfit(Com$新保单数目,Com$加班时间) 
$coefficients 
 Intercept X 
0.118129074 0.003585132 

2.lm函数

  运行程序：

1. fm=lm(Com$加班时间~Com$新保单数目,data=Com) 
2. summary(fm)

  运行结果：

> fm=lm(Com$加班时间~Com$新保单数目,data=Com) 
> summary(fm) 
 
Call: 
lm(formula = Com$加班时间 ~ Com$新保单数目, data = Com) 
 
Residuals: 
 Min 1Q Median 3Q Max 
-0.83899 -0.33483 0.07842 0.37228 0.52594 
 
Coefficients: 
 Estimate Std. Error t value Pr(>|t|) 
(Intercept) 0.1181291 0.3551477 0.333 0.748 
Com$新保单数目 0.0035851 0.0004214 8.509 2.79e-05 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
 
Residual standard error: 0.48 on 8 degrees of freedom 
Multiple R-squared: 0.9005, Adjusted R-squared: 0.8881 
F-statistic: 72.4 on 1 and 8 DF, p-value: 2.795e-05 

  根据运行结果可以得到回归方程显著性情况（见表2），以及回归系数显著性（见表3）

表2 回归方程显著性表

表3 回归系数显著性表

3.公式原理法

  根据公式：
$${\hat{\beta }}_1=\frac{{\sum }_{i=1}^{10}{x}_i{y}_i}{{\sum }_{i=1}^{10}{x}_i^2-n{\bar{x}}^2}$$

$${\hat{\beta }}_0=\bar{y}-{\hat{\beta }}_1\bar{x}$$

  运行程序：

1. Com1=mean(Com$新保单数目) #平均新增保单数 
2. Com2=mean(Com$加班时间) #平均加班时间 
3. a=(sum(np_number*overtime)-10*Com1*Com2)/(sum(np_number*np_number)-10*Com1*Com
1) 
4. b=Com2-a*Com1 
5. a 
6. b

  运行结果：

> a 
[1] 0.003585132 
> b 
[1] 0.1181291

  综上，回归方程为：$y=0.1181+0.036x$。

4.2 回归方程标准误差

求解方法：

  1.由表2可以看出，回归方程标准误差$\hat{\sigma }=0.48$。

  2.公式原理法：

  根据$\hat{\sigma }$求解公式：
$${\hat{\sigma }}^2=\frac{1}{n-2}\sum _{i=1}^{10}{e}_i^2=\frac{1}{n-2}\sum _{i=1}^{10}(y_i-y_i{)}^2$$

$$\hat{\sigma }=\sqrt{{\hat{\sigma }}^2}$$

  运行程序：

1. overtime1=a*np_number+b #拟合值 
2. var=sum((overtime-overtime1)^2)/(10-2) #方差 
3. var 
4. std=sqrt(var) #标准差 
5. std 

  运行结果：

> var 
[1] 0.2304223 
> std 
[1] 0.4800233

  综上，回归方程标准误差$\hat{\sigma }=0.48$。

4.3 ${\hat{\beta }}_0$和${\hat{\beta }}_1$的置信度为 95%的区间估计。

求解方法：

1.confint（）函数

1. confint(fm) #参数区间估计

  运行结果：

> confint(fm) #参数区间估计 
 2.5 % 97.5 % 
(Intercept) -0.700843004 0.937101152 
Com$新保单数目 0.002613486 0.004556779

  根据结果，可得两个参数置信度为 95%的参数区间估计（见表4）。

表4 置信度为 95%的参数区间估计

  由表4可得${\hat{\beta }}_0\in [-0.700843004,0.937101152]$和${\hat{\beta }}_1\in [0.002613486,0.004556779]$。

2.公式原理法

  由于
$${\hat{\beta }}_1\sim N({\beta }_1,\frac{{\sigma }^2}{L_{x}x})$$
  可得：
$$t=\frac{{\hat{\beta }}_1-{\beta }_1}{\sqrt{{\hat{\sigma }}^2/L_{xx}}}=\frac{({\hat{\beta }}_1-{\beta }_1{)}^2\sqrt{L_{xx}}}{\hat{\sigma }}$$
  服从自由度为n-2（8）的 t 分布，因而：
$$P({\hat{\beta }}_1-t_{\alpha /2}\frac{L_{xx}}{\hat{\sigma }}＜{\beta }_1＜{\hat{\beta }}_1+t_{\alpha /2}\frac{L_{xx}}{\hat{\sigma }})=1-\alpha$$
  即,${\beta }_1$的置信度为|$1-\alpha$的置信区间为:
$$({\hat{\beta }}_1-t_{\alpha /2}\frac{L_{xx}}{\hat{\sigma }},{\hat{\beta }}_1+t_{\alpha /2}\frac{L_{xx}}{\hat{\sigma }})$$
  其中：
$$t_{0.05/2}(8)=t_{0.025}(8)=2.3060$$

$$L_{xx}=\sum _{i=1}^{10}{x}_i^2-n{\bar{x}}^2$$

  由于：
$${\hat{\beta }}_0\sim N({\beta }_0,(\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}}){\sigma }^2)$$
  可得：
$$t=\frac{{\hat{\beta }}_0-{\beta }_0}{\sqrt{(\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}}){\hat{\sigma }}^2}}=\frac{{\hat{\beta }}_0-{\beta }_0}{\hat{\sigma }\sqrt{(\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}})}}$$
  服从自由度为n-2（8）的t分布，因而：
$$P({\hat{\beta }}_0-\frac{{\hat{\beta }}_0-{\beta }_0}{\hat{\sigma }\sqrt{(\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}})}}＜{\beta }_0＜{\hat{\beta }}_0+\frac{{\hat{\beta }}_0-{\beta }_0}{\hat{\sigma }\sqrt{(\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}})}})=1-\alpha$$
  即，${\beta }_0$的$1-\alpha$的置信区间为：
$$({\hat{\beta }}_0-\frac{{\hat{\beta }}_0-{\beta }_0}{\hat{\sigma }\sqrt{(\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}})}},{\hat{\beta }}_0+\frac{{\hat{\beta }}_0-{\beta }_0}{\hat{\sigma }\sqrt{(\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}})}})$$
  其中：
$$t_{0.05/2}(8)=t_{0.025}(8)=2.3060$$

$$L_{xx}=\sum _{i=1}^{10}{x}_i^2-n{\bar{x}}^2$$

  运行程序：

1. Lxx=sum(np_number^2)-10*(Com1^2) #均方离差 
2. c1=a-2.3060*std/sqrt(Lxx) #常数项置信下限 
3. d1=a+std/sqrt(Lxx) #常数项置信上限 
4. c2=b-2.3060*std*sqrt(1/10+Com1^2/Lxx) #参数项置信上限 
5. d2=b+2.3060*std*sqrt(1/10+Com1^2/Lxx) #参数项置信上限 
6. c1 
7. d1 8. c2 
9. d2

  运行结果：

> c1 
[1] 0.002613487 
> d1 
[1] 0.004006488 
> c2 
[1] -0.7008415 
> d2 
[1] 0.9370997

  综上，${\hat{\beta }}_0$和${\hat{\beta }}_1$的置信度为95%的置信区间分别为： ${\hat{\beta }}_0\in [-0.701,0.937]$和${\hat{\beta }}_1\in [0.003,0.004]$。

4.4 $x$和$y$的决定系数

求解方法：

  1.由表2可以看出，$x$和$y$的决定系数为：$R^2=0.9005$

  2.公式法：
$$R^2=\frac{SSR}{SST}=\frac{{\sum }_{i=1}^n({\hat{y}}_i-\bar{y}{)}^2}{{\sum }_{i=1}^n(y_i-\bar{y}{)}^2}$$
运行程序：

1. R_squared=sum((overtime1-Com2)^2)/sum((overtime-Com2)^2) 
2. R_squared 

  运行结果：

> R_squared 
[1] 0.9004924

  综上，$x$和$y$的决定系数为 0.9005。

4.5 方差分析

求解：

  运行程序：

1. anova(fm) #方差分析

  运行结果：

> anova(fm) #方差分析 
Analysis of Variance Table 
 
Response: Com$加班时间 
 Df Sum Sq Mean Sq F value Pr(>F) 
Com$新保单数目 1 16.6816 16.6816 72.396 2.795e-05 *** 
Residuals 8 1.8434 0.2304 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

  方差分析结果见表5所示。

表5 方差分析表

  结果显示：p=0.00002795<0.001,说明在 0.001 的显著性水平下签发的新保单 数 x 对每周加班 时间有显著影响。

4.6 回归系数${\beta }_1$的显著性检验

求解方法：

  1.由4.1节 中 summary（）函数结果可以看出，在 0.001 的显著性水平下，回 归系数 ${\beta }_1$ 呈显著正相关。

  2.公式原理法

  根据t检验统计量：
$$t=\frac{{\hat{\beta }}_1}{\sqrt{\hat{\sigma }/Lxx}}=\frac{{\hat{\beta }}_1\sqrt{Lxx}}{\sqrt{\hat{\sigma }}}$$
  其中：
$${\hat{\sigma }}^2,a=\frac{1}{10-2}\sum _{i-1}^{10}{e}_i^2=\frac{1}{10-2}\sum _{i-1}^{10}(y_i-{\hat{y}}_i){)}^2$$
  运行程序：

1. t<-a*sqrt(Lxx)/std #计算带估参数 t 值 
2. p<-pt(t,8,lower.tail=FALSE) #计算 P（t>t 值）概率 
3. p1<-2*p #计算 P（|t|>|t 值|），即 P 值概率 
4. t 
5. p 
6. p1

  运行结果：

> t 
[1] 8.508575 
> p 
[1] 1.397387e-05 
> p1 
[1] 2.794774e-05

  结果显示：回归系数 ${\beta }_1$ 的显著性检验 p=0.00002795<0.001,说明在 0.001 的显 著性水平下签发的新保单数 x 对每周加班时间有显著影响。

4.7 相关系数的显著性检验

求解方法：

1.cor.test 函数计算相关系数的显著性检验

  运行程序：

1. tr1<-cor.test(np_number,overtime) #x与y的相关系数的显著性检验 
2. tr1

  运行结果：

 Pearson's product-moment correlation 
 
data: np_number and overtime 
t = 8.5086, df = 8, p-value = 2.795e-05 
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval: 
 0.7932921 0.9881624 
sample estimates:  cor 
0.9489428

  结果显示：x 与 y 的相关系数为 0.9489428，表明呈显著正相关，且 p=0.00002795<0.05，通过在 0.001 的显著性水平下的显著性检验。

2.公式法：

  根据 Pearson 计算公式：
$$r=\frac{s_{xy}}{\sqrt{s_X^2{x}_y^2}}=\frac{l_{xx}}{l_{xx}{l}_{xy}}=\frac{\sum (x_i-\bar{x})(y_i-\bar{y})}{\sqrt{(\sum x_i-\bar{x}{)}^2\sum y_i-\bar{y}{)}^2}}$$
  t检验统计量：
$$tr=\frac{r-0}{\frac{1-r^2}{10-2}}$$
  运行程序：

1. lxy<-function(x,y){ 
2. n=length(x); 
3. sum(x*y)-sum(x)*sum(y)/n 
4. } #建立离均差乘积和函数 
5. lxy(np_number,np_number) #x 的离均差平方和 
6. lxy(overtime,overtime) #y 的离均差平方和 
7. lxy(np_number,overtime) #x 与 y 的离均差平方和 
8. r<-lxy(np_number,overtime)/sqrt(lxy(np_number,np_number)*lxy(overtime,overtime
)) #x 与 y 的 Pearson 相关系数 
9. n<-length(np_number) #向量的长度 
10. tr2<-r/sqrt((1-r^2)/(n-2)) 
11. p<-pt(tr2,8,lower.tail=FALSE) #计算 P（t>t 值）概率 
12. p1<-2*p #计算 P（|t|>|t 值|），即 P 值概率 
13. r #相关系数 
14. tr2 #相关系数 t 值 
15. p1 #P 值

  运行结果：

> r #相关系数 
[1] 0.9489428 
> tr2 #相关系数 t 值 
[1] 8.508575 
> p1 #P 值 
[1] 2.794774e-05

  结果表明：x 与 y 的相关系数为 0.9489428，表明呈显著正相关，且p=0.00002795<0.05，通过在 0.001 的显著性水平下的显著性检验。

4.8 对回归分析做残差图并做相应的分析

  运行程序：

1. plot(fm,which=1) #绘制回归方程残差图

  运行结果：

图3 回归方程残差图

  结果分析：通过回归方程残差图可以看出，几乎所有点 e=0 附近波动，且不 呈任何趋势，说明选用的线性回归模型较为合适。

4.9 预测下一周签发新保单 $x_0$ =1000 张时需要的加班时间

求解方法：

1.Predict（）函数

  运行程序：

1. x<-c(825,215,1070,550,480,920,1350,325,670,1215) #读取新保单数目 
2. y<-c(3.5,1.0,4.0,2.0,1.0,3.0,4.5,1.5,3.0,5.0) #读取加班时间 
3. Com<-data.frame(x,y) #合并数据 
4. fm=lm(y~x,data=Com) #最小二乘法求回归方程 
5. x0<-data.frame(x = 1000) #输入新值 1000 并存储为数据框 
6. pre0<-predict(fm,x0) #计算预测值 
7. pre0 #预测值

  运行结果：

> pre0 
 1 
3.703262

2.公式法

  通过回归方程公式：
$$\hat{y}={\beta }_{1}x+{\beta }_0$$
  运行程序：

1. y0 <-a*1000+b #x0=1000 时，y 预测值 
2. y0

  运行结果：

> y0 y0 
[1] 3.703262

  综上，结果显示：下一周签发新保单$x_0$ =1000 张时需要的加班时间 3.703 小 时。

4.10 给出 $y_0$ 的置信度为 95%的精确预测区间和近似预测区间

  求解： $y_0$的置信水平为 0.05 的置信区间精确为：
$$\hat{y}=\pm t_{0.05/2}(8-2)\sqrt{1+h_{00}}\hat{\sigma }$$
  其中：
$$h_{00}=\frac{1}{10}+\frac{(x_0-\bar{x}{)}^2}{Lxx}$$
   $y_0$ 的置信水平为 0.05 的置信区间近似为：
$${\hat{y}}_0\pm 2\hat{\sigma }$$
  运行程序：

1. ypred<-predict(fm,x0,interval = "prediction",level = 0.95) 
2. pre1<-pre0+2*std #计算近似置信上限 
3. pre2<-pre0-2*std #计算近似置信下限 4. ypred #预测值及精确预测区间 
5. pre1 #近似置信上限 
6. pre2 #近似置信上限

  运行结果：

> ypred #预测值及精确预测区间 
 fit lwr upr 
1 3.703262 2.51949 4.887033 
> pre1 #近似置信上限 
 1 
4.663308 
> pre2 #近似置信上限 
 1 
2.743215

  结果见表6所示。

表6 y0置信度为95%的精确及近似区间

  由表 6 得知， $y_0$的置信度为 95%精确预测区间 $y_0\in$(2.52,4.89) ；近似预测区间 $y_0\in$(2.74,4.66) 。

4.11 给出$E(y_0)$的置信度为95%的区间估计

求解：

  运行程序：

> yconf #预测值 E（y0）及精确预测区间 
 fit lwr upr 
1 3.703262 3.283728 4.122795 

  运行结果：

> yconf #预测值 E（y0）及精确预测区间 
 fit lwr upr 
1 3.703262 3.283728 4.122795 

  结果见表7所示：

表7 E(y0）置信度为95%的精确及近似区间

  由表 7得知， $E(y_0)$ 的置信度为 95%的区间估计： $E(y_0)\in$ (3.28, 4.12) 。