《精通机器学习：基于R 第二版》学习笔记
setwd("C:/Users/Admin/Documents/R")
library(pacman)
p_load(dplyr,ggplot2)
1、单变量线性回归
# 怀俄明州蛇河流域的水量数据
snake <- tribble(~id,~x,~y,
1, 23.1, 10.5,
2, 32.8, 16.7,
3, 31.8, 18.2,
4, 32.0, 17.0,
5, 30.4, 16.3,
6, 24.0, 10.5,
7, 39.5, 23.1,
8, 24.2, 12.4,
9, 52.5, 24.9,
10, 37.9, 22.8,
11, 30.5, 14.1,
12, 25.1, 12.9,
13, 12.4, 8.8,
14, 35.1, 17.4,
15, 31.5, 14.9,
16, 21.1, 10.5,
17, 27.6, 16.1)
# 更改有意义的变量名
names(snake) <- make.names(c("id","content","yield"))
head(snake)
## # A tibble: 6 x 3
## id content yield
## <dbl> <dbl> <dbl>
## 1 1 23.1 10.5
## 2 2 32.8 16.7
## 3 3 31.8 18.2
## 4 4 32 17
## 5 5 30.4 16.3
## 6 6 24 10.5
1.1 散点图，可以看到前后有两个明显的离群点
ggplot(snake,aes(content,yield)) +
geom_point() +
geom_smooth(method = "lm",se=F) +
labs(x="water content id snow",y="water yield") +
theme_bw()
1.2 线性回归
yield.fit <- lm(yield ~ content,snake)
summary(yield.fit)
##
## Call:
## lm(formula = yield ~ content, data = snake)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1793 -1.5149 -0.3624 1.6276 3.1973
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.72538 1.54882 0.468 0.646
## content 0.49808 0.04952 10.058 4.63e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.743 on 15 degrees of freedom
## Multiple R-squared: 0.8709, Adjusted R-squared: 0.8623
## F-statistic: 101.2 on 1 and 15 DF, p-value: 4.632e-08
理论上，Multiple R-squared 的取值范围在0和1之间，用来表示X和Y的相关程度。在本例中，它的意义是水量87%的方差可以被降雪含水量解释。顺便说一句，R平方项就是[X, Y]的相关系数的平方。
线性回归必须通过假设检验，其中的假设可以总结如下：
线性：预测变量与响应变量之间的关系是线性的。如果线性关系不能清晰呈现，可以对变量X或Y进行数据转换（对数转换、多项式转换、指数转换等）以解决问题。
误差不相关：在时间序列和面板数据中，如果误差是相关的，那么就有可能建立一个非常不规范的模型。
同方差性：误差是正态分布的，并具有相同的方差。这意味着对于不同的输入值，误差的方差是个固定值。如果违背了这个假设，参数估计就有可能产生偏差，导致对显著性的统计检验结果过高或者过低，从而得到错误的结论。这种情况就称为异方差性。
非共线性：两个预测变量之间不存在线性关系，也就是说，特征之间不应该存在相关性。同样地，共线性也会导致估计偏差。
不存在异常值：异常值会严重影响参数估计。理想情况下，必须在使用线性回归拟合模型之前就除去异常值。
1.3 假设检验
par(mfrow=c(2,2))
plot(yield.fit)
正态Q-Q图（Normal Q-Q）可以帮助我们确定残差是否服从正态分布。
分位数-分位数图（Scale-Location）表示一个变量的分位数对应于另一个变量的分位数画出的图，从图中可以看出有离群点（第7、9、10个观测），这可能违反假设。
残差杠杆图（Residuals vs Leverage）可以告诉我们哪个观测值（如果有）会对模型造成过度影响，换句话说，是否存在我们应该关注的异常值。鉴别强影响点的统计量是库克距离，一般认为，如果这个统计量的值大于1（第9个观测），就需要进行更深入的检查。
离群点如果靠删除来解决的话，删除后也可能引入其他的离群点，因为删除观测值，可能导致其他的观测值因为库克距离大于1而落到正常区间之外。
1.4 使用qqplot检查置信区间，确定残差服从正态分布
car::qqPlot(yield.fit)
## [1] 7 10
2、多变量线性回归
p_load(alr3,corrplot)
data(water)
str(water)
## 'data.frame': 43 obs. of 8 variables:
## $ Year : int 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 ...
## $ APMAM : num 9.13 5.28 4.2 4.6 7.15 9.7 5.02 6.7 10.5 9.1 ...
## $ APSAB : num 3.58 4.82 3.77 4.46 4.99 5.65 1.45 7.44 5.85 6.13 ...
## $ APSLAKE: num 3.91 5.2 3.67 3.93 4.88 4.91 1.77 6.51 3.38 4.08 ...
## $ OPBPC : num 4.1 7.55 9.52 11.14 16.34 ...
## $ OPRC : num 7.43 11.11 12.2 15.15 20.05 ...
## $ OPSLAKE: num 6.47 10.26 11.35 11.13 22.81 ...
## $ BSAAM : int 54235 67567 66161 68094 107080 67594 65356 67909 92715 70024 ...
2.1 相关系数又称Pearson’s r，可以用来测量两个变量之间线性相关性的强度和方向
socal.water <- water[,-1]
socal.water %>% cor() %>% corrplot(method = "number")
或者直接画图观察相关性：
socal.water %>% pairs(~ .,data = .)
2.2 构建模型与模型评价
对于一个数据集，你先用前向逐步回归，然后再用后向逐步回归，可能会得到两个完全矛盾的模型。最重要的一点是，逐步回归会使回归系数发生偏离，换句话说，会使回归系数的值过大，置信区间过窄。
对于特征选择，最优子集回归是逐步回归的一个可接受的替代方案。在最优子集回归中，算法使用所有可能的特征组合来拟合模型，所以，如果有3个特征，将生成7个模型。然后和逐步回归一样，分析者需要应用自己的判断和统计分析来选择最优模型，模型选择就是后面工作的关键。
正如你猜想的那样，如果数据集有多个特征，工作量就会非常大。当特征数多于观测数时（p大于n），这个方法的效果就不会好。
2.3 使用所有特征构建线性模型
socal.water <- water[,-1]
fit <- lm(BSAAM ~ .,data = socal.water)
summary(fit)
##
## Call:
## lm(formula = BSAAM ~ ., data = socal.water)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12690 -4936 -1424 4173 18542
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15944.67 4099.80 3.889 0.000416 ***
## APMAM -12.77 708.89 -0.018 0.985725
## APSAB -664.41 1522.89 -0.436 0.665237
## APSLAKE 2270.68 1341.29 1.693 0.099112 .
## OPBPC 69.70 461.69 0.151 0.880839
## OPRC 1916.45 641.36 2.988 0.005031 **
## OPSLAKE 2211.58 752.69 2.938 0.005729 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7557 on 36 degrees of freedom
## Multiple R-squared: 0.9248, Adjusted R-squared: 0.9123
## F-statistic: 73.82 on 6 and 36 DF, p-value: < 2.2e-16
可以看到，OPRC 和 OPSLAKE 这两个参数具有显著的p值。有趣的是，OPBPC并不显著，尽管它与响应变量高度相关。简言之，当我们控制其他OP开头的特征时，OPBPC 无法对预测方差提供任何有意义的解释。这就是说，模型中存在 OPRC 和 OPSLAKE 时，特征 OPBPC 从统计学角度来看没有任何作用。
2.4 使用最优子集法
sub.fit <- leaps::regsubsets(BSAAM ~ .,data=socal.water)
best.summary <- summary(sub.fit)
names(best.summary)
## [1] "which" "rsq" "rss" "adjr2" "cp" "bic" "outmat" "obj"
which.min(best.summary$rss)
## [1] 6
以上代码告诉我们，有6个特征的模型具有最小的RSS。本应如此，因为它有最多的输入，输入越多，RSS越小。请注意，增加特征必然会减少RSS！而且必然会增加R方。
2.5 比较CP和BIC
par(mfrow=c(1,2))
plot(best.summary$cp,xlab = "number of features",ylab="cp")
plot(sub.fit,scale = "Cp")
在左侧的图中，可以看出有3个特征的模型具有最小的Cp值。在右侧的图中，显示了能给出最小Cp的特征组合。这张图应该这么看：先在Y轴的最高点找到最小的Cp值，此处是1.2；然后向右在X轴上找到与之对应的色块。通过这张图，我们可以看到这个具有最小Cp值的模型中的3个特征是 APSLAKE 、OPRC 和 OPSLAKE 。
2.6 Cp与BIC和修正R方的比较
which.min(best.summary$bic)
## [1] 3
which.max(best.summary$adjr2)
## [1] 3
可以看出，在本例中，BIC、修正R方和Cp选择的最优模型是一致的。
2.7 假设检验
best.fit <- lm(BSAAM ~ APSLAKE + OPRC + OPSLAKE,data = socal.water)
summary(best.fit)
##
## Call:
## lm(formula = BSAAM ~ APSLAKE + OPRC + OPSLAKE, data = socal.water)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12964 -5140 -1252 4446 18649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15424.6 3638.4 4.239 0.000133 ***
## APSLAKE 1712.5 500.5 3.421 0.001475 **
## OPRC 1797.5 567.8 3.166 0.002998 **
## OPSLAKE 2389.8 447.1 5.346 4.19e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7284 on 39 degrees of freedom
## Multiple R-squared: 0.9244, Adjusted R-squared: 0.9185
## F-statistic: 158.9 on 3 and 39 DF, p-value: < 2.2e-16
2.8诊断图
par(mfrow=c(2,2))
plot(best.fit)
通过这4张图，我们完全可以认为，残差具有固定的方差并且服从正态分布。杠杆图中也没有什么需要我们进一步处理的异常。
2.9 多重共线性
car::vif(best.fit)
## APSLAKE OPRC OPSLAKE
## 1.011499 6.452569 6.444748
一般认为， VIF （方差膨胀因子）值超过 5 （有些人认为是 10 ）就说明存在严重的共线性。所以 OPRC 和 OPSLAKE 存在潜在的共线性问题（ VIF 值大于 5 ）。
从下图也可以看出，两个变量之间高度相关：
ggplot(socal.water,aes(OPRC,OPSLAKE)) +
geom_point() +
theme_bw()
看一下最优子集法中生成的修正 R 方的值就会发现，APSLAKE 和 OPSLAKE 组成的两变量模型的值为 0.90，加入OPRC 之后仅有一个微不足道的提升，到了0.92。
best.summary$adjr2
## [1] 0.8777515 0.9001619 0.9185369 0.9168706 0.9146772 0.9123079
fit.2 <- lm(BSAAM ~ APSLAKE + OPSLAKE,data = socal.water)
summary(fit.2)
##
## Call:
## lm(formula = BSAAM ~ APSLAKE + OPSLAKE, data = socal.water)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13335.8 -5893.2 -171.8 4219.5 19500.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19144.9 3812.0 5.022 1.1e-05 ***
## APSLAKE 1768.8 553.7 3.194 0.00273 **
## OPSLAKE 3689.5 196.0 18.829 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8063 on 40 degrees of freedom
## Multiple R-squared: 0.9049, Adjusted R-squared: 0.9002
## F-statistic: 190.3 on 2 and 40 DF, p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(fit.2)
car::vif(fit.2)
## APSLAKE OPSLAKE
## 1.010221 1.010221
模型是显著的，诊断图中也没发现什么问题，VIF小于5，共线性问题应该得到了解决。
2.10 同方差检验
BP检验：原假设是误差方差为 0，对应的备择假设是误差方差不为0。
lmtest::bptest(fit.2)
##
## studentized Breusch-Pagan test
##
## data: fit.2
## BP = 0.0046205, df = 2, p-value = 0.9977
我们没有证据拒绝认为“误差方差为 0 ”的原假设，因为p值 = 0.9977 。检验结果中， BP = 0.0046是一个卡方值。
2.11 使用该模型生成的预测值与实际值的散点图
tibble(predict=fit.2$fitted.values,actual=socal.water$BSAAM) %>% ggplot(aes(predict,actual)) +
geom_point() +
geom_smooth(method = lm) +
theme_bw()
2.12 留一法交叉验证检测预测误差平方和
MPV::PRESS(best.fit)
## [1] 2426757258
MPV::PRESS(fit.2)
## [1] 2992801411
如果仅参考这个统计量，我们应该选择模型 best.fit。但我还是认为更简约的模型更好。
2.13 交互项
当一个特征的预测效果依赖于另一个特征的值时，这两个特征就是交互作用的。
p_load(MASS)
data(Boston)
str(Boston)
## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
lstat 表示低社会经济地位家庭百分比
age表示房龄年数
在 lm() 函数中使用 feature1*feature2 ，将两个特征及其交互项加入模型。
value.fit <- lm(medv ~ lstat * age,data = Boston)
summary(value.fit)
##
## Call:
## lm(formula = medv ~ lstat * age, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.806 -4.045 -1.333 2.085 27.552
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.0885359 1.4698355 24.553 < 2e-16 ***
## lstat -1.3921168 0.1674555 -8.313 8.78e-16 ***
## age -0.0007209 0.0198792 -0.036 0.9711
## lstat:age 0.0041560 0.0018518 2.244 0.0252 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.149 on 502 degrees of freedom
## Multiple R-squared: 0.5557, Adjusted R-squared: 0.5531
## F-statistic: 209.3 on 3 and 502 DF, p-value: < 2.2e-16
检查输出可以知道，社会经济地位是个具有高预测性的特征，而房龄则不是。但这两个变量具有显著的交互作用，可以对房屋价值进行正向解释。
