Anscombe

1973年， Francis Anscombe创造了一个叫anscombe的数据集。这个数据集含有四个分组(set)。每个分组之间有x和y两个坐标。通过对每组数据的分析，我们发现每组数据中的统计结果都是相似的，但是他们的图形是完全不同的。

• 这个数据集告诉我们，在进行数据分析的时候，除了查看主要的结果参数之外也要认为的观察其数据的分布。

• 另外相关分析一定是要基于一定的相关基础上进行的分析。不能天马行空的做相关分析。

# Compute properties of Anscombe
Anscombe <- list()
for (i in as.character(1:4)) {
    Anscombe[[i]] <- anscombe %>% select(matches(i)) %>% mutate(set = i)
    names(Anscombe[[i]]) <- c("x", "y", "set")
}
Anscombe <- bind_rows(Anscombe)
Anscombe %>%
  group_by(set) %>%
  summarize(N = n(), mean(x), sd(x), mean(y), sd(y), cor(x,y))

## # A tibble: 4 x 7
##   set       N `mean(x)` `sd(x)` `mean(y)` `sd(y)` `cor(x, y)`
##   <chr> <int>     <dbl>   <dbl>     <dbl>   <dbl>       <dbl>
## 1 1        11         9    3.32      7.50    2.03       0.816
## 2 2        11         9    3.32      7.50    2.03       0.816
## 3 3        11         9    3.32      7.5     2.03       0.816
## 4 4        11         9    3.32      7.50    2.03       0.817

### plot scatter plot
ggplot(Anscombe, aes(x, y)) + geom_point() + 
    geom_smooth(method = "lm", se = F) + 
    facet_wrap(.~set, ncol = 2)

线性回归的简单可视化

简单的线性回归我们可以通过geom_smooth(method = "lm")进行可视化。同时也可以通过geom_abline来截距(intercept)和系数(hgt)

其中

ggplot(data = bdims, aes(x = hgt, y = wgt)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)

回归诊断

在得到预测模型之后，我们要评判模型是否成立。最基本的是通过plot来进行预测判断的。

mod <- lm(wgt ~ hgt, data = bdims)
par(mfrow = c(2,2))
plot(mod)

library(car)

## Loading required package: carData

## 
## Attaching package: 'carData'

## The following object is masked _by_ '.GlobalEnv':
## 
##     Anscombe

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:openintro':
## 
##     densityPlot

qqPlot(mod)

## [1] 124 474

durbinWatsonTest(mod)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.0519505      1.894424   0.192
##  Alternative hypothesis: rho != 0

crPlots(mod)

ncvTest(mod)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.1361065, Df = 1, p = 0.71218

spreadLevelPlot(mod)

## 
## Suggested power transformation:  0.6747609

library(gvlma)
gvlma(mod)

## 
## Call:
## lm(formula = wgt ~ hgt, data = bdims)
## 
## Coefficients:
## (Intercept)          hgt  
##    -105.011        1.018  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = mod) 
## 
##                      Value   p-value                   Decision
## Global Stat        95.4846 0.000e+00 Assumptions NOT satisfied!
## Skewness           59.6610 1.132e-14 Assumptions NOT satisfied!
## Kurtosis           34.0772 5.297e-09 Assumptions NOT satisfied!
## Link Function       0.4014 5.264e-01    Assumptions acceptable.
## Heteroscedasticity  1.3451 2.461e-01    Assumptions acceptable.

数据转换

dat <- read.table("~/dat.txt", header = T)
hist(dat$WBC)

f2 <- powerTransform(dat$WBC)
hist(bcPower(dat$WBC, f2$lambda)) ##最近转换

hist(bcPower(dat$WBC, f2$roundlam)) ##近似lambda值转换

boxTidwell(prestige ~ income + education, ~ type + poly(women, 2), data=Prestige)

##           MLE of lambda Score Statistic (z)  Pr(>|z|)    
## income         -0.34763             -4.4824 7.381e-06 ***
## education       1.25383              0.2170    0.8282    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## iterations =  8

线性回归统计量

• 进行线性回归的基本参数是lm其介绍format和一个data。返回的结果是包含所有结果的list。我们可以通过summary查看汇总的结果；coef查看截距和系数；confint查看95%可信区间;fitted.values查看每个点的拟合值(根据模型预测到的y值)；residuals查看残差(真实的y到模型的距离)。

