R の相関: cor()、cor.test()、相関行列

2 変数間の相関: cor()

相関は次のことを問います。一方の変数が上がるとき、もう一方も上がる傾向があるか（正）、下がるか（負）、それとも独自に動くか（ゼロ付近）。R では 1 回の呼び出しで済みます。

print(cor(mtcars$wt, mtcars$mpg))

答えは約 −0.87 です。重い車ほど燃費が悪く、その関係は強いということです。この 1 つの数値が Pearson の相関係数で、一般に r と書かれます。

r の読み方

係数は必ず −1 と +1 の間に収まります。符号 が方向を、大きさ が強さを表します。

| |r| | 一般的な解釈 | | --- | --- | | 0.0 – 0.2 | ほぼ無視できる | | 0.2 – 0.4 | 弱い | | 0.4 – 0.6 | 中程度 | | 0.6 – 0.8 | 強い | | 0.8 – 1.0 | 非常に強い |

これらの区分は法律ではなく会話の糸口として扱ってください。物理学では r = 0.6 は期待外れですが、心理学ではキャリアの目玉になります。身につけておくべき性質が 2 つあります。r には単位がない（重量をトンで測っても kg で測っても mpg との r は同じです。r は標準化された値で計算されるからです）こと、そして r は 線形 の関連しか測らないこと。完全な U 字型の関係でも r ≈ 0 になり得ます。

そして言っておかねばならない一文があります。相関は因果ではありません。 アイスクリームの売上と溺死者数は年間の各月にわたって強く相関しますが、それはアイスクリームが人を溺れさせるからではなく、夏が両方を押し上げるからです。相関は 2 つの変数が連動していることを教えますが、なぜ かについては沈黙します。x が y を動かしているのかもしれず、y が x を動かしているのかもしれず、第 3 の何か（この場合は季節）が両方を動かしているのかもしれません。そのどれかを決めるには、より大きな r ではなく、実験や慎重な因果推論が必要です。

Spearman と Kendall: Pearson が適切でないとき

Pearson は生の値で計算するため、外れ値に敏感で、曲がった関係には目が届きません。method 引数で順位ベースの代替手法に切り替えられます。

x <- c(1, 2, 3, 4, 5, 6, 7, 8)
y <- c(1, 4, 9, 16, 25, 36, 49, 300)  # increasing, with a wild last value

print(cor(x, y))                        # Pearson - dragged by the outlier
print(cor(x, y, method = "spearman"))   # ranks: perfectly monotonic
print(cor(x, y, method = "kendall"))

Spearman はすべての値をその順位に置き換え、順位に対して Pearson を計算します。y は x が増えるたびに必ず増えるので、順位はぴったり揃い、Spearman はちょうど 1 を報告します。外れ値の 大きさ は問題でなくなり、その 位置 だけが効きます。データが順序尺度（アンケートの段階評価）のとき、大きく歪んでいるとき、あるいは関係が単調だが直線的でないときは Spearman を選びましょう。Kendall は同順・逆順の組から似た問いに答えます。小標本ではより頑健ですが計算が遅く、一般的な既定は Spearman です。

相関行列

多数の変数にわたって関係を見渡すには、cor() に複数の数値列を渡します。

m <- cor(mtcars[, c("mpg", "wt", "hp", "disp")])
print(round(m, 2))

すべての変数を他のすべての変数と突き合わせた結果が、対角に 1（すべては自分自身と完全に相関します）を、その対角線を挟んで鏡像を伴って表示されます。小数第 2 位に丸めることは、聞こえる以上に重要です。丸めない行列は数字の壁であり、行列の目的は 見渡す ことだからです。ここで見渡すと、mpg は 3 つすべてと負に相関し（重く、強力で、排気量の大きい車ほど燃料を食う）、wt、hp、disp は互いに強く正に相関していることが分かります。「大きな車」変数のクラスタであり、線形回帰と多重共線性の頭痛に進むときに効いてきます。

欠測値: use 引数

欠測データがあるとき、cor() の既定は推測ではなく NA を返すことです。

d <- mtcars[, c("mpg", "wt", "hp")]
d$wt[c(3, 7)] <- NA  # knock two values out

print(cor(d)[1, 2])
print(cor(d, use = "complete.obs")[1, 2])
print(cor(d, use = "pairwise.complete.obs")[1, 2])

• use = "complete.obs" は NA を ひとつでも 含む行をすべて削除し、生き残った行から行列全体を計算します。一貫性はありますが無駄が多く、wt の欠測が mpg–hp の組からもその行を取り除いてしまいます。
• use = "pairwise.complete.obs" は、その組 が揃っているすべての行から各セルを計算します。より多くのデータを残せますが、セルごとに異なる部分集合に基づくため、ごくまれに内部的に整合しない行列が生じることがあります。

NA が数個散らばっている程度ならどちらでも構いません。ただし、どちらを使ったかは明記しましょう。

有意なのか? cor.test()

cor() は数値をくれますが、それがノイズでありうるかについては何も言いません。cor.test() が仮説検定を実行します。

print(cor.test(mtcars$wt, mtcars$mpg))

出力をブロックごとに見ていきましょう。

• t = −9.56, df = 30 — 検定統計量です。帰無仮説は真の相関がゼロであるというもので、観測された r は自由度 n − 2（32 台 − 2）の t 統計量に変換されます。
• p-value = 1.29e-10 — 真の相関がゼロだったとして、32 個の標本でこれほどゼロから離れた r が見られる確率は約 0.0000000001 です。関連が本物であることの圧倒的な証拠です。ただし、p 値は r がゼロと異なるかどうかを語るのであって、関係が 大きい か 因果的 かを語るものではないことを忘れないでください。
• 95 percent confidence interval: −0.93 to −0.74 — 真の相関としてありうる範囲です。p 値より有用なことがよくあります。この区間の楽観的な端でさえ強い負の相関です。
• sample estimates: cor = −0.87 — cor() が返したのと同じ数値です。

小標本にはここで特別な敬意が必要です。n = 10 では ±0.5 の相関が偶然に驚くほど頻繁に現れ、幅の広い信頼区間がそのことを教えてくれます。p 値だけでなく区間を報告しましょう。

目で見る: 必ず図を描く

相関係数は関係全体を 1 つの数値に圧縮するもので、その圧縮は曲がりやクラスタ、あるいは 1 点がすべてを担っている状況を隠しかねません。どんな r も信じる前に散布図を見てください。

plot(mtcars$wt, mtcars$mpg)              # one pair
pairs(mtcars[, c("mpg", "wt", "hp", "disp")])  # every pair in the matrix

pairs() は相関行列に対応する散布図のグリッドを描きます。数値が思ったとおりの意味かを確かめる最速の方法です。洗練されたヒートマップ形式の行列グラフィックスには、corrplot パッケージ（install.packages("corrplot") の後に corrplot(cor(m))）が標準的な道具です。

この章のまとめ

• cor(x, y) は Pearson の r を返します。符号が方向、大きさが強さで、常に [−1, 1] の範囲、単位はありません。
• 相関は 線形 の共変動を測るもので、因果については何も語りません。潜んでいる第 3 の変数は常に候補です。
• 順位を使う method = "spearman" は、順序データ、外れ値、単調だが曲がった関係に向きます。
• 数値列に対する cor(df) が行列を返します。round(, 2) し、欠測があるときは use = 引数に注意しましょう。
• cor.test(x, y) は p 値と信頼区間を加えます。区間を報告しましょう。
• 数値を信じる前に、必ず散布図を見ましょう。

次は、問いが「連動しているか?」から「このグループの平均はあちらと違うか?」へと鋭くなるとき — t 検定です。