阅读 26 SEO

R语言数据分析,从入门到自闭(图缓慢补充中……)

写在前面

本文章参考了多个国内外教程,在这里就不列出来了。这么写的原因还是想告诉大家:本文章并非原创,感谢各位前辈的分享。

加载需要的数据包

本文章用到的数据包中,“tidyverse”,“skimr”,“FactoMineR”,“factoextra”,“pheatmap”是必须的。但是为了更好的查看数据,所以又加入了“GGally”,“patchwork”,"ggstatsplot和"“ggpubr”。后几个包并不是必须的,但会让你的数据可视化更方便。各种包按自己的好恶来,例如有人就极度讨厌“ggpubr”。最后“pheatmap”虽然已经被作者放弃了,但是我觉得他搞的“ComplexHeatmap”实在是有点走火入魔……可能他也意识到这个问题,现在给了个CompleHeatmap::pheatmap的选项(还能说什么,绝了)。什么你问knitr干嘛的?你猜……

knitr::opts_chunk$set(echo = TRUE, warning = FALSE)

library('tidyverse')
library('ggpubr')
library('patchwork')
library('skimr')
library('GGally')
library('lme4')
library("FactoMineR")
library("factoextra")
library('ComplexHeatmap')
library('ggstatsplot')
library('agricolae')
library('car')
library('vip')
library('onewaytests')
library('jmv')

读取数据

这里rm命令清空存在环境中的所有变量,避免先前环境中的变量对接下来的操作带来影响。


rm(list = ls())

为了方便大家测试,这里使用了R语言自带数据集iris,如果你想用mtcars或者别的请随意。比如想看汽车各种特性对油耗的影响就可以用mtcars。skim可以很华丽的展示你的数据结构。当然,没啥别的用途了。


head(iris)
data <- iris
skim(data %>% group_by(Species))

skim_variable Species n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Sepal.Length setosa 0 1 5.01 0.35 4.3 4.80 5.00 5.20 5.8 ▃▃▇▅▁
Sepal.Length versicolor 0 1 5.94 0.52 4.9 5.60 5.90 6.30 7.0 ▂▇▆▃▃
Sepal.Length virginica 0 1 6.59 0.64 4.9 6.23 6.50 6.90 7.9 ▁▃▇▃▂
Sepal.Width setosa 0 1 3.43 0.38 2.3 3.20 3.40 3.68 4.4 ▁▃▇▅▂
Sepal.Width versicolor 0 1 2.77 0.31 2.0 2.52 2.80 3.00 3.4 ▁▅▆▇▂
Sepal.Width virginica 0 1 2.97 0.32 2.2 2.80 3.00 3.18 3.8 ▂▆▇▅▁
Petal.Length setosa 0 1 1.46 0.17 1.0 1.40 1.50 1.58 1.9 ▁▃▇▃▁
Petal.Length versicolor 0 1 4.26 0.47 3.0 4.00 4.35 4.60 5.1 ▂▂▇▇▆
Petal.Length virginica 0 1 5.55 0.55 4.5 5.10 5.55 5.88 6.9 ▃▇▇▃▂
Petal.Width setosa 0 1 0.25 0.11 0.1 0.20 0.20 0.30 0.6 ▇▂▂▁▁
Petal.Width versicolor 0 1 1.33 0.20 1.0 1.20 1.30 1.50 1.8 ▅▇▃▆▁
Petal.Width virginica 0 1 2.03 0.27 1.4 1.80 2.00 2.30 2.5 ▂▇▆▅▇

随时随地的可视化

这是我认为R语言跟spss和其他软件比最大的优势了。是的,在Rstudio中你可以随时随地的可视化,无限制的切片数据可视化。相比于单纯的统计分析,我认为视觉往往来的更准确更直接。话不多说让我们开始吧。

首先,我们在R语言里进行一些传统的方法。我们进行一个切片求平均值。这里我们使用了tidyverse套件(dplyr)。其中的 %>% 是通道符号,他的含义是将前面的参数传入到下一个命令中作为第一个参数(可以用.代替)。统计各个组的平均值和标准偏差,过去大家都在用spss得到这些,统计分析,最后获得结果,大家都满意了,在R里你同样可以做到。


data_summary <- data %>% group_by(Species) %>% summarise_each(funs(mean,sd),
                                                           Sepal.Length, Sepal.Width, 
                                                           Petal.Length, Petal.Width)
data_summary