• 通过bromm::augment可以查看相关的结果

mod <- lm(wgt ~ hgt, data = bdims)
summary(mod)

## 
## Call:
## lm(formula = wgt ~ hgt, data = bdims)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.743  -6.402  -1.231   5.059  41.103 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -105.01125    7.53941  -13.93   <2e-16 ***
## hgt            1.01762    0.04399   23.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.308 on 505 degrees of freedom
## Multiple R-squared:  0.5145, Adjusted R-squared:  0.5136 
## F-statistic: 535.2 on 1 and 505 DF,  p-value: < 2.2e-16

coef(mod)

## (Intercept)         hgt 
## -105.011254    1.017617

head(data.frame(fitted.values = fitted.values(mod), residuals = residuals(mod)))

##   fitted.values  residuals
## 1      72.05406  -6.454065
## 2      73.37697  -1.576967
## 3      91.89759 -11.197592
## 4      84.77427 -12.174274
## 5      85.48661  -6.686606
## 6      79.68619  -4.886191

broom::augment(mod)

## # A tibble: 507 x 9
##      wgt   hgt .fitted .se.fit  .resid    .hat .sigma   .cooksd .std.resid
##    <dbl> <dbl>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>     <dbl>      <dbl>
##  1  65.6  174     72.1   0.432  -6.45  0.00215   9.31 0.000520     -0.694 
##  2  71.8  175.    73.4   0.452  -1.58  0.00236   9.32 0.0000340    -0.170 
##  3  80.7  194.    91.9   1.07  -11.2   0.0131    9.30 0.00976      -1.21  
##  4  72.6  186.    84.8   0.792 -12.2   0.00724   9.30 0.00628      -1.31  
##  5  78.8  187.    85.5   0.818  -6.69  0.00773   9.31 0.00203      -0.721 
##  6  74.8  182.    79.7   0.615  -4.89  0.00437   9.31 0.000607     -0.526 
##  7  86.4  184     82.2   0.700   4.17  0.00566   9.32 0.000575      0.449 
##  8  78.4  184.    82.7   0.718  -4.34  0.00596   9.32 0.000655     -0.468 
##  9  62    175     73.1   0.447 -11.1   0.00230   9.30 0.00164      -1.19  
## 10  81.6  184     82.2   0.700  -0.630 0.00566   9.32 0.0000131    -0.0679
## # … with 497 more rows

• 线性回归一般就是均值拟合。因此对于拟合值的均值和具体值的均值其实式样的。另外的话，对于残差的均值也是无限接近于0的。

# Mean of weights equal to mean of fitted values?
mean(bdims$wgt) == mean(fitted.values(mod))

## [1] TRUE

# Mean of the residuals
mean(residuals(mod))

## [1] -1.266971e-15

进行预测

• 我们在预测好模型之后，需要用预测的模型在新的数据集中进行预测。以便评估其预测能力。一般来说我们可以通过predict来得到预测的y值。同样的也是可以通过broom::augment(mod, newdata)来获得新的预测值。

ben <- data.frame(wgt = 74.8, hgt = 182.8)
predict(mod, ben)

##        1 
## 81.00909

broom::augment(mod, newdata = ben)

## # A tibble: 1 x 4
##     wgt   hgt .fitted .se.fit
##   <dbl> <dbl>   <dbl>   <dbl>
## 1  74.8  183.    81.0   0.659

残差标准误(Residual standard error)

如果想要查看拟合的模型和真实数据之间的差异，可以是通过RMSE来评价。

• 如果一个RMSE为4.67则代表在模型拟合的时候，预测到的y值和真实的y值之间存在4.67%的差异。

• RMSE越小代表其拟合的越好。

• 一个模型的RMSE是通过残方差/自由度得到的。这个结果可以在summary中看到

RMSE <- sqrt(sum(residuals(mod)^2) / df.residual(mod)); RMSE

## [1] 9.30804

R平方(R squared)

RMSE可以用来评估模型拟合的值和实际y之间的差异是多少。而R squared则代表了在y的变异中，x的解释度。

• 例如如果美国各县贫困率(y)和高中毕业率(x)的线性回归模型的R squared为0.464。 那么可以说明，美国各县贫困率变异率的46.4％可以用高中毕业率来解释。

