By

4、异常观测值

  • 一个全面的回归分析要覆盖对异常值的分析，包括离群点、高杠杆值点和强影响点。

4.1、离群点

1
2
3
4
5
6
7
8
> #离群点：模型预测效果不佳的观测点。它们通常有很大的、或正或负的残差。
> #outlierTest()函数可以求得求得最大标准化残差绝对值Bonferroni（显著水平检验使用预设显著水平的1/2）调整后的p值
> outlierTest(fit)
       rstudent unadjusted p-value Bonferonni p
Nevada 3.542929         0.00095088     0.047544
> #可以看到Nevada被判定为离群点（p=0.048）。注意，该函数只是根据单个最大（或
> #正或负）残差值的显著性来判断是否有离群点。若不显著，则说明数据集中没有离群点；若显著，
> #则你必须删除该离群点，然后再检验是否还有其他离群点存在。

4.2、高杠杆值点

1
2
3
4
5
6
7
8
9
10
11
#高杠杆值点：是与其他预测变量有关的离群点，它们是由许多异常的预测变量值组合起来的，与响应变量值没有关系
#可根据帽子统计量（hat statistic）判断。对于一个给定的数据集，帽子均值为p/n，其中p 是模型估计的参数数目
#（包含截距项）， n 是样本量。一般来说，若观测点的帽子值大于帽子均值的2或3倍，即可以认定为高杠杆值点。
hat.plot<-function(fit){
  p<-length(coefficients(fit))# 输出模型参数数量，包含截距项
  n<-length(fitted(fit))#预测值的数量
  plot(hatvalues(fit),main="Index Plot of Hat Values")
  abline(h=c(2,3)*p/n,col="red",lty=2)
  identify(1:n,hatvalues(fit),names(hatvalues(fit)))
}
hat.plot(fit)

4.3、强影响点

1
2
3
4
5
6
7
8
9
#强影响点：对模型参数估计值影响有些比例失衡的点，例如，若移除模型的一个观测点时模型会发生巨大的改变，
#那么你就需要检测一下数据中是否存在强影响点了。
#方法：Cook距离，或称D统计量，以及变量添加图（added variable plot）
#Cook’s D值大于4/(n-k-1)，则表明它是强影响点，其中n 为样本量大小， k 是预测变量数目。
cutoff<-4/(nrow(states)-length(fit$coefficients)-2)
plot(fit,which = 4,cook.levels = cutoff)
abline(h=cutoff,lty=2,col="red")
#注意，虽然该图对搜寻强影响点很有用，但我逐渐发现以1为分割点更具一般性。
#若设定D=1为判别标准，则数据集中没有点看起来像是强影响点。-作者

1
2
3
4
#所谓变量添加图，即对于每个预测变量Xk，绘制Xk 在其他k-1个预测变量上回归的残差值相对于
#响应变量在其他k-1个预测变量上回归的残差值的关系图。 
library(car)
avPlots(fit, ask=FALSE, onepage=TRUE, id.method="identify")

1
2
3
4
5
6
7
#influencePlot 影响图
influencePlot(fit, id.method = "identify", main="Influence Plot",
              sub="Circle size is proportional to Cook's distance")
#纵坐标超过+2或小于2的州可被认为是离群点，水平轴超过0.2或0.3
#的州有高杠杆值（通常为预测值的组合）。圆圈大小与影响成比例，圆圈很大
#的点可能是对模型参数的估计造成的不成比例影响的强影响点

5、改进措施

5.1、删除观测点

  • 删除那些离群点或者强影响点
  • 谨慎！异常点可能是未知规律的表现！

5.2、变量变换

变换一般使用Y^λ替代Y
powerTransform()通过λ的极大似然估计来正太化变量X^λ

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
> summary(powerTransform(states$Murder))
bcPower Transformation to Normality 
              Est.Power Std.Err. Wald Lower Bound Wald Upper Bound
states$Murder    0.6055   0.2639           0.0884           1.1227
Likelihood ratio tests about transformation parameters
                           LRT df       pval