## # A tibble: 3 x 9
##   Species Sepal.Length_me~ Sepal.Width_mean Petal.Length_me~ Petal.Width_mean
##   <fct>              <dbl>            <dbl>            <dbl>            <dbl>
## 1 setosa              5.01             3.43             1.46            0.246
## 2 versic~             5.94             2.77             4.26            1.33 
## 3 virgin~             6.59             2.97             5.55            2.03 
## # ... with 4 more variables: Sepal.Length_sd <dbl>, Sepal.Width_sd <dbl>,
## #   Petal.Length_sd <dbl>, Petal.Width_sd <dbl>

当然你也可以这样……


data_res <- data %>%  pivot_longer(col = -Species, names_to = 'Name', values_to = 'Value') %>% 
  group_by(Species,Name) %>% 
  summarise(mean = mean(Value), sd = mean(Value))

data_res

## # A tibble: 12 x 4
## # Groups:   Species [3]
##    Species    Name          mean    sd
##    <fct>      <chr>        <dbl> <dbl>
##  1 setosa     Petal.Length 1.46  1.46 
##  2 setosa     Petal.Width  0.246 0.246
##  3 setosa     Sepal.Length 5.01  5.01 
##  4 setosa     Sepal.Width  3.43  3.43 
##  5 versicolor Petal.Length 4.26  4.26 
##  6 versicolor Petal.Width  1.33  1.33 
##  7 versicolor Sepal.Length 5.94  5.94 
##  8 versicolor Sepal.Width  2.77  2.77 
##  9 virginica  Petal.Length 5.55  5.55 
## 10 virginica  Petal.Width  2.03  2.03 
## 11 virginica  Sepal.Length 6.59  6.59 
## 12 virginica  Sepal.Width  2.97  2.97

有了平均值和标准差,我们自然能画出第一个柱状图。


cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- ggplot(data_summary,aes(x = Species, y = Sepal.Length_mean, fill = Species)) + 
  geom_bar(stat =  'identity', position = 'dodge', width = 0.5)  + 
  geom_errorbar(aes(ymin = Sepal.Length_mean - Sepal.Length_sd, 
                    ymax = Sepal.Length_mean + Sepal.Length_sd), width = 0.25) + 
  scale_fill_manual(values = cbPalette) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme_pubr()
p1

重复四次,你能获得四个柱状图


compare_group <- list(c('setosa','versicolor'),c('versicolor','virginica'),c('setosa','virginica'))

p1 <- ggplot(data_summary,aes(x = Species, y = Sepal.Length_mean, fill = Species)) + 
  geom_bar(stat =  'identity', position = 'dodge', width = 0.5)  + 
  geom_errorbar(aes(ymin = Sepal.Length_mean - Sepal.Length_sd, 
                    ymax = Sepal.Length_mean + Sepal.Length_sd), width = 0.25) + 
  scale_fill_manual(values = cbPalette) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme_pubr()

p2 <- ggplot(data_summary,aes(x = Species, y = Sepal.Width_mean, fill = Species)) + 
  geom_bar(stat =  'identity', position = 'dodge', width = 0.5)  + 
  geom_errorbar(aes(ymin = Sepal.Width_mean - Sepal.Width_sd, 
                    ymax = Sepal.Width_mean + Sepal.Width_sd), width = 0.25) + 
  scale_fill_manual(values = cbPalette) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme_pubr()

p3 <- ggplot(data_summary,aes(x = Species, y = Petal.Length_mean, fill = Species)) + 
  geom_bar(stat =  'identity', position = 'dodge', width = 0.5)  + 
  geom_errorbar(aes(ymin = Petal.Length_mean - Petal.Length_sd, 
                    ymax = Petal.Length_mean + Petal.Length_sd), width = 0.25) + 
  scale_fill_manual(values = cbPalette) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme_pubr()

p4 <- ggplot(data_summary,aes(x = Species, y = Petal.Width_mean, fill = Species)) + 
  geom_bar(stat =  'identity', position = 'dodge', width = 0.5)  + 
  geom_errorbar(aes(ymin = Petal.Width_mean - Petal.Width_sd, 
                    ymax = Petal.Width_mean + Petal.Width_sd), width = 0.25) + 
  scale_fill_manual(values = cbPalette) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme_pubr()

p1+p2+p3+p4+plot_layout(guides = 'collect')&theme(legend.position = 'bottom')

