Correlation and Regression
目录
- 写在前面
- 数据的可视化
- 相关分析
- Anscombe
- 简单的线性回归
- 线性回归的简单可视化
- 回归诊断
- 标准方式
- car包的使用
- 线性模型假设的综合验证
- 线性回归统计
- 进行预测
- 模型拟合评价
- 残差标准误(Residual standard error)
- R平方(R squared)
- 异常点
- 模型比较
- 嵌套模型比较
- AIC比较
- 最优模型选择
- 逐步回归
- 全子集回归
这里介绍的主要是datacamp中关于相关和回归以及R语言之战中关于回归的综合学习笔记
如果要研究两个连续性变量之间的关系，一般来说就是使用相关分析和回归分析。回归分析的话，最基础的还是线性回归。在datacamp中主要介绍的还是线性回归。
- 加载所需要的数据集
library(openintro)
## Please visit openintro.org for free statistics materials
##
## Attaching package: 'openintro'
## The following objects are masked from 'package:datasets':
##
## cars, trees
library(HistData)
library(tidyverse)
## ─ Attaching packages ────────────────── tidyverse 1.2.1 ─
## ✔ ggplot2 3.1.0 ✔ purrr 0.3.2
## ✔ tibble 2.1.1 ✔ dplyr 0.8.0.1
## ✔ tidyr 0.8.3 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ─ Conflicts ─────────────────── tidyverse_conflicts() ─
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
data(anscombe)
data(mammals)
data(mlbBat10)
data(possum)
data(ncbirths)
data(bdims)
data(GaltonFamilies)
数据的可视化
在进行分析之前是首先对数据进行基本的了解。对于两个连续性变量如果要最基本的可视化就是散点图。进一步的，可以通过cut函数对散点图进行切割转换为箱式图。 - 通过可视化除了观察数据的分布，也可以看到异常值。 - cut函数接受breaks参数可以吧连续性变量分成多个部分进行箱式图可视化
p1 <- ggplot(ncbirths, aes(weeks, weight)) + geom_point()
p2 <- ggplot(data = ncbirths,
aes(x = cut(weeks, breaks = 5), y = weight)) +
geom_boxplot()
cowplot::plot_grid(p1, p2, ncol = 2)
## Warning: Removed 2 rows containing missing values (geom_point).
!upload-images.jianshu.io/upload_images/10631927-784333f01cef80db.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
- 可视化的实话如何数据跨度特别大的时候，可以通过转换坐标轴开可视化数据
p1 <- ggplot(mammals, aes(BodyWt, BrainWt)) + geom_point()
p2 <- ggplot(data = mammals, aes(x = BodyWt, y = BrainWt)) +
geom_point() +
scale_x_log10() + scale_y_log10()
cowplot::plot_grid(p1, p2, ncol = 2)
image.png
相关分析
-
通过图观察完两个数据的分布之后，我们需要在数据上分析两个数据是否相关。这个时候就要用到相关分析了。基本的相关分析的参数是
cor。其接受两个相关分析的数据。另外的话，相关分析对于NA的处理是很保守的。如果数据中含有NA的话，则会返回NA。如果需要分析的时候需要自动去掉NA，可以使用use参数。只需要设置为use = "pairwise.complete.obs"即可。 -
两个连续性数据相关程度可以通过相关系数来查看。对于相关系数的解释是：一般是一个-1到1的数值。其中正负代表方向。数值的绝对值代表相关程度。越大相关性越高。另外相关分析得到的相关系数只是线性的相关程度。
# Compute correlation
ncbirths %>%
summarize(N = n(), r = cor(weeks, weight),
r_nona = cor(weeks, weight, use = "pairwise.complete.obs"))
## N r r_nona
## 1 1000 NA 0.6701013
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)
image.png
简单的线性回归
简单的线性回归模型我们通过lm函数就可以得到。
线性回归的简单可视化
简单的线性回归我们可以通过geom_smooth(method = "lm")进行可视化。同时也可以通过geom_abline来截距(intercept)和系数(hgt)
其中
ggplot(data = bdims, aes(x = hgt, y = wgt)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)
image.png
回归诊断
在得到预测模型之后，我们要评判模型是否成立。最基本的是通过plot来进行预测判断的。
标准方式
mod <- lm(wgt ~ hgt, data = bdims)
par(mfrow = c(2,2))
plot(mod)
image.png
模型诊断主要是分为四个方面：正态性；独立性；线性；方差齐性。结合上图
- 左上图：观察数据是否是线性，如果是线性则会发现残差和拟合值之间的分布不存在任何的关系。如上图就是存在线性关系
- 右上图：观察数据是否服从正态分布，一般所有的点都在对角线上可以说明是服从正态分布。
- 左下图：观察数据是否是方差齐性，如果数据满足假设，那么数据会随着水平线随机分布。
- �右下图：观察数据是是否存在异常点。具体异常点的介绍可以看下面。
car包的使用
但是上图其实有些数据显示的并不是很好。比如我们想要观察的独立性。所以我们可以使用car包的函数就进行更好的数据发
- 正态性：与基础包中的plot()函数相比，qqPlot()函数提供了更为精确的正态假设检验方法
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)
image.png
## [1] 124 474
- 独立性：car包提供了一个可做
Durbin-Watson检验的函数，能够检测误差的序列相关性。结果中如果p-value> 0.代表样本之间存在独立性。
durbinWatsonTest(mod)
## lag Autocorrelation D-W Statistic p-value
## 1 0.0519505 1.894424 0.206
## Alternative hypothesis: rho != 0
- 线性：通过成分残差图(component plus residual plot)也称偏残差图(partial residual plot)，你可以看看因变量与自变量之间是否呈非线性关系，也可以看看是否有不同于已设定线性模型的系统偏差。
crPlots(mod)
image.png
- 方差齐性：car包提供了两个有用的函数，可以判断误差方差是否恒定。
ncvTest函数生成一个计分 检验，零假设为误差方差不变。因此如果p-value> 0.05则方差是齐的。另外还有一个函数为spreadLevelPlot是创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值 与拟合值的关系。
ncvTest(mod)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.1361065, Df = 1, p = 0.71218
spreadLevelPlot(mod)
image.png
##
## Suggested power transformation: 0.6747609
线性模型假设的综合验证
通过使用gvlma包中的gvlma函数可以对线性模型假设进行综合检验。 - 其中结果当中Global Stat代表总体评价。如果都通过则p > 0.05。如果没有。则可以通过上面的car包判断那些没有通过。
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.
线性回归统计
-
进行线性回归的基本参数是
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
异常点
异常点通常分为三种:
-
异常点：那些距离回归直线很远，通常这些点的
残差都会很大。一般来说
- 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
- 高杠杆点：这样的点拥有很极端的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
- 强影响点：有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/(nk 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")
image.png
-
综合评估:利用
car::influencePlot函数，你还可以将离群点、杠杆值和强影响点的 信息整合到一幅图形中
- 纵坐标超过+2或小于2的州可被认为是离群点，水平轴超过0.2或0.3的有高杠杆值(通常为预测值的组合)。圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点
influencePlot(mod)
image.png
## 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")
image.png
通过矫正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")
image.png
基于Mallows Cp统计量的四个最佳模型.越好的模型离截距项和斜率均为1的直线越近.图形表明，你可以选择这几个模型，其余可能的模型都可以不予考虑:含Population和Illiteracy的双变量模型;含Population、Illiteracy和Frost的三变量模型，或Population、Illiteracy和Income的三变量模型(它们在图形上重叠了，不易分辨);含Population、Illiteracy、Income和Frost的四变量模型。
PS:大部分情况中，全子集回归要优于逐步回归，因为考虑了更多模型。但是，当有大量预测变量时，全子集回归会很慢。一般来说，变量自动选择应该被看做是对模型选择的一种辅助方法，而不是直接方法。拟合效果佳而没有意义的模型对你毫无帮助，主题背景知识的理解才能最终指引你获得理想的模型
网友评论