相関の概念をさらにもう一歩深めるために。-ケンドールのタウを考えてみる-

昨日のエントリーで相関係数の有意性を確かめるためにどうすればいいのかということを書きました。
<a href="http://tomoshige-n.hatenablog.com/entry/2014/09/26/221108">相関係数の有意性を確かめる方法について -相関係数について1歩踏み込む- - Data Science by R and Python</a>
相関係数の有意性を確かめる方法について -相関係数について1歩踏み込む- - Data Science by R and Python

ですが、最後に書いたようにこの指標の有効性が検証できるのは「2つの変数が正規分布」に従うときに限られるということも最後に書きました。変数で正規分布に従わない変数はたくさんあります。身近な変数では「体重」は正規分布には従わないとされています。裾の重い分布、下の図のようになります。

f:id:tomoshige_n:20140927181754g:plain
URLBody Weight Distribution of Adult CT Scan Patients

この図、リンクを見ればわかる通りですが、右に引きづられています。この結果、分布は正規分布しないことがわかります。このようなデータ、変数同士の相関はどのようにして確認すればよいかという点を説明してみようと思います。

注意）以降の議論で使う「相関」というのは、これまでのような「直線的関係性」を意味していません。相関＝関連性を意味するものとします。

ケンドールのタウ（ $\tau$ ）

ケンドールの $\tau$ という指標（ケンドールの順位相関係数）を取り上げます。この指標は変数X, Yの単調性の指標です。単調性というのは「直線的」な増加や減少を指していません。「Xが増加すれば、Yは増加する」「Yが減少すれば、Xも減少する」などの傾向を指します。それが「直線的な増加」や「直線的な減少」ではなくてもよいのです。指数的、対数的な関係であっても、単調でありさえすれば良いという便利な指標です。

ただし、どのような関係で変動するのかはわかりませんし、さらに「増加、減少」などが周期的に起こるような場合にも検出ができません。この点に付いては注意する必要があります。

しかしながら、この指標の優れているところは、相関係数とは異なり、2つの変数X,Yに特別な分布（正規分布など）を仮定しないという点です。このような手法を「ノンパラメトリック」と呼びます。

ケンドールの $\tau$ の詳しい説明

まず、ケンドールの $\tau$ の定義を書くと、 $(X_1, Y_1), (X_2, Y_2)$ を同一の2変量分布に従う独立な組であるとします。このとき、ケンドールのタウは次のように定義されます。

$\tau = P[sgn{(X_1 - X_2)(Y_1 - Y_2)} = 1] - P[sgn{(X_1 - X_2)(Y_1 - Y_2)} = -1]$

さて、sgn{x}というのはxが正であれば1を返し、負であれば-1となるような関数です。この点を踏まえると、組の間に $sgn{(X_1 - X_2)(Y_1 - Y_2)} = 1$ が成立すれば、増加傾向があり、逆に-1が成立するときには、減少の傾向があるということがわかります。これが定義です。

しかしながら、母集団のことは私たちにはわからないので、次のようなケンドールのタウの不偏推定量を用いてこの値をサンプルから推定するということを行います。なので、利用する際にはこちらを使います。

$K = \frac{(符号が一致する数) - (符号が一致しない数)}{n(n-1)/2}$

この値は実際、期待値を取ると上で定義した $\tau$ になることが示されます（簡単なので、興味のある人はやってみてください）。

このケンドールのタウという指標は、相関係数とにたような性質を持っています。例えば、変数XとYが独立であるとき、ケンドールの $\tau$ = 0が成立することが示されます。もちろん、逆は成立しません。相関が0 -> XとYは独立は成立しません。

ただし、ケンドールのタウが1となるような場合は相関係数の場合と状況が異なります。先ほども説明した通り、ケンドールのタウは「単調性」に関する指標なので、XとY=X**2のタウの値は1になります。一方で、相関係数は1にはなりません。

また、データのノイズが大きいと、ケンドールのタウは検出力が大幅に低下します。この点は注意が必要になります。

実際に実行する方法

これを実際に使うためのコマンドが以下になります。データについては、ニューヨークの気象データがRに初期装備されていたので、これを使います。まずは、データを図示するところから始めます。

head(airquality)
    Ozone Solar.R Wind Temp Month Day
1      41     190  7.4   67     5   1
2      36     118  8.0   72     5   2
3      12     149 12.6   74     5   3
4      18     313 11.5   62     5   4
5      NA      NA 14.3   56     5   5

#必要な変数は最初の4つ
df = airquality[,1:4]

