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

简介: 一元线性回归模型及诊断(原理+实例+代码)

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 最小二乘法回归

  设回归方程为:

image.png

求解方法:

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.公式原理法

  根据公式:

image.png

  运行程序:

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.036 x

4.2 回归方程标准误差

求解方法:

  1.由表2可以看出,回归方程标准误差σ ^ = 0.48 \hat\sigma=0.48σ^=0.48

  2.公式原理法:

image.png

  运行程序:

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

  综上,回归方程标准误差σ ^ = 0.48

4.3 β ^ 0 β ^ 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%的参数区间估计

   image.png

2.公式原理法

image.png

image.png

image.png

  运行程序:

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

   image.png

4.4 x y 的决定系数

求解方法:

image.png

运行程序:

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

  运行结果:

> R_squared 
[1] 0.9004924

  综上,x xxy yy的决定系数为 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 回归系数β 1 \beta_1β1的显著性检验

求解方法:

  1.由4.1节 中 summary()函数结果可以看出,在 0.001 的显著性水平下,回 归系数 β 1 呈显著正相关。

  2.公式原理法

image.png

  运行程序:

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

  结果显示:回归系数 β 1 \beta_1β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.公式法:

   image.png

  运行程序:

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.公式法

  通过回归方程公式:

y^=β1x+β0

  运行程序:

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

  运行结果:

> y0 y0 
[1] 3.703262

  综上,结果显示:下一周签发新保单x 0 x_0x0 =1000 张时需要的加班时间 3.703 小 时。

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

image.png

  运行程序:

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%精确预测区间 y0(2.52,4.89) ;近似预测区间 y0(2.74,4.66) 。

4.11 给出E(y0)的置信度为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(y0) 的置信度为 95%的区间估计:E(y0) (3.28, 4.12) 。


相关文章
|
20天前
|
Python
探索LightGBM:异常值处理与鲁棒建模
探索LightGBM:异常值处理与鲁棒建模【2月更文挑战第2天】
65 0
|
7月前
|
数据可视化 数据挖掘
【因果推断】Day01- 实用计量方法图解与概述
【因果推断】Day01- 实用计量方法图解与概述
128 2
|
20天前
|
机器学习/深度学习 数据采集 监控
机器学习-特征选择:如何使用递归特征消除算法自动筛选出最优特征?
机器学习-特征选择:如何使用递归特征消除算法自动筛选出最优特征?
190 0
|
20天前
|
机器学习/深度学习 算法 数据可视化
Python用KNN(K-近邻)回归、分类、异常值检测预测房价、最优K值选取、误差评估可视化
Python用KNN(K-近邻)回归、分类、异常值检测预测房价、最优K值选取、误差评估可视化
|
20天前
|
机器学习/深度学习 运维 算法
【视频】检测异常值的4种方法和R语言时间序列分解异常检测
【视频】检测异常值的4种方法和R语言时间序列分解异常检测
|
20天前
|
数据可视化 Python
PYTHON贝叶斯推断计算:用BETA先验分布推断概率和可视化案例
PYTHON贝叶斯推断计算:用BETA先验分布推断概率和可视化案例
|
20天前
|
数据可视化
结构方程模型 SEM 多元回归和模型诊断分析学生测试成绩数据与可视化
结构方程模型 SEM 多元回归和模型诊断分析学生测试成绩数据与可视化
|
20天前
|
机器学习/深度学习 数据可视化 算法
支持向量回归SVR拟合、预测回归数据和可视化准确性检查实例
支持向量回归SVR拟合、预测回归数据和可视化准确性检查实例
|
20天前
|
机器学习/深度学习 数据可视化 算法
Python支持向量回归SVR拟合、预测回归数据和可视化准确性检查实例
Python支持向量回归SVR拟合、预测回归数据和可视化准确性检查实例
|
20天前
|
存储 监控 算法
R语言贝叶斯非参数模型:密度估计、非参数化随机效应META分析心肌梗死数据
R语言贝叶斯非参数模型:密度估计、非参数化随机效应META分析心肌梗死数据