LR test, lambda = (0) 5.665991  1 0.01729694
LR test, lambda = (1) 2.122763  1 0.14512456
#结果表明可以使用Murder^0.6正则化Murder，但是λ=1的假设北邮被拒绝，所有没有强有力的证据表明需要变换
> #boxTidwell() 通过预测变量幂数的最大似然估计来改善线性关系
> boxTidwell(Murder~Population+Illiteracy,data=states)
           Score Statistic   p-value MLE of lambda
Population      -0.3228003 0.7468465     0.8693882
Illiteracy       0.6193814 0.5356651     1.3581188
iterations =  19 
#使用变换Population^0.87和Illiteracy^1.36能够大大改善线性关系。但是对Population
#（p=0.75）和Illiteracy（p=0.54）的计分检验又表明变量并不需要变换。

  • 谨慎对待变量变换:如果你不能证明A，那就证明B，假装它就是A。（对于统计学家来说，这很滑稽好笑。）如果你变换了变量，你的解释必须基于变换后的变量，而不是初始变量。若变换得没有意义，应该避免。如，你怎样解释自杀意念的频率与抑郁程度的立方根间的关系呢？

    5.3、增删变量

  • 删除变量在处理多重共线性时是一种非常重要的方法。做预测多重共线性没有影响，但是影响对变量的解释，解决通常删除共线性变量，或者使用岭回归

    5.4、尝试其他方法

6、选择最佳模型

6.1、模型比较

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
> #anova()函数可以比较两个嵌套模型（一些项完全包含在另一个模型中）的拟合优度。
> fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
> fit2<-lm(Murder~Population+Illiteracy,data=states)
> anova(fit2,fit1)
Analysis of Variance Table
Model 1: Murder ~ Population + Illiteracy
Model 2: Murder ~ Population + Illiteracy + Income + Frost
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     47 289.25                           
2     45 289.17  2  0.078505 0.0061 0.9939
> #函数同时还对是否应该添加Income和Frost到线性模型中进行了检验。
> #检验不显著,所以可以剔除
> 
> #AIC（Akaike Information Criterion，赤池信息准则）也可以用来比较模型，它考虑了模型的
> #统计拟合度以及用来拟合的参数数目。 AIC值越小越好，说明较少的参数获得了足够的拟合度。
> fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
> fit2<-lm(Murder~Population+Illiteracy,data=states)
> AIC(fit1,fit2)
     df      AIC
fit1  6 241.6429
fit2  4 237.6565

6.2、变量选择

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#逐步回归法：变量每次进入一个，但是每一步中，变量都会被重新评价，对模型没有贡献的变量将会被删除，
#预测变量可能会被添加、删除好几次，直到获得最优模型为止。
library(MASS)
fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
stepAIC(fit1)#默认使用向后逐步回归
Start:  AIC=97.75
Murder ~ Population + Illiteracy + Income + Frost
             Df Sum of Sq    RSS     AIC
- Frost       1     0.021 289.19  95.753
- Income      1     0.057 289.22  95.759
<none>                    289.17  97.749
- Population  1    39.238 328.41 102.111
- Illiteracy  1   144.264 433.43 115.986
Step:  AIC=95.75
Murder ~ Population + Illiteracy + Income
             Df Sum of Sq    RSS     AIC
- Income      1     0.057 289.25  93.763
<none>                    289.19  95.753
- Population  1    43.658 332.85 100.783
- Illiteracy  1   236.196 525.38 123.605
Step:  AIC=93.76
Murder ~ Population + Illiteracy
             Df Sum of Sq    RSS     AIC
<none>                    289.25  93.763
- Population  1    48.517 337.76  99.516
- Illiteracy  1   299.646 588.89 127.311
Call:
lm(formula = Murder ~ Population + Illiteracy, data = states)
Coefficients:
(Intercept)   Population   Illiteracy  
  1.6515497    0.0002242    4.0807366
