回归诊断

1.样本是否符合正态分布假设？

2.是否存在离群值导致模型发生较大误差？

3.线性模型是否合理？

4.误差是否满足独立性、等方差、正态分布等假设条件？

5.是否存在多重共线性

正态分布检验：函数shapiro.test()

P>0.05,正态分布

例题1

Anscomber数据

表给出的四组人造数据，每组数据集由11对点(x_i,y_i)组成，拟合于简单线性模型



试分析四组数据是否通过回归方程的检验，并用图形分析每组数据的基本情况。

解：输入数据，作回归分析（程序名：exam0611.R）

> Anscombe<-data.frame(
+ X=c(10.0, 8.0, 13.0, 9.0, 11.0, 14.0, 6.0, 4.0, 12.0, 7.0, 5.0),
+ Y1=c(8.04,6.95, 7.58,8.81,8.33,9.96,7.24,4.26,10.84,4.82,5.68),
+ Y2=c(9.14,8.14, 8.74,8.77,9.26,8.10,6.13,3.10, 9.13,7.26,4.74),
+ Y3=c(7.46,6.77,12.74,7.11,7.81,8.84,6.08,5.39, 8.15,6.44,5.73),
+ X4=c(rep(8,7), 19, rep(8,3)),
+ Y4=c(6.58,5.76,7.71,8.84,8.47,7.04,5.25,12.50, 5.56,7.91,6.89)
+ )
> summary(lm(Y1~X, data=Anscombe))

Call:
lm(formula = Y1 ~ X, data = Anscombe)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.92127 -0.45577 -0.04136  0.70941  1.83882 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.0001     1.1247   2.667  0.02573 * 
X             0.5001     0.1179   4.241  0.00217 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared:  0.6665,    Adjusted R-squared:  0.6295 
F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217

> summary(lm(Y2~X, data=Anscombe))

Call:
lm(formula = Y2 ~ X, data = Anscombe)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9009 -0.7609  0.1291  0.9491  1.2691 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)    3.001      1.125   2.667  0.02576 * 
X              0.500      0.118   4.239  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared:  0.6662,    Adjusted R-squared:  0.6292 
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179

> summary(lm(Y3~X, data=Anscombe))

Call:
lm(formula = Y3 ~ X, data = Anscombe)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1586 -0.6159 -0.2325  0.1510  3.2407 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.0075     1.1244   2.675  0.02542 * 
X             0.4994     0.1179   4.237  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared:  0.666, Adjusted R-squared:  0.6289 
F-statistic: 17.95 on 1 and 9 DF,  p-value: 0.002185

> summary(lm(Y4~X4, data=Anscombe))

Call:
lm(formula = Y4 ~ X4, data = Anscombe)

Residuals:
   Min     1Q Median     3Q    Max 
-1.751 -0.831  0.000  0.809  1.839 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.0017     1.1239   2.671  0.02559 * 
X4            0.4999     0.1178   4.243  0.00216 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared:  0.6667,    Adjusted R-squared:  0.6297 
F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165

系数	估计值	标准差	t-值	P-值
	3.0	1.125	2.67	0.026
	0.5	0.118	4.24	0.0022
方程	=1.24	=0.667	F=17.99	P=0.002

这四组数据的计算结果由上表所示（最多有0.01的误差）。从表所列结果，说明这四组数据全部能通过模型检验和方程的系数检验。由于每个数据集得到的各种统计量的值是相同的，因此，可能会认为每个数据集合对于线性模型会同等的适用，但事实并非如此。

我们画出四组数据的散点图和相应的回归直线。

> attach(Anscombe)
> par(mfrow=c(2,2))
> plot(Y1~X);abline(lm(Y1~X))
> plot(Y2~X);abline(lm(Y2~X))
> plot(Y3~X);abline(lm(Y3~X))
> plot(Y4~X4);abline(lm(Y4~X4))

第一个数据集合，如果简单线性回归模型合适的话，这就是我们所期望看到的数据集合。第二个数据集合，它给出一个不同的结论，即基于简单线性回归分析是不正确的，而一条光滑曲线，可能是二次多项式，可以以较小的剩余变异拟合数据。

> lm2.sol<-lm(Y2~X+I(X^2));summary(lm2.sol)

Call:
lm(formula = Y2 ~ X + I(X^2))

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0013287 -0.0011888 -0.0006294  0.0008741  0.0023776 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -5.9957343  0.0043299   -1385   <2e-16 ***
X            2.7808392  0.0010401    2674   <2e-16 ***
I(X^2)      -0.1267133  0.0000571   -2219   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001672 on 8 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 7.378e+06 on 2 and 8 DF,  p-value: < 2.2e-16

因此，回归方程为



更合理（见下图）

> x<-seq(min(X),max(X),len=200)
> y<-predict(lm2.sol,data.frame(X=x))
> plot(Y2~X);lines(x,y)

第三个数据集合表示，简单回归的描述对于大部分数据是正确的，但一个样本距离拟合回归直线太远，这称为异常值问题。很可能需要从数据集合中删除那个与其他数据不匹配的数据样本，回归需要根据剩下的10个样本重新拟合。

> i<-1:11;Y31<-Y3[i!=3];X3<-X[i!=3]
> lm3.sol<-lm(Y31~X3);summary(lm3.sol)

Call:
lm(formula = Y31 ~ X3)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0060173 -0.0012121 -0.0010173 -0.0008225  0.0140693 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.0106277  0.0057115   702.2   <2e-16 ***
X3          0.3450433  0.0006262   551.0   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.006019 on 8 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 3.036e+05 on 1 and 8 DF,  p-value: < 2.2e-16

得到的线性回归方程为



下图绘出修正后的直线方程

> plot(Y3~X);abline(lm3.sol)

