By

文章目录
  1. 4、异常观测值
    1. 4.1、离群点
    2. 4.2、高杠杆值点
    3. 4.3、强影响点
  2. 5、改进措施
    1. 5.1、删除观测点
    2. 5.2、变量变换
    3. 5.3、增删变量
    4. 5.4、尝试其他方法
  3. 6、选择最佳模型
    1. 6.1、模型比较
    2. 6.2、变量选择
  4. 7、深入分析
    1. 7.1、交叉验证
    2. 7.2、相对重要性

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

文章目录
  1. 4、异常观测值
    1. 4.1、离群点
    2. 4.2、高杠杆值点
    3. 4.3、强影响点
  2. 5、改进措施
    1. 5.1、删除观测点
    2. 5.2、变量变换
    3. 5.3、增删变量
    4. 5.4、尝试其他方法
  3. 6、选择最佳模型
    1. 6.1、模型比较
    2. 6.2、变量选择
  4. 7、深入分析
    1. 7.1、交叉验证
    2. 7.2、相对重要性