无法运行ANOVA来比较随机森林模型。

Question

我正在使用tidymodels来拟合多个随机森林模型。然后我按照这个教程进行比较模型结果。但是我遇到了错误：Error in

 UseMethod("anova") : 
  no applicable method for 'anova' applied to an object of class "ranger"

作为一个例子：

set.seed(123)
iris <- iris %>% mutate(
  is_versicolor = ifelse(Species == "versicolor", "versicolor", "not_versicolor")) %>%
  mutate(is_versicolor = factor(is_versicolor, levels = c("versicolor", "not_versicolor")))

iris_split <- initial_split(iris, strata = is_versicolor, prop = 0.8)
iris_train <- training(iris_split)
iris_test  <- testing(iris_split)

rec_normal <- recipe(is_versicolor ~ Petal.Width + Species, data = iris_train)
rec_interaction <- rec_normal %>% 
  step_interact(~ Petal.Width:starts_with("Species"))

iris_model <- rand_forest() %>% set_engine("ranger") %>% set_mode("classification")

# normal workflow
iris_wf <- workflow() %>% 
  add_model(iris_model) %>% 
  add_recipe(rec_normal)

# interaction workflow
iris_wf_interaction <- iris_wf %>% 
  update_recipe(rec_interaction)

# fit models
iris_normal_lf <- last_fit(iris_wf, split = iris_split)
iris_inter_lf <- last_fit(iris_wf_interaction, split = iris_split)

normalmodel <- iris_normal_lf %>% extract_fit_engine()
intermodel  <- iris_inter_lf %>% extract_fit_engine()

anova(normalmodel, intermodel) %>% tidy()

我该如何运行ANOVA或类似的模型比较，以查看哪个模型更好？

- Adam_G

Answer 1

只需使用您的代码，并适应Julia Silge关于workflowsets的博客：

使用workflowsets预测#TidyTuesday巨型南瓜重量

由于ranger不支持ANOVA，因此请生成折叠以进行重新采样：

set. Seed(234)
iris_folds <- vfold_cv(iris_train)
iris_folds

将您的配方组合成一个工作流集：

iris_set <-
  workflow_set(
    list(rec_normal, rec_interaction),
    list(iris_model),
    cross = TRUE
  )

iris_set

设置并行处理：

doParallel::registerDoParallel()
set. Seed(2021)

使用交叉验证进行拟合：

iris_rs <-
  workflow_map(
    iris_set,
    "fit_resamples",
    resamples = iris_folds
  )

autoplot(iris_rs)

这个图表通常会回答你如何比较模型的问题。

由于"species"在两个配方公式的右侧，并且响应变量"is_versicolor"是由species计算得出的，所以这些模型是完全准确的。

完成输出：

collect_metrics(iris_rs)

final_fit <-
  extract_workflow(iris_rs, "recipe_2_rand_forest") %>%
  fit(iris_train)

游侠模型中没有整洁的人。

如果您在代码中进行更改：

rec_normal <- recipe(is_versicolor ~ Sepal.Length + Sepal.Width, data = iris_train)
rec_interaction <- recipe(is_versicolor ~ Petal.Width + Petal.Length, data = iris_train)

你可以玩得开心！

希望这能对Adam有所帮助。像你一样，我也在学习美妙的Tidymodels，并期待着您的评论。 :-)

- Isaiah

---

谢谢！这很有帮助。那么我如何将它们转换为类似于 P 值的东西？ - Adam_G

使用类似性能指标这样的东西并不会给你一个p值；相反，它会给你一个度量标准上的平均值和方差，以便比较两个模型配置。如果你想要更像ANOVA的东西，可以尝试在重采样上使用lme4::lmer()。查看tune_race_anova()的内部以了解可能的方法。 - Julia Silge

限制性的，但你可以将你的引擎设置为GLM：iris_model <- logistic_reg() |> set_engine("glm") |> set_mode("classification") - Isaiah

谢谢@JuliaSilge！（顺便说一句，我很喜欢你在NormConf上的演讲；正试图在这里获得灵感 :-)） - Adam_G

Answer 2

您可以使用 aov 函数比较随机森林模型的准确性。首先，您可以使用 collect_metrics 收集准确性，并将其保存在数据框中以运行带有 aov 的模型来获取结果。这是一个可再现的示例：