1
2
3
4
5
6
7
#全子集回归
library(leaps)
leaps<-regsubsets(Murder~Population+Illiteracy+Income+Frost,data=states,nbest=4)
plot(leaps, scale="adjr2")
library(car)
subsets(leaps,statistic = "cp", main="Cp Plot for All Subsets Regression")
abline(1,1,lty=2,col="red")

  • 可以看出，双变量对应的模型能够得到最多的R平方增量
  • 基于Mallows Cp统计量的四个最佳模型。越好的模型离截距项和斜率均为1的直线越近。图中可以得出4个不错的模型（P-Il,P-Il-F,P-Il-In,P-Il-In-F）

7、深入分析

7.1、交叉验证

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#交叉验证
shrinkage<-function(fit,k=10){
  require(bootstrap)
  
  theta.fit<-function(x,y){lsfit(x,y)}
  theta.predict<-function(fit,x){cbind(1,x)%*%fit$coef}
  
  x<-fit$model[,2:ncol(fit$model)]
  y<-fit$model[,1]
  
  results<-crossval(x,y,theta.fit,theta.predict,ngroup=k)
  r2<-cor(y,fit$fitted.values)^2
  r2cv<-cor(y,results$cv.fit)^2
  cat("Original R-square=",r2,"\n")
  cat(k,"Fold Cross-Validated R-square=",r2cv,"\n")
  cat("Change=",r2-r2cv,"\n")
}
fit1<-lm(Murder~Population+Income+Illiteracy+Frost,data=states)
library(bootstrap)
shrinkage(fit1)
> Original R-square= 0.5669502 
> 10 Fold Cross-Validated R-square= 0.4174945 
> Change= 0.1494558 
fit2<-lm(Murder~Population+Illiteracy,data=states)
library(bootstrap)
shrinkage(fit2)
> Original R-square= 0.5668327 
> 10 Fold Cross-Validated R-square= 0.5193488 
> Change= 0.04748383 
#双变量比全变量模型交叉验证后R方减少的更少，说明其泛化能力强，R平方减少得越少，预测则越精确

7.2、相对重要性

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
> #相对重要性
> #scale将数据标准化为均值为0、标准差为1的数据集，这样用R回归即可获得标准化的回归系数。
> #评价预测变量相对重要性的方法一直在涌现。最简单的莫过于比较标准化的回归系数，它表
> #示当其他预测变量不变时，该预测变量一个标准差的变化可引起的响应变量的预期变化（以标准
> #差单位度量）。
> zstates <- as.data.frame(scale(states))
> zfit<-lm(Murder~Population+Income+Illiteracy+Frost,data=zstates)
> coef(zfit)
  (Intercept)    Population        Income    Illiteracy         Frost 
-2.054026e-16  2.705095e-01  1.072372e-02  6.840496e-01  8.185407e-03 
> #当其他因素不变时，文盲率一个标准差的变化将增加0.68个标准差的谋杀率。
> #根据标准化的回归系数，我们可认为Illiteracy是最重要的预测变量，而Frost是最不重要的。
#相对权重法：它是对所有可能子模型添加一个预测变量引起的R平方平均增加量的一个近似值
relweights <- function(fit, ...) {
  R <- cor(fit$model)
  nvar <- ncol(R)
  rxx <- R[2:nvar, 2:nvar]
  rxy <- R[2:nvar, 1]
  svd <- eigen(rxx)
  evec <- svd$vectors
  ev <- svd$values
  delta <- diag(sqrt(ev))
  
  # correlations between original predictors and new orthogonal variables
  lambda <- evec %*% delta %*% t(evec)
  lambdasq <- lambda^2
  
  # regression coefficients of Y on orthogonal variables
  beta <- solve(lambda) %*% rxy
  rsquare <- colSums(beta^2)
  rawwgt <- lambdasq %*% beta^2
  import <- (rawwgt/rsquare) * 100
  lbls <- names(fit$model[2:nvar])
  rownames(import) <- lbls
  colnames(import) <- "Weights"
  
  # plot results
  barplot(t(import), names.arg = lbls, ylab = "% of R-Square", 
          xlab = "Predictor Variables", main = "Relative Importance of Predictor Variables", 
          sub = paste("R-Square = ", round(rsquare, digits = 3)), 
          ...)
  return(import)
}
fit <- lm(Murder ~ Population + Illiteracy + Income + 
            Frost, data = states)
