從相關(guān)矩陣到網(wǎng)絡(luò)圖一次把R的相關(guān)分析與可視化說透先說個我遇見的事兒。某次項目匯報對方興沖沖拉出一張布滿星星的相關(guān)矩陣熱圖指著上面一片“***”說“你看這些變量全都顯著相關(guān)項目有戲?！蔽颐榱艘谎凵Ⅻc圖心里直嘆氣那分明是一個變量被極端值拉出來的“虛假直線”幾個所謂顯著相關(guān)換了穩(wěn)健方法立馬變不顯著。這不是個例——太多人把“統(tǒng)計分析”做成了“跑個函數(shù)出個圖”卻忽略了相關(guān)分析真正見功夫的地方數(shù)據(jù)形態(tài)、系數(shù)選型、缺失值處理、可視化后的正確解讀。這篇文章就圍繞“用R做相關(guān)分析及其制圖可視化”展開我會先從相關(guān)分析的底層邏輯講起幫你搞明白r值和p值到底在說什么再比較Pearson、Spearman、Kendall三種系數(shù)怎么選接著用R的mtcars數(shù)據(jù)集完整走一遍數(shù)據(jù)準備、相關(guān)計算、顯著性檢驗的實操流程然后重點展示從散點圖、相關(guān)矩陣圖到熱圖和網(wǎng)絡(luò)圖的四級可視化最后聊幾個實戰(zhàn)中繞不開的坑包括多重比較校正、異常值影響、置信區(qū)間解讀這些不踩一遍不會重視的問題。無論你是剛上手R的初學(xué)者還是從SPSS遷移過來的熟練工這篇都能幫你把相關(guān)分析這步做得更扎實。1. 想清楚相關(guān)分析在解決什么問題再打開RStudio1.1 散點圖先行先看形態(tài)再算系數(shù)我在帶新人做數(shù)據(jù)分析時定的第一條規(guī)矩是任何相關(guān)分析之前必須先把兩兩變量的散點圖畫出來。原因很簡單——相關(guān)系數(shù)是一組數(shù)據(jù)經(jīng)過復(fù)雜計算后壓縮出來的單個數(shù)字它一旦被算出來原始數(shù)據(jù)的形態(tài)就看不到了。這就是著名的Anscombe四重奏告訴我們的道理四組數(shù)據(jù)計算出來的相關(guān)系數(shù)幾乎完全一樣約0.816但把它們畫成散點圖一組是漂亮的線性關(guān)系一組是明顯的曲線關(guān)系一組是被一個離群點統(tǒng)治的“偽直線”一組則是一堆散亂點外加一個杠桿點支撐起整個趨勢。如果你只報告相關(guān)系數(shù)這四組數(shù)據(jù)在你眼里就“一模一樣”但這顯然是荒謬的。所以實操習慣應(yīng)該是這樣先pairs()掃一眼再用cor()算數(shù)。R里面畫基礎(chǔ)散點圖矩陣一句話就夠pairs(mtcars[, c(mpg, disp, hp, wt)], pch 19, col rgb(0, 0, 0, 0.6))看到對角線之外的小圖你會很快建立對數(shù)據(jù)整體面貌的判斷哪些變量看起來線性、哪些明顯彎、哪些全是離群點。這個步驟花不了10秒但能幫你免掉后面一堆“假顯著”的麻煩。1.2 r值的幾何直覺標準化之后就是夾角余弦很多人背過公式但沒建立直覺。相關(guān)系數(shù)的數(shù)學(xué)定義是協(xié)方差除以兩個變量標準差的乘積$$r \frac{Cov(x,y)}{s_x s_y}$$概念上可以理解為先把x和y都做標準化減去均值除以標準差然后求這兩個標準化向量的內(nèi)積再除以n-1。學(xué)過線性代數(shù)的朋友會發(fā)現(xiàn)這本質(zhì)上就是兩個向量夾角的余弦值——所以r永遠落在[-1, 1]區(qū)間內(nèi)。當兩個變量的變化方向完全一致時標準化后的向量重合夾角為零r等于1完全相反時r等于-1毫無關(guān)系時兩個向量正交r接近0。這個幾何直覺非常有用。它提醒你r衡量的是“線性關(guān)系”的強度和方向不是任意的關(guān)聯(lián)。兩個變量即使存在完美的二次曲線關(guān)系r也可能算出來是0。這也是我在第一節(jié)強調(diào)散點圖先行的原因——只有看到真實的數(shù)據(jù)形態(tài)你才知道用線性相關(guān)系數(shù)是不是合適的選擇。1.3 p值不是效應(yīng)量解讀相關(guān)分析的第二支柱算出r之后R會順手給你一個p值。它的邏輯是這樣的在原假設(shè)“兩個變量的總體相關(guān)系數(shù)等于0”的前提下根據(jù)樣本量和r的大小計算出一個t統(tǒng)計量然后看這個t值有多“極端”。p值小說明當前樣本數(shù)據(jù)與“總體無關(guān)”的假設(shè)不太兼容于是你很自然地認為變量之間確實存在某種相關(guān)。但是這個定義有兩處必須強調(diào)。第一p值受樣本量影響極大樣本量5000時r0.03也會有漂亮的p0.001樣本量15時r0.6可能也達不到顯著。第二p值衡量的是“是否顯著”r值衡量的是“關(guān)系多強”這兩個維度缺一不可。實際寫報告時我的習慣是總報告置信區(qū)間cor.test(mtcars$mpg, mtcars$wt, method pearson)cor.test輸出的置信區(qū)間比單獨的p值信息量大得多。比如mpg和wt的相關(guān)系數(shù)它可能是-0.8795%置信區(qū)間為[-0.94, -0.74]這意味著在大量重復(fù)抽樣中總體相關(guān)的可能范圍都不會越過0結(jié)論相當穩(wěn)固。如果置信區(qū)間跨過0就算點估計是-0.4你也應(yīng)該把這個結(jié)果視作“方向未定”。2. 別上來就cor()三種相關(guān)系數(shù)的選型邏輯用R做相關(guān)分析時最常見的懶人寫法是這樣cor(data)默認算的是Pearson積差相關(guān)。很多教程也默認這個選擇但套用一句話說默認參數(shù)只代表“軟件默認”不代表“統(tǒng)計正確”。實際工作中你得先回答一個問題我的數(shù)據(jù)配得上Pearson嗎2.1 Pearson連續(xù)、線性、干凈的數(shù)據(jù)才配得上它Pearson相關(guān)系數(shù)計算時直接使用原始數(shù)值所以它對三個前提敏感變量是等距或比率尺度連續(xù)型變量之間的關(guān)系是線性的數(shù)據(jù)中沒有極端異常值。醫(yī)療研究里的身高體重、生物學(xué)里的基因表達量適度清洗后、工業(yè)場景里的溫度和壓力這些連續(xù)且大致呈線性關(guān)系的數(shù)據(jù)用Pearson都很自然。它的優(yōu)點也很明顯統(tǒng)計檢驗的效能最高在同樣樣本量下最容易檢出真實存在的相關(guān)。而且在變量近似正態(tài)分布時它是最佳線性無偏估計意義上的最優(yōu)選擇。但如果數(shù)據(jù)里有一個明顯的離群點Pearson可能被單個點拉扯得面目全非這一點我放在第5章專門講。2.2 Spearman把數(shù)據(jù)變成排名再說話Spearman的做法是先把每個變量的原始值轉(zhuǎn)換成秩也就是排名然后對秩計算Pearson相關(guān)。這樣做的好處是它只關(guān)心單調(diào)關(guān)系不要求線性對異常值的敏感度大大降低因為極端值再多也只占一個排名位置。四六級成績與績點的關(guān)系、用戶評分與購買頻次的關(guān)系、各種李克特量表數(shù)據(jù)這類“序數(shù)型”或“分布偏得厲害”的數(shù)據(jù)Spearman都是更穩(wěn)妥的選擇。代碼上和Pearson的區(qū)別只是一個參數(shù)cor.test(mtcars$mpg, mtcars$hp, method spearman)實測下來你會發(fā)現(xiàn)當數(shù)據(jù)比較“乖”時Spearman和Pearson的結(jié)果通常非常接近一旦數(shù)據(jù)出現(xiàn)非線性單調(diào)關(guān)系或異常值Spearman往往更能反映真實趨勢。2.3 Kendall小樣本和同秩數(shù)據(jù)多的穩(wěn)健派Kendall的τ系數(shù)思路更樸素對兩個變量的每一對觀測值判斷這對值在兩個變量上的變化方向是否一致然后統(tǒng)計一致性比例。它不依賴于數(shù)值本身只依賴于“對”的方向所以對異常值幾乎“免疫”在小樣本或同秩tie很多的數(shù)據(jù)里表現(xiàn)更穩(wěn)。它的數(shù)值通常比Pearson和Spearman略小但這不是錯誤而是定義方式不同。比如你要分析12個樣本之間兩種測量方法的一致性或者調(diào)查問卷得來的幾十個有序分類數(shù)據(jù)我建議優(yōu)先考慮Kendall。它的計算速度在大樣本下會慢一些但樣本量通常撐不到那個級別。2.4 三種系數(shù)怎么選一張表說清楚系數(shù)類型數(shù)據(jù)尺度關(guān)系形態(tài)要求對異常值敏感度典型場景Pearson積差相關(guān)連續(xù)型等距/等比線性關(guān)系高生理指標、物理測量、表達量數(shù)據(jù)Spearman秩相關(guān)連續(xù)型或有序型單調(diào)關(guān)系即可低問卷量表、偏態(tài)數(shù)據(jù)、序數(shù)數(shù)據(jù)Kendall秩相關(guān)任意有序數(shù)據(jù)單調(diào)關(guān)系即可極低小樣本、同秩多的數(shù)據(jù)、一致性評價我在實際項目中采用的原則很樸素畫完散點圖覺得形態(tài)還線性且沒有明顯異常值用Pearson只要形態(tài)彎曲或分布偏態(tài)立刻切Spearman或Kendall不用糾結(jié)太多。因為分析服務(wù)的是判斷不是形式。3. 從原始數(shù)據(jù)到相關(guān)矩陣R語言里的標準動作3.1 選數(shù)據(jù)別把不該算的變量一把梭塞進矩陣分析的第一步不是寫cor()而是想清楚“哪些變量應(yīng)該被放進相關(guān)分析”。我見過太多人把編號列、ID列、單位不同的一堆指標全扔進cor()里結(jié)果輸出一片亂碼。正確做法是先確定分析目標再構(gòu)造一個只含分析所需變量的數(shù)據(jù)框。用內(nèi)置數(shù)據(jù)集mtcars演示假設(shè)我們要研究汽車油耗mpg和哪些車輛屬性關(guān)聯(lián)最緊密df - mtcars[, c(mpg, disp, hp, wt, qsec)]這幾列分別是油耗、排量、馬力、車重、四分之一英里加速用時。檢查一下數(shù)據(jù)結(jié)構(gòu)確認沒有亂七八糟的因子或字符串混在里面str(df) summary(df)str確認類型summary幫你快速掃一遍取值范圍為后面的異常值風險做心理準備。3.2 缺失值怎么處理complete.obs、pairwise.complete.obs和rcorr的默認行為真實數(shù)據(jù)幾乎沒有完整的所以這個坑必繞不開。R的cor()函數(shù)里控制缺失值的參數(shù)是use常見的取值有三種everything默認只要數(shù)據(jù)里有NA結(jié)果直接輸出NA最老實但也最沒用。complete.obs只保留所有變量都非缺失的整行樣本即所謂的列刪法listwise deletion。數(shù)據(jù)缺失較多時樣本量會被砍掉一大截造成浪費。pairwise.complete.obs逐對計算時只使用這對變量都非缺失的樣本即所謂的成對刪法pairwise deletion。每個相關(guān)系數(shù)可能基于不同數(shù)量的樣本但能最大限度利用現(xiàn)有數(shù)據(jù)。我看到不少分析報告直接默認cor(x)結(jié)果輸出一片NA然后人傻在那里。如果你確認缺失值比例可以接受我通常建議用pairwise.complete.obs做探索性分析如果樣本量足夠且缺失模式是隨機的用complete.obs做正式報告也可以。還有第三種更穩(wěn)妥的姿態(tài)用Hmisc包中的rcorr()它內(nèi)部也是成對刪除而且能同時返回相關(guān)系數(shù)矩陣和p值矩陣library(Hmisc) res - rcorr(as.matrix(df), type pearson) res$r # 相關(guān)系數(shù)矩陣 res$P # 顯著性p值矩陣 res$n # 各對樣本的有效樣本量rcorr需要矩陣輸入所以先as.matrix。輸出后你會同時拿到r矩陣和P矩陣這對后續(xù)篩選顯著相關(guān)非常方便。3.3 單個關(guān)系細看cor.test是必做動作相關(guān)矩陣適合全局掃視但如果你對某對變量特別感興趣比如“油耗和車重到底關(guān)系多強”那必須用cor.test細看cor.test(df$mpg, df$wt, method pearson)它會給出r值、t統(tǒng)計量、自由度、p值以及95%置信區(qū)間。實操時我不會只看p值而是重點看置信區(qū)間上下界。比如區(qū)間是[-0.94, -0.74]這個關(guān)系相當穩(wěn)定如果區(qū)間是[-0.85, 0.12]那你得承認雖然點估計看著有點嚇人但證據(jù)其實撐不住“確定存在相關(guān)”的結(jié)論。很多人報告相關(guān)矩陣時只寫r和p這不算錯但丟了區(qū)間就是丟了信息量。3.4 批量篩選從幾十個變量里快速找到重點關(guān)系當變量多了以后手工一對對看cor.test不現(xiàn)實我習慣在rcorr的結(jié)果上直接做篩選r_mat - res$r p_mat - res$P sig_mask - p_mat 0.05然后用sig_mask去過濾r_mat只保留顯著的相關(guān)系數(shù)再結(jié)合5.1的多重比較校正做進一步篩查。這個流程跑下來一次性從20個變量里挑出5個核心關(guān)系非常高效。4. 制圖可視化的三層遞進散點圖、相關(guān)矩陣圖、網(wǎng)絡(luò)圖4.1 基礎(chǔ)圖形不能省pairs和ggplot2的散點圖矩陣相關(guān)分析的可視化我從不用一張熱圖打發(fā)全部需求。第一件該做的事是散點圖矩陣它能讓你看到所有兩兩關(guān)系的形態(tài)?；A(chǔ)pairs()畫得快但不好看要想拿得出手我推薦用ggplot2手動構(gòu)造一個精修版library(ggplot2) library(GGally) df_plot - df ggpairs(df_plot, upper list(continuous wrap(cor, size 4)), lower list(continuous smooth))GGally包的ggpairs會在對角線放密度曲線或直方圖上三角放相關(guān)系數(shù)下三角放散點圖和平滑擬合線。血管里的東西一目了然mpg與wt的散點圖是一條清晰的負斜率直線而mpg與qsec的散點圖有點彎這時就該警覺Pearson的r值可能會低估真實關(guān)系的強度。4.2 corrplot一張相關(guān)矩陣圖把格局拉開當變量數(shù)量在5到15個之間時corrplot包畫的圓形相關(guān)矩陣圖非常經(jīng)典。它用圓圈大小表示|r|的大小用顏色表示正負信息密度和美觀度都不錯library(corrplot) cor_matrix - cor(df, use complete.obs, method pearson) corrplot(cor_matrix, method color, type upper, order hclust, addCoef.col black, tl.col black, tl.srt 45, number.cex 0.8)這里的幾個參數(shù)值得解釋method color用顏色方塊表示數(shù)值既有直觀視覺沖擊又能看清大小。type upper只保留上三角因為相關(guān)矩陣是對稱的全畫出來徒增視覺負擔。當然如果就是要展示對稱之美或者強調(diào)某個變量的關(guān)系用full也無妨。order hclust按層次聚類結(jié)果的順序重排變量讓相關(guān)性相近的變量靠在一起矩陣看起來有聚類的結(jié)構(gòu)。這在變量多時尤其重要否則只是一團混亂的顏色。圖上每一個格子里面我還加了相關(guān)系數(shù)數(shù)值addCoef.col這樣看圖的人不用再猜顏色塊到底代表0.6還是0.7。你可以把這個圖直接導(dǎo)出成PDF或PNG放進報告外觀完全夠格。4.3 熱圖加聚類ggcorrplot更現(xiàn)代的做法如果團隊里的成員更習慣看類似Python seaborn風格的熱圖R這邊我用ggcorrplot替代library(ggcorrplot) p - ggcorrplot(cor_matrix, hc.order TRUE, type lower, lab TRUE, lab_size 3, colors c(#6D9EC1, white, #E46726), outline.color gray60) print(p)hc.order自動做聚類順序排列typelower取左下三角。這樣一來變量會自動“物以類聚”相關(guān)性強的一組在視覺上擠成一個深色塊相關(guān)性弱的分隔出白色區(qū)域。這張圖我做業(yè)務(wù)匯報時幾乎總會放一張因為聽眾不用思考“格子坐標是什么”只看顏色聚集就能讀懂“誰和誰有關(guān)系”。4.4 再進階一步用igraph把相關(guān)矩陣畫成網(wǎng)絡(luò)圖當變量數(shù)量突破20個相關(guān)矩陣再精巧也是密密麻麻一片這時候網(wǎng)絡(luò)圖反而是更好的表達工具。網(wǎng)絡(luò)圖的思路很直接每個變量是節(jié)點相關(guān)系數(shù)絕對值超過某個閾值的變量之間連一條邊邊顏色區(qū)分正負相關(guān)節(jié)點大小反映它在整個網(wǎng)絡(luò)中的連接強度。library(igraph) cor_matrix_threshold - cor_matrix cor_matrix_threshold[abs(cor_matrix_threshold) 0.6] - 0 net - graph_from_adjacency_matrix(cor_matrix_threshold, mode undirected, weighted TRUE, diag FALSE) E(net)$color - ifelse(E(net)$weight 0, #E46726, #6D9EC1) E(net)$width - abs(E(net)$weight) * 3 V(net)$size - degree(net) * 8 set.seed(123) plot(net, layout layout_with_fr(net))需要注意閾值不是亂定的。0.6是我在探索階段常用的經(jīng)驗值但你應(yīng)該根據(jù)樣本量和業(yè)務(wù)背景調(diào)整。閾值太高網(wǎng)絡(luò)里節(jié)點會孤立得七零八落閾值太低全網(wǎng)連成一團等于沒畫。我通常會在0.4到0.7之間掃一遍挑一個既能看清主要群落、節(jié)點又不至于全部“單飛”的閾值。網(wǎng)絡(luò)圖的最大價值不是替代相關(guān)矩陣而是幫你找到“變量群落”。比如你可能會發(fā)現(xiàn)某幾個業(yè)務(wù)指標形成了一個密集連接的小團體另幾個指標孤立在角落——這背后的業(yè)務(wù)含義往往比單個r值有價值得多。5. 實戰(zhàn)里的坑與解多重比較、異常值、置信區(qū)間一個都不能少5.1 多重比較一張矩陣幾十個“顯著”真沒這么簡單假設(shè)你有10個變量pairwise相關(guān)對數(shù)就有45對。即便總體中完全不存在任何相關(guān)在0.05的顯著性水平下平均也會有約2.25對顯著。變量到了20個190對里憑空冒出來的“顯著結(jié)果”接近10個。這就是做相關(guān)矩陣最容易踩的大坑把每一對的p值當成“獨立檢驗”來看。要解決它思路和方差分析后的兩兩比較一致——對p值做多重比較校正p_vec - p_mat[upper.tri(p_mat)] p_adj - p.adjust(p_vec, method BH)BHBenjamini-Hochberg方法控制錯誤發(fā)現(xiàn)率在探索性分析里是常用且不太“暴力”的選擇。如果你想更保守可以用Bonferronimethodbonferroni它會直接拿顯著性水平除以檢驗次數(shù)代價是可能會漏掉很多真實存在的弱相關(guān)。我個人的選擇是探索階段用BH正式報告前用更嚴格的方法復(fù)核一遍。不過有一點也得說清楚嚴格校正后沒通過不代表關(guān)系一定不存在更可能是功效不足你可以把它標記為“待驗證”而不是“無相關(guān)”。5.2 異常值一個點毀掉一條線異常值對Pearson相關(guān)的影響是災(zāi)難性的。比如10個點原本呈現(xiàn)明顯的負相關(guān)r-0.3只要在右上方塞進一個“離群點”r可能瞬間變成正0.4。這個數(shù)字反轉(zhuǎn)的魔術(shù)我在很多項目里見過。應(yīng)對方法有三層。第一層是“可視化發(fā)現(xiàn)”散點圖矩陣上掃一圈看看有沒有孤獨的點飄在主體數(shù)據(jù)之外。第二層是“穩(wěn)健替換”用Winsorize方法把極端值壓縮到合理分位數(shù)然后用壓縮后的數(shù)據(jù)再做相關(guān)或者直接用Spearman相關(guān)系數(shù)秩相關(guān)天然降權(quán)極端值。第三層是“分而治之”把數(shù)據(jù)按某個關(guān)鍵維度分組比如不同地區(qū)分開做相關(guān)分析看每組內(nèi)的相關(guān)方向和強度是否一致。不分組只看整體有時候會得到“平均效應(yīng)”掩蓋真實關(guān)系的結(jié)論。5.3 置信區(qū)間點估計只說了一半的故事有人在報告里寫“r0.8p0.001”聽起來關(guān)系強又顯著但實際上可能樣本量只有10置信區(qū)間是[-0.1, 0.96]從強正相關(guān)到零相關(guān)都在合理范圍內(nèi)?？尚艈嵴f實話證據(jù)相當弱。所以在正式報告里相關(guān)系數(shù)至少要附帶置信區(qū)間或者明確報告樣本量。R里cor.test給別人看它天然輸出置信區(qū)間批量計算的話可以用corrplot包里的cor.mtest或者自寫bootstrap函數(shù)library(boot) boot_cor - function(data, indices) { d - data[indices, ] return(cor(d$mpg, d$wt, method pearson)) } set.seed(123) b - boot(df, boot_cor, R 1000) boot.ci(b, type perc)這個bootstrap置信區(qū)間不依賴正態(tài)假設(shè)對異常值也更寬容。我通常在變量關(guān)系不太干凈時用它替代cor.test自帶的近似區(qū)間。5.4 報告呈現(xiàn)數(shù)字背后要讓人看得懂最后聊一下報告習慣。我見過一個反面的極端有人把四五十個變量的相關(guān)矩陣整頁貼在PPT上完全沒有任何標注讀者根本不知道重點在哪。建議的做法是只放重點變量或重點關(guān)系最多不超過10個核心結(jié)果。r值保留兩位小數(shù)p值寫成具體數(shù)字p0.012或者用“*p0.05, **p0.01”標注同時在校注里說明沒有做多重比較校正或采用哪種校正法。把結(jié)論寫成一句話比如“油耗與車重呈顯著負相關(guān)r-0.8795%CI[-0.94,-0.74]說明車重越大油耗越高”。這句話比十個星星都有力。圖和表搭配使用散點圖矩陣證明數(shù)據(jù)形態(tài)corrplot展示全局格局最后報告里只點名核心結(jié)論。我每次跑相關(guān)分析無論多急都要給自己留三樣?xùn)|西一張散點圖矩陣確認線性形態(tài)、一張帶顯著性標注的相關(guān)矩陣圖觀察全局格局、一張經(jīng)過多重比較校正后的篩選結(jié)果表。先看圖再算數(shù)后下結(jié)論。這不是流程儀式感而是實打?qū)崕臀以陧椖坷锉苊膺^好幾次“拿假相關(guān)講故事”的尷尬。如果你之前是直接cor()然后截個圖放進報告下次不妨試試這個思路——你的結(jié)論會穩(wěn)很多報告也不那么容易被同行問倒。