【相関係数・順位相関係数】2変数の関係を把握しよう

こんにちは、たくろー(@takuro_data )です！

前回は、データのバラツキを表す分位数や分散について学習しました。

【分位数・分散・標準偏差・外れ値・変動係数】平均とあわせて確認しよう

前回までは、１変数の記述統計を扱ってきました。今回は２変数の相関についてです。２つの変数間に想定される直線的な関係を相関といいます。

要点

ピアソンの積率相関係数は、一般的に使われる相関係数
スピアマンの順位相関係数は、順位データに対して計算されたピアソンの積率相関係数
ケンドールの順位相関係数は、xについての順位とyについての順位が一致しているかどうかに着目して、相関の程度を測る

ピアソンの積率相関係数

一般的に使われる相関係数です。-1から1の値をとり、相関の程度を表します。変数$x$と変数$y$の相関係数は以下で計算します。$Cov$は共分散、$s$は標準偏差です。

$$r=\frac{Cov(x,y)}{s_{x}s_{y}}$$

前回同様、Rのサンプルデータwomen（30-39歳のアメリカ人の女性の平均の身長と体重データ）を使って算出します。

#データ読み込み
#単位変換：インチからセンチメートル、ポンドからキログラム
data("women")
women$height <- round(women$height / 0.39370,1)
women$weight <- round(women$weight / 2.2046,1)
#グラフで確認
plot(women$weight,women$height)

#ピアソンの積率相関係数
cor(women$height,women$weight, method="pearson")

#共分散から計算
corvar <- var(women$height, women$weight)
sd_h <- sd(women$height)
sd_w <- sd(women$weight)
corvar / (sd_h * sd_w)

> cor(women$height,women$weight, method="pearson")
[1] 0.9952961
> corvar / (sd_h * sd_w)
[1] 0.9952961

スピアマンの順位相関係数

順位のデータしか利用できない場合や、２つの変数間に曲線的な関係が想定される場合に、順位相関係数を用います。スピアマンの順位相関係数は、順位データに対して計算されたピアソンの積率相関係数です。

Rのサンプルデータcarsを使って算出します。carsは1920年代に記録された自動車の速度と停止距離のデータです。

#データ取得
data("cars")
plot(cars$speed,cars$dist)

#相関係数算出
cor(cars$speed,cars$dist, method="pearson")
cor(cars$speed,cars$dist, method="spearman")

#順位から計算
rank_s <- rank(cars$speed)
rank_d <- rank(cars$dist)
cor(rank_s, rank_d, method="pearson")

> cor(cars$speed,cars$dist, method="pearson")
[1] 0.8068949
> cor(cars$speed,cars$dist, method="spearman")
[1] 0.8303568
> cor(rank_s, rank_d, method="pearson")
[1] 0.8303568

ケンドールの順位相関係数

xについての順位とyについての順位が一致しているかどうかに着目して、相関の程度を測る指標です。

XとYのケンドールの順位相関係数を求めるためには、2対ずつ全ての組み合わせを取り出し比較します。X_i > X_jかつY_i > Y_jという同方向のケース数P、X_i > X_jかつY_i < Y_jという逆方向のケース数Q、X_i = X_jの数をT_x、Y_i = Y_jの数をT_yとして、以下の式で計算します。

$$\tau = \frac{P-Q}{\sqrt{_{n}C_{2}-T_{x}}\sqrt{_{n}C_{2}-T_{y}}}$$

Rにいいサンプルデータがなかったので、Jリーグの順位データを作成して算出してみます。2017年と2018年のリーグ順位を使い、順位の相関をみてみました（昇格・降格チームの順位は適当に埋めています）。年によって好調、不調の波があるチームがあったので、順位の相関は低かったです。

#データ作成
j2018 <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18)
j2017 <- c(1,15,2,11,7,13,3,14,10,9,12,5,16,8,17,6,4,18)
team_name <- c("Kawasaki","Hiroshima","Kashima","Sapporo","Urawa","FTokyo","COsaka","Shimizu","GOsaka","Kobe","Sendai","Yokohama","Shonan","Tosu","Nagoya","Iwata","Kashiwa","Nagasaki")
jrank <- cbind(j2017,j2018)
row.names(jrank) <- team_name
plot(j2017,j2018,type="n",xlim=c(18,1),ylim=c(18,1))
text(j2017,j2018,team_name)

#相関係数算出
cor(jrank[,1],jrank[,2], method="spearman")
cor(jrank[,1],jrank[,2], method="kendall")

#順位から計算
#関数定義
hikaku <- function(x1,x2,y1,y2){
  if(((x1<x2)&&(y1<y2))||((x1>x2)&&(y1>y2))){
    1
  } else if(((x1<x2)&&(y1>y2))||(x1>x2)&&(y1<y2)) {
    2
  } else{
    3
  }
}
kumiawase2 <- function(n){
  n*(n-1)/2
}
#順位行列を作成して計算
n <- length(j2018)
rankmat <- matrix(0,n-1,n-1)
for(i in 1:(n-1)){
  for(j in (i+1):n){
    rankmat[j-1,i] <- hikaku(j2017[i],j2017[j],j2018[i],j2018[j])
  }
}
dimnames(rankmat) <- list(team_name[2:n],team_name[1:(n-1)])
(sum(rankmat==1)-sum(rankmat==2))/kumiawase2(n)

> cor(jrank[,1],jrank[,2], method="spearman")
[1] 0.2837977
> cor(jrank[,1],jrank[,2], method="kendall")
[1] 0.1895425
> (sum(rankmat==1)-sum(rankmat==2))/kumiawase2(n)
[1] 0.1895425