#散布図とヒストグラムを書く
pairs3 <- function (data) {
  oldpar <- par(no.readonly = TRUE) # 現在のグラフィックスパラメータ退避
  on.exit(par(oldpar)) # (関数がエラー中断しても)パラメータ復帰
  panel.hist <- function(x, ...) # ヒストグラムを描くための関数を定義
  { usr <- par("usr") # 現在のユーザー領域座標情報を得る
    on.exit(par(usr)) # 関数終了時に usr パラメータ復帰
    par(usr = c(usr[1:2], 0, 1.5) ) # ユーザー領域座標情報を変更
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
  }
  pairs(data, panel=panel.smooth,
        cex = 1.0, pch = 3, bg="light blue",
        diag.panel=panel.hist, cex.labels = 2, font.labels=2)
}

#散布図 + HISTOGRAM
pairs3(df)

f:id:tomoshige_n:20140927195254p:plain

この図からわかるように、Ozoneと、Solarは正規分布に近くなさそうです。この2つを前回のエントリーの関数でqqplot, histogram重ね合わせた視覚的確認、さらにコルモゴロフ・スミルノフ検定を行ってみます。すると、以下の図と結果が得られます（結果はOzoneのみ確認します）。

f:id:tomoshige_n:20140927195457p:plain

> ks.test(na.omit(df[,1]),pnorm,mean=mean(na.omit(df[,1])),sd=sd(na.omit(df[,1])))

	One-sample Kolmogorov-Smirnov test

data:  na.omit(df[, 1])
D = 0.148, p-value = 0.01243
alternative hypothesis: two-sided

 警告メッセージ: 
In ks.test(na.omit(df[, 1]), pnorm, mean = mean(na.omit(df[, 1])),  :
   コルモゴロフ・スミノフ検定において、タイは現れるべきではありません

ここで、警告メッセージが出されます。が、これは詳しくやり始めると大変なので無視しておきます。ks.testは値の中に「タイ（同じ数値）」があっては本当はいけないんです。

ただ、結果からもわかるように、5%有意で正規分布であるという帰無仮説を棄却することができ、かつ視覚的確認からも、qqnormの裾の方がかなり外れていることがわかります。よって、正規分布と仮定するのは適切ではないという結論を得ることができました。

変数が正規分布に従わないときの検定

このような場合、例えばOzoneとTempの関係性の数値を確認したいと思って「相関係数」を計算しても、この同時分布は前回も書きましたが「2変量正規分布」ではありませんので、検定をおこなうことはできません。なので、前回の手法で有意かどうかを確認することはできないということになります。

そこで、ケンドールのタウの登場です。ケンドールのタウを実際に計算すると以下のような値が得られます。

> cor(df$Ozone,df$Temp,use="pairwise",method="kendall")
[1] 0.5862988

そして、この統計量が有意であるかどうかを調べるためには、前回同様検定を行う必要があります。このとき「検定」のためには以下の統計量を使います（この統計量は、漸近的に標準正規分布に従います。またKは上で求めたケンドールのタウです）。

$Z_{k} = \frac{K}{\sqrt{2(2n+1)/9n(n-1)}} \sim Normal(0,1)$

この値を、実際に計算して検定するための関数が次の関数です。

> cor.test(df$Ozone,df$Temp,use="pairwise",method="kendall")

	Kendall's rank correlation tau

data:  df$Ozone and df$Temp
z = 9.1599, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.5862988

以上の結果より、「z = 9.1599, p-value < 2.2e-16」であるから、帰無仮説「ケンドールのタウは0」という仮説が成立する可能性は極めて低く棄却することができ、対立仮説「ケンドールのタウは0ではない」を採用することができます。よって、オゾン（Ozone）レベルと、気温（Temp）の間には「関連性」があることが認められました。

まとめ

さて、今回はケンドールのタウを説明してきました。おそらく、これまでの「相関係数」の考え方と随分離れており、取っ付きにくい方もいるかもしれません。ただ、正規分布しないデータに対して、正規分布だと決めつけて分析を行うのは、やはり危険ですし、少し難しくても別の手法を引き出してくる方が良いと思います。特に、ケンドールのタウはノンパラメトリックですから、ある意味どんな分布に従うデータに対しても使うことができます。

ただ、直感的に理解しづらい部分もありましたので、次回は「相関係数」とよく似ているけど、「ノンパラメトリック」な考え方である「スピアマンの順位相関係数」について書こうと思います。読んでくださってありがとうございました。