データ

# u, v, wのカラムデータを持つデータフレームを作成
d <- data.frame(
  u = c(7, 4, 6, 2, 9, 3, 8, 1, 6, 3),
  v = c(3, 6, 5, 8, 3, 4, 2, 9, 2, 6),
  w = c(6, 3, 4, 9, 1, 8, 4, 2, 5, 7))

library(DT)
datatable(d)

相関分析

相関係数

x <- c(1, 2, 3, 5)
y <- c(0, 3, 3, 6)

# cor関数で2変量の相関係数を計算する。
cor(x, y)
## [1] 0.9561829

変数u vs. 変数v

cor(d$u, d$v)
## [1] -0.851136
cor(d$u, d$v, method = 'spearman') # スピアマン順位相関
## [1] -0.836927
cor(d$u, d$v, method = 'kendall')  # ケンドール順位相関
## [1] -0.7059312

相関行列

r <- cor(d)
r
##            u          v          w
## u  1.0000000 -0.8511360 -0.4151543
## v -0.8511360  1.0000000  0.1015166
## w -0.4151543  0.1015166  1.0000000
library(kableExtra)
kable(round(r, 2), caption = '相関表') |> kable_classic('striped', full_width = F)
相関表
u v w
u 1.00 -0.85 -0.42
v -0.85 1.00 0.10
w -0.42 0.10 1.00

相関分析図

psych パッケージを利用したグラフ

library(psych)
pairs.panels(d)

cor.plot(d)

corrplot パッケージを利用したグラフ

library(corrplot)
corrplot.mixed(r, lower = 'ellips', upper = 'number')

無相関検定

検定統計量\(T\)は\(H_0\) のもと自由度 \(n-2\)の\(t\)分布に従う
\[T=\frac{|r|\sqrt{n-2}}{\sqrt{1-r^2}}\sim t(n-2)\] ここで，\(r\)：標本相関係数，\(n\)：サンプルサイズ

無相関検定 U vs. V

(uv <- cor.test(d$u, d$v)) # 無相関検定
## 
##  Pearson's product-moment correlation
## 
## data:  d$u and d$v
## t = -4.586, df = 8, p-value = 0.001788
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9641022 -0.4772863
## sample estimates:
##       cor 
## -0.851136
# p値による有意性判定
ifelse(uv$p.value < 0.05, '有意', '有意でない')
## [1] "有意"

無相関検定 U vs. W

(uw <- cor.test(d$u, d$w)) # 無相関検定
## 
##  Pearson's product-moment correlation
## 
## data:  d$u and d$w
## t = -1.2907, df = 8, p-value = 0.2329
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8282759  0.2903733
## sample estimates:
##        cor 
## -0.4151543
# p値による有意性判定
ifelse(uw$p.value < 0.05, '有意', '有意でない')
## [1] "有意でない"

無相関検定 V vs. W

(vw <- cor.test(d$v, d$w)) # 無相関検定
## 
##  Pearson's product-moment correlation
## 
## data:  d$v and d$w
## t = 0.28862, df = 8, p-value = 0.7802
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5641701  0.6872176
## sample estimates:
##       cor 
## 0.1015166
# p値による有意性判定
ifelse(vw$p.value < 0.05, '有意', '有意でない')
## [1] "有意でない"

Hmiscパッケージ利用

相関行列に加え無相関検定の\(p\)値も同時に計算されるで便利。

library(Hmisc)

rcorr(as.matrix(d))
##       u     v     w
## u  1.00 -0.85 -0.42
## v -0.85  1.00  0.10
## w -0.42  0.10  1.00
## 
## n= 10 
## 
## 
## P
##   u      v      w     
## u        0.0018 0.2329
## v 0.0018        0.7802
## w 0.2329 0.7802

インタラクティブグラフ

無相関検定で有意でないものは「☓」マークが付くので分かりやすい。

#install.packages("ggcorrplot")
library(ggcorrplot)
library(plotly)
cor(d) |> ggcorrplot(lab = T, hc.order = T, outline.color = "white", p.mat = cor_pmat(d)) |> ggplotly() |>
layout(font  = list(size = 11, color = 'blue', family = 'UD Digi Kyokasho NK-R'),
       title = '主タイトル',
       xaxis = list(title = 'x軸カテゴリラベル'),
       yaxis = list(title = 'y軸カテゴリラベル'))

R 演習課題（相関分析「擬似相関・無相関検定」の回の課題）

次のデータを用いて相関分析を行え。

※相関分析１の課題は講義スライドに書いてあるので注意。

データ

相関係数\(\rho = 0.2\)，サンプルサイズ\(n = 100\)の２次元正規乱数

library(MASS)

n   <- 100 # サンプルサイズ
rho <- 0.2 # 母相関係数（ρ）

set.seed(5)

# 多変量正規乱数を生成する関数（MASS::mvrnorm）
m <- mvrnorm(n = n, mu = c(10, 5),
             Sigma = rbind(c(1, rho),
                           c(rho, 1)))

d <- as.data.frame(m) # データフレーム型に変換

colnames(d) <- c('u', 'v')

library(DT)
datatable(round(d, 2))

Python

データ

import pandas as pd

d0 = {'u': [7, 4, 6, 2, 9, 3, 8, 1, 6, 3],
      'v': [3, 6, 5, 8, 3, 4, 2, 9, 2, 6],
      'w': [6, 3, 4, 9, 1, 8, 4, 2, 5, 7]}
      
d = pd.DataFrame(d0)

相関分析

相関係数

変数u vs. 変数v

import numpy as np

r = np.corrcoef(d['u'], d['v'])
r[0, 1]
## -0.8511360272327114

相関行列

cor = d.corr()
cor
##           u         v         w
## u  1.000000 -0.851136 -0.415154
## v -0.851136  1.000000  0.101517
## w -0.415154  0.101517  1.000000
d.corr(method = 'spearman') # スピアマン順位相関
##           u         v         w
## u  1.000000 -0.836927 -0.382265
## v -0.836927  1.000000  0.070553
## w -0.382265  0.070553  1.000000
d.corr(method = 'kendall')  # ケンドール順位相関
##           u         v         w
## u  1.000000 -0.705931 -0.321860
## v -0.705931  1.000000  0.069786
## w -0.321860  0.069786  1.000000

相関図

import seaborn as sb

sb.heatmap(cor, cmap = 'RdBu_r', annot = True)

無相関検定

from scipy.stats import pearsonr
rho, p = pearsonr(d['u'], d['v'])
print("相関係数：", rho.round(2))
## 相関係数： -0.85
print("p値：", p)
## p値： 0.001788003147831374