皮爾森相關係數(Pearson correlation coefficient)

套路34: 皮爾森相關係數 (Pearson correlation coefficient)

1. 使用時機: 皮爾森相關係數是用以反映兩組變數之間關係密切程度的統計指標。
2. 分析類型: 母數分析(parametric analysis)。直接使用資料數值算統計叫parametric方法，把資料排序之後用排序的名次算統計叫non-parametric方法。
3. 皮爾森相關係數前提假設: 兩組變數之資料均為常態分布(normal distribution)或接近常態分布。
4. 範例資料1: 某研究餵食量對斑馬魚(zebrafish)存活的影響。在恆溫下，投入飼料X (mg)，斑馬魚存活比例Y的資料如下：
X (mg)        
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
0.55
Y (存活比例)
1
0.95
0.95
0.9
0.85
0.7
0.65
0.6
0.55
0.42
5. 畫圖看資料分佈:
第一步: 輸入建立資料，儲存在變數名稱為dat的data frame中。
  v1 <- c(0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55)
  v2 <- c(1, 0.95, 0.95, 0.9, 0.85, 0.7, 0.65, 0.6, 0.55, 0.42)
  dat <- data.frame(v1, v2)
第二步: 安裝ggplot2程式套件。
第三步: 呼叫ggplot2程式套件備用。
  library(ggplot2)
第四步: 使用函數ggplot畫x-y散佈圖及回歸線。
  ggplot(dat, aes(x = v1, y = v2)) +
    geom_point(shape = 1) +    # 畫空心圓
    geom_smooth(method = lm)   # 加回歸線
  # 灰色區域是回歸線的95%信賴區間。
6. 檢查資料是否符合常態分佈(normal distribution):
第一步: 輸入建立資料。
  v1 <- c(0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55)
  v2 <- c(1, 0.95, 0.95, 0.9, 0.85, 0.7, 0.65, 0.6, 0.55, 0.42)
第二步: 使用基本模組(base)函數shapiro.test代入v1及v2檢查資料是否符合常態分佈。
  shapiro.test(v1)
  shapiro.test(v2)
第三步: 判讀結果。
        Shapiro-Wilk normality test
data:  v1
W = 0.97016, p-value = 0.8924
        Shapiro-Wilk normality test
data:  v2
W = 0.92593, p-value = 0.4091
  # 如計算結果p-value < 0.05，H0:資料為常態分佈，不成立。
  # 如計算結果p-value > 0.05，H0:資料為常態分佈，成立。
7. 使用R計算皮爾森相關係數，方法一:
第一步: 輸入建立資料。
  v1 <- c(0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55)
  v2 <- c(1, 0.95, 0.95, 0.9, 0.85, 0.7, 0.65, 0.6, 0.55, 0.42)
第二步: 使用基本模組(base)函數cor代入v1及v2計算皮爾森相關係數。
  cor(v1, v2, method = "pearson")
[1] -0.9811093  # 計算結果皮爾森相關係數 = -0.9811093。
  # 無提供p值。
第三步: 檢定皮爾森相關係數是否為零 (X、Y無關): H0: ρ = 0，HA: ρ ≠ 0。
使用基本(base)模組函數cor.test代入v1及v2計算皮爾森相關係數並進行假設檢定。
  cor.test(v1, v2, alternative = "two.sided", method = "pearson")
第四步: 判讀結果。
         Pearson's product-moment correlation
data:  v1 and v2
t = -14.344, df = 8, p-value = 5.446e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.995675 -0.919468
sample estimates:  cor: -0.9811093   # 計算結果皮爾森相關係數 = -0.9811093。
  # 如計算結果p-value < 0.05，H0: ρ = 0，不成立，X、Y有關。
  # 如計算結果p-value > 0.05，H0: ρ = 0，成立，X、Y無關。
  # 如果要檢定: H0: ρ ≥ 0 & HA: ρ < 0或H0: ρ > 0 & HA: ρ ≤ 0，alternative = "less"。
  # 如果要檢定: H0: ρ ≤ 0 & HA: ρ > 0或H0: ρ < 0 & HA: ρ ≥ 0，alternative = "greater"。
