【DS】R的整洁交叉验证教程

笔者邀请您，先思考：

1 超参数如何调整和优化？

介绍

这篇文章将使用{tidymodels}软件包集合中的几个软件包，即{recipes}，{rsample}和{parsnip}以整洁的方式来训练随机森林。 我还将使用{mlrMBO}来调整随机森林的超参数。

设置

让我们加载所需要的包：

1library("tidyverse")
2library("tidymodels")
3library("parsnip")
4library("brotools")
5library("mlbench")

加载包含在mlrbench包中的数据:

1data("BostonHousing2")

我将训练一个随机森林来预测房价，这是cmedv列:

1head(BostonHousing2)

仅保留相关列

1boston <- BostonHousing2 %>% 
2    select(-medv, -town, -lon, -lat) %>% 
3    rename(price = cmedv)

我移除town, lat和lon列，因为tract列中包含的信息已经足够了。

为了训练和评估模型的性能，我将数据分成两部分。 一个数据集，我称之为训练集，将在下面进一步分为两个。 我不会触摸第二个数据集，测试集，直到最后。

1train_test_split <- initial_split(boston, prop = 0.9)
2housing_train <- training(train_test_split)
3housing_test <- testing(train_test_split)

我想训练一个随机森林来预测房屋的价格，但是随机森林有所谓的超参数，这些参数是无法从数据中估算或学习的参数。相反，这些参数必须由分析师选择。为了选择它们，您可以使用看起来表现良好的文献中的值（如在宏观经济学中完成）或者您可以进一步将训练集拆分为两个，创建超参数网格，在一部分数据中对网格中的所有值训练模型，并在第二部分数据上对比模型的预测。然后，您将坚持使用性能最佳的模型，例如，具有最低RMSE的模型。问题是，您无法仅使用一个值来估计RMSE的真实值。这就像你想通过从人口中抽取一个单一的观察来估计人口的高度。你需要更多的观察。为了给出一组超参数的RMSE的真实值，而不是做一次拆分，我会做30次，然后计算平均RMSE，这意味着我对为超参数的每个值的组合训练30个模型有兴趣。

首先，让我们使用{rsample}包中的mc_cv()函数再次分割训练数据。该函数实现蒙特卡罗交叉验证:

1validation_data <- mc_cv(housing_train, prop = 0.9, times = 30)
2

validation_data是什么样的？

1validation_data

让我们进一步往下看：

1validation_data$splits[[1]]

第一个值是第一个数据集的行数，第二个值是第二个数据集的行数，第三个值是原始数据集的行数。

我们应该如何称呼这两个新数据集？ {rsample}的作者Max Kuhn谈到了分析和评估集：

现在，为了继续，我需要预处理数据。 我将分三步完成。 第一步和第二步用于标准化数值变量，第三步将字符和因子变量转换为虚拟变量。这是必要的，因为我将训练一个随机森林，它不能直接处理因子变量。 让我们定义一个配方来做到这一点，并从预处理测试集开始。 我在配方周围编写了一个包装函数，因为我需要将这个配方应用于各种数据集：

1simple_recipe <- function(dataset){
2    recipe(price ~ ., data = dataset) %>%
3        step_center(all_numeric()) %>%
4        step_scale(all_numeric()) %>%
5        step_dummy(all_nominal())
6}

定义配方后，我可以使用prep（）函数，该函数从数据中估算处理数据所需的参数。 例如，对于中心化，prep（）估计平均值，然后从变量中减去平均值。 使用bake（），然后将估算值应用于数据：

1testing_rec <- prep(simple_recipe(housing_test), testing = housing_test)
2test_data <- bake(testing_rec, newdata = housing_test)
3

在使用prep（）和bake（）之前分割数据很重要，因为如果没有，您将使用prep（）步骤中测试集的观察结果，从而将测试集中的知识引入训练数据中。 这称为数据泄漏，必须避免。这就是为什么有必要首先将训练数据分成分析集和评估集，然后分别预处理这些数据集。 但是，validation_data对象现在不能与recipe（）一起使用，因为它不是数据框。 不用担心，我只需编写一个函数，从validation_data对象中提取分析集和评估集，应用预处理，训练模型，并返回RMSE。 这将是一个重要的函数，是分析的核心。