最后一个数据集合与上述三个不同，没有足够的信息来对拟合模型作出判断。斜率参数的估计值很大程度上由的值来决定。如果第8号样本被删除，则不能估计。因此，我们无法相信这个一个综合分析，它对单个样本如此依赖。

例题2

Forbes数据

在十九世纪四、五十年代，苏格兰物理学家James D.Forbes，试图通过水的沸点来估计海拔高度。他知道通过气压计测得的大气压可用于得到海拔高度，高度越高，气压越低。在这里讨论的实验中，他研究了气压和沸点之间的关系，由于在当时，运输精密的气压计相当困难，这引起了他研究此问题的兴趣，测量沸点将给旅行者提供一个快速估计高度的方法。

Forbes在阿尔卑斯山及苏格兰收集数据，选定地点后，他装起仪器，测量气压及沸点，气压单位采用水银柱高度，并根据测量时周围气温与标准气温之间的差异校准气压，沸点用华氏温度来表示。我们从他1857年的论文中选取了个地方的数据，见下表。

在阿尔卑斯山及苏格兰的17个地方沸点（℉）及大气压（英寸汞柱）的Forbes数据

分析过程：

Forbes的理论认为，在观测值范围内，沸点和气压值的对数成一直线。由此，取10作为对数的底数。事实上，统计分析与对数的底是没有关系的。由于气压对数据值变化不大，最小值为1.318，而最大的为1.478，因此将所有气压的对数值乘以100，如表中第5列。这将不改变分析的主要性质的同时，避免研究非常小的数字。

求解过程：

输入数据，画出散点图（程序名：exam0804.R）

> X <- matrix(c(
+ 194.5, 20.79, 1.3179, 131.79,
+ 194.3, 20.79, 1.3179, 131.79,
+ 197.9, 22.40, 1.3502, 135.02,
+ 198.4, 22.67, 1.3555, 135.55,
+ 199.4, 23.15, 1.3646, 136.46,
+ 199.9, 23.35, 1.3683, 136.83,
+ 200.9, 23.89, 1.3782, 137.82,
+ 201.1, 23.99, 1.3800, 138.00,
+ 201.4, 24.02, 1.3806, 138.06,
+ 201.3, 24.01, 1.3805, 138.05,
+ 203.6, 25.14, 1.4004, 140.04,
+ 204.6, 26.57, 1.4244, 142.44,
+ 209.5, 28.49, 1.4547, 145.47,
+ 208.6, 27.76, 1.4434, 144.34,
+ 210.7, 29.04, 1.4630, 146.30,
+ 211.9, 29.88, 1.4754, 147.54,
+ 212.2, 30.06, 1.4780, 147.80),
+ ncol=4, byrow=T,
+ dimnames = list(1:17, c("F", "h", "log", "log100")))
> forbes<-as.data.frame(X)
> plot(forbes$F,forbes$log100)

Forbes数据的散点图的总的印象是，这些点基本上，但并不精确地，落在一条直线上，作回归分析

> lm.sol<-lm(log100~F,data=forbes)
> summary(lm.sol)

Call:
lm(formula = log100 ~ F, data = forbes)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.32261 -0.14530 -0.06750  0.02111  1.35924 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -42.13087    3.33895  -12.62 2.17e-09 ***
F             0.89546    0.01645   54.45  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3789 on 15 degrees of freedom
Multiple R-squared:  0.995, Adjusted R-squared:  0.9946 
F-statistic:  2965 on 1 and 15 DF,  p-value: < 2.2e-16

由计算结果得到：



对应于两个系数的P-值均，是非常显著的。

关于方程的检验，残差标准差，。相关系数的平方，，关于F-分布的P-值，也是非常显著的。

该模型能过t检验和F检验，因此，回归方程为



我们将得到的直线方程画在散点图上

> abline(lm.sol)

得到散点图和相应的回归直线，如图所示

下面分析残差，称



为回归方程的残差

在R中，函数residuals()计算回归方程残差。计算残差，并画出关于残差的散点图

> y.res<-residuals(lm.sol);plot(y.res)
> text(12,y.res[12],labels=12,adj=1.2)

其中text(12,y.res[12],labels=12,adj=1.2)是将第12号残差点标出。

从图中可以看到，第12个样本点可能会有问题，它比其他的样本点的残差大得多，因为其他点的残差的绝对值都小于0.35，而此点残差的绝对值均为1.3，因此，这个点可能不正确，或者模型的差假设不正确，或者是不是常数，等等。总之，需要对这个问题进行分析。

这里作简单处理，在，数据汇总，去掉第12号样本点

> i<-1:17;forbes12<-as.data.frame(X[i!=12,])
> lm12<-lm(log100~F, data=forbes12)
> summary(lm12)

Call:
lm(formula = log100 ~ F, data = forbes12)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.21175 -0.06194  0.01590  0.09077  0.13042 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -41.30180    1.00038  -41.29 5.01e-16 ***
F             0.89096    0.00493  180.73  < 2e-16 ***

在去掉第12号样本点后，回归方程的系数没有太大变化，但系数标准差和残差的标准差有很大的变化，减少了约3倍左右，相关系数也有提高。

在计算完回归模型后，计算其残差，并用shapiro.test()函数作残差的正态性检验

> y.res<-residuals(lm.sol)
> shapiro.test(y.res)

    Shapiro-Wilk normality test

data:  y.res
W = 0.54654, p-value = 3.302e-06

因此，残差不满足正态性假设。

在去掉第12号样本之后，再对所得到的回归模型的残差进行正态性检验

> y12.res<-residuals(lm12)
> shapiro.test(y12.res)

    Shapiro-Wilk normality test

data:  y12.res
W = 0.92215, p-value = 0.1827

能通过正态性检验，因此，去掉12号样本点是合理的。