寫(xiě)在前面
本文章參考了多個(gè)國(guó)內(nèi)外教程，在這里就不列出來(lái)了。這么寫(xiě)的原因還是想告訴大家：本文章并非原創(chuàng)，感謝各位前輩的分享。
加載需要的數(shù)據(jù)包
本文章用到的數(shù)據(jù)包中，“tidyverse”，“skimr”，“FactoMineR”，“factoextra”，“pheatmap”是必須的。但是為了更好的查看數(shù)據(jù)，所以又加入了“GGally”，“patchwork”,"ggstatsplot和"“ggpubr”。后幾個(gè)包并不是必須的，但會(huì)讓你的數(shù)據(jù)可視化更方便。各種包按自己的好惡來(lái)，例如有人就極度討厭“ggpubr”。最后“pheatmap”雖然已經(jīng)被作者放棄了，但是我覺(jué)得他搞的“ComplexHeatmap”實(shí)在是有點(diǎn)走火入魔……可能他也意識(shí)到這個(gè)問(wèn)題，現(xiàn)在給了個(gè)CompleHeatmap::pheatmap的選項(xiàng)（還能說(shuō)什么，絕了）。什么你問(wèn)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')
讀取數(shù)據(jù)
這里rm命令清空存在環(huán)境中的所有變量，避免先前環(huán)境中的變量對(duì)接下來(lái)的操作帶來(lái)影響。
rm(list = ls())
為了方便大家測(cè)試，這里使用了R語(yǔ)言自帶數(shù)據(jù)集iris，如果你想用mtcars或者別的請(qǐng)隨意。比如想看汽車各種特性對(duì)油耗的影響就可以用mtcars。skim可以很華麗的展示你的數(shù)據(jù)結(jié)構(gòu)。當(dāng)然，沒(méi)啥別的用途了。
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 | ▂▇▆▅▇ |
隨時(shí)隨地的可視化
這是我認(rèn)為R語(yǔ)言跟spss和其他軟件比最大的優(yōu)勢(shì)了。是的，在Rstudio中你可以隨時(shí)隨地的可視化，無(wú)限制的切片數(shù)據(jù)可視化。相比于單純的統(tǒng)計(jì)分析，我認(rèn)為視覺(jué)往往來(lái)的更準(zhǔn)確更直接。話不多說(shuō)讓我們開(kāi)始吧。
首先，我們?cè)赗語(yǔ)言里進(jìn)行一些傳統(tǒng)的方法。我們進(jìn)行一個(gè)切片求平均值。這里我們使用了tidyverse套件（dplyr）。其中的 %>% 是通道符號(hào)，他的含義是將前面的參數(shù)傳入到下一個(gè)命令中作為第一個(gè)參數(shù)（可以用.代替）。統(tǒng)計(jì)各個(gè)組的平均值和標(biāo)準(zhǔn)偏差，過(guò)去大家都在用spss得到這些，統(tǒng)計(jì)分析，最后獲得結(jié)果，大家都滿意了，在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>
當(dāng)然你也可以這樣……
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
有了平均值和標(biāo)準(zhǔn)差，我們自然能畫(huà)出第一個(gè)柱狀圖。
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
重復(fù)四次，你能獲得四個(gè)柱狀圖
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')
然后你可以做統(tǒng)計(jì)分析。別忘了檢查各個(gè)品種的變量是否滿足正態(tài)分布，會(huì)發(fā)現(xiàn)正態(tài)性檢驗(yàn)還好。除了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
方差齊性檢驗(yàn)。發(fā)現(xiàn)各組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
可以換成箱線圖，再加上散點(diǎn)，標(biāo)記好顯著性。因?yàn)橹鞍l(fā)現(xiàn)有點(diǎn)不滿足正態(tài)和方差齊性，所以這里選擇了非參數(shù)檢驗(yàn)。
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()
如果你覺(jué)得wilcox.test還不夠穩(wěn)健，你可以試試Welch's ANOVA檢驗(yàn)和Games-Howell事后檢驗(yàn)
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)
篇幅原因這里只放一部分結(jié)果
##
## 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語(yǔ)言對(duì)數(shù)據(jù)進(jìn)行可視化
到了這里，事情解決了，對(duì)么？很傳統(tǒng)的方法，得到了想知道的一切。數(shù)據(jù)有差異性，完美。但其實(shí)，R語(yǔ)言能夠做的更多。我更愿意稱R語(yǔ)言為可視化統(tǒng)計(jì)檢驗(yàn)方法。相比于傳統(tǒng)方法的先檢驗(yàn)差異性再作圖，R語(yǔ)言是先做可視化再做具體的方差分析。例如，我們可以一行代碼查看數(shù)據(jù)的關(guān)系。這里面包含了數(shù)據(jù)的相關(guān)散點(diǎn)圖，直方圖，相關(guān)系數(shù)和箱線圖。依靠這些你可以更好的提出你需要的科學(xué)假設(shè)。
ggpairs(data, mapping = aes(color = Species)) + theme_bw()
通過(guò)直方圖可以很容易看到數(shù)據(jù)的分布情況。是否滿足正態(tài)還有方差齊性是不是變得更加具象了？
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
跑個(gè)qq圖，配合置信區(qū)間。
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")
然后是數(shù)據(jù)結(jié)果的分組可視化。
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')
還可以換個(gè)炫酷的效果。這里用的Games_Howell test。不喜歡可以換回t檢驗(yàn)，大家可以自己試試看。
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查看各個(gè)數(shù)據(jù)在四個(gè)維度上的表現(xiàn)，這里我因?yàn)橥瑫r(shí)裝了兩個(gè)包所以，沖突了。所以需要申明到底用的是哪個(gè)包的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)
最后根據(jù)四個(gè)petal和sepal的長(zhǎng)寬數(shù)據(jù)，我們做個(gè)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建模
到這里，我們體驗(yàn)了大部分統(tǒng)計(jì)可視化的內(nèi)容。那么R還能做什么？既然我們知道所謂統(tǒng)計(jì)分析本身就是建模，那么我們能否用更復(fù)雜的模型對(duì)petal和sepal的長(zhǎng)寬等因子進(jìn)行重要性檢驗(yàn)?zāi)?？所以又有了這部分，利用tidymodels的建模，并根據(jù)模型對(duì)數(shù)據(jù)分類。
library('tidymodels')
library('caret')
#固定隨機(jī)數(shù)，方便復(fù)現(xiàn)
set.seed(2020)
下面是整個(gè)建模的流程
#將數(shù)據(jù)分成3份兩份作為訓(xùn)練數(shù)據(jù)，一份作為測(cè)試準(zhǔn)確性數(shù)據(jù)。
data_split <- initial_split(data,prop=.66)
data_train <- training(data_split)
data_test <-testing(data_split)
#bootstrap創(chuàng)造一個(gè)數(shù)據(jù)用來(lái)tuning模型
Spec_boost <- bootstraps(data_train, times = 30)
#設(shè)定模型數(shù)據(jù)轉(zhuǎn)換，在這里可以做中心化，標(biāo)準(zhǔn)化，數(shù)據(jù)合并等一系列操作。最后傳入的prep記錄所有的操作。另外在這里Species~.意思是Species作為分類結(jié)果，其他幾列作為predictor預(yù)測(cè)分類。
data_rec <-
recipe(Species ~.,data = data_train) %>%
prep()
#查看訓(xùn)練的數(shù)據(jù)，這里用juice很形象的說(shuō)你的數(shù)據(jù)是榨出來(lái)的。但這里沒(méi)對(duì)表型數(shù)據(jù)做過(guò)多調(diào)整。
juice(data_rec)
#設(shè)定模型，這里mtry和min_n使用tune函數(shù)調(diào)整。模式選擇分類，隨機(jī)森林的引擎選擇randomForest，另外還有一個(gè)引擎叫ranger。具體到后面會(huì)有細(xì)微差別。
rf_model<-rand_forest(mtry=tune(),min_n = tune())%>%
set_mode("classification")%>%
set_engine("randomForest")
#設(shè)定工作流，主要是使用哪個(gè)模型以及使用哪個(gè)數(shù)據(jù)方便后期做具體的調(diào)參。
rf_wflow <-
workflow() %>%
add_recipe(data_rec)%>%
add_model(rf_model)
#激活平行線程
doParallel::registerDoParallel()
#使用tune_grid函數(shù)和boostrap的數(shù)據(jù)進(jìn)行大致的調(diào)參。
rf_results <-
rf_wflow %>%
tune_grid(resamples = Spec_boost)
#參數(shù)打分可視化
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")
#選擇最佳模型參數(shù)
so_best <-
rf_results %>%
select_best(metric = "roc_auc")
#構(gòu)建模型
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()
#設(shè)定參數(shù)
rf_grid <- grid_regular(
mtry(range = c(1,4)),
min_n(range = c(10,30)),
levels = 10
)
#調(diào)參結(jié)果
regular_res <- rf_wflow %>%
tune_grid(
resamples = Spec_boost,
grid = rf_grid
)
#結(jié)果可視化
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")
#調(diào)參后通過(guò)roc_auc指標(biāo)選擇最好的模型參數(shù)
so_best <-
regular_res %>%
select_best(metric = "roc_auc")
#使用最佳參數(shù)去構(gòu)建模型
rf_final_fit<-rf_wflow%>%
finalize_workflow(so_best)%>%
fit(data = data_train)
#設(shè)定好數(shù)據(jù)調(diào)整的方法（還記得recipes么）
data_rec3 <-rf_final_fit%>%
pull_workflow_prepped_recipe()
#用同樣方法去變換test數(shù)據(jù)集，然后用模型擬合test數(shù)據(jù)集，輸出結(jié)果。
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
#查看參數(shù)重要性
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()