但在此之前，让我们运行一个简单的线性回归作为基准。 对于线性回归，我不会使用任何CV，所以让我们预处理训练集：

1trainlm_rec <- prep(simple_recipe(housing_train), testing = housing_train)
2trainlm_data <- bake(trainlm_rec, newdata = housing_train)
3linreg_model <- lm(price ~ ., data = trainlm_data)
4broom::augment(linreg_model, newdata = test_data) %>% 
5    rmse(price, .fitted)
6

broom :: augment（）将预测以新列.fitted添加到test_data中。 我不会在随机森林中使用这个技巧，因为我将使用{ranger}的随机森林没有augment（）方法。 我会自己将预测添加到数据中。 我不会在随机森林中使用这个技巧，因为我将使用{ranger}的随机森林没有augment（）方法。 我会自己将预测添加到数据中。

好的，现在让我们回到随机森林并编写这个强大的函数：

 1my_rf <- function(mtry, trees, split, id){
 2
 3    analysis_set <- analysis(split)
 4
 5    analysis_prep <- prep(simple_recipe(analysis_set), training = analysis_set)
 6
 7    analysis_processed <- bake(analysis_prep, newdata = analysis_set)
 8
 9    model <- rand_forest(mtry = mtry, trees = trees) %>%
10        set_engine("ranger", importance = 'impurity') %>%
11        fit(price ~ ., data = analysis_processed)
12
13    assessment_set <- assessment(split)
14
15    assessment_prep <- prep(simple_recipe(assessment_set), testing = assessment_set)
16
17    assessment_processed <- bake(assessment_prep, newdata = assessment_set)
18
19    tibble::tibble("id" = id,
20        "truth" = assessment_processed$price,
21        "prediction" = unlist(predict(model, new_data = assessment_processed)))
22}
23

rand_forest（）函数可从{parsnip}包中获得。 该软件包为许多其他机器学习包提供了统一的接口。 这意味着不必学习range（）和randomForest（）的语法，而且你可以简单地使用rand_forest（）函数并将引擎参数更改为你想要的那个（ranger，randomForest等）。

我们试试这个函数：

1results_example <- map2_df(.x = validation_data$splits,
2                           .y = validation_data$id,
3                           ~my_rf(mtry = 3, trees = 200, split = .x, id = .y))
4

1head(results_example)

我现在可以在mtry = 3和trees = 200时计算RMSE：

1results_example %>%
2    group_by(id) %>%
3    rmse(truth, prediction) %>%
4    summarise(mean_rmse = mean(.estimate)) %>%
5    pull

随机森林的RMSE已经低于线性回归。 现在的目标是通过调整mtry和trees超参数来降低此RMSE。 为此，我将使用{mlrMBO}包中实现的贝叶斯优化方法。

贝叶斯超参数优化

我将重用上面的代码，并定义一个函数，该函数执行从预处理到返回我希望通过调整超参数以最小化RMSE的度量：

 1tuning <- function(param, validation_data){
 2
 3    mtry <- param[1]
 4    trees <- param[2]
 5
 6    results <- purrr::map2_df(.x = validation_data$splits,
 7                       .y = validation_data$id,
 8                       ~my_rf(mtry = mtry, trees = trees, split = .x, id = .y))
 9
10    results %>%
11        group_by(id) %>%
12        rmse(truth, prediction) %>%
13        summarise(mean_rmse = mean(.estimate)) %>%
14        pull
15}
16

这正是之前的代码，但它现在返回RMSE。 让我们尝试使用之前的值的函数：

1tuning(c(3, 200), validation_data)

让我们画出mtry=3和trees从200到300的RMSE值，这需要一些时间，因为我需要100次评估这个成本函数。 如果评估函数是简单的，我可以通过改变mtry的值来制作3D图，但是如果评估函数的简单，我会进行详尽的网格搜索以找到超参数而不是使用贝叶斯优化。