然后你可以做统计分析。别忘了检查各个品种的变量是否满足正态分布,会发现正态性检验还好。除了setosa的petal.width不是很满足。


data %>% group_by(Species) %>% 
  summarise(statistic = shapiro.test(Sepal.Length)$statistic,
            p.value = shapiro.test(Sepal.Length)$p.value)

data %>% group_by(Species) %>% 
  summarise(statistic = shapiro.test(Sepal.Width)$statistic,
            p.value = shapiro.test(Sepal.Width)$p.value)

data %>% group_by(Species) %>% 
  summarise(statistic = shapiro.test(Petal.Length)$statistic,
            p.value = shapiro.test(Petal.Length)$p.value)

data %>% group_by(Species) %>% 
  summarise(statistic = shapiro.test(Petal.Width)$statistic,
            p.value = shapiro.test(Petal.Width)$p.value)

## # A tibble: 3 x 3
##   Species    statistic p.value
##   <fct>          <dbl>   <dbl>
## 1 setosa         0.978   0.460
## 2 versicolor     0.978   0.465
## 3 virginica      0.971   0.258

## # A tibble: 3 x 3
##   Species    statistic p.value
##   <fct>          <dbl>   <dbl>
## 1 setosa         0.972   0.272
## 2 versicolor     0.974   0.338
## 3 virginica      0.967   0.181

## # A tibble: 3 x 3
##   Species    statistic p.value
##   <fct>          <dbl>   <dbl>
## 1 setosa         0.955  0.0548
## 2 versicolor     0.966  0.158 
## 3 virginica      0.962  0.110

## # A tibble: 3 x 3
##   Species    statistic     p.value
##   <fct>          <dbl>       <dbl>
## 1 setosa         0.800 0.000000866
## 2 versicolor     0.948 0.0273     
## 3 virginica      0.960 0.0870

方差齐性检验。发现各组petal.length和petal.width方差不齐。


leveneTest(Sepal.Length ~ Species, data=data)

leveneTest(Sepal.Width ~ Species, data=data)

leveneTest(Petal.Length ~ Species, data=data)

leveneTest(Petal.Width ~ Species, data=data)

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value   Pr(>F)   
## group   2  6.3527 0.002259 **
##       147                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.5902 0.5555
##       147

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group   2   19.48 3.129e-08 ***
##       147                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group   2  19.892 2.261e-08 ***
##       147                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

可以换成箱线图,再加上散点,标记好显著性。因为之前发现有点不满足正态和方差齐性,所以这里选择了非参数检验。


data %>% pivot_longer(col = -Species, names_to = 'Name', values_to = 'Value') %>% 
  ggplot(aes(x = Species, y = Value, color = Species)) + 
  geom_boxplot( ) + 
  geom_jitter(width = 0.2, height = 0.5, size = .1) + 
  scale_color_manual(values = cbPalette) + 
  scale_y_continuous(expand = c(0, 0),
                     limits = c(0, 12)) +
  facet_wrap(~Name) + 
  stat_compare_means(comparisons = compare_group, label = 'p.signif',method = 'wilcox.test') + theme_pubr()

如果你觉得wilcox.test还不够稳健,你可以试试Welch's ANOVA检验和Games-Howell事后检验


anovaOneW(deps=Sepal.Length,group=Species,data=data,desc=T,
          descPlot = T,norm=T,qq=T,eqv=T,phMeanDif = T,
          phMethod= "gamesHowell",phTest = T,phFlag=T)

anovaOneW(deps=Sepal.Width,group=Species,data=data,desc=T,
          descPlot = T,norm=T,qq=T,eqv=T,phMeanDif = T,
          phMethod= "gamesHowell",phTest = T,phFlag=T)

anovaOneW(deps=Petal.Length,group=Species,data=data,desc=T,
          descPlot = T,norm=T,qq=T,eqv=T,phMeanDif = T,
          phMethod= "gamesHowell",phTest = T,phFlag=T)

anovaOneW(deps=Petal.Width,group=Species,data=data,desc=T,
          descPlot = T,norm=T,qq=T,eqv=T,phMeanDif = T,
          phMethod= "gamesHowell",phTest = T,phFlag=T)

篇幅原因这里只放一部分结果

