厚缊 業(yè)余的R語言可視化重度患者 個人博客:houyun.xyz 郵箱:houyunhuang@163.com 轉(zhuǎn)載本文(包括長期轉(zhuǎn)載賬號)必須聯(lián)系厚缊授權(quán) 相關(guān)系數(shù)矩陣可視化已經(jīng)至少有兩個版本的實現(xiàn)了,魏太云基于base繪圖系統(tǒng)寫了corrplot 包,應(yīng)該說是相關(guān)這個小領(lǐng)域中最精美的包了,使用簡單,樣式豐富,只能用驚艷來形容。Kassambara的ggcorrplot 基于ggplot2 重寫了corrplot ,實現(xiàn)了corrplot 中絕大多數(shù)的功能,但僅支持“square”和“circle”的繪圖標(biāo)記,樣式有些單調(diào),不過整個ggcorrplot 包的代碼大概300行,想學(xué)習(xí)用ggplot2 來自定義繪圖函數(shù),看這個包的源代碼很不錯。還有部分功能相似的corrr 包(在寫ggcor 之前完全沒有看過這個包,寫完之后發(fā)現(xiàn)在相關(guān)系數(shù)矩陣變data.frame方面驚人的相似),這個包主要在數(shù)據(jù)相關(guān)系數(shù)提取、轉(zhuǎn)換上做了很多的工作,在可視化上稍顯不足。ggcor 的核心是為相關(guān)性分析、數(shù)據(jù)提取、轉(zhuǎn)換、可視化提供一整套解決方案,目前的功能大概完成了70%,后續(xù)會根據(jù)實際需要繼續(xù)擴展。這里要介紹的ggcor 是corrplot 的有一種實現(xiàn),在吸收借鑒(或者說是全般)corrplot 的基礎(chǔ)上,略有提升,使用上會更靈活簡單。這里不得不說下整個開發(fā)過程,一方面是為了感謝來自五湖四海的朋友提的許許多多的專業(yè)意見,另一方面是深感歉意,由于我考慮不周,ggcor (第一次上傳github測試大約在十天前,名字還是ggcorrr )存在各種各樣的小bug,也有很多限制,這次為了升級和簡化,不僅名字變了,里面的函數(shù)也變得完全不一樣,若你們以前寫了ggcorrr 的代碼,現(xiàn)在基本不可以重用了。盡管這一版本的更新我考慮了很多應(yīng)用場景,綜合了方方面面的意見,但是bug和限制還是在所難免,希望朋友們繼續(xù)提意見,提需求。當(dāng)然,這一次有兩點是可以保證的:一是包名不會變了,再變都成精了;二是包內(nèi)的主體函數(shù)和參數(shù)不會變了,而且即使要變,我也會考慮兼容性。安裝if(!require(devtools)) install.packages('devtools') if(!require(ggcor)) devtools::install_github('houyunhuang/ggcor')
數(shù)據(jù)預(yù)處理函數(shù)ggplot2 要求的數(shù)據(jù)格式是data.frame ,要把相關(guān)系數(shù)矩陣處理成理想的數(shù)據(jù)格式需要一系列的操作。ggcor 提供了兩個主要的函數(shù),一個是as_cor_tbl() 函數(shù),另一個是fortify_cor() 函數(shù),兩者結(jié)合應(yīng)該能滿足絕大多數(shù)的應(yīng)用場景需求。as_cor_tbl() 是更底層的函數(shù),fortify_cor() 本質(zhì)上調(diào)用了as_cor_tbl() 來得到最終的數(shù)據(jù)結(jié)果。這個函數(shù)適用于已經(jīng)知道,或者需要用其它更特殊的函數(shù)(非stats::cor() )來處理得到系數(shù)的情況,常用的參數(shù)是前三個。x —— 相關(guān)系數(shù)矩陣(或者數(shù)據(jù)框),矩陣行名和列名是必要的,若沒有或者缺失值會自動補全名字,行名以“Y”開頭,附上遞增的整數(shù)序列,列名以“X”開頭,附上附上遞增的整數(shù)序列。type —— 相關(guān)系數(shù)矩陣圖樣式,“upper”截斷下三角,“l(fā)ower”截斷上三角。show.diag —— 相關(guān)系數(shù)矩陣圖中是否包含對角線,僅對對稱矩陣有效。p —— 相關(guān)系數(shù)檢驗p值矩陣(或者數(shù)據(jù)框),必須與x 一一對應(yīng)。low —— 相關(guān)系數(shù)置信區(qū)間下界矩陣(或者數(shù)據(jù)框),必須與x 一一對應(yīng)。upp —— 相關(guān)系數(shù)置信區(qū)間上界矩陣(或者數(shù)據(jù)框),必須與x 一一對應(yīng)。cluster.type —— 是否對相關(guān)系數(shù)矩陣進行重新排序,“none”表示不重排,“all”表示行列均重排,“row”表示對行重排,“col”則只對列重排。... —— 其它傳遞給matrix_order() 函數(shù)的參數(shù)。
library(ggcor) ## function(x, ## type = c('full', 'upper', 'lower'), ## show.diag = TRUE, ## p = NULL, ## low = NULL, ## upp = NULL, ## cluster.type = c('none', 'all', 'row', 'col'), ## ...)
corr <- cor(mtcars) df <- as_cor_tbl(corr) df ## return a tibble
## # A tibble: 121 x 3 ## x y r ## * <int> <int> <dbl> ## 1 1 11 1 ## 2 1 10 -0.852 ## 3 1 9 -0.848 ## 4 1 8 -0.776 ## 5 1 7 0.681 ## # … with 116 more rows ## Extra attributes: ## xname: mpg cyl disp hp drat wt qsec vs am gear carb ## yname: carb gear am vs qsec wt drat hp disp cyl mpg ## show.diag: TRUE
fortify_cor() 主要適用于處理原始數(shù)據(jù)表,即調(diào)用cor() 求相關(guān)系數(shù),cor 函數(shù)對數(shù)據(jù)按列進行兩兩相關(guān)性計算,默認使用pearson 方法,當(dāng)然理論解讀中提前的spearman 和kendall 方法也都支持。調(diào)用(若需要)cor.test() 函數(shù)進行統(tǒng)計檢驗,并返回ggcor 需要的數(shù)據(jù)類。x —— 原數(shù)據(jù)矩陣(或者數(shù)據(jù)框),列名是必要的,若沒有或者缺失值會自動補全名字,列名以“X”開頭,附上附上遞增的整數(shù)序列。y —— 原數(shù)據(jù)矩陣(或者數(shù)據(jù)框),列名是必要的,若沒有或者缺失值會自動補全名字,列名以“X”開頭,附上附上遞增的整數(shù)序列。當(dāng)y 不為空(NULL)時,相關(guān)系數(shù)是x 中的每一列和y 中的每一列的相關(guān)性。type —— 相關(guān)系數(shù)矩陣圖樣式,“upper”截斷下三角,“l(fā)ower”截斷上三角。show.diag —— 相關(guān)系數(shù)矩陣圖中是否包含對角線,僅對對稱矩陣有效。cor.test —— 邏輯值,是否進行相關(guān)性檢驗。cor.test.alt —— 相關(guān)性檢驗備擇假設(shè),詳細請查看cor.test() 幫助。cor.test.method —— 相關(guān)性檢驗方法,詳細請查看cor.test() 幫助。cluster.type —— 是否對相關(guān)系數(shù)矩陣進行重新排序,“none”表示不重排,“all”表示行列均重排,“row”表示對行重排,“col”則只對列重排。cluster.method —— 當(dāng)cluster.order 為“HC”(默認)時算法,詳細請查看ggcor::matrix_order() 。... —— 其它傳遞給cor() 函數(shù)的參數(shù)。
## function( ## x, ## y = NULL, ## type = c('full', 'upper', 'lower'), ## show.diag = FALSE, ## cor.test = FALSE, ## cor.test.alt = 'two.sided', ## cor.test.method = 'pearson', ## cluster.type = c('none', 'all', 'row', 'col'), ## cluster.method = 'HC', ## ... )
df01 <- fortify_cor(mtcars, cor.test = TRUE, cluster.type = 'all') df01
## # A tibble: 121 x 6 ## x y r p low upp ## * <int> <int> <dbl> <dbl> <dbl> <dbl> ## 1 1 11 1 7.44e-232 1.000 1.000 ## 2 1 10 0.750 7.83e- 7 0.543 0.871 ## 3 1 9 0.428 1.46e- 2 0.0927 0.676 ## 4 1 8 0.527 1.94e- 3 0.218 0.740 ## 5 1 7 0.395 2.53e- 2 0.0537 0.654 ## # … with 116 more rows ## Extra attributes: ## xname: carb hp wt cyl disp drat am gear qsec mpg vs ## yname: vs mpg qsec gear am drat disp cyl wt hp carb ## show.diag: TRUE
as_cor_tbl() 和fortify_cor() 返回值均是cor_tbl 類的數(shù)據(jù)框(準確的說是tibble ),并包含其它額外特殊屬性,要充分利用ggcor 包中一系列工具函數(shù)帶來的便捷性,必須調(diào)用(手動比較麻煩)這兩個函數(shù)預(yù)處理數(shù)據(jù)。初始化繪圖函數(shù)ggcor 封裝了一個基本的初始化函數(shù)ggcor() ,用來處理數(shù)據(jù)、繪圖類型、背景、坐標(biāo)軸、顏色映射、圖例等??瓷先?/span>ggcor() 函數(shù)非常復(fù)雜,但若仔細觀察參數(shù)命名規(guī)則和分類,又是非常簡單的。繪圖區(qū)網(wǎng)格線參數(shù)名字均以grid.* 開頭;坐標(biāo)軸相關(guān)參數(shù)均以axis.* 開頭;圖例相關(guān)(主要是colourbar)的參數(shù)均以legend.* 開頭;還有panel.backgroud 、xlim 、ylim 均是常見的參數(shù),panel.backgroud 參數(shù)用來設(shè)置繪圖區(qū)背景顏色,xlim 、ylim 則是設(shè)置x/y軸的范圍。把這幾大塊參數(shù)去掉后,剩下的參數(shù)只有6個了,瞬間清爽了不少。要完全理解ggcor() 的作用原理及相關(guān)參數(shù)的設(shè)置,需要先講講ggcor() 內(nèi)部構(gòu)成。ggcor() 本質(zhì)上是調(diào)用了ggplot() 來初始化,然后根據(jù)相關(guān)系數(shù)圖的樣式添加了一些輔助的圖層。x 、y 、mapping 、is.cor 、show.diag 和... 參數(shù)均和數(shù)據(jù)預(yù)處理和映射相關(guān)。x 可以是cor_tbl 、矩陣(數(shù)據(jù)框)。當(dāng)為cor_tbl 時直接作為data 參數(shù)傳遞給ggplot() ;為矩陣(數(shù)據(jù)框)時,若是(is.cor = TRUE )相關(guān)系數(shù)矩陣(數(shù)據(jù)框)時,調(diào)用as_cor_tbl() 函數(shù)處理成cor_tbl ,若不是(is.cor = FALSE )相關(guān)系數(shù)矩陣(數(shù)據(jù)框)時,調(diào)用fortify_cor() 函數(shù)處理成cor_tbl ,此時x 、y 、show.diag 和... 均傳遞給as_cor_tbl() 或者fortify_cor() 。mapping 對應(yīng)ggplot() 中的mapping 參數(shù),當(dāng)為空(默認)時,根據(jù)cor_tbl 中的變量情況添加,基礎(chǔ)形式是aes(x = x, y = y, r = r, fill = r) 。若cor_tbl 包含“p”(進行了相關(guān)系數(shù)顯著性檢驗),則最基礎(chǔ)形式基礎(chǔ)上額外添加p = p ,若檢驗方法(cor.test.method = 'pearson' ),再加上low = low 和upp = upp 。
fill.* 均是fill顏色映射函數(shù)相關(guān)的參數(shù)。- 若
fill.scale.add 為FALSE 不添加顏色映射函數(shù)。若為TRUE (默認),則會在初始化中自動添加顏色映射函數(shù)。 fill.colours 是顏色字符向量,用來控制填充色(fill)的顏色映射,若為空(默認),使用corrplot 默認配色,可以通過ggcor:::.defualt.colors 查看具體顏色。fill.bin 是否對連續(xù)性顏色進行分組(默認FALSE )。當(dāng)fill.scale.add = TRUE 時,若fill.bin = TRUE ,ggcor() 的顏色映射函數(shù)是使用scale_fill_steps2n() ,若fill.bin = FALSE ,ggcor() 的顏色映射函數(shù)是使用scale_fill_gradient2n() 。
- 背景網(wǎng)格線是通過
geom_tile() (目前來看,geom_tile() 會存在一些難以處理的缺點,以后可能會添加專用的網(wǎng)格線圖層函數(shù)),panel.backgroud (fill)、grid.* 系列參數(shù)傳遞給geom_tile() 。 legend.position 傳遞給theme() 中對應(yīng)的參數(shù),用來控制圖例位置,其它legend.* 開頭參數(shù)傳遞給guide_colourbar() 或者guide_colorsteps() 。legend.breaks 用來控制圖例顏色棒標(biāo)簽顯示位置,legend.labels 是對應(yīng)的標(biāo)簽。若fill.bin = TRUE ,legend.breaks 也是圖例顏色棒切割分組的位置。coord.fixed 邏輯值,為真xlim 、ylim 傳遞給coord.fixed() 函數(shù),為假傳遞給coord_cartesian() 函數(shù)。在ggcor 包中,相關(guān)系數(shù)矩陣若是n * m的矩陣,那么第i行對應(yīng)的坐標(biāo)點(即as_cor_tbl() 返回結(jié)果中的y)為n-i(為了和表格呈現(xiàn)樣式一致,行方向翻轉(zhuǎn)了),第j列對應(yīng)的坐標(biāo)點(即as_cor_tbl() 返回結(jié)果中的x)為j,得到第(i, j)個數(shù)據(jù)點所在的方格坐標(biāo)為(xmin = j-0.5, xmax = j+0.5, ymin = n-i-0.5, ymax = n-i+0.5)。從這個規(guī)律我們不難算出默認的xlim 、ylim 取值范圍。注意,對cor_tbl 數(shù)據(jù)框,或者說是ggcor 包中函數(shù)操作數(shù)據(jù),均不會改變每個數(shù)據(jù)單元格在圖中的坐標(biāo)位置。
ggcor() 初始化之后,本質(zhì)上返回的是ggplot 對象,若是想改變默認設(shè)置,可以按照ggplot2 的相應(yīng)的函數(shù)和設(shè)置方法去調(diào)整。## function( ## x, ## y = NULL, ## mapping = NULL, ## is.cor = FALSE, ## show.diag = TRUE, ## fill.scale.add = TRUE, ## fill.colours = NULL, ## fill.bin = FALSE, ## panel.backgroud = NA, ## grid.colour = 'grey50', ## grid.size = 0.25, ## grid.linetype = 'solid', ## axis.x.position = c('auto', 'bottom', 'top'), ## axis.y.position = c('auto', 'left', 'right'), ## fill.colours = NULL, ## legend.title = 'correalation', ## legend.position = 'auto', ## legend.breaks = NULL, ## legend.labels = NULL, ## coor.fixed = TRUE, ## xlim = NULL, ## ylim = NULL, ## ...) 我寫ggcor() 函數(shù)花的時間最長,耗費的精力最多,但到現(xiàn)在仍然是半成品的感覺,有些問題直到現(xiàn)在我也沒找到相對完美的解決方案。盡管如此,對于新手,我還是建議調(diào)用ggcor() 來進行初始化,若自己去研究各種圖層函數(shù),折騰很多細節(jié),一天也難得出一幅圖,對于心里的打擊比較大??磶讉€初始化之后的效果。ggcor(mtcars)
ggcor(mtcars, type = 'lower')
df02 <- fortify_cor(mtcars, type = 'upper') ggcor(df02, panel.backgroud = '#66C2A5')
圖層函數(shù)ggcor 提供了定制的geom_square() 、geom_circle2() 、geom_ellipse2() 、geom_pie2() 、geom_colour() 、geom_confbox() 、geom_num() 、geom_mark() 、geom_cross() 9個ggplot2 圖層函數(shù),可以根據(jù)需要進行疊加。除了ggplot2 中一般化的參數(shù)(x 、y 、fill 、colour 、size 等)最常用參數(shù)r 、p 、low 、upp 、num 、r0 、sig.thres 、sig.level 、mark 等。r —— 相關(guān)系數(shù),適用于geom_square() 、geom_circle2() 、geom_ellipse2() 、geom_pie2() 、geom_confbox() 、geom_mark() (統(tǒng)計顯著性標(biāo)記的系數(shù),與 1.234??類似)。這些參數(shù)之所以都設(shè)置為“r”,主要是因為在相關(guān)系數(shù)可視化中基本都映射為相關(guān)系數(shù),統(tǒng)一命名可以減少一些參數(shù)記憶,方便使用。p —— 相關(guān)系數(shù)檢驗P值,適用于geom_mark() 、geom_cross() ,結(jié)合sig.thres 等參數(shù)來根據(jù)顯著性水平做一些輔助標(biāo)記。r 、low 、upp —— 適用于geom_confbox() ,三個參數(shù)均是必須的,low 、r 、upp 分別確定置信區(qū)間盒子的下端、中間和上端線條位置。num —— 適用于geom_num() ,數(shù)值變量。這個函數(shù)封裝于geom_text() ,做了一點點優(yōu)化,以后可能會刪除。
r0 ,外接圓半徑縮放系數(shù),適用于geom_square() 、geom_circle2() 、geom_ellipse2() 、geom_pie2() 函數(shù)。該參數(shù)的主要意義是處理圖形覆蓋問題,當(dāng)在每個單元格畫半徑為0.5的方塊、圓等圖標(biāo)時,會相互覆蓋掉背景網(wǎng)格線,影響視覺效果。該參數(shù)默認值是0.48。sig.thres 統(tǒng)計顯著性臨界值,用來過濾非顯著的數(shù)據(jù)行。適用于geom_mark() 、geom_cross() 。sig.level 、mark 適用于geom_mark() ,前者為統(tǒng)計顯著性水平向量(如c(0.001, 0.01, 0.05) ),后者為對應(yīng)的標(biāo)記符號(c('***', '**', '*') )。統(tǒng)計顯著性水平向量不要求一定要按順序排列,只要求和標(biāo)記符號一一對應(yīng)就行。
ggcor(mtcars) + geom_square()
ggcor(mtcars, type = 'upper') + geom_circle2()
ggcor(mtcars, type = 'lower', show.diag = TRUE) + geom_ellipse2()
ggcor(mtcars, type = 'full', cluster.type = 'all') + geom_pie2()
ggcor(mtcars, cluster.type = 'all') + geom_colour() + geom_num(aes(num = r), colour = 'grey90', size = 3.5)
ggcor(mtcars, type = 'full', cor.test = TRUE) + geom_confbox()
ggcor(mtcars, type = 'full', cor.test = TRUE, cluster.type = 'all') + geom_colour() + geom_cross()
ggcor(mtcars, type = 'full', cor.test = TRUE) + geom_square() + geom_cross() 其它圖層函數(shù)的使用都比較符合直覺,需要設(shè)置的地方也很少,而geom_mark() 會涉及一些其它問題,這里拉出來說說。很多時候,我們并不關(guān)心不具備統(tǒng)計顯著性的相關(guān)系數(shù),也不需要在圖中顯示,這時需要設(shè)置sig.thres ,即要過濾的顯著性臨界值。ggcor(mtcars, type = 'full', cor.test = TRUE, cluster.type = 'all') + geom_raster() + geom_mark(sig.thres = 0.05, size = 3, colour = 'grey90') 若是只要統(tǒng)計顯著性標(biāo)記的'*'號,不要系數(shù)怎么整?在geom_mark() 中也很簡單,直接設(shè)置r 就好。ggcor(mtcars, type = 'full', cor.test = TRUE, cluster.type = 'all') + geom_raster() + geom_mark(r = NA, sig.thres = 0.05, size = 5, colour = 'grey90')
那要改變星號標(biāo)記規(guī)則,只要小于0.05的標(biāo)記一顆星,其它什么都不標(biāo)記呢?注意:因為星號在文本中顯示在偏上的位置,若不設(shè)置vjust 參數(shù),看上去縱向會不居中。ggcor(mtcars, type = 'full', cor.test = TRUE, cluster.type = 'all') + geom_raster() + geom_mark(r = NA, sig.level = 0.05, mark = '*', vjust = 0.65, size = 6, colour = 'grey90')
非對稱相關(guān)系數(shù)矩陣非對稱相關(guān)系數(shù)矩陣和非對稱矩陣是有細微的區(qū)別的,前者表示行列代表不同的變量集合,相互之間的順序可以打亂。所以,有時候要分析兩個表中每個變量之間的相關(guān)性,此時得到的結(jié)果就是非對稱的相關(guān)系數(shù)矩陣。library(vegan) # 使用vegan包所帶的數(shù)據(jù)集 data(varechem) data(varespec) df03 <- fortify_cor(x = varechem, y = varespec[ , 1:30], cluster.type = 'col') ggcor(df03) + geom_colour()
df04 <- fortify_cor(x = varespec[ , 1:30], y = varechem, cor.test = TRUE) ggcor(df04) + geom_square() + geom_cross(size = 0.2)# 數(shù)據(jù)集不好,沒幾個顯著的
上下三角不一樣怎么畫?ggcor 提供了很多輔助函數(shù)來對cor_tbl 數(shù)據(jù)進行過濾,函數(shù)命名規(guī)則上都以get_* 開頭。get_lower_data() —— 獲取相關(guān)系數(shù)矩陣下三角所在行,僅支持對稱的相關(guān)系數(shù)矩陣。get_upper_data() —— 獲取相關(guān)系數(shù)矩陣上三角所在行,僅支持對稱的相關(guān)系數(shù)矩陣。get_diag_data() —— 獲取相關(guān)系數(shù)矩陣對角線所在行,僅支持對稱的相關(guān)系數(shù)矩陣。get_diag_tri() —— 刪除相關(guān)系數(shù)矩陣對角線所在行,僅支持對稱的相關(guān)系數(shù)矩陣。get_data() —— 是以上四個函數(shù)的重新包裝,主要在畫圖時使用,稍后通過例子詳細說明。
df05 <- fortify_cor(x = varechem, cor.test = TRUE, cluster.type = 'all') ggcor(df05) + geom_circle2()
df05_lower <- get_lower_data(df05, show.diag = FALSE) ggcor(df05_lower) + geom_circle2()
ggcor(df05) + geom_pie2(data = get_data(type = 'upper', show.diag = FALSE)) + geom_ellipse2(data = get_data(type = 'lower', show.diag = TRUE))
ggcor(df05) + geom_segment(aes(x = x - 0.5, y = y + 0.5, xend = x + 0.5, yend = y - 0.5), data = get_data(type = 'diag'), size = 0.5, colour = 'grey60') + geom_colour(data = get_data(type = 'upper', show.diag = FALSE)) + geom_mark(data = get_data(type = 'upper', show.diag = FALSE), size = 3) + geom_circle2(data = get_data(r >= 0.5, type = 'lower', show.diag = FALSE), r = 0.8, fill = '#66C2A5') + geom_num(aes(num = r), data = get_data(type = 'lower', show.diag = FALSE), size = 3)
列名放在對角線?
ggcor(mtcars, cor.test = TRUE, cluster.type = 'all') + geom_confbox(data = get_data(type = 'upper', show.diag = FALSE)) + geom_num(aes(num = r), data = get_data(type = 'lower', show.diag = FALSE), size = 3.5) + add_diaglab(size = 4.56) + remove_axis()
想對顏色分組?
很多情況下,連續(xù)性顏色棒并不是很好分區(qū)每個單元格對應(yīng)的數(shù)值區(qū)間,這時根據(jù)相關(guān)系數(shù)大小對顏色進行分組可能更適合。說來非常巧合,一直頭疼這個問題難以解決,就在前不久,ggraph 包的作者Thomas Lin向ggplot2 包提交了一組新的顏色映射函數(shù)(scale_fill/colour_steps*() ),這個問題迎刃而解,非常幸運。ggcor() 函數(shù)有參數(shù)fill.binned ,默認為FALSE ,設(shè)置為TRUE 就會根據(jù)相關(guān)系數(shù)大小對顏色分組。若要控制分組的數(shù)量和區(qū)間,可以通過legend.breaks 來設(shè)置。ggcor(mtcars, fill.bin = TRUE) + geom_square() # 默認分組
ggcor(mtcars, fill.bin = TRUE, legend.breaks = seq(-1, 1, length.out = 11)) + geom_square() #指定分組,0.2為一個區(qū)間
col <- col <- ggcor:::.default_colors col[6] <- '#F2F2F2' ggcor(mtcars, cluster.type = 'all', fill.colours = col, fill.bin = T, legend.breaks = c(-1, -0.8, -0.5, 0.5, 0.8, 1)) + geom_colour()
玩點花活
這部分內(nèi)容要在線下載表情,很多時候會因為網(wǎng)絡(luò)問題下載失敗。不給圖了。library(ggimage) emoji <- c('1f004', '1f0cf', '1f170', '1f171', '1f17e', '1f17f', '1f18e', '1f191', '1f192', '1f193', '1f194', '1f195', '1f196', '1f197') ggcor(df05) + geom_pokemon(aes(image=ifelse(r > 0.5, 'pikachu', 'tauros')), data = get_data(type = 'lower', show.diag = FALSE)) + geom_emoji(aes(image = ifelse(p <= 0.05, '1f600', '1f622')), data = get_data(type = 'upper', show.diag = FALSE)) + geom_emoji(aes(image = emoji), data = get_data(type = 'diag'))
ggcor(df05) + geom_pokemon(aes(image=ifelse(r > 0.5, 'pikachu', 'tauros')), data = get_data(type = 'lower', show.diag = FALSE)) + geom_colour(data = get_data(type = 'upper', show.diag = FALSE)) + geom_shade(data = get_data(type = 'upper', show.diag = FALSE), sign = -1, size = 0.1) + geom_emoji(aes(image = emoji), data = get_data(type = 'diag'))
mantel 檢驗組合圖
mantel 檢驗(Mantel test 是對兩個矩陣相關(guān)關(guān)系的檢驗)的組合圖已經(jīng)十分流行了,用各種工具做的都有。大概5月份的時候,我基于corrplot 模擬重現(xiàn)了那幅圖,直到現(xiàn)在每周都有人詢問我相關(guān)實現(xiàn)的問題,我基本都是回答說等新方案,因為那個實現(xiàn)很復(fù)雜,沒有基本的R知識,很難替換成自己的數(shù)據(jù)。8月底,我在ggcorrr (ggcor 前身)開發(fā)了gghpcc ,幾乎可以滿足可視化的要求,但是在優(yōu)化ggcor 的過程中,我才恍然大悟,mantel檢驗本質(zhì)上也是相關(guān)性分析,何不統(tǒng)一整合進ggcor 呢,就這樣,gghpcc 被干掉了。進行案例展示之前,我先解釋一個關(guān)鍵性問題,這個問題在目前框架下沒有很好的辦法去自動化解決方案,那就是坐標(biāo)軸范圍,需要手動設(shè)置。前文有提及,這里再重復(fù)一次,ggcor() 初始化的默認坐標(biāo)范圍是和correlation matrix的行列數(shù)相關(guān)的,若行數(shù)為n,列數(shù)為m,x軸范圍是c(0.5, m + 0.5),y軸的范圍是c(0.5, n + 0.5)。數(shù)據(jù)預(yù)處理函數(shù)ggcor 提供了mantel檢驗的封裝函數(shù)fortify_mantel() ,支持vegan 包中的mantel() 、mantel.partial() 和ade4 包中的mantel.randtest() 、mantel.rtest() 函數(shù),差別上說mantel.partial() 是偏mantel檢驗(有控制變量),其它三個是mantel檢驗,當(dāng)不使用并行計算時,mantel.randtest() 速度最快(底層是C語言),mantel.rtest() 最慢,純粹R代碼實現(xiàn)。spec 是物種數(shù)據(jù),支持列表(list,非data.frame)或者數(shù)據(jù)框(data.frame)。- 若是數(shù)據(jù)框(最常見的情況),數(shù)據(jù)框中的每一列是一個物種(OTU),每行是一個樣本,可以通過
spec.select 參數(shù)來指定哪些列構(gòu)成一個群落。例如若spec 中1-4列是spec01群落,5-12是spec02群落,spec.select = list(spec01 = 1:4, spec02 = 5:12) (當(dāng)然,你也可以(最好不)不指定群落名稱,這是名字自動命名)。還有一種情況(設(shè)置spec.group 參數(shù)的情況)后面單獨解釋。
env 是環(huán)境數(shù)據(jù),支持列表(list,非data.frame)或者數(shù)據(jù)框(data.frame),env 中的每個元素對應(yīng)一個環(huán)境變量(當(dāng)然,若是列表,也可以支持多個環(huán)境變量組合成一個環(huán)境因素的情況)。還有一種情況(設(shè)置env.group 參數(shù)的情況)后面單獨解釋。env.ctrl 是環(huán)境控制變量,支持列表(list,非data.frame)或者數(shù)據(jù)框(data.frame)。需要注意,當(dāng)env.ctrl 非列表時,每次計算的控制環(huán)境是相同的,若需要分別設(shè)置不同的控制環(huán)境,需要通過列表手動設(shè)置。還有一種情況(設(shè)置env.ctrl.group 參數(shù)的情況)后面單獨解釋。mantel.fun mantel函數(shù),“mantel”、“mantel.randtest”、“mantel.rtest”和“mantel.partial”。spec.dist.fun 物種距離函數(shù),vegdist 或者dist 。env.dist.fun 環(huán)境距離函數(shù),vegdist 或者dist 。注意,當(dāng)前的情況下,控制環(huán)境使用環(huán)境距離函數(shù)。spec.dist.method 物種距離計算方法,參數(shù)默認是“bray”。env.dist.method 環(huán)境距離計算方法,參數(shù)默認是“euclidean”。... 其他傳遞給mantel.fun 的參數(shù)。 group相關(guān)的參數(shù)是為了處理需要根據(jù)樣本進行分組的情況,比如我A、B、C三個不同的樣本分組,物種、環(huán)境和控制環(huán)境(均必須為數(shù)據(jù)框)同樣如此,可以通過向量索引(和樣本量等長)來指定分組。這時,每個樣本組的物種只和對應(yīng)樣本組的環(huán)境列表的每個元素進行mantel測試。注意:以上情況均需要設(shè)置is.pair = TRUE 。## function(spec, ## env, ## env.ctrl = NULL, ## mantel.fun = 'mantel', ## is.pair = FALSE, ## spec.select = NULL, # a list of index vector ## spec.group = NULL, ## env.group = NULL, ## env.ctrl.group = NULL, ## spec.dist.fun = 'vegdist', ## env.dist.fun = 'vegdist', ## spec.dist.method = 'bray', ## env.dist.method = 'euclidean', ## ...)
案例一:按列設(shè)置物種群落library(vegan) # 使用vegan包所帶的數(shù)據(jù)集 data(varechem) data(varespec) fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20))
## # A tibble: 56 x 4 ## spec env r p ## <chr> <chr> <dbl> <dbl> ## 1 spec01 N 0.100 0.129 ## 2 spec02 N 0.254 0.011 ## 3 spec03 N 0.0596 0.238 ## 4 spec04 N 0.182 0.038 ## 5 spec01 P 0.165 0.028 ## 6 spec02 P 0.112 0.129 ## 7 spec03 P -0.0240 0.573 ## 8 spec04 P 0.201 0.02 ## 9 spec01 K 0.0481 0.295 ## 10 spec02 K 0.307 0.008 ## # … with 46 more rows
案例二:按樣本分組group <- rep(LETTERS[1:3], 8) fortify_mantel(varespec[ , 38:43], varechem, spec.group = group, env.group = group, is.pair = TRUE, mantel.fun = 'mantel.randtest')
## # A tibble: 42 x 4 ## spec env r p ## <chr> <chr> <dbl> <dbl> ## 1 A N -0.334 0.884 ## 2 B N 0.212 0.167 ## 3 C N 0.241 0.223 ## 4 A P -0.295 0.888 ## 5 B P 0.250 0.099 ## 6 C P 0.312 0.14 ## 7 A K 0.0909 0.31 ## 8 B K 0.310 0.101 ## 9 C K 0.0720 0.33 ## 10 A Ca -0.365 0.942 ## # … with 32 more rows
案例三:偏mantel測試fortify_mantel(varespec, varechem[ , 1:10], env.ctrl = varechem[ , 11:14], spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest')
## # A tibble: 40 x 4 ## spec env r p ## <chr> <chr> <dbl> <dbl> ## 1 spec01 N 0.100 0.134 ## 2 spec02 N 0.254 0.017 ## 3 spec03 N 0.0596 0.243 ## 4 spec04 N 0.182 0.039 ## 5 spec01 P 0.165 0.031 ## 6 spec02 P 0.112 0.122 ## 7 spec03 P -0.0240 0.582 ## 8 spec04 P 0.201 0.009 ## 9 spec01 K 0.0481 0.271 ## 10 spec02 K 0.307 0.007 ## # … with 30 more rows
mantel檢驗可視化mantel檢驗既然是相關(guān)性分析的一種特殊情況,那么相關(guān)性分析的可視化也同樣適用于mantel檢驗。當(dāng)然,由于fortify_mantel() 返回數(shù)據(jù)類和fortify_cor() 返回的數(shù)據(jù)類并不一樣,需要經(jīng)過一道轉(zhuǎn)換程序。mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20)) df06 <- as_cor_tbl(mantel) ggcor(df06) + geom_pie2() + geom_cross()
mantel檢驗組合圖
mantel組合圖是與相關(guān)性分析高度整合的,依賴于相關(guān)性分析函數(shù),換句話說mantel組合圖只是在相關(guān)性分析圖的基礎(chǔ)上額外疊加了一個圖層。核心函數(shù)是add_link() ,是個泛型函數(shù)(后面另說)。參數(shù)暫時不解釋了,看例子。實在抱歉,add_mantel() 又被我干掉了。corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest') ggcor(corr, xlim = c(-5, 14.5)) + add_link(mantel, diag.label = TRUE) + add_diaglab(angle = 45) + geom_square() + remove_axis('y')
corr <- fortify_cor(varechem, type = 'upper', show.diag = FALSE, cor.test = TRUE, cluster.type = 'all') ggcor(corr, xlim = c(-5, 14.5)) + add_link(mantel, diag.label = TRUE) + add_diaglab(angle = 45) + geom_pie2() + remove_axis('y')
corr <- fortify_cor(varechem, type = 'lower', show.diag = FALSE, cor.test = TRUE, cluster.type = 'all') ggcor(corr, xlim = c(0.5, 20)) + add_link(mantel, diag.label = TRUE) + add_diaglab(angle = 45) + geom_ellipse2() + remove_axis('y')
## 僅是測試效果,沒有實際意義 corr <- fortify_cor(varechem, varechem[ , 1:7], type = 'full', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest', nrepet = 2000) extra.params <- extra_params(group.label = text_params(size = 6), link.params = link_params(group.point.hjust = 2)) ggcor(corr, axis.y.position = 'left', legend.position = 'left', xlim = c(0.5, 14.5)) + add_link(mantel, extra.params = extra.params) + geom_circle2()
group <- rep(LETTERS[1:3], 8)
corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec[ , 38:43], varechem, spec.group = group, env.group = group, is.pair = TRUE, mantel.fun = 'mantel.randtest')
ggcor(corr, xlim = c(-5, 14.5)) + add_link(mantel, diag.label = TRUE) + add_diaglab(angle = 45) + geom_colour() + geom_shade(sign = -1) + remove_axis('y')
corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest') mantel <- dplyr::filter(mantel, p <= 0.05)
ggcor(corr, xlim = c(-5, 14.5)) + add_link(mantel, diag.label = TRUE, legend.drop = TRUE) + add_diaglab(angle = 45) + geom_square() + remove_axis('y')
## 個人覺得很丑
只要簡單的相關(guān)性?corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') corr01 <- fortify_cor(varechem, varespec[ , 38:39], type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') mantel <- fortify_mantel(varespec, varechem, spec.select = list(spec01 = 22:25, spec02 = 1:4, spec03 = 38:43, spec04 = 15:20), mantel.fun = 'mantel.randtest')
ggcor(corr, xlim = c(-5, 14.5)) + add_link(x = corr01, diag.label = TRUE, link.line.colours = c('#E31A1C', '#33A02C')) + add_diaglab(angle = 45) + geom_square() + remove_axis('y')
ggcor(corr, xlim = c(-5, 14.5)) + add_link(x = corr01, mapping = aes(size = abs(r)), diag.label = TRUE) + add_diaglab(angle = 45) + geom_square() + scale_size_continuous(limits = c(0, 1), range = c(0.25, 3)) + guides(size = guide_legend(title = 'abs r', override.aes = list(colour = 'grey35'), order = 1)) + remove_axis('y')
更靈活的方案若要自己更靈活的處理,或者是要用這個組合圖展示其它類型的數(shù)據(jù),add_link() 支持自定義的數(shù)據(jù)框做參數(shù)。第一個參數(shù)df 需要一個數(shù)據(jù)框,包含x和y列,x列類似于mantel檢驗中的物種群落(或者是樣本組),y類似于mantel檢驗中的環(huán)境變量。在內(nèi)部,會自動用df 中的y和相關(guān)系數(shù)矩陣的行名進行匹配確定坐標(biāo)位置。corr <- fortify_cor(varechem, type = 'upper', show.diag = TRUE, cor.test = TRUE, cluster.type = 'all') df <- data.frame(x = rep(LETTERS[1:3], 14), y = rep(cor_tbl_yname(corr), 3), r = runif(42, -1, 1), p = runif(42, 0, 0.5), stringsAsFactors = FALSE) ggcor(corr, xlim = c(-5, 14.5)) + add_link(df, diag.label = TRUE, colour = 'red') + add_diaglab(angle = 45) + geom_square() + remove_axis('y')
相關(guān)性網(wǎng)絡(luò)圖這塊內(nèi)容不會整合在ggcor 包里面,但是利用ggcor 里面的函數(shù)很容易導(dǎo)出相關(guān)性分析數(shù)據(jù)供其它函數(shù)使用。df07 <- fortify_cor(varespec, type = 'upper', show.diag = FALSE, cor.test = TRUE, keep.name = TRUE) df07
## # A tibble: 946 x 6 ## x y r p low upp ## <fct> <fct> <dbl> <dbl> <dbl> <dbl> ## 1 Empenigr Callvulg -0.262 0.217 -0.602 0.159 ## 2 Rhodtome Callvulg -0.0608 0.778 -0.453 0.351 ## 3 Rhodtome Empenigr 0.529 0.00781 0.160 0.769 ## 4 Vaccmyrt Callvulg -0.0925 0.667 -0.478 0.323 ## 5 Vaccmyrt Empenigr 0.205 0.337 -0.216 0.562 ## # … with 941 more rows ## Extra attributes: ## xname: Callvulg Empenigr Rhodtome Vaccmyrt Vaccviti Pinusylv ... ## yname: Cladphyl Claddefo Cladcerv Icmaeric Peltapht Stersp N ... ## show.diag: FALSE
## make graph library(ggraph) library(tidygraph) graph <- as_tbl_graph(df07) ggraph(graph, layout = 'linear', circular = TRUE) + geom_edge_arc(aes(colour = r, alpha = p)) + scale_edge_alpha_continuous(range = c(1, 0.1)) + coord_fixed() #顏色是相關(guān)性,線條濃淡是統(tǒng)計檢驗P值
寫在最后的話
后續(xù)除了優(yōu)化現(xiàn)有的功能,最重要的一個模塊是寫一個交互界面,能讓不懂R的人,或者不喜歡敲代碼的人也能輕松出圖。當(dāng)然,任何好用的包,都是無數(shù)的人共同協(xié)作、共同努力的結(jié)果。我希望能有更多的人參與測試,更多的人提供一些與相關(guān)性分析相關(guān)的新功能。ggcor 還有很多的小技巧沒有來得及介紹,全部留在下一期吧。高顏值免費在線繪圖R統(tǒng)計和作圖
|