1、关于线性回归#
1.1 公式理解#
由于实际问题很少遇到单变量线性回归，所以更常见的表示为通用线性模型：
$$
y = \beta_0 + \beta_1x_1 + \beta_2x_2 + …+\beta_kx_k + \epsilon
$$
（1）β0表示截距，即所有预测变量取0时的值；
（2）β1等表示斜率，在其它预测变量不变时，每个预测变量每变动一个单位时输出变量的变化值。
（3）ε表示模型的未知误差，可以理解为残差。
如上公式，线性回归常用普通最小二乘法(OLS)寻找最佳的截距与斜率值的组合，最小化残差(预测值与真实值的差异)平方和。
对于分类变量（k类），会选择1类作为参考水平，同时生成(k-1)个虚拟变量，进而得到(k-1)个斜率值。此时的斜率值表示，该类的预测变量结果相对于参考水平的变化倍数。
1.2 注意点#
（1）线性回归的优势在于模型的可解释性（相对于随机森林、SVM等黑匣子）。即每个预测变量都有斜率值以表示对预测变量的影响；并且不同预测变量的影响是累加的。
（2）对于线性回归结果需要检查拟合值与真实值的残差是否满足两个条件。
- 正态分布：残差值是否符合以0为均值的正态分布。
- 同方差：输出变量的方差不会随着预测拟合值的变化而变化，如下图a的结果是符合的。
如果不满足上面两个条件，提示可能不适合使用通用线性回归，使得模型性能降低，可以尝试之后会学习的通用线性模型。
（3）对于回归问题的性能度量指标
- 平均绝对误差MAE：残差绝对值的均值；
- 均方误差MSE：残差平方和的均值；
- 均方根误差RMSE：MSE的平方根。
2、mlr建模#
2.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
|
data(Ozone, package = "mlbench")
ozoneTib <- Ozone
colnames(ozoneTib) <- c("Month", "Date", "Day", "Ozone", "Press_height",
"Wind", "Humid", "Temp_Sand", "Temp_Monte",
"Inv_height", "Press_grad", "Inv_temp", "Visib")
str(ozoneTib)
# 'data.frame': 366 obs. of 13 variables:
# $ Month : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ Date : Factor w/ 31 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
# $ Day : Factor w/ 7 levels "1","2","3","4",..: 4 5 6 7 1 2 3 4 5 6 ...
# $ Ozone : num 3 3 3 5 5 6 4 4 6 7 ...
# $ Press_height: num 5480 5660 5710 5700 5760 5720 5790 5790 5700 5700 ...
# $ Wind : num 8 6 4 3 3 4 6 3 3 3 ...
# $ Humid : num 20 NA 28 37 51 69 19 25 73 59 ...
# $ Temp_Sand : num NA 38 40 45 54 35 45 55 41 44 ...
# $ Temp_Monte : num NA NA NA NA 45.3 ...
# $ Inv_height : num 5000 NA 2693 590 1450 ...
# $ Press_grad : num -15 -14 -25 -24 25 15 -33 -28 23 -2 ...
# $ Inv_temp : num 30.6 NA 47.7 55 57 ...
# $ Visib : num 200 300 250 100 60 60 100 250 120 120 ...
### 如上`Ozone`列表示每天的臭氧水平，其余是其它相关记录数据。
### 从数据角度来看，还存在一些问题，需要清洗一下：
#(1)Ozone列本身就有NA缺失值，需要去除该记录
#(2)Month等列为因子形，根据后续目的需要转为数值类型
ozoneClean <- mutate_all(ozoneTib, as.numeric) %>%
dplyr::filter(is.na(Ozone) == FALSE)
|
注意到其它列也有缺失值。如果要求比较严苛，可以删去所有包含NA值的观测记录。这里没有删除，因为后面会使用imputation的方式填补缺失值。
2.2 确定预测目标与训练方法#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
|
##(1)首先需要解决NA缺失值的问题
#查看提供有哪些填补方法
?imputations
#如下采用了机器学习(决策树)方法预测缺失值
imputeMethod <- imputeLearner("regr.rpart")
ozoneImp <- impute(as.data.frame(ozoneClean),
classes = list(numeric = imputeMethod))
ozoneImp$data #可以看到已经填补了缺失值的数据
##(2)然后就可以确定预测目标：根据每天的空气相关记录数据，预测那一天的臭氧水平
ozoneTask <- makeRegrTask(data = ozoneImp$data, target = "Ozone")
##(3)使用线性回归模型
lin <- makeLearner("regr.lm")
|
2.3 特征筛选建模#
这一步是去除不具备预测价值的干扰变量，提高模型精度。有如下两种思路：
2.3.1 TopK个预测变量#
这种思路是评价每一个预测变量与待拟合变量的"关系"，有很多评价指标。
1
2
|
#提供有33种评价指标，默认为 randomForestSRC_importance
listFilterMethods()
|
然后根据指标结果，选取最有价值的若干个预测变量。具体的数量，可设置为超参数，选取最优的。
1
2
3
4
5
6
7
8
9
10
11
12
|
##(1)计算每个预测变量与待预测变量的关系指标
如下采用默认randomForestSRC_importance指标，可设置method改为其它指标
filterVals <- generateFilterValuesData(ozoneTask)
head(filterVals$data)
# name type filter value
# 1: Month numeric randomForestSRC_importance 1.36212556
# 2: Date numeric randomForestSRC_importance -0.15006233
# 3: Day numeric randomForestSRC_importance -0.01314694
# 4: Press_height numeric randomForestSRC_importance 0.79176752
# 5: Wind numeric randomForestSRC_importance -0.09471144
# 6: Humid numeric randomForestSRC_importance 2.09659900
plotFilterValues(filterVals) + theme_bw()
|
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
|
##(2)定义超参数，遍历TopK个可能的模型结果，选取最优K值
filterWrapper <- makeFilterWrapper(learner = lin, fw.method = "randomForestSRC_importance")
lmParamSpace <- makeParamSet(
makeIntegerParam("fw.abs", lower = 1, upper = 12)
)
gridSearch <- makeTuneControlGrid() #gridSearch
kFold <- makeResampleDesc("CV", iters = 10)
tunedFeats <- tuneParams(filterWrapper, task = ozoneTask, resampling = kFold,
par.set = lmParamSpace, control = gridSearch)
tunedFeats
# Tune result:
# Op. pars: fw.abs=10
# mse.test.mean=20.7037511
##如上结果表明，保留Top10个预测变量的模型性能最好(mse最小)
##(3)根据上述结果，选取Top10预测变量建模
filteredTask <- filterFeatures(ozoneTask, fval = filterVals,
abs = unlist(tunedFeats$x))
filteredModel <- train(lin, filteredTask)
getLearnerModel(filteredModel) %>% summary()
# Call:
# stats::lm(formula = f, data = d)
#
# Residuals:
# Min 1Q Median 3Q Max
# -14.0584 -2.9348 -0.4092 2.7411 13.7259
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 38.0035915 27.9394230 1.360 0.174638
# Month -0.2819539 0.0781255 -3.609 0.000352 ***
# Day 0.0048944 0.1194186 0.041 0.967331
# Press_height -0.0097813 0.0052630 -1.858 0.063938 .
# Humid 0.0743308 0.0182334 4.077 5.66e-05 ***
# Temp_Sand 0.2182749 0.0508613 4.292 2.30e-05 ***
# Temp_Monte 0.2468545 0.0734535 3.361 0.000863 ***
# Inv_height -0.0004474 0.0002719 -1.646 0.100720
# Press_grad -0.0001682 0.0108708 -0.015 0.987666
# Inv_temp 0.0257259 0.0765361 0.336 0.736976
# Visib -0.0056351 0.0035464 -1.589 0.112970
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 4.479 on 350 degrees of freedom
# Multiple R-squared: 0.6887, Adjusted R-squared: 0.6798
# F-statistic: 77.43 on 10 and 350 DF, p-value: < 2.2e-16
|
2.3.2 遍历特征组合#
相对于上一种方法考虑每个预测变量的价值，这种思路是希望选出性能最好的预测变量组合。
具体搜索最优预测变量组合，有如下方式
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
?FeatSelControl
##穷举法：比较所有可能的预测变量组合
FeatSelControlExhaustive()
##随机法：随机搜索若干个预测变量组合
FeatSelControlRandom()
##顺序搜索
FeatSelControlSequential()
###具体有4种思路
####向前搜索(sfs)：从空模型依次添加可提高性能的特征，直到不能改进为止；
####浮动向前(sffs)：从空模型，在每个步骤添加/删除对模型性能最优改进的特征，直到不再有助于改善模型性能位置
####向后搜索(sbs)：从全特征模型依次删除可提高模型性能的特征，直到不能改进为止；
####浮动向后(sfbs)：从全特征模型模型，在每个步骤添加/删除对模型性能最优改进的特征，直到不再有助于改善模型性能位置
##....
|
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
|
##(1) 搜索、遍历最佳参数组合
featSelControl <- makeFeatSelControlSequential(method = "sffs")
kFold <- makeResampleDesc("CV", iters = 10)
selFeats <- selectFeatures(learner = lin, task = ozoneTask,
resampling = kFold, control = featSelControl)
selFeats
# FeatSel result:
# Features (8): Month, Press_height, Wind, Humid, Temp_Sand, Te...
# mse.test.mean=20.1119257
selFeats$x
# [1] "Month" "Press_height" "Wind" "Humid" "Temp_Sand" "Temp_Monte"
# [7] "Inv_height" "Visib"
##如上结果表明，保留Month等8个预测变量的模型性能最好(mse最小)
##(2) 根据上述结果，选取其中8个预测变量建模
ozoneSelFeat <- ozoneImp$data[, c("Ozone", selFeats$x)]
ozoneSelFeatTask <- makeRegrTask(data = ozoneSelFeat, target = "Ozone")
wrapperModel <- train(lin, ozoneSelFeatTask)
getLearnerModel(wrapperModel) %>% summary()
# Call:
# stats::lm(formula = f, data = d)
#
# Residuals:
# Min 1Q Median 3Q Max
# -13.9342 -2.9498 -0.2835 2.7223 13.8290
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 41.7966703 27.8005624 1.503 0.133620
# Month -0.2966591 0.0782718 -3.790 0.000177 ***
# Press_height -0.0103530 0.0051610 -2.006 0.045620 *
# Wind -0.1225207 0.1285928 -0.953 0.341355
# Humid 0.0764341 0.0149821 5.102 5.52e-07 ***
# Temp_Sand 0.2270547 0.0433972 5.232 2.89e-07 ***
# Temp_Monte 0.2665338 0.0636193 4.190 3.54e-05 ***
# Inv_height -0.0004740 0.0001854 -2.557 0.010987 *
# Visib -0.0052262 0.0035579 -1.469 0.142749
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 4.462 on 352 degrees of freedom
# Multiple R-squared: 0.6894, Adjusted R-squared: 0.6823
# F-statistic: 97.65 on 8 and 352 DF, p-value: < 2.2e-16
|
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
|
#先封装缺失值模拟
imputeMethod <- imputeLearner("regr.rpart")
imputeWrapper <- makeImputeWrapper(lin,
classes = list(numeric = imputeMethod))
#再封装特征选择
featSelControl <- makeFeatSelControlSequential(method = "sfbs")
kFold <- makeResampleDesc("CV", iters = 10)
featSelWrapper <- makeFeatSelWrapper(learner = imputeWrapper,
resampling = kFold, #内部循环
control = featSelControl)
library(parallel)
library(parallelMap)
ozoneTaskWithNAs <- makeRegrTask(data = ozoneClean, target = "Ozone")
kFold3 <- makeResampleDesc("CV", iters = 3)
parallelStartSocket(cpus = detectCores())
lmCV <- resample(featSelWrapper, ozoneTaskWithNAs, resampling = kFold3) #~ 1.5 min
parallelStop()
lmCV
# Resample Result
# Task: ozoneClean
# Learner: regr.lm.imputed.featsel
# Aggr perf: mse.test.mean=21.2937631
# Runtime: 52.349
|
2.5 模型诊断#
1
2
3
4
5
6
|
wrapperModelData <- getLearnerModel(wrapperModel)
summary(wrapperModelData)
par(mfrow = c(2, 2))
plot(wrapperModelData)
par(mfrow = c(1, 1))
|
如下图左边的两个子图用于检查残差的同方差性，右上QQ图用于检查残差的正态分布。
具体来看残差随着预测值的增加而变多，表明存在异方差关系；而QQ图表明残差符合正太分布。