## 
##  ONE-WAY ANOVA
## 
##  One-Way ANOVA (Welch's)                                       
##  ------------------------------------------------------------- 
##                    F           df1    df2         p            
##  ------------------------------------------------------------- 
##    Sepal.Length    138.9083      2    92.21115    < .0000001   
##  ------------------------------------------------------------- 
## 
## 
##  Group Descriptives                                                          
##  --------------------------------------------------------------------------- 
##                    Species       N     Mean        SD           SE           
##  --------------------------------------------------------------------------- 
##    Sepal.Length    setosa        50    5.006000    0.3524897    0.04984957   
##                    versicolor    50    5.936000    0.5161711    0.07299762   
##                    virginica     50    6.588000    0.6358796    0.08992695   
##  --------------------------------------------------------------------------- 
## 
## 
##  ASSUMPTION CHECKS
## 
##  Normality Test (Shapiro-Wilk)              
##  ------------------------------------------ 
##                    W            p           
##  ------------------------------------------ 
##    Sepal.Length    0.9878974    0.2188639   
##  ------------------------------------------ 
##    Note. A low p-value suggests a
##    violation of the assumption of
##    normality
## 
## 
##  Homogeneity of Variances Test (Levene's)                
##  ------------------------------------------------------- 
##                    F           df1    df2    p           
##  ------------------------------------------------------- 
##    Sepal.Length    7.381092      2    147    0.0008818   
##  ------------------------------------------------------- 
## 
## 
##  POST HOC TESTS
## 
##  Games-Howell Post-Hoc Test – Sepal.Length                                  
##  -------------------------------------------------------------------------- 
##                                     setosa       versicolor    virginica    
##  -------------------------------------------------------------------------- 
##    setosa        Mean difference            —    -0.9300000    -1.5820000   
##                  t-value                    —     -10.52099    -15.386196   
##                  df                         —      86.53800      76.51587   
##                  p-value                    —    < .0000001    < .0000001   
##                                                                             
##    versicolor    Mean difference                          —    -0.6520000   
##                  t-value                                  —     -5.629165   
##                  df                                       —      94.02549   
##                  p-value                                  —     0.0000006   
##                                                                             
##    virginica     Mean difference                                        —   
##                  t-value                                                —   
##                  df                                                     —   
##                  p-value                                                —   
##  -------------------------------------------------------------------------- 
##    Note. * p < .05, ** p < .01, *** p < .001

利用R语言对数据进行可视化

到了这里,事情解决了,对么?很传统的方法,得到了想知道的一切。数据有差异性,完美。但其实,R语言能够做的更多。我更愿意称R语言为可视化统计检验方法。相比于传统方法的先检验差异性再作图,R语言是先做可视化再做具体的方差分析。例如,我们可以一行代码查看数据的关系。这里面包含了数据的相关散点图,直方图,相关系数和箱线图。依靠这些你可以更好的提出你需要的科学假设。


ggpairs(data, mapping = aes(color = Species)) + theme_bw()

通过直方图可以很容易看到数据的分布情况。是否满足正态还有方差齐性是不是变得更加具象了?


p_SL <- data %>% gghistogram(data = ., x= 'Sepal.Length',  add = 'mean', fill = 'Species', palette = cbPalette,add_density = TRUE, legend = 'right') 

p_SW <- data %>% gghistogram(data = ., x= 'Sepal.Width',  add = 'mean', fill = 'Species', palette = cbPalette,add_density = TRUE, legend = 'right') 

p_PL <- data %>% gghistogram(data = ., x= 'Petal.Length',  add = 'mean', fill = 'Species', palette = cbPalette,add_density = TRUE, legend = 'right') 

p_PW <- data %>% gghistogram(data = ., x= 'Petal.Width',  add = 'mean', fill = 'Species', palette = cbPalette,add_density = TRUE, legend = 'none' ) 

(p_SL + p_SW +  p_PL +  guide_area() + plot_layout(guides = 'collect')) /p_PW 

跑个qq图,配合置信区间。



qqPlot(data$Sepal.Length ,main="qq plot", groups = c("setosa","versicolor","virginica"), col="blue", col.lines="red")

qqPlot(data$Sepal.Width ,main="qq plot", groups = c("setosa","versicolor","virginica"), col="blue", col.lines="red")

qqPlot(data$Petal.Length ,main="qq plot", groups = c("setosa","versicolor","virginica"), col="blue", col.lines="red")

qqPlot(data$Petal.Width ,main="qq plot", groups = c("setosa","versicolor","virginica"), col="blue", col.lines="red")


然后是数据结果的分组可视化。