1plot_points <- crossing("mtry" = 3, "trees" = seq(200, 300))
2
3plot_data <- plot_points %>% 
4    mutate(value = map_dbl(seq(200, 300), ~tuning(c(3, .), validation_data)))

1plot_data %>% 
2    ggplot(aes(y = value, x = trees)) + 
3    geom_line(colour = "#82518c") + 
4    theme_blog() +
5    ggtitle("RMSE for mtry = 3")

对于mtry = 3，最小值似乎在255左右。最小化的函数根本不是平滑的。

我现在按照arxiv论文中的代码来运行优化。 我想我得到了论文的要点，但我还不了解一切。 目前，我目前仍在试验该库，但据我所知，一个更简单的模型，称为代理模型，用于寻找有希望的点，并在这些点评估函数的价值。 这似乎与Gourieroux，Monfort，Renault中描述的间接推理方法有些相似（在精神上）。

让我们首先加载包并创建要优化的函数：

1library("mlrMBO")
2fn <- makeSingleObjectiveFunction(name = "tuning",
3                                 fn = tuning,
4                                 par.set = makeParamSet(makeIntegerParam("x1", lower = 3, upper = 8),
5                                                        makeIntegerParam("x2", lower = 50, upper = 500)))

此函数基于我之前定义的函数。 要优化的参数也被定义为它们的界限。 我在在3到8搜索mtry,在50到500搜索trees。

现在是我没有完全理解的部分。

1# Create initial random Latin Hypercube Design of 10 points
2library(lhs)# for randomLHS
3des <- generateDesign(n = 5L * 2L, getParamSet(fn), fun = randomLHS)

我认为这意味着这10个点是用于启动整个过程的点。 我不明白为什么他们必须从超立方体中采样，但没关系。 然后我选择代理模型，随机森林，并预测标准误差。 在这里，我也不太明白为什么标准错误可以是一个选项。

1# Specify kriging model with standard error estimation
2surrogate <- makeLearner("regr.ranger", predict.type = "se", keep.inbag = TRUE)

这里我定义了一些选项:

1# Set general controls
2ctrl <- makeMBOControl()
3ctrl <- setMBOControlTermination(ctrl, iters = 10L)
4ctrl <- setMBOControlInfill(ctrl, crit = makeMBOInfillCritEI())

这是优化部分:

1# Start optimization
2result <- mbo(fn, des, surrogate, ctrl, more.args = list("validation_data" = validation_data))
3result

因此推荐的参数mtry为6，trees为381。 RMSE的值低于之前，等于0.393。 现在让我们用这个值训练训练数据上的随机森林。 首先，我预处理训练数据：

1training_rec <- prep(simple_recipe(housing_train), testing = housing_train)
2
3train_data <- bake(training_rec, newdata = housing_train)
4

现在让我们训练我们的最终模型并预测价格:

1final_model <- rand_forest(mtry = 6, trees = 381) %>%
2        set_engine("ranger", importance = 'impurity') %>%
3        fit(price ~ ., data = train_data)
4
5price_predict <- predict(final_model, new_data = select(test_data, -price))

让我们将数据转换回来，并将预测的价格与真实的价格进行直观比较：

1cbind(price_predict * sd(housing_train$price) + mean(housing_train$price), 
2      housing_test$price)

现在我们来计算RMSE:

1tibble::tibble("truth" = test_data$price,
2        "prediction" = unlist(price_predict)) %>% 
3    rmse(truth, prediction)

非常好。

希望你喜欢！ 如果您发现此博客文章有用，您可能想在Twitter上关注我的博客帖子更新并给我买浓咖啡。

原文链接：
https://www.brodrigues.co/blog/2018-11-25-tidy_cv/

内容推荐

• 如何阅读论文？

• 论文管理工具，我用Zotero

• R机器学习介绍第一部分

• R机器学习介绍第二部分

• R中你应该学习的7种可视化

• R中用线性回归进行预测建模

• 使用RMarkdown沟通结果

• 使用LIME探索模型

• RStudio1.2新功能介绍
  

数据人网：数据人学习，交流和分享的平台，诚邀您创造和分享数据知识，共建和共享数据智库。