8. 使用R計算皮爾森相關係數，方法二:
第一步: 輸入建立資料。
  v1 <- c(0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55)
  v2 <- c(1, 0.95, 0.95, 0.9, 0.85, 0.7, 0.65, 0.6, 0.55, 0.42)
第二步: 安裝fBasics程式套件。
第三步: 呼叫fBasics程式套件備用。
  library(fBasics)
第四步: 閱讀fBasics模組的pearsonTest函數的使用說明。
  help(pearsonTest)
第五步: 使用fBasics模組的pearsonTest函數代入v1及v2計算皮爾森相關係數。
  pearsonTest(v1, v2)
第六步: 判讀結果。
Title:
 Pearson's Correlation Test
Test Results:
  PARAMETER:
    Degrees of Freedom: 8
  SAMPLE ESTIMATES:
     Correlation: -0.9811
  STATISTIC:
     t: -14.3445
   P VALUE:
     Alternative   Two-Sided: 5.446e-07
     Alternative   Less: 2.723e-07
     Alternative   Greater: 1
   CONFIDENCE INTERVAL:
     Two-Sided: -0.9957, -0.9195
     Less: -1, -0.936
     Greater: -0.9945, 1
  # 如計算結果p < 0.05，H0: ρ = 0，不成立，X、Y有關。
  # 如計算結果p > 0.05，H0: ρ = 0，成立，X、Y無關。
  # 如果要檢定: H0: ρ ≥ 0 & HA: ρ < 0或H0: ρ > 0 & HA: ρ ≤ 0，alternative = "less"。
  # 如果要檢定: H0: ρ ≤ 0 & HA: ρ > 0或H0: ρ < 0 & HA: ρ ≥ 0，alternative = "greater"。
9. 檢定二組樣本數不同之皮爾森相關係數統計上是否有差: H0: ρ1 = ρ2，HA: ρ1 ≠ ρ2。
第一步: 輸入建立資料。
  v1 <- c(0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55)
  v2 <- c(1, 0.95, 0.95, 0.9, 0.85, 0.7, 0.65, 0.6, 0.55, 0.41)
  v3 <- c(0.1, 0.16, 0.21, 0.26, 0.31, 0.37, 0.41, 0.47, 0.51, 0.56, 0.6, 0.64, 0.69, 0.74)
  v4 <- c(0.99, 0.94, 0.93, 0.9, 0.86, 0.74, 0.66, 0.6, 0.55, 0.5, 0.45, 0.4, 0.35, 0.3)
第二步: 使用基本(base)模組函數cor，代入v1、v2及v3、v4算出兩組皮爾森相關係數r1、r2。
  r1 <- cor(v1, v2, method = "pearson")
  r2 <- cor(v3, v4, method = "pearson")
第三步: 安裝cocor程式套件。
第四步: 呼叫cocor程式套件備用。
  library(cocor)
第五步: 閱讀cocor程式套件的cocor.indep.groups函數使用說明。
  help(cocor.indep.groups)
第六步: 使用cocor程式套件的cocor.indep.groups函數，代入相關係數r1、r2及樣本數(n1 = 10, n2 = 14)，進行假設檢定。
  cocor.indep.groups(r1, r2, 10, 14, alternative = "two.sided", test = "all", alpha = 0.05, conf.level = 0.95,   
  null.value = 0, data.name = NULL, var.labels = NULL, return.htest = FALSE)
第七步: 判讀結果。
  Results of a comparison of two correlations based on independent groups
Comparison between r1.jk = -0.9801 and r2.hm = -0.9925
Difference: r1.jk - r2.hm = 0.0124
Group sizes: n1 = 10, n2 = 14
Null hypothesis: r1.jk is equal to r2.hm
Alternative hypothesis: r1.jk is not equal to r2.hm (two-sided)
Alpha: 0.05
fisher1925: Fisher's z (1925)
  z = 1.0161, p-value = 0.3096   # p-value > 0.05，H0: ρ1 = ρ2成立，二相關係數統計上沒差。
  Null hypothesis retained