compare_group <- list(c('setosa','versicolor'),c('versicolor','virginica'),c('setosa','virginica'))

p1 <- ggboxplot(data, x = 'Species', y = 'Sepal.Length', add = c('jitter','mean_se'), color = 'Species', palette = cbPalette) + 
  stat_compare_means(comparisons = compare_group, label = 'p.signif',method = 'wilcox.test')

p2 <- ggboxplot(data, x = 'Species', y = 'Sepal.Width', add = c('jitter','mean_se'), color = 'Species', palette = cbPalette) + 
  stat_compare_means(comparisons = compare_group, label = 'p.signif',method = 'wilcox.test')

p3 <- ggboxplot(data, x = 'Species', y = 'Petal.Length', add = c('jitter','mean_se'), color = 'Species', palette = cbPalette) + 
  stat_compare_means(comparisons = compare_group, label = 'p.signif',method = 'wilcox.test')

p4 <- ggboxplot(data, x = 'Species', y = 'Petal.Width', add = c('jitter','mean_se'), color = 'Species', palette = cbPalette) + 
  stat_compare_means(comparisons = compare_group, label = 'p.signif',method = 'wilcox.test')

p1 + p2 + p3 + p4 + plot_layout(guides = 'collect') & theme(legend.position = 'bottom')


还可以换个炫酷的效果。这里用的Games_Howell test。不喜欢可以换回t检验,大家可以自己试试看。


p1 <- data %>% group_by(Species) %>% ggstatsplot::ggbetweenstats(x = Species, y = Sepal.Length, nboot = 10, messages = FALSE)

p2 <- data %>% group_by(Species) %>% ggstatsplot::ggbetweenstats(x = Species, y = Sepal.Width, nboot = 10, messages = FALSE)

p3 <- data %>% group_by(Species) %>% ggstatsplot::ggbetweenstats(x = Species, y = Petal.Length, nboot = 10, messages = FALSE)

p4 <- data %>% group_by(Species) %>% ggstatsplot::ggbetweenstats(x = Species, y = Petal.Width, nboot = 10, messages = FALSE)

p1 + p2 + p3 + p4

最后我们可以用heatmap查看各个数据在四个维度上的表现,这里我因为同时装了两个包所以,冲突了。所以需要申明到底用的是哪个包的pheatmap。


set.seed(2020)
anno <- data.frame(Species = data$Species)
row.names(anno) <- row.names(data)

ComplexHeatmap::pheatmap(data[,-5], border_color = gpar(col = "black", lty = 2), cluster_rows = T, cluster_cols = T, show_rownames = F, show_colnames = T, annotation_row = anno)

最后根据四个petal和sepal的长宽数据,我们做个pca分类图。


data.pca <- PCA(data[,-5], graph = FALSE)

fviz_pca_ind(data.pca,
             geom.ind = "point",
             col.ind = data$Species,
             addEllipses = TRUE,
             ellipse.type = 'convex',
             legend.title = "Groups"
) + theme_minimal()
row.names(data) <- paste0("row_", seq(nrow(data)))

res <- factoextra::hcut(data[,-5], k = 3, stand = TRUE)

fviz_dend(res, rect = TRUE, cex = 0.5,
          k_colors = c("#00AFBB","#2E9FDF", "#E7B800", "#FC4E07"))

Tidymodel建模

到这里,我们体验了大部分统计可视化的内容。那么R还能做什么?既然我们知道所谓统计分析本身就是建模,那么我们能否用更复杂的模型对petal和sepal的长宽等因子进行重要性检验呢?所以又有了这部分,利用tidymodels的建模,并根据模型对数据分类。


library('tidymodels')
library('caret')
#固定随机数,方便复现
set.seed(2020)

下面是整个建模的流程


#将数据分成3份两份作为训练数据,一份作为测试准确性数据。
data_split <- initial_split(data,prop=.66)
data_train <- training(data_split)
data_test <-testing(data_split) 

#bootstrap创造一个数据用来tuning模型
Spec_boost <- bootstraps(data_train, times = 30)

#设定模型数据转换,在这里可以做中心化,标准化,数据合并等一系列操作。最后传入的prep记录所有的操作。另外在这里Species~.意思是Species作为分类结果,其他几列作为predictor预测分类。
data_rec <-
  recipe(Species ~.,data = data_train) %>%
  prep()