library(tidymodels)
set.seed(123)
iris <- iris %>% mutate(
  is_versicolor = ifelse(Species == "versicolor", "versicolor", "not_versicolor")) %>%
  mutate(is_versicolor = factor(is_versicolor, levels = c("versicolor", "not_versicolor")))

iris_split <- initial_split(iris, strata = is_versicolor, prop = 0.8)
iris_train <- training(iris_split)
iris_test  <- testing(iris_split)

rec_normal <- recipe(is_versicolor ~ Petal.Width + Species, data = iris_train)
rec_interaction <- rec_normal %>% 
  step_interact(~ Petal.Width:starts_with("Species"))

iris_model <- rand_forest() %>% set_engine("ranger") %>% set_mode("classification")

# normal workflow
iris_wf <- workflow() %>% 
  add_model(iris_model) %>% 
  add_recipe(rec_normal)

# interaction workflow
iris_wf_interaction <- iris_wf %>% 
  update_recipe(rec_interaction)

# fit models
iris_normal_lf <- last_fit(iris_wf, split = iris_split)
iris_inter_lf <- last_fit(iris_wf_interaction, split = iris_split)
#> ! train/test split: preprocessor 1/1: Categorical variables used in `step_interact` should probably be avoided...

normalmodel <- iris_normal_lf %>% extract_fit_engine()
intermodel  <- iris_inter_lf %>% extract_fit_engine()

# Check confusion matrix
iris_normal_lf %>%
  collect_predictions() %>% 
  conf_mat(is_versicolor, .pred_class) 
#>                 Truth
#> Prediction       versicolor not_versicolor
#>   versicolor             10              0
#>   not_versicolor          0             20

iris_inter_lf %>%
  collect_predictions() %>% 
  conf_mat(is_versicolor, .pred_class) 
#>                 Truth
#> Prediction       versicolor not_versicolor
#>   versicolor             10              0
#>   not_versicolor          0             20

# Extract accuracy of models and create dataset
acc_normalmodel <- iris_normal_lf %>% collect_metrics() %>% select(.estimate) %>% slice(1)
acc_intermodel <- iris_normal_lf %>% collect_metrics() %>% select(.estimate) %>% slice(1)
results = data.frame(model = c("normalmodel", "intermodel"),
                     accuracy = c(acc_normalmodel$.estimate, acc_intermodel$.estimate))

# perform ANOVA on the classification accuracy
aov_results <- aov(accuracy ~ model, data = results)
summary(aov_results)
#>             Df   Sum Sq  Mean Sq
#> model        1 4.93e-32 4.93e-32

^{使用reprex v2.0.2于2022年12月15日创建}

正如你所看到的，结果并未显示出P值，因为自由度太低 (为什么在R中不能从此ANOVA获取P值)

---

你还可以在两个模型的预测结果上使用 aov 并比较它们的性能。以下是一个可复现的例子：

# Get predictions of both models for not_versicolor
normalmodel_pred<-as.data.frame(normalmodel$predictions)$not_versicolor
intermodel_pred<-as.data.frame(intermodel$predictions)$not_versicolor

summary(aov(normalmodel_pred~intermodel_pred))
#>                  Df Sum Sq Mean Sq F value Pr(>F)    
#> intermodel_pred   1 25.032  25.032    9392 <2e-16 ***
#> Residuals       118  0.314   0.003                   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{由 reprex v2.0.2 创建于2022年12月17日}

正如您所看到的，p值小于0.05，这表明模型之间存在差异，如果您查看预测的概率，这是正确的。

---

关于ANOVA的更多信息，请查看以下内容：

• 在R中理解ANOVA第7章

- Quinten

---

你能展示一下怎么做吗？不幸的是，如果没有这个，这并不能真正回答问题。 - Adam_G

@Adam_G，我添加了一些代码来对两个模型的预测执行aov，这将返回一个显著性。 - Quinten

该模型是一个分类模型，因此原始预测将被转换为not_versicolor或versicolor。两个模型的分类准确性相同，均为100％。您可以使用 round (normalmodel [["predictions"]] [,1]) - intermodel [["predictions"]] [,1] 来查看这一点，结果会全部显示为零。作为分类模型，这两个模型之间没有区别。 - Isaiah

