完整详细的回归分析实例R语言实现（含数据代码）

先做声明，代码中表示的只是按照理论计算的过程，并非R语言的内置函数,大家可以用R语言内置的回归函数做对比,帮助理解学习。本书引用书籍为何晓群的《应用回归分析》

问题2.15

完整代码如下

setwd("E:/AllClass/junior1/RegressionAnalysis/unit2")#设置文件路径 #保留原路径setwd("C:/User/10854/documents") #以下利用理论方法使用一元回归模型 #导入数据 data<-read.csv("2-7.csv")#书本2.15，表数据2-7 x<-data[,1] y<-data[,2] n<-length(x) split.screen(c(1,3)) screen(1) plot(x,y,pch=16) title(main="数据散点图") #求均值与回归变量lxx，lyy，lxy meanx<-mean(x) meany<-mean(y) lxx<-sum((x-meanx)^2) lyy<-sum((y-meany)^2) lxy<-sum((x-meanx)*(y-meany)) #回归系数估计 beta_1<-lxy/lxx#beta_1 beta_0<-meany-beta_1*meanx#beta_0 screen(2) plot(x,y,pch=16) points(x,beta_0+beta_1*x,type="l") title(main="回归图") #预测值与平方和 y_hat<-beta_0+beta_1*x sse<-sum((y_hat-y)^2)#残差平方和 ssr<-sum((y_hat-meany)^2)#回归平方和 sst<-ssr+sse#总离差平方和 #回归误差ε的方差sigma估计 sigma_hat<-sqrt(1/(n-2)*sse) #对bet_0、beta_1的95%区间估计 alpha<-0.05 #beta_0,beta_1的分布标准差 sd.beta_0<-sqrt((1/n+(meanx^2)/lxx))*sigma_hat sd.beta_1<-sqrt(sigma_hat^2/lxx) beta_1l<-beta_1-qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信下限 beta_1u<-beta_1+qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信上限 beta_0l<-beta_0-qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信下限 beta_0u<-beta_0+qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信上限 #remark #qt是求出置信度1-α对应的统计量值t（1-α） #dt是求出统计量对应的置信度值，即p值(这里用不上t分布) #dt返回概率密度，pt返回分布函数，qt返回分位数函数，rt生成随机数。 #qf\df都是同理，对应的是F分布 #计算xy决定系数 R<-ssr/sst #回归方程的显著性检验 #法一方差分析F检验 f<-(ssr/1)/(sse/(n-2)) p1<-pf(f,1,n-2)#F为统计量、1为第一个自由度，n-2为第二个自由度 #法二回归系数的beta_1的t检验 t1<-beta_1/sd.beta_1#t统计量 p2<-pt(t1,n-2) #法三相关系数r的t检验 r<-lxy/(sqrt(lxx*lyy)) t2<-sqrt(n-2)*r/sqrt(1-r^2); p3<-pt(t2,n-2) #p值 #残差图 screen(3) e<-y_hat-y #残差 n<-length(e) sigma_u<-seq(2*sigma_hat,2*sigma_hat,length.out=n) #残差2σ原则 sigma_l<-seq(-2*sigma_hat,-2*sigma_hat,length.out=n) plot(x,e,pch=16,ylim=c(5,-5)) points(x,sigma_u,type="l") #画2σ上下界 points(x,sigma_l,type="l") title(main="残差图") #预测广告费用为1000万元时，销售达多少 x0<-1000 y0<-beta_0+beta_1*x0 #因变量新值得95%置信区间 h00<-1/n+((x0-meanx)^2)/lxx y0_l<-y0-qt(1-alpha/2,n-2)*sqrt(1+h00)*sigma_hat#预测新值得置信下限 y0_u<-y0+qt(1-alpha/2,n-2)*sqrt(1+h00)*sigma_hat#预测新值得置信上限 #近似置信区间 y0_l_ <- y0-2*sigma_hat y0_u_ <- y0+2*sigma_hat #因变量新值得平均值的95%置信区间 y0_l_E<-y0-qt(1-alpha/2,n-2)*sqrt(h00)*sigma_hat#预测新值均值得置信下限 y0_u_E<-y0+qt(1-alpha/2,n-2)*sqrt(h00)*sigma_hat#预测新值均值得置信上限 #以下利用R函数回归 

(1)画散点图

1.1问题求解

1.1.1输入

#导入数据 data<-read.csv("2-7.csv")#书本2.15，表数据2-7 x<-data[,1] y<-data[,2] n<-length(x) # split.screen(c(1,3))  # screen(1) plot(x,y,pch = 16) title(main="数据散点图") 

