1、关于GAM非线性回归

(1) n阶多项式

  • 如前所说，线性回归的假设是每个预测变量与输出变量之间为线性相关。即类似 y = ax + b。

  • 当预测变量与输出变量之间为非线性相关，即呈曲线特征时，可尝试使用高阶多项式进行拟合。 $$ y = \beta_0 + \beta_1x + \beta_2x^2 + …+\beta_nx^n + \varepsilon $$ 正交多项式介绍及应用- 优化与算法- 博客园

  • 一般拟合高阶多项式时，除了最高阶n外，还会包括所有低阶项(1,2,3… n-1)次幂。目的是为了避免最值点必须处于x=0的位置。

  • 样条曲线是分段的多项式函数：将预测变量分为若干区域，每个区域单独拟合一个多项式曲线。相邻的多项式之间的连接成为knot。

(2) GAM

  • GAM会尝试将每个预测变量与输出变量关系拟合为平滑曲线，通常是多个样条曲线的组合，其中每一个样条曲线称为基函数。

    • 对于预测变量Xi，平滑函数可以表示为如下。a表示每个基函数的权重。

    $$ f(x_i) = a_1b_1(x_i)+a_2b_2(x_i)+…+a_nb_n(x_n) $$

  • 所以GAM广义加性模型公式表示为： $$ y = \beta_0 + f_1(x_1) + f_2(x_2) + … + f_n(x_n) + \varepsilon $$

上一节提到的通用线性模型公式为： $$ y = \beta_0 + \beta_1x_1 + \beta_2x_2 + … + \beta_nx_n + \varepsilon $$ 可以这么理解：通用线性模型是广义加性模型的特例。

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)使用gam非线性回归模型
gam <- makeLearner("regr.gamboost")

2.3 特征筛选建模

参考上一节，有两种筛选思路。如下选取第二种筛选最佳组合的方式

 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

##(1) 搜索、遍历最佳参数组合
featSelControl <- makeFeatSelControlSequential(method = "sffs")
kFold <- makeResampleDesc("CV", iters = 10)
selFeats <- selectFeatures(learner = gam, task = ozoneTask, 
                           resampling = kFold, control = featSelControl)
selFeats
# FeatSel result:
# Features (7): Month, Day, Press_height, Wind, Temp_Monte, Pre...
# mse.test.mean=15.0966815
selFeats$x
# [1] "Month"        "Day"          "Press_height" "Wind"         "Temp_Monte"  
# [6] "Press_grad"   "Visib"
##如上结果表明，保留Month等7个预测变量的模型性能最好(mse最小)

##(2) 根据上述结果，选取其中7个预测变量建模,并可视化
ozoneSelFeat <- ozoneImp$data[, c("Ozone", selFeats$x)]
ozoneSelFeatTask <- makeRegrTask(data = ozoneSelFeat, target = "Ozone")
gamModel <- train(gam, ozoneSelFeatTask)
gamModelData <- getLearnerModel(gamModel)
summary(gamModelData)

par(mfrow = c(3, 3))
plot(gamModelData, type = "l")
plot(gamModelData$predict(), resid(gamModelData))
qqnorm(resid(gamModelData))
qqline(resid(gamModelData))
par(mfrow = c(1, 1))
image-20220428222622566

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(gam, 
                                   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 = 10)
parallelStartSocket(cpus = detectCores())
gamCV <- resample(featSelWrapper, ozoneTaskWithNAs, resampling = kFold3) #~ 1.5 min

parallelStop()
gamCV
# Resample Result
# Task: ozoneClean
# Learner: regr.gamboost.imputed.featsel
# Aggr perf: mse.test.mean=15.9090278
# Runtime: 264.97