嗨@Isaiah，那是真的，但每个模型的分类概率可能不同，对吧？ - Quinten

rec_normal <- recipe(is_versicolor ~ Sepal.Width, data = iris_train) rec_alternative <- recipe(is_versicolor ~ Sepal.Length, data = iris_train) 这是一对更实用的模型，因为分类准确性并不是100%。 - Isaiah

Answer 3

使用不同的模型对进行配对，并基于重复样本的分类准确性比较模型。可以轻松扩展到其他指标。

library(dplyr)
library(tibble)
library(ggplot2)
library(tidyr)
library(rsample)
library(recipes)
library(parsnip)
library(workflows)
library(tune)
library(yardstick)
library(workflowsets)

set.seed(123)
iris <- iris %>% mutate(
  is_versicolor = ifelse(Species == "versicolor", "versicolor", "not_versicolor")) %>%
  mutate(is_versicolor = factor(is_versicolor, levels = c("versicolor", "not_versicolor")))

iris_split <- initial_split(iris, strata = is_versicolor, prop = 0.8)
iris_train <- training(iris_split)
iris_test  <- testing(iris_split)

# replacing normal and interaction recipes with models
# that give less than 100% accuracy.
rec_normal <- recipe(is_versicolor ~ Sepal.Width, data = iris_train)
rec_alternative <- recipe(is_versicolor ~ Sepal.Length, data = iris_train)

iris_model <- rand_forest() %>% set_engine("ranger") %>% set_mode("classification")

# Create folds
set.seed(234)
iris_folds <- vfold_cv(iris_train)
iris_folds

# Combine models into set
iris_set <-
  workflow_set(
    list(rec_normal, rec_alternative),
    list(iris_model),
    cross = TRUE
  )

doParallel::registerDoParallel()
set.seed(2021)

# fit models
iris_rs <-
  workflow_map(
    iris_set,
    "fit_resamples",
    resamples = iris_folds
  )

# Visualise model performance
autoplot(iris_rs)

# Extract resample accuracies
model_1_rs <- iris_rs[1,][[4]][[1]]$.metrics
model_2_rs <- iris_rs[2,][[4]][[1]]$.metrics
model_acc <- tibble(model_1 = NA, model_2 = NA)
for (i in 1:10) {
  model_acc[i, 1] <- model_1_rs[[i]][[".estimate"]][1]
  model_acc[i, 2] <- model_2_rs[[i]][[".estimate"]][1]
}
model_acc <- model_acc |> pivot_longer(cols = everything(), names_to = "model", values_to = "acc")

# Do ANOVA
aov_results <- aov(acc ~ model, data = model_acc)
summary(aov_results)
ggplot(data = model_acc, aes(fill = model)) +
  geom_density(aes(x = acc, alpha = 0.2)) +
  labs(x = "accuracy")

给出 p 值：

> summary(aov_results)

            Df Sum Sq Mean Sq F value Pr(>F)

model        1 0.0281 0.02813   1.378  0.256

Residuals   18 0.3674 0.02041

从不同的角度观察模型准确性的p值： 首先可视化变化：

model_acc |> ggplot(aes(x = model, y = acc)) +
  geom_boxplot() +
  labs(y = 'accuracy')

然后计算一个测试统计量：

observed_statistic <- model_acc %>%
  specify(acc ~ model) %>%
  calculate(stat = "diff in means", order = c("model_1", "model_2"))

observed_statistic

然后进行分布的模拟：

null_dist_2_sample <- model_acc %>%
  specify(acc ~ model) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means" ,order = c("model_1", "model_2"))

并绘制：

null_dist_2_sample %>%
  visualize() + 
  shade_p_value(observed_statistic,
                direction = "two-sided") +
  labs(x = "test statistic")

并获取p值：

p_value_2_sample <- null_dist_2_sample %>%
  get_p_value(obs_stat = observed_statistic,
              direction = "two-sided")

p_value_2_sample

# A tibble: 1 × 1
  p_value
    <dbl>
1   0.228

这与aov的p值几乎相同。 请注意，由于两个模型的准确性接近，p值较高。

- Isaiah

---

添加了第二个p值的计算。 - Isaiah