zou2007: Zou's (2007) confidence interval
  95% confidence interval for r1.jk - r2.hm: -0.0103 0.0774
  Null hypothesis retained (Interval includes 0)
  # 如計算結果p-value < 0.05，H0: ρ1 = ρ2，不成立，二相關係數統計上有差。
  # 如計算結果p-value > 0.05，H0: ρ1 = ρ2，成立，二相關係數統計上沒差。
  # 如果要檢定: H0: ρ1 ≥ ρ2 & HA: ρ1 < ρ2或H0: ρ1 > ρ2 & HA: ρ1 ≤ ρ2，alternative = "less"。
  # 如果要檢定: H0: ρ1 ≤ ρ2 & HA: ρ1 > ρ2或H0: ρ1 < ρ2 & HA: ρ1 ≥ ρ2，alternative = "greater"。
10. 範例資料2: Bacteria_Genera.txt為空白分隔純文字資料檔，內含不同實驗樣本中細菌組成比例。
第一步: 輸入建立資料，使用read.table函數直接從檔案讀入資料或將資料貼到函數中如下所示，然後將資料放入變數dataset_g。
dataset_g <- read.table(header = T, text = "
  Acidipila Aciditerrimonas Acinetobacter Amaricoccus Aminicenantes Angustibacter
A1 0 0 0 0 0.046488287 0
A2 0.019110994 0 0 0.027027027 0.019110994 0
A3 0 0 0 0 0.016568013 0
B1 0.091030514 0.016349563 0 0 0 0.032699126
B2 0.09476955 0 0.032037955 0.027745683 0 0
B3 0.097683083 0 0.017834409 0 0 0.017834409
C1 0.090243422 0.018420861 0 0 0.026051032 0
C2 0.08234428 0 0.023290479 0.016468856 0 0.016468856
C3 0.06573422 0 0.017568209 0.017568209 0 0
D1 0.122666303 0.025039154 0 0 0 0
D2 0.099660276 0.022011273 0.026958193 0.022011273 0 0
D3 0.082080737 0.013874177 0 0 0 0
E1 0.05484085 0.019389168 0 0.027420425 0 0
E2 0.096534949 0.015457963 0.026773978 0.037864122 0 0
E3 0.072001213 0 0 0 0 0.033036422")
  # 資料為空白分隔的純文字。
  #  header = T資料有標題。
  # 資料第一列Acidipila前有一空白。
  # 變數名稱(dataset_g) 是使用者自定。
11. 相關係數的應用I: 相關係數矩陣(correlation matrix)，顯示多變數資料中任二變數間的相關性。
第一步:使用as.matrix函數將dataset_g資料轉換成matrix資料型態(data type)，放入變數m。
  m <- as.matrix(dataset_g)
  # 變數名稱(m)是使用者自定。
第二步: 安裝Hmisc程式套件。
第三步: 呼叫Hmisc程式套件備用。
  library(Hmisc)
第四步: 閱讀Hmisc程式套件的rcorr函數使用說明。
  help(rcorr)
第五步: 執行rcorr函數代入變數m資料，計算任二變數資料之相關係數矩陣及相對應之p值矩陣。
  rcorr(m, type = "pearson")
第六步: 判讀結果。
12. 相關係數的應用II: 相關係數矩陣(correlation matrix)，顯示多變數資料中任二變數間的相關性。
第一步:使用前述cor函數代入dataset_g計算任二變數資料相關係數放入變數M。
  M <- cor(dataset_g, method = "pearson")
# 變數名稱(M)是使用者自定。
第二步: 安裝corrplot程式套件。
第三步: 呼叫corrplot程式套件備用。
  library(corrplot)
第四步: 閱讀corrplot函數使用說明。
  help(corrplot)
第五步: 執行corrplot函數代入M計算並畫出相關矩陣。
  corrplot(M, type = "upper", method = "ellipse", order = "hclust", tl.col = "black")
第六步: 儲存圖檔。
  # 右上往左下斜(藍色): 正相關。左上往右下斜(紅色): 負相關。橢圓越細長相關係數越接近+1或-1。
13. 相關係數的應用III: 數值分佈矩陣(scatter plot matrix)及相關係數矩陣，顯示多變數資料中任二變數間的數值分佈及相關性。
第一步: 使用as.matrix函數將資料dataset_g轉換成matrix資料型態(data type)，放入變數m。
  m <- as.matrix(dataset_g)
  # 變數名稱(m)是使用者自定。
第二步: 安裝PerformanceAnalytics程式套件。
第三步: 呼叫PerformanceAnalytics程式套件備用。
  library(PerformanceAnalytics)
第四步: 閱讀PerformanceAnalytics程式套件chart.Correlation函數的使用說明。
  help(chart.Correlation)
第五步: 執行PerformanceAnalytics程式套件的chart.Correlation函數代入變數m資料，計算任二變數資料之數值分佈矩陣及相關係數矩陣。
  chart.Correlation(m, histogram = TRUE, method = "pearson", pch = 19)
第六步: 儲存圖檔。
  # 上圖中每個變量的分佈顯示在對角線上。
  # 在對角線的底部(下三角): 顯示具有迴歸線的雙變量(X-Y)散佈圖。
  # 在對角線的頂端(上三角): 相關係數加上作為恆星的顯著性水準。
  # 每個顯著性級別都與一個符號相關聯: p值(0: "***", 0.001: "**", 0.01: "*", 0.05: ".", 0.1: " ")。
14. 相關係數的應用IV: 使用相關係數矩陣作變數群聚分析(cluster analysis)及相關係數熱圖(heatmap)。
第一步: 使用基本模組(base)函數cor代入dataset_g計算相關係數矩陣，放入變數res。
  res <- cor(dataset_g, method = "pearson")
  # 變數名稱(res)是使用者自定。
第二步: 使用基本模組(base)函數colorRampPalette選擇顏色。
  col <- colorRampPalette(c("blue", "white", "red"))(20)
  # 變數名稱(col)是使用者自定。
第三步: 使用基本模組(base)函數heatmap代入res及col進行群聚分析並繪製熱圖。
  heatmap(x = res, col = col, symm = TRUE)
第四步: 儲存圖檔。
  # 上圖中紅色為正相關(相關係數正值)，藍色為反相關(相關係數負值)，白色為無關(相關係數近0)。

來勁了嗎? 想知道更多?? 補充資料(連結):
1. 關於Karl Pearson (https://en.wikipedia.org/wiki/Karl_Pearson)
2. 關於correlation coefficient (https://en.wikipedia.org/wiki/Correlation_coefficient)
3. 皮爾森相關係數計算公式 (https://en.wikipedia.org/wiki/Pearson_correlation_coefficient)
4. 關於R基礎，R繪圖及統計快速入門:
   a. R Tutorial: https://www.tutorialspoint.com/r/index.htm
   b. Cookbook for R: http://www.cookbook-r.com/
   c. Quick-R: https://www.statmethods.net/
   d. Statistical tools for high-throughput data analysis (STHDA): http://www.sthda.com/english/
e. The Handbook of Biological Statistics: http://www.biostathandbook.com/
f. An R Companion for the Handbook of Biological Statistics: http://rcompanion.org/rcompanion/index.html
5. 關於cocor Diedenhofen B, Musch J. 2015. cocor: a comprehensive solution for the statistical comparison of correlations. PLoS One. 10(3):e0121945. doi: 10.1371/journal.pone.0121945.
6. Zar, JH. 2010. Biostatistical Analysis, Fifth Edition, Pearson.

留言

張貼留言

這個網誌中的熱門文章

統計不球人 目錄 (Table of Contents)

用 R 做生醫資料科學   (Biomedical Data Science using R) 如何選擇統計方法 1 如何選擇統計方法  2 ***************************   第一段 : 起手式   ****************************** 套路 0: 安裝 R 不球人 套路 1: R 算統計不球人 ， 使用 R 函數 套路 2: 使用 ggplot2 繪圖 (Graphics using ggplot2) 套路 3: 敘述性統計 (Descriptive Statistics) 套路 4: 如何檢查資料是否接近常態分布 (Check Normality using R) 套路 5: 機率分布 (Probability Distributions) ***************************   第二段 : 比較平均值、中數   *************************** 套路 6: 單一樣本 Z 檢定 (One-Sample Z test) 套路 7: 單一樣本 t 檢定 (One-Sample t test) 套路 8: 單一樣本 Wilcoxon signed-rank 檢定 ( One-Sample Wilcoxon signed-rank test ， non-parametric ) 套路 9: 單一樣本符號檢定 (One Sample Sign test) 套路 10: 兩組獨立樣本 Z 檢定 (Two-Sample Z Test) 套路 11: 兩組變異數相同獨立樣本 t 檢定 ( t Test for Two-Samples with Equal Variances) 套路 12: 兩組變異數不同樣本 t 檢定 ( t Test for Two-Samples with Unequal Variances ， Welch Two Sample t test) 套路 13: 成對樣本 t 檢定 (Paired Sample t test) 套路 14: 兩組獨立樣本曼恩 - 惠尼 U 檢定 (...
閱讀完整內容

如何選擇統計方法 1

如何選擇統計方法 1 第一步 : 選擇統計方法可 依照統計方法用途的分類 來選擇統計方法 。 第二步 : 要了解確認統計方法的用途可看他的虛無假設 (Null hypothesis) H 0 。 以下為統計方法用途分類及本站網頁聯結 。 類別 1. 對資料進行基本分析 ， 以了解資料的大概情況 (summary) ， 例如平均值 、 中位數 、 四分位數 、 有無離群值 (outlier) 及是否為常態分布等 … 敘述性統計 (Descriptive Statistics) 如何檢查資料是否接近常態分布 (Normality Test using R) 類別 2. 比較平均值 (mean) 或中位數 (median) 2.1. 單一樣本 ( 與期望值或特定值相比 ) 單一樣本 Z 檢定 (One-Sample Z test) 單一樣本 t 檢定 (One-Sample t test) 單一樣本 Wilcoxon signed-rank 檢定 ( One-Sample Wilcoxon signed-rank test) 單一樣本符號檢定 (One Sample Sign test) 2.2. 兩組獨立樣本 ( 兩組獨立樣本相比 ) 兩組獨立樣本 Z 檢定 (Two-Sample Z Test) 兩組變異數相同樣本 t 檢定 (t Test for Two-Samples with Equal Variances) 兩組變異數不同樣本 t 檢定 (t Test for Two-Samples with Unequal Variances ， Welch Two Sample t-test) 兩組獨立樣本曼恩 - 惠尼 U 檢定 (Two-Sample Mann-Whitney U Test) 比較兩組獨立次序資料 (Two-Sample Test for Ordinal Data) 兩組獨立樣本 Mood's 中位數檢定 (Two-Sample Mood's Median Test) 2.3. 兩組成對樣本 ( 兩組成對樣本相比 ) 成對樣本 t 檢定 (Pai...
閱讀完整內容

單因子多樣本中位數差異檢定 (Kruskal-Wallis test)

圖片
套路 21: 單因子多樣本中位數差異檢定 (Kruskal-Wallis test) 變異數分析適用連續資料之差異分析，若自變項有一個，就是單因子變異數分析。 多樣本中位數差異檢定是對應 單因子變異數分析的無母數分析方法。當數據不符合單因子變異數分析的前提假設 ( 相同變異數，數據為常態分布 ) 時，可選擇 多樣本中位數差異檢定來做為替代分析方法 。 ** 注意 ** ” 比較變異數 ” (comparing variances) 和 ” 變異數分析 ” (analysis of variance) 不同。變異數分析是多組資料比較平均值或中位數。統計假設檢定檢定什麼 ? 看 H 0 。例如多組獨立樣本假設檢定 H 0 : k 組樣本中位數相同， H A : 至少有一組中位數不同 ， 是檢定多組資料的中位數是否相同。假設相等時為雙尾 (two-tailed test) 檢定。 1. 使用時機 : 用於比較多組樣本中位數 (median) 。 2. 分析類型 : 無母數分析 (non-parametric analysis) 。 直接使用資料數值算統計叫 parametric 方法，把資料排序之後用排序的名次算統計叫 non-parametric 方法。 3. 假設前提 : 無。 4. 範例資料 : 咪路調查餵食不同飼料的肉雞體重 (g) ， 資料如下： 飼料 1 飼料 2 飼料 3 飼料 4 61.8 78.8 70.5 60.3 65.1 79.5 72.6 63.8 61.7 76.0 71.7 64.1 63.3 73.4 72.0 61.4 77.3 71.1 60.9 不同飼料是否效果不同 ? H 0 : 餵食不同飼料的肉雞體重 ( 中位數 ) 相同。 H A : 餵食不同飼料的肉雞體重 ...
閱讀完整內容