#查看训练的数据,这里用juice很形象的说你的数据是榨出来的。但这里没对表型数据做过多调整。
juice(data_rec)

#设定模型,这里mtry和min_n使用tune函数调整。模式选择分类,随机森林的引擎选择randomForest,另外还有一个引擎叫ranger。具体到后面会有细微差别。
rf_model<-rand_forest(mtry=tune(),min_n = tune())%>%
           set_mode("classification")%>%
           set_engine("randomForest")

#设定工作流,主要是使用哪个模型以及使用哪个数据方便后期做具体的调参。
rf_wflow <-
    workflow() %>%
    add_recipe(data_rec)%>%
    add_model(rf_model) 



#激活平行线程
doParallel::registerDoParallel()

#使用tune_grid函数和boostrap的数据进行大致的调参。
rf_results <-
  rf_wflow %>% 
  tune_grid(resamples = Spec_boost)

#参数打分可视化
rf_results %>%
  collect_metrics() %>%
  filter(.metric == "roc_auc") %>%
  select(mean, min_n, mtry) %>%
  pivot_longer(min_n:mtry,
    values_to = "value",
    names_to = "parameter"
  ) %>%
  ggplot(aes(value, mean, color = parameter)) +
  geom_point(show.legend = FALSE) +
  facet_wrap(~parameter, scales = "free_x") +
  labs(x = NULL, y = "AUC")
#选择最佳模型参数
so_best <-
  rf_results %>% 
    select_best(metric = "roc_auc")

#构建模型
rf_fit <- rand_forest(mode = 'classification',
                       mtry = so_best$mtry,
                       min_n = so_best$min_n) %>% 
  set_engine("randomForest") %>% 
  fit( Species ~ ., data = data_train)

plot(rf_fit$fit)
#各因子重要性可视化
Importance <- varImp(rf_fit$fit)
Importance$Feature <- row.names(Importance)
Importance

##                   Overall      Feature
## Sepal.Length  0.007138445 Sepal.Length
## Sepal.Width   0.000000000  Sepal.Width
## Petal.Length 28.186545520 Petal.Length
## Petal.Width  31.611573933  Petal.Width
ggdotchart(Importance,x = 'Feature', y = 'Overall', color = 'Feature', 
           palette = 'mpg', sorting = "descending",
           font.label = list(color = "white", size = 9, vjust = 0.5),add.params = list(color = "Feature"), add = "segments", rotate = TRUE, group = "Feature", dot.size = 6,
          ggtheme = theme_pubr())

#激活平行线程
doParallel::registerDoParallel()

#设定参数
rf_grid <- grid_regular(
  mtry(range = c(1,4)),
  min_n(range = c(10,30)),
  levels = 10
)

#调参结果
regular_res <- rf_wflow %>% 
  tune_grid(
  resamples = Spec_boost,
  grid = rf_grid
)

#结果可视化
regular_res%>%
  collect_metrics() %>%
  filter(.metric == "roc_auc") %>%
  mutate(mtry = factor(mtry)) %>%
  ggplot(aes(min_n, mean, color = mtry)) +
  geom_line(alpha = 0.5, size = 1.5) +
  geom_point() +
  labs(y = "AUC")
#调参后通过roc_auc指标选择最好的模型参数
so_best <-
  regular_res %>% 
    select_best(metric = "roc_auc")

#使用最佳参数去构建模型
rf_final_fit<-rf_wflow%>%
              finalize_workflow(so_best)%>%
              fit(data = data_train)

#设定好数据调整的方法(还记得recipes么)
data_rec3 <-rf_final_fit%>% 
pull_workflow_prepped_recipe()  

#用同样方法去变换test数据集,然后用模型拟合test数据集,输出结果。
rf_final_fit%>%
        pull_workflow_fit()%>%
        predict( new_data = bake(data_rec3, data_test))%>% 
        bind_cols(data_test, .)%>%
        metrics(truth = Species, estimate =.pred_class)
## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy multiclass     0.96 
## 2 kap      multiclass     0.940
#查看参数重要性
vi(rf_final_fit$fit$fit)
vip(rf_final_fit$fit$fit, mapping = aes(fill = row.names(rf_final_fit$fit$fit$fit$importance))) + theme_pubr()

作者:邵扬_Barnett

原文链接:https://www.jianshu.com/p/43fa029048ba

文章分类
后端
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容, 请发送邮件至 gxwowoo@163.com 举报,一经查实,本站将立刻删除。
相关推荐