relweights(fit, col = "lightgrey")
             Weights
Population 14.723401
Illiteracy 59.000195
Income      5.488962
Frost      20.787442

By

回归

1、回归的多面性

  • 简单线性：用一个量化的解释变量预测一个量化的响应变量
  • 多项式：用一个量化的解释变量预测一个量化的响应变量，模型的关系是n阶多项式
  • 多元线性：用两个或多个量化的解释变量预测一个量化的响应变量
  • 多变量：用一个或多个解释变量预测多个响应变量
  • Logistic：用一个或多个解释变量预测一个类别型响应变量
  • 泊松：用一个或多个解释变量预测一个代表频数的响应变量
  • Cox比例风险：用一个或多个解释变量预测一个事件（死亡、失败或旧病复发）发生的时间
  • 时间序列：对误差项相关的时间序列数据建模
  • 非线性：用一个或多个量化的解释变量预测一个量化的响应变量，不过模型是非线性的
  • 非参数：用一个或多个量化的解释变量预测一个量化的响应变量，模型的形式源自数据形式，不事先设定
  • 稳健：用一个或多个量化的解释变量预测一个量化的响应变量，能抵御强影响点的干扰

    2、最小二乘回归


    为了能够恰当地解释OLS模型的系数，数据必须满足以下统计假设。
  • 正态性: 对于固定的自变量值，因变量值成正态分布。
  • 独立性: Yi值之间相互独立。
  • 线性: 因变量与自变量之间为线性相关。
  • 同方差性: 因变量的方差不随自变量的水平不同而变化。也可称作不变方差，但是说同方差性感觉上更犀利。

2.1 使用lm()拟合回归模型

1
myfit <- lm(formula, data)

常用的参数和解释

对拟合线性模型非常有用的其他函数:

2.2 简单线性回归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
> #OLS回归
> 
> #简单线性回归 女性身高体重数据
> 
> fit<-lm(weight~height, data=women)
> 
> summary(fit)
Call:
lm(formula = weight ~ height, data = women)
Residuals:
    Min      1Q  Median      3Q     Max 
-1.7333 -1.1333 -0.3833  0.7417  3.1167 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
height        3.45000    0.09114   37.85 1.09e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared:  0.991,	Adjusted R-squared:  0.9903 
F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14
> 
> women$weight
 [1] 115 117 120 123 126 129 132 135 139 142 146 150 154 159 164
> 
> fitted(fit)
       1        2        3        4        5        6        7        8        9       10       11 
112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 133.2833 136.7333 140.1833 143.6333 147.0833 
      12       13       14       15 
150.5333 153.9833 157.4333 160.8833 
> 
> residuals(fit)
          1           2           3           4           5           6           7           8 
 2.41666667  0.96666667  0.51666667  0.06666667 -0.38333333 -0.83333333 -1.28333333 -1.73333333 
          9          10          11          12          13          14          15 
-1.18333333 -1.63333333 -1.08333333 -0.53333333  0.01666667  1.56666667  3.11666667 
> 
> plot(women)
> abline(fit)
# 此模型的得到的方程式是：weight= -87.51667+3.45*height

2.3 多项式回归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
> #增加一个二次项提高精度
> 
> fit2<-lm(weight~height+I(height^2),data=women)
> 
> #此处I()的含义是height^2表示为一个常规的R表达式，因为^在在表达式中有特殊含义
> 
> summary(fit2)
Call:
lm(formula = weight ~ height + I(height^2), data = women)
Residuals:
     Min       1Q   Median       3Q      Max 