• 一个模型的R squared主要是通过1 - 残差方差/y方差得到的。这个结果在sumamry这也可以看到

bdims_tidy <- broom::augment(mod)
bdims_tidy %>%
  summarize(var_y = var(wgt), var_e = var(.resid)) %>%
  mutate(R_squared = 1 - var_e/var_y)

## # A tibble: 1 x 3
##   var_y var_e R_squared
##   <dbl> <dbl>     <dbl>
## 1  178.  86.5     0.515

异常点

异常点通常分为三种:

1. 异常点：那些距离回归直线很远，通常这些点的残差都会很大。一般来说

• Q-Q图，落在置信区间带外的点即可被认为是离群点。
• 另外一个粗糙的判断准则:标准化残差值大于2或者小于-2的可能是离群点。
• car包中的outlierTest可以计算离群值

outlierTest(model = mod)

##     rstudent unadjusted p-value Bonferonni p
## 474 4.505712         8.2305e-06    0.0041729
## 124 4.435050         1.1309e-05    0.0057337

2. 高杠杆点：这样的点拥有很极端的x值。比如其他样本点的x值都只有几十，而有一个样本点的x值超过了100，这时候这个点将会极大地影响回归模型。就像物理中的力矩一样，这种样本点有很高的leverage，所以叫做leverage points。有的leverage points可以很好地融入模型中，不会对模型造成很大的影响，通常也叫做good leverage points。有时，这种good leverage points证明了模型的普适性；但是它们同时会增大 R^2 ，使我们对模型过度自信。所以good leverage points也并不是完全没有弊端。

• 高杠杆值的观测点可通过帽子统计量(hat statistic)判断。对于一个给定的数据集，帽子均值为p/n，其中p是模型估计的参数数目(包含截距项)，n是样本量。一般来说，若观测点的帽子值大于帽子均值的2或3倍。在broom包中的.hat列则为帽子统计值。同时hatvalues(mod)也可以得到帽子值。

hat_avg <- length(coef(mod))/length(fitted(mod)); hat_avg

## [1] 0.003944773

mod %>% broom::augment() %>% select(wgt, hgt, .fitted, .resid, .hat) %>% arrange(desc(.hat)) %>% head()

## # A tibble: 6 x 5
##     wgt   hgt .fitted .resid   .hat
##   <dbl> <dbl>   <dbl>  <dbl>  <dbl>
## 1  85.5  198.    96.6 -11.1  0.0182
## 2  90.9  197.    95.6  -4.66 0.0170
## 3  49.8  147.    44.8   5.02 0.0148
## 4  80.7  194.    91.9 -11.2  0.0131
## 5  95.9  193     91.4   4.51 0.0126
## 6  44.8  150.    47.1  -2.32 0.0124

3. 强影响点：有good leverage points自然有bad leverage points。既是leverage points又是outliers的点就被称为bad leverage points，也就是influential points。它们的存在极大地影响了模型的可靠性，因为它们会把回归直线向自己的方向“拉扯”。

• Cook距离，一般来说，Cook’s D值大于4/(n-k-1)，则表明它是强影响点，其中n 为样本量大小，k是预测变量数目。通过broom::augment中的.cooksd代表其值。Cook’s D图有助于鉴别强影响点，但是并不提供关于这些点如何影响模型的信息。变量添加 图弥补了这个缺陷。 PS:虽然该图对搜寻强影响点很有用，逐渐发现以1为分割 点比4/(nk 1)更具一般性。若设定D=1为判别标准，则数据集中没有点看起来像是强影响点。

cutoff <- 4/(nrow(bdims) - length(mod$coefficients) - 2)
plot(mod, which = 4, cook.levels = cutoff)
abline(h = cutoff, lty = 2, col = "red")

4. 综合评估:利用car::influencePlot函数，你还可以将离群点、杠杆值和强影响点的 信息整合到一幅图形中

• 纵坐标超过+2或小于2的州可被认为是离群点，水平轴超过0.2或0.3的有高杠杆值(通常为预测值的组合)。圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点

influencePlot(mod)

