小男孩‘自慰网亚洲一区二区,亚洲一级在线播放毛片,亚洲中文字幕av每天更新,黄aⅴ永久免费无码,91成人午夜在线精品,色网站免费在线观看,亚洲欧洲wwwww在线观看

分享

ggcor |相關(guān)系數(shù)矩陣可視化

 生物_醫(yī)藥_科研 2019-11-01
厚缊 業(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ù)擴展。

這里要介紹的ggcorcorrplot的有一種實現(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)然理論解讀中提前的spearmankendall方法也都支持。調(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 = lowupp = upp。
  • fill.*均是fill顏色映射函數(shù)相關(guān)的參數(shù)。
    • fill.scale.addFALSE不添加顏色映射函數(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ù)(xy、fillcolour、size等)最常用參數(shù)r、p、low、upp、num、r0、sig.thres、sig.level、mark等。
  • 可映射參數(shù)
    • 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、lowupp —— 適用于geom_confbox(),三個參數(shù)均是必須的,lowr、upp 分別確定置信區(qū)間盒子的下端、中間和上端線條位置。
    • num —— 適用于geom_num(),數(shù)值變量。這個函數(shù)封裝于geom_text(),做了一點點優(yōu)化,以后可能會刪除。
  • 非映射參數(shù)
    • 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.levelmark 適用于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月底,我在ggcorrrggcor前身)開發(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)。
    • 若是列表,列表中每個元素構(gòu)成一個群落;
    • 若是數(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)計和作圖

    本站是提供個人知識管理的網(wǎng)絡(luò)存儲空間,所有內(nèi)容均由用戶發(fā)布,不代表本站觀點。請注意甄別內(nèi)容中的聯(lián)系方式、誘導(dǎo)購買等信息,謹防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容,請點擊一鍵舉報。
    轉(zhuǎn)藏 分享 獻花(0

    0條評論

    發(fā)表

    請遵守用戶 評論公約

    類似文章 更多