-0.50941 -0.29611 -0.00941  0.28615  0.59706 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 261.87818   25.19677  10.393 2.36e-07 ***
height       -7.34832    0.77769  -9.449 6.58e-07 ***
I(height^2)   0.08306    0.00598  13.891 9.32e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3841 on 12 degrees of freedom
Multiple R-squared:  0.9995,	Adjusted R-squared:  0.9994 
F-statistic: 1.139e+04 on 2 and 12 DF,  p-value: < 2.2e-16
> 
> women$weight
 [1] 115 117 120 123 126 129 132 135 139 142 146 150 154 159 164
> 
> fitted(fit2)
       1        2        3        4        5        6        7        8        9       10       11       12       13 
115.1029 117.4731 120.0094 122.7118 125.5804 128.6151 131.8159 135.1828 138.7159 142.4151 146.2804 150.3118 154.5094 
      14       15 
158.8731 163.4029 
> 
> residuals(fit2)
           1            2            3            4            5            6            7            8            9 
-0.102941176 -0.473109244 -0.009405301  0.288170653  0.419618617  0.384938591  0.184130575 -0.182805430  0.284130575 
          10           11           12           13           14           15 
-0.415061409 -0.280381383 -0.311829347 -0.509405301  0.126890756  0.597058824 
> 
> plot(women)
> lines(women$height,fitted(fit2))

1
2
3
4
5
6
7
8
9
10
11
#car包中的scatterplot方便的绘制二元关系图
library(car)
scatterplot(weight~height,
            data=women,
            spread=FALSE,#spread=FALSE选项删除了残差正负均方根在平滑
            #曲线上的展开和非对称信息
            pch=19,
            main="Women Age 30-39",
            xlab="Height (inches)",
            ylab="Weight (lbs.)"
            )

2.4 多元线性回归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
> #多项式回归
> #观察谋杀率 人口 文盲率 平均收入 结霜天数关系
> states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])
> #观察二元关系
> cor(states)
               Murder Population Illiteracy     Income      Frost
Murder      1.0000000  0.3436428  0.7029752 -0.2300776 -0.5388834
Population  0.3436428  1.0000000  0.1076224  0.2082276 -0.3321525
Illiteracy  0.7029752  0.1076224  1.0000000 -0.4370752 -0.6719470
Income     -0.2300776  0.2082276 -0.4370752  1.0000000  0.2262822
Frost      -0.5388834 -0.3321525 -0.6719470  0.2262822  1.0000000
> 
> library(car)
> 
> #大致观察数据图示关系(见下图)
> scatterplotMatrix(states, spread = FALSE, lty = 2, main="Scatter Plot Matrix")
> 
> #回归
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
> 
> summary(fit)
Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
    data = states)
Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
Population  2.237e-04  9.052e-05   2.471   0.0173 *  
Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
Income      6.442e-05  6.837e-04   0.094   0.9253    
Frost       5.813e-04  1.005e-02   0.058   0.9541    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared:  0.567,	Adjusted R-squared:  0.5285 
F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08
> 
> #当预测变量不止一个时，回归系数的含义为，一个预测变量增加一个单位，其他预测变量保
> #持不变时，因变量将要增加的数量。
> #例如本例中，文盲率的回归系数为4.14，表示控制人口、收入和温度不变时，文盲率上升1%，
> #谋杀率将会上升4.14%，它的系数在p<0.001的水平下显著不为0。

2.5 有交互项的多元线性回归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
> #有交互项的回归  y=x-y+xy
> fit <- lm(mpg ~ hp + wt +hp:wt, data=mtcars)
> summary(fit)
Call:
lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars)
Residuals:
    Min      1Q  Median      3Q     Max 