##        StudRes         Hat       CookD
## 80  -0.5046944 0.017018035 0.002208169
## 124  4.4350503 0.002961811 0.028173878
## 127 -1.2017305 0.018199677 0.013373432
## 349  2.6569840 0.010944356 0.038595550
## 474  4.5057124 0.002788117 0.027335737

嵌套模型比较

使用anova函数比较两个嵌套模型的拟合优度。所谓嵌套模型就是一个模型的所有参数位于另外一个模型当中。

States <- as.data.frame(state.x77) %>% select(Population:Illiteracy, Murder, Frost)
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = States)
fit2 <- lm(Murder ~ Population + Illiteracy, data = States)
anova(fit1, fit2)

## Analysis of Variance Table
## 
## Model 1: Murder ~ Population + Illiteracy + Income + Frost
## Model 2: Murder ~ Population + Illiteracy
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     45 289.17                           
## 2     47 289.25 -2 -0.078505 0.0061 0.9939

通过结果发现P = 0.9939，说明增加两个参数对于模型没有影响。

AIC比较

AIC(Akaike Information Criterion，赤池信息准则)，可以用来比较模型。AIC值越小的模型要优先选择，它说明模型用较少的参数 获得了足够的拟合度。

AIC(fit1, fit2)

##      df      AIC
## fit1  6 241.6429
## fit2  4 237.6565

结果上来看fit2的AIC更小，说明模型更好

逐步回归

逐步回归中，模型会一次添加或者删除一个变量，直到达到某个判停准则为止。分为逐步向前；逐步向后以及向前向后逐步回归。在基础包当中可以使用step可是实现。另外的话在MASS包的stepAIC可以根据AIC更好的实现逐步回归，通过direction决定方向

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked _by_ '.GlobalEnv':
## 
##     mammals

## The following object is masked from 'package:dplyr':
## 
##     select

## The following objects are masked from 'package:openintro':
## 
##     housing, mammals

#### direction接受：backward,forward和both(默认)
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

全子集回归

逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但是不能保证模型就是最佳模 型，因为不是每一个可能的模型都被评价了。全子集回归所做的就是把模型中所有的组合全部评估一遍。使用leaps::regsubsets可以来进行全子集回归。在函数中的nsets可以设置展示结果的方式。如果不设置则全部展示所有结果。如果例如：nsets = 2则代表先展示两个最佳的单预测变量模型，然后展示两个最佳的双预测变量模型，以此类推，直到包含 所有的预测变量。

• 结果刚中会返回R平方；调整R平方和Mallows Cp统计量。 R平方含义是预测变量解释响应变量的程度;调整R平方与之类似，但考虑了模型的参数数 目。R平方总会随着变量数目的增加而增加。当与样本量相比，预测变量数目很大时，容易导致 3 过拟合。R平方很可能会丢失数据的偶然变异信息，而调整R平方则提供了更为真实的R平方估计。 另外，Mallows Cp统计量也用来作为逐步回归的判停规则。广泛研究表明，对于一个好的模型， 它的Cp统计量非常接近于模型的参数数目(包括截距项)。
• 可视化方面：可以通过plot来展现。同时也可以通过car::subsets来展示

library(leaps)
leap <- regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data = States, nbest = 4)
plot(leap, scale = "adjr2")

通过矫正R平方来查看结果。结果当中越高的模型代表模型越好。

library(car)
subsets(leap, statistic = "cp",legend = FALSE)

##            Abbreviation
## Population            P
## Illiteracy           Il
## Income               In
## Frost                 F

abline(1,1,lty = 2, col = "red")

基于Mallows Cp统计量的四个最佳模型.越好的模型离截距项和斜率均为1的直线越近.图形表明，你可以选择这几个模型，其余可能的模型都可以不予考虑:含Population和Illiteracy的双变量模型;含Population、Illiteracy和Frost的三变量模型，或Population、Illiteracy和Income的三变量模型(它们在图形上重叠了，不易分辨);含Population、Illiteracy、Income和Frost的四变量模型。

PS:大部分情况中，全子集回归要优于逐步回归，因为考虑了更多模型。但是，当有大量预测变量时，全子集回归会很慢。一般来说，变量自动选择应该被看做是对模型选择的一种辅助方法，而不是直接方法。拟合效果佳而没有意义的模型对你毫无帮助，主题背景知识的理解才能最终指引你获得理想的模型。