【Rによる統計学 No.3】変数の相関

こんにちは、タクロウ(@takuro_109)です!

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

【Rによる統計学 No.2】データのバラツキ

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

要点
  1. ピアソンの積率相関係数は、一般的に使われる相関係数
  2. スピアマンの順位相関係数は、順位データに対して計算されたピアソンの積率相関係数
  3. ケンドールの順位相関係数は、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

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

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

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対ずつ全ての組み合わせを取り出し比較します。Xi > XjかつYi > Yjという同方向のケース数P、Xi > XjかつYi < Yjという逆方向のケース数Q、Xi = Xjの数をTx、Yi = Yjの数をTyとして、以下の式で計算します。

$$\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

さいごに

今回は、変数の相関について学習しました。あまり順位相関係数を使ってこなかったのですが、散布図をみて曲線的な関係がありそうな場合は、順位相関係数を使っていこうと思いました。普段は関数コマンドに頼っているので、計算式から改めてプログラミングしてみると理解が深まりました。次回は確率分布についてです。

以上、お読みいただきありがとうございました。