-3.0632 -1.6491 -0.7362  1.4211  4.5513 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 49.80842    3.60516  13.816 5.01e-14 ***
hp          -0.12010    0.02470  -4.863 4.04e-05 ***
wt          -8.21662    1.26971  -6.471 5.20e-07 ***
hp:wt        0.02785    0.00742   3.753 0.000811 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.153 on 28 degrees of freedom
Multiple R-squared:  0.8848,	Adjusted R-squared:  0.8724 
F-statistic: 71.66 on 3 and 28 DF,  p-value: 2.981e-13
> 
> #若两个预测
> #变量的交互项显著，说明响应变量与其中一个预测变量的关系依赖于另外一个预测变量的水平。
> #因此此例说明，每加仑汽油行驶英里数与汽车马力的关系依车重不同而不同。
> 
> #通过effects包中的effect()函数，你可以用图形展示交互项的结果。
> library(effects)
> plot(effect("hp:wt", fit, xlevels=list(wt=c(2.2,3.2,4.2))),multiline=TRUE)


从图中可以很清晰地看出，随着车重的增加，马力与每加仑汽油行驶英里数的关系减弱了。
当wt = 4.2时，直线几乎是水平的，表明随着hp的增加， mpg不会发生改变。

3、回归诊断

  • 对回归结果的信念，都只建立在你的数据满足统计假设的前提之上。

    3.1 标准方法

    1
    2
    3
    4
    5
    #回归诊断
    fit<- lm(weight~height,data=women)
    #将plot()函数绘制的四幅图形组合在一个大的2×2的图中。
    par(mfrow=c(2,2))
    plot(fit)

  • 正态性：当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正
    态分布。正态Q-Q图（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概率
    图。若满足正态假设，那么图上的点应该落在呈45度角的直线上；若不是如此，那么就
    违反了正态性的假设。

  • 独立性：由样本的数据收集去确保

  • 线性：若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。
    换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“残差图与拟合图”
    （Residuals vs Fitted，左上）中可以清楚的看到一个曲线关系，这暗示着你可能需要对回
    归模型加上一个二次项。

  • 同方差性 ：若满足不变方差假设，那么在位置尺度图（Scale-Location Graph，左下）中，
    水平线周围的点应该随机分布。该图似乎满足此假设。

  • 最后一幅“残差与杠杆图”（Residuals vs Leverage，右下）提供了你可能关注的单个观测点的信息。从图形可以鉴别出离群点、高杠杆值点和强影响点。

看下二次拟合的图：

1
2
3
4
5
6
7
fit<- lm(weight~height+I(height^2),data=women)
#将plot()函数绘制的四幅图形组合在一个大的2×2的图中。
par(mfrow=c(2,2))
plot(fit)
#这第二组图表明多项式回归拟合效果比较理想，基本符合了线性假设、残差正态性（除了
#观测点13）和同方差性（残差方差不变）。观测点15看起来像是强影响点（根据是它有较大的
#Cook距离值），删除它将会影响参数的估计。事实上，删除观测点13和15，模型会拟合得会更好。

3.2 改进的方法

car包提供多种函数用以评价回归模型：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
qqPlot() 分位数比较图
durbinWatsonTest() 对误差自相关性做Durbin-Watson检验
crPlots() 成分与残差图
ncvTest() 对非恒定的误差方差做得分检验
spreadLevelPlot() 分散水平检验
outlierTest() Bonferroni离群点检验
avPlots() 添加的变量图形
inluencePlot() 回归影响图
scatterplot() 增强的散点图
scatterplotMatrix() 增强的散点图矩阵
vif() 方差膨胀因子

  • 正态性：qqPlot()
    1
    2
    3
    4
    5
    6
    7
    #它画出了在n-p-1个自由度的t分布下的学生化残差（studentized residual，也称学生化删除残差或折叠
    #化残差）图形，其中n是样本大小， p是回归参数的数目（包括截距项）。
    library(car)
    fit<-lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
    par(mfrow=c(1,1))
    qqPlot(fit,labels=row.names(states),id.method = "identify",simulate=TRUE, main="Q-Q Plot")
    #当simulate=TRUE时， 95%的置信区间将会用参数自助法

另一种和可视化误差的方法

