本次筆記內(nèi)容：
- 用
ggpubr的stat_compare_means()給boxplot加上significant labels- 編寫(xiě)函數(shù)選擇P值小于某一值的pairs，并只在圖上標(biāo)注這些pairs
- 多圖使用
ggpubr的ggarrange()和annotate_figure()整理，注意aes_string()的使用- 把一個(gè)List的plots:
list(plot1, plot2, plot3...)拼在一起
對(duì)手頭的數(shù)據(jù)，你需要對(duì)每個(gè)變量按照固定分組畫(huà)一個(gè)boxplot，并加上Kruskal-Wallis多組檢驗(yàn)結(jié)果及post-hoc兩兩比較結(jié)果，其中兩兩比較的結(jié)果只需要有顯著性的。以此畫(huà)多個(gè)圖后將他們“拼圖”成一個(gè)。
以iris數(shù)據(jù)集為例，根據(jù)以下代碼用ggplot2繪制了"Sepal.Width"變量的基本箱線圖，并用ggpubr包的stat_compare_means()加上了significant labels。這里是加上了所有組合的labels。combn()用于獲取所有不重復(fù)的兩兩組合pairs。注意comparisons = 這里需要傳遞一個(gè)list(), list中包含了以c()儲(chǔ)存的pairs，你也可以自己指定pairs是什么，比方說(shuō)comparisons = list(c(a,b), c(b,c))
p <- ggplot(iris, aes(x = Species, y = Sepal.Width, fill = Species)) +
geom_boxplot(position = position_dodge(0.8)) +
geom_point(position = position_jitterdodge()) +
scale_fill_aaas() +
theme_classic(base_size = 16)+
labs(x = "", y = 'Sepal.Width') +
stat_compare_means() +
stat_compare_means(comparisons = combn(levels(iris$Species), 2, simplify =FALSE))
在iris數(shù)據(jù)集中組間都很顯著，但是如果是很多組在Kruskal-Wallis(或者ANOVA)檢驗(yàn)后的兩兩比較(post-hoc)，可能會(huì)出現(xiàn)boxplot上有一堆label, 圖被壓縮的很丑。而且主要是想把顯著的label標(biāo)注上去。
首先stat_compare_means()是通過(guò)pairwise.wilcox.test做post-hoc兩兩比較的：
> test <- pairwise.wilcox.test(iris$Sepal.Width, iris$Species, p.adjust.method = "none")
> test
Pairwise comparisons using Wilcoxon rank sum test
data: iris$Sepal.Width and iris$Species
setosa versicolor
versicolor 2.1e-13 -
virginica 7.1e-09 0.0046
P value adjustment method: none
于是可以寫(xiě)個(gè)函數(shù)選擇p < "你想要一個(gè)cutoff值" 的pairs, 為了防止因?yàn)闆](méi)有p值小于你想要的cutoff值而報(bào)錯(cuò)，在沒(méi)有的時(shí)候返回一個(gè)空l(shuí)ist()。
which_pair_to_use <- function(data, variable, group, cutoff = 0.1, p.adjust_method = "none") {
pair_p <- pairwise.wilcox.test(data[[variable]], data[[group]], p.adjust.method = p.adjust_method)
ind <- as.matrix(which(pair_p$p.value < cutoff, arr.ind = TRUE))
if (nrow(ind) > 0) {
compairs <- list()
for (i in 1:nrow(ind)){
pair <- c(rownames(pair_p$p.value)[ind[i,"row"]], colnames(pair_p$p.value)[ind[i,"col"]])
compairs[[i]] <- pair
}
return(compairs)
} else {
return(list())}
}
注意由于是把變量的名稱字符串傳遞給上述函數(shù)，所以參數(shù)需要用引號(hào)括起來(lái)，cutoff =一般來(lái)說(shuō)是0.05，這里由于p值都太小了所以我用一個(gè)比較小的值來(lái)舉例。 為了可以畫(huà)多個(gè)圖方便，可以簡(jiǎn)單寫(xiě)成一個(gè)函數(shù)。注意在把ggplot2寫(xiě)進(jìn)函數(shù)的時(shí)候，給aes()傳參數(shù)會(huì)出現(xiàn)問(wèn)題。因?yàn)?code>aes()中并不是引號(hào)括起來(lái)的字符串。這時(shí)需要使用aes_string()
fun_to_plot <- function(data, group, variable) {
p <- ggplot(data, aes_string(x = group, y = variable, fill = group)) +
geom_boxplot(position = position_dodge(0.8)) +
geom_point(position = position_jitterdodge()) +
scale_fill_aaas() +
theme_classic(base_size = 16)+
labs(x = "", y = variable) +
stat_compare_means() +
stat_compare_means(comparisons = which_pair_to_use(data,
variable,
group,
cutoff = 0.0005))
return(p)
}
fun_to_plot(iris, "Species", "Sepal.Width")
#得到下圖
接下來(lái)把Iris中Sepal的width和length畫(huà)在同一個(gè)圖里，并且共用同一個(gè)legend。并且給合并的這個(gè)圖加上大title。用ggpubr的ggarrange()來(lái)組合兩個(gè)圖，用annotate_figure()來(lái)給組合好的圖加以修飾。
annotate_figure(
ggarrange(fun_to_plot(iris,"Species", "Sepal.Width"),
fun_to_plot(iris,"Species","Sepal.Length"),
common.legend = TRUE, legend = "right"),
top = text_grob("Sepal_species", size = 18)
)
另外如果你用函數(shù)畫(huà)出了lists of plots: 即 plots_list = list(plot1, plot2, plot3...)這樣的形式?？梢杂?code>grid.arrange(grobs = plots_list, ncol = 7)把他們拼在一起。ncol = 7即一行放7個(gè)圖。
除了使用ggpubr來(lái)給圖添加顯著label以外，還有ggsignif包。（可以了解一下。我如果用到了會(huì)在這里加上）
參考鏈接：
annotate_figure() , ggarrange()
ggplot2傳參數(shù)：aes_string()
arrange a list of plot
ggsignif參考鏈接