1.1.2输出

(2)$x$与$y$之间是否大致呈线性关系

由第一问的散点图，大致正相关，并且呈线性关系

(3)用最小二乘估计求回归方程

3.1问题分析

$$L_{xx}={\Sigma }_{i=1}^n(x_i-\bar{x}{)}^2$$

$$L_{yy}={\Sigma }_{i=1}^n(y_i-\bar{y}{)}^2$$

$$L_{xy}={\Sigma }_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})$$

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

3.2问题求解

3.2.1输入

#求均值与回归变量lxx，lyy，lxy meanx<-mean(x) meany<-mean(y) lxx<-sum((x-meanx)^2) lyy<-sum((y-meany)^2) lxy<-sum((x-meanx)*(y-meany)) #回归系数估计 beta_1<-lxy/lxx#beta_1 beta_0<-meany-beta_1*meanx#beta_0 # screen(2) plot(x,y,pch=16) #散点 points(x,beta_0+beta_1*x,type="l") #回归线 title(main="回归图") 

3.2.2输出

得出各参数

回归图：

(4)求回归标准误差 $\hat{\sigma }$

4.1问题分析

4.2问题求解

4.2.1输入

sse<-sum((y_hat-y)^2)#残差平方和 #回归误差ε的方差sigma估计 sigma_hat<-sqrt(1/(n-2)*sse) 

4.2.2输出

$SSE=1.843,\hat{\sigma }=0.48$

(5)给出${\hat{\beta }}_0$与${\hat{\beta }}_1$的置信度为95%的区间估计

5.1问题分析

${\hat{\beta }}_1$与${\hat{\beta }}_0$的分布为
$${\hat{\beta }}_1\sim N\left({\beta }_1,\frac{{\sigma }^2}{L_{xx}}\right),{\hat{\beta }}_0\sim N\left({\beta }_0,\left[\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}}\right]{\sigma }^2\right)$$
由${\hat{\beta }}_1$分布构造了服从自由度为$n-2$的$t$分布统计量
$$t=\frac{({\hat{\beta }}_1-{\beta }_1)\sqrt{L_{xx}}}{\hat{\sigma }}$$
因而
$$P\left(\left|\frac{({\hat{\beta }}_1-{\beta }_1)\sqrt{L_{x}x}}{\hat{\sigma }}\right|>t_{\alpha /2}\left(n-2\right)\right)$$
得到${\beta }_1$的置信度为$1-\alpha$的置信区间为($\alpha =0.05$)
$$\left({\hat{\beta }}_1-t_{\alpha /2}\frac{\hat{\sigma }}{\sqrt{L_{xx}}},{\hat{\beta }}_1+t_{\alpha /2}\frac{\hat{\sigma }}{\sqrt{L_{xx}}}\right)$$
对${\hat{\beta }}_0$同理得置信度为$1-\alpha$的置信区间为($\alpha =0.05$)
$$\left({\hat{\beta }}_1-t_{\alpha /2}\sqrt{\left[\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}}\right]}\hat{\sigma },{\hat{\beta }}_1+t_{\alpha /2}\sqrt{\left[\frac{1}{n}+\frac{{\bar{x}}^2}{L_{xx}}\right]}\hat{\sigma }\right)$$

5.2问题求解

5.2.1输入

alpha<-0.05 #置信度 #beta_0,beta_1的分布标准差 sd.beta_0<-sqrt((1/n+(meanx^2)/lxx))*sigma_hat sd.beta_1<-sqrt(sigma_hat^2/lxx) beta_1l<-beta_1-qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信下限 beta_1u<-beta_1+qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信上限 beta_0l<-beta_0-qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信下限 beta_0u<-beta_0+qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信上限 

5.2.2输出

得到${\beta }_1$的置信度为0.05的置信区间为$\left[0.0026,0.0046\right]$

得到${\beta }_1$的置信度为0.05的置信区间为[-0.7,0.937]

(6) $x$与$y$的决定系数

6.1问题分析

$$SSR=\Sigma (\hat{y}-\bar{y}{)}^2$$

$$SST=SSR+SSE$$

6.2问题求解

6.2.1输入

sse<-sum((y_hat-y)^2)#残差平方和 ssr<-sum((y_hat-meany)^2)#回归平方和 sst<-ssr+sse#总离差平方和 #计算xy决定系数 R<-ssr/sst 

6.2.2输出

$SSE=1.84,SSR=16.68,SST=18.525,r^2=0.9$

(7)对回归方程做方差分析

7.1问题求解

7.1.1输入

f<-(ssr/1)/(sse/(n-2)) p1<-pf(f,1,n-2)#F为统计量、1为第一个自由度，n-2为第二个自由度 #1-p1为p值 

7.2.2输出

一元线性回归方差分析表如下：

(8)做回归系数${\beta }_1$的显著性检验

8.1问题分析

${\hat{\beta }}_1$的分布为
$${\hat{\beta }}_1\sim N\left({\beta }_1,\frac{{\sigma }^2}{L_{xx}}\right)$$
假设检验原假设为$H_0:{\beta }_1=0$ ,构造检验统计量$t$服从自由度为 $n-2$ 的 $t$ 分布
$$t=\frac{{\hat{\beta }}_1\sqrt{L_{xx}}}{\hat{\sigma }}$$
并计算对应的 $p$ 值

8.2问题求解

8.2.1输入

#回归系数的beta_1的t检验 t1<-beta_1/sd.beta_1 #t统计量 p2<-pt(t1,n-2) #（1-p2）为p值 

8.2.2输出

$t=8.50,p=1.40\ast 10^{-5}$

(9)做相关系数的显著性检验

9.1问题分析

相关系数为$r$
$$r=\frac{L_{xy}}{\sqrt{L_{xx}{L}_{yy}}}={\hat{\beta }}_1\sqrt{\frac{L_{xx}}{L_{yy}}}$$
构造检验统计量
$$t=\frac{\sqrt{n-2}r}{\sqrt{1-r^2}}\sim t(n-2)$$
当$|t|>t_{\alpha /2}(n-2)$ 时，认为简单回归系数显著不为零

9.2问题求解

9.2.1输入

#相关系数r的t检验 r<-lxy/(sqrt(lxx*lyy)) t2<-sqrt(n-2)*r/sqrt(1-r^2); p3<-pt(t2,n-2) #（1-p3）为p值 

9.2.2输出

$t=8.50,p=1.40\ast 10^{-5}$

(10)对回归方程作残差图并做相应的分析

10.1问题求解

10.1.1输入

#残差图 # screen(3) e<-y_hat-y #残差 n<-length(e) sigma_u<-seq(2*sigma_hat,2*sigma_hat,length.out=n) #残差2σ原则 sigma_l<-seq(-2*sigma_hat,-2*sigma_hat,length.out=n) plot(x,e,ylim=c(5,-5)) points(x,sigma_u,type="l") #画2σ上下界 points(x,sigma_l,type="l") title(main="残差图") 

10.1.2输出

10.2问题分析

计算出残差后，及自变量x为横轴，以残差为总纵轴画散点图得到残差图，从残差图上看出，产茶时围绕着 $e=0$ 随机波动，并且波动范围在方差估计$\hat{\sigma }$的两倍。所以可以判定模型的基本假定是满足的。

(11)该公司预计下一周签发新保单$x_0=1000$张，需要的加班时间时多少？

11.1问题求解

11.1.1输入

x0<-1000 y0<-beta_0+beta_1*x0 

11.1.2输出

$y=3.7$

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

12.1问题分析

可以计算得到 ${\hat{y}}_0$ 的分布为
$${\hat{y}}_0\sim N\left({\beta }_0+{\beta }_1{x}_0,(\frac{1}{n}+\frac{(x_0-\bar{x}{)}^2}{L_{xx}}){\sigma }^2\right)$$
所以枢轴量为
$$y_0-{\hat{y}}_0\sim N\left(0,(1+h_{00}){\sigma }^2\right),h_{00}=(\frac{1}{n}+\frac{(x_0-\bar{x}{)}^2}{L_{xx}})$$
统计量为
$$t=\frac{y_0-\widehat{y_0}}{\sqrt{1+h_{00}}\hat{\sigma }}\sim t(n-2)$$
得精确置信区间为
$$\widehat{y_0}\pm t_{\alpha /2}(n-2)\sqrt{1+h_{00}}\hat{\sigma }$$
近似预测区间为
$${\hat{y}}_0\pm 2\hat{\sigma }$$

12.2问题求解

12.2.1输入

#因变量新值95%置信区间 h00<-1/n+((x0-meanx)^2)/lxx y0_l<-y0-qt(1-alpha/2,n-2)*sqrt(1+h00)*sigma_hat#预测新值得置信下限 y0_u<-y0+qt(1-alpha/2,n-2)*sqrt(1+h00)*sigma_hat#预测新值得置信上限 

12.2.2输出

精确置信区间为[2.519,4.887]

近似预测区间为[2.743,4.663]

所以两区间比较接近，可以用近似区间估计

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

13.1问题分析

根据 ${\hat{y}}_0$ 构造枢轴量（含$E(y_0)$）
$$\widehat{y_0}-E（y_0）\sim N\left(0,h_{00}{\sigma }^2\right),h_{00}=(\frac{1}{n}+\frac{(x_0-\bar{x}{)}^2}{L_{xx}})$$
统计量为
$$t=\frac{y_0-\widehat{y_0}}{\sqrt{h_{00}}\hat{\sigma }}\sim t(n-2)$$
得置信水平95%得精确置信区间为
$$\widehat{y_0}\pm t_{\alpha /2}(n-2)\sqrt{h_{00}}\hat{\sigma }$$

13.2问题求解

13.2.1输入

#因变量新值得平均值的95%置信区间 y0_l_E<-y0-qt(1-alpha/2,n-2)*sqrt(h00)*sigma_hat#预测新值均值得置信下限 y0_u_E<-y0+qt(1-alpha/2,n-2)*sqrt(h00)*sigma_hat#预测新值均值得置信上限 

13.2.2输出

$E(y_0)$的95%置信区间为[3.283,4.122]

例题2.1(火灾)

（1）问题求解

根据上一题的程序，我们对例题2.1直接求解

1.1输入

setwd("E:/AllClass/junior1/RegressionAnalysis/unit2")#设置文件路径 #保留原路径setwd("C:/User/10854/documents") #以下利用理论方法使用一元回归模型 #导入数据 data<-read.csv("2-1.csv")#书本火灾题目，表数据2-1 x<-data[,1];y<-data[,2] n<-length(x) split.screen(c(1,3)) screen(1) plot(x,y) title(main="数据散点图") #求均值与回归变量lxx，lyy，lxy meanx<-mean(x);meany<-mean(y) lxx<-sum((x-meanx)^2) lyy<-sum((y-meany)^2) lxy<-sum((x-meanx)*(y-meany)) #回归系数估计 beta_1<-lxy/lxx#beta_1 beta_0<-meany-beta_1*meanx#beta_0 screen(2) plot(x,y) points(x,beta_0+beta_1*x,type="l") title(main="回归图") #预测值与平方和 y_hat<-beta_0+beta_1*x sse<-sum((y_hat-y)^2)#残差平方和 ssr<-sum((y_hat-meany)^2)#回归平方和 sst<-ssr+sse#总离差平方和 #回归误差ε的方差sigma估计 sigma_hat<-sqrt(1/(n-2)*sse) #对bet_0、beta_1的95%区间估计 alpha<-0.05 #beta_0,beta_1的分布标准差 sd.beta_0<-sqrt((1/n+(meanx^2)/lxx))*sigma_hat sd.beta_1<-sqrt(sigma_hat^2/lxx) beta_1l<-beta_1-qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信下限 beta_1u<-beta_1+qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信上限 beta_0l<-beta_0-qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信下限 beta_0u<-beta_0+qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信上限 #计算xy决定系数 R<-ssr/sst #回归方程的显著性检验 #方差分析F检验 f<-(ssr/1)/(sse/(n-2)) p1<-pf(f,1,n-2)#F为统计量、1为第一个自由度，n-2为第二个自由度 #回归系数的beta_1的t检验 t1<-beta_1/sd.beta_1#t统计量 p2<-pt(t1,n-2) #残差图 screen(3) e<-y_hat-y n<-length(e) sigma_u<-seq(2*sigma_hat,2*sigma_hat,length.out=n)#残差2σ原则 sigma_l<-seq(-2*sigma_hat,-2*sigma_hat,length.out=n) plot(x,e,ylim=c(12,-12)) points(x,sigma_u,type="l") points(x,sigma_l,type="l") title(main="残差图") 

1.2输出

回归方程的 ${\hat{\beta }}_1$的检验中，$p$值为$6.239\ast 10^{-9}<0.05$

数据的散点图，回归图，残差图如下

（2）结果分析

因为决定系数为0.92，具有比较强的线性相关性，且残差均在 $\pm 2\hat{\sigma }$ 内波动，而线性回归系数${\hat{\beta }}_1$的检验通过，所以可以认为火灾发生地点与最近的消防站距离和火灾损失呈线性关系，符合模型$y_i=10.28+4.91x_i$

数据与源代码链接