1
2
3
4
5
6
7
8
9
10
residplot<-function(fit,nbreaks=10){
  z<-rstudent(fit)
  hist(z,breaks=nbreaks,freq = FALSE, xlab="Studentized Residual",main="Distribution of Errors")
  rug(jitter(z), col = "brown")
  curve(dnorm(x,mean=mean(z),sd=sd(z)),
        add=TRUE, col="blue",lwd=2)
  lines(density(z)$x,density(z)$y,col="red",lwd=2,lty=2)
  legend("topright",legend = c("Normal Curve","Kernel Density Curve"),lty=1:2,col=c("blue","red"),cex=.7)
}
residplot(fit)

  • 独立性

    1
    2
    3
    4
    5
    6
    7
    8
    9
    > #误差的独立性
    > durbinWatsonTest(fit)
     lag Autocorrelation D-W Statistic p-value
       1      -0.2006929      2.317691   0.258
     Alternative hypothesis: rho != 0
    > #p值不显著说明无相关性，误差项之间独立，滞后项（lag=1）表明数据集中
    > #每个数据都是与其后一个数据进行比较的。该检验适用于时间独立的数据，对于非聚集型的数据
    > #并不适用。注意， durbinWatsonTest()函数使用自助法（参见第12章）来导出p值，如果添加
    > #了选项simulate=FALSE，则每次运行测试时获得的结果都将略有不同。

  • 线性

    1
    2
    3
    4
    crPlots(fit)
    #若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分，
    #那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代
    #替X），或用其他回归变体形式而不是线性回归。

  • 同方差性
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    > #线性
    > crPlots(fit)
    > #若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分，
    > #那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代
    > #替X），或用其他回归变体形式而不是线性回归。
    > 
    > #同方差性
    > #ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。
    > #若检验显著，则说明存在异方差性（误差方差不恒定）
    > ncvTest(fit)
    Non-constant Variance Score Test 
    Variance formula: ~ fitted.values 
    Chisquare = 1.746514    Df = 1     p = 0.1863156 
    > p=0.18不显著接受零假设
    > #spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。
    > spreadLevelPlot(fit)
    Suggested power transformation:  1.209626 
    #你将会看到一个非水平的曲线。代码结果建议幂次变换（suggested power transformation）的含义
    #是，经过p次幂（Y^p）变换，非恒定的误差方差将会平稳。
    #建议幂次转换为0.5，在回归等式中用 Y 代替Y，可能会使模型满足同方差性。若建议幂次为0，
    #则使用对数变换。对于当前例子，异方差性很不明显，因此建议幂次接近1（不需要进行变换）。

3.3 线性模型假设的综合验证

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
> library(gvlma)
> gvmodel<-gvlma(fit)
> summary(gvmodel)
Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
    data = states)
Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
Population  2.237e-04  9.052e-05   2.471   0.0173 *  
Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
Income      6.442e-05  6.837e-04   0.094   0.9253    
Frost       5.813e-04  1.005e-02   0.058   0.9541    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared:  0.567,	Adjusted R-squared:  0.5285 
F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 
Call:
 gvlma(x = fit) 
                    Value p-value                Decision
Global Stat        2.7728  0.5965 Assumptions acceptable.
Skewness           1.5374  0.2150 Assumptions acceptable.
Kurtosis           0.6376  0.4246 Assumptions acceptable.
Link Function      0.1154  0.7341 Assumptions acceptable.
Heteroscedasticity 0.4824  0.4873 Assumptions acceptable.
#从输出项（Global Stat中的文字栏）我们可以看到数据满足OLS回归模型所有的统计假设
#（p=0.5965597） 。若Decision下的文字表明违反了假设条件（比如p<0.05） ，你可以使用前几节讨论
#的方法来判断哪些假设没有被满足。

3.4 多重共线性

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
> #多重共线性
> #例如 握力~年龄+DOB(Date Of Birth)
> #F检验显著， p<0.001。但是当你观察DOB和年龄的回归系数时，却
> #发现它们都不显著（也就是说无法证明它们与握力相关）
> #原因是DOB与年龄在四舍五入后相关性极大。回归系数测量的是当其他预测变量不变时，某
> #个预测变量对响应变量的影响。那么此处就相当于假定年龄不变，然后测量握力与年龄的关系，
> #这种问题就称作多重共线性（multicollinearity）。它会导致模型参数的置信区间过大，使单个系数
> #解释起来很困难。
> #使用VIF度量（Variance Inflation Factor，方差膨胀因子），
> #VIF的平方根表示变量回归参数的置信区间能膨胀为与模型无关的预测变量的程度（因此而得名）
> library(car)
> vif(fit)
Population Illiteracy     Income      Frost 
  1.245282   2.165848   1.345822   2.082547 
> sqrt(vif(fit))>2
Population Illiteracy     Income      Frost 
     FALSE      FALSE      FALSE      FALSE

By

7个设计原则(前五个重要SOLID)

单一职责原则（SRP：Single responsibility principle）

  • 一个类，只有一个引起它变化的原因。应该只有一个职责。每一个职责都是变化的一个轴线，如果一个类有一个以上的职责，这些职责就耦合在了一起。这会导致脆弱的设计。当一个职责发生变化时，可能会影响其它的职责。另外，多个职责耦合在一起，会影响复用性。例如：要实现逻辑和界面的分离。

    开闭原则（OCP；Open Close Principle）

  • 开闭原则就是说对扩展开放，对修改关闭。在程序需要进行拓展的时候，不能去修改原有的代码，实现一个热插拔的效果。所以一句话概括就是：为了使程序的扩展性好，易于维护和升级。想要达到这样的效果，我们需要使用接口和抽象类。

    里氏代换原则（LSP：Liskov Substitution Principle）

  • 里氏代换原则(Liskov Substitution Principle LSP)面向对象设计的基本原则之一。 里氏代换原则中说，任何基类可以出现的地方，子类一定可以出现。 LSP是继承复用的基石，只有当衍生类可以替换掉基类，软件单位的功能不受到影响时，基类才能真正被复用，而衍生类也能够在基类的基础上增加新的行为。里氏代换原则是对“开-闭”原则的补充。实现“开-闭”原则的关键步骤就是抽象化。而基类与子类的继承关系就是抽象化的具体实现，所以里氏代换原则是对实现抽象化的具体步骤的规范。

    接口隔离原则（ISP：Interface Segregation Principle）

  • 这个原则的意思是：使用多个隔离的接口，比使用单个接口要好。还是一个降低类之间的耦合度的意思，从这儿我们看出，其实设计模式就是一个软件的设计思想，从大型软件架构出发，为了升级和维护方便。所以上文中多次出现：降低依赖，降低耦合。

    依赖倒转原则（DIP：Dependence Inversion Principle）

  • 所谓依赖倒置原则（Dependence Inversion Principle）就是要依赖于抽象，不要依赖于具体。简单的说就是要求对抽象进行编程，不要对实现进行编程，这样就降低了客户与实现模块间的耦合。
  • 实现开闭原则的关键是抽象化，并且从抽象化导出具体化实现，如果说开闭原则是面向对象设计的目标的话，那么依赖倒转原则就是面向对象设计的主要手段。

    合成复用原则（CRP：Composite Reuse Principle）

  • 合成复用原则就是指在一个新的对象里通过关联关系（包括组合关系和聚合关系）来使用一些已有的对象，使之成为新对象的一部分；新对象通过委派调用已有对象的方法达到复用其已有功能的目的。简言之：要尽量使用组合/聚合关系，少用继承。

    迪米特法则（最少知道原则）（DP：Demeter Principle）

  • 为什么叫最少知道原则，就是说：一个实体应当尽量少的与其他实体之间发生相互作用，使得系统功能模块相对独立。也就是说一个软件实体应当尽可能少的与其他实体发生相互作用。这样，当一个模块修改时，就会尽量少的影响其他的模块，扩展会相对容易，这是对软件实体之间通信的限制，它要求限制软件实体之间通信的宽度和深度。

设置模式之间的关系