R言語で相関係数を出力する例を紹介します。
相関係数の出力だけでなく、相関係数と散布図の関係についてもみていきます。
実際には2つのデータ間に明らかな関係があるのに相関係数が0に近く相関がほとんどないと解釈できてしまう例なども触れていきます。
ピアソンの積率相関係数については以下の記事を参照。
-
【統計学】相関係数・ピアソンの積率相関係数
相関係数とピアソンの積率相関係数について解説する。 相関係数およびピアソンの積率相関係数の定義を与え、その性質や幾何学的解釈などをみていく。 相関係数の検定については以下の記事を参照されたい。 相関係 ...
続きを見る
また、相関係数の検定については以下を参照。
-
【統計学】相関係数の検定・無相関性の検定
相関係数の検定を解説する。 相関係数の検定の検定統計量や棄却域の導出について解説する。 相関係数の分布の導出を行い、検定統計量をどのように構成すればよいかみていく。 相関係数については以下の記事を参照 ...
続きを見る
-
【R言語】相関係数の検定・無相関性の検定 関数cor.testの使い方
R言語で相関係数の検定を行う方法を紹介します。 この記事では、t分布に基づく無相関性の検定を実行する関数やその使い方について見ていきます。 フィッシャーのz変換を用いた相関係数の検定、差の検定について ...
続きを見る
標本相関係数（ピアソンの積率相関係数）
よくデータ解析で耳にする「相関係数」とはピアソンの積率相関係数のことを指します。相関係数の検定などではまた意味合いが変わりますが、2つのデータ間の相関を計算したい場合は次の標本相関係数（ピアソンの積率相関係数）を用います。\begin{align} r_{x, y} = \cfrac{\sum_{i = 1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i = 1}^n(x_i - \bar{x})^2 \sum_{i = 1}^n(y_i - \bar{y})^2}}.\end{align}\((x_1, y_1), (x_2, y_2), \ldots, (x_n , y_n)\)は2変数の標本です。
上の標本相関係数が2つのデータ間の相関の尺度として用いられる理由として、次のように標本相関係数\(r_{x, y}\)は\(x\)と\(y\)の線形相関の強さを表すからです。
- r_{x, y}が正であるときは、2つのデータ\(x\)と\(y\)の間に正の相関があり、\(x\)が大きいとき\(y\)も大きくなる。
- 逆に\(r_{x, y}\)が負であるときは、2つのデータ\(x\)と\(y\)の間に負の相関があり、\(x\)が大きくなると\(y\)は小さくなる。
このような理由から、相関の尺度としてよく用いられます。
標本相関係数の計算例
Rで標本相関係数を算出する例を以下にまとめました。標本相関係数がどのようなものなのか分かりやすいように、2変量正規分布（山みたいな分布）の確率密度のプロットと、そこから発生させた乱数の標本相関係数と散布図を載せています。
また、実際の解析にも応用できるように、Rの標準データセットを用いた標本相関係数の計算例やプロットも解説しています。
2変量正規分布からの乱数
2変量正規分布から乱数を発生させて、その標本相関係数および散布図をプロットするには次のように実行します。
パッケージmvtnormの関数dmvnormを用いることで多変量正規分布の確率密度が得られます。
次を実行することで次のパラメータの多変量正規分布の確率密度を描くことが可能です。\begin{align}\boldsymbol{\mu} &= \begin{pmatrix}0 \\ 0\end{pmatrix},\\ \boldsymbol{\Sigma} &= \begin{pmatrix}1 & \rho \\ \rho & 1\end{pmatrix}.\end{align}関数outerでxとy軸を-3から3まで0.1での動かしたときの確率密度dmvnormを計算し確率密度を描いています。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | library(mvtnorm) n <- 20 mu <- c(0, 0) rhos <- c(0.01, 0.5, 0.8, 0.99) mkMatrix <- function(rho) { return(matrix(c(1, rho, rho, 1), ncol = 2)) } #確率密度 par(mfrow = c(2, 2)) par(mar = c(4, 3, 1, 1)) par(oma = c(0, 0, 0, 0)) par(mgp = c(2, 1, 0)) for (rho in rhos) { Sigma <- mkMatrix(rho) x <- seq(-3, 3, 0.1) y <- x f <- function(u, v) { c(dmvnorm(matrix(c(u, v), ncol = 2), mu, Sigma)) } density <- outer(x, y, f) persp(x, y, density, theta = 0, phi = 60, expand = 0.5, xlim = c(-3, 3), ylim = c(-3, 3), sub = paste0("ρ=", rho)) } |
確率密度
左上から\(\rho = 0.1, 0.5, 0.8, 0.99\)の多変量正規分布の確率密度関数がプロットされます。相関係数\(\rho\)が小さくなるにつれて円上に分布し、高くなるにつれて\(y = x\)の直線の周りに分布することが分かります。
次に標本相関係数と散布図の描き方についてみていきます。
標本相関係数は次の関数を用いることで計算することができます。
| x | numeric型のベクトルまたは行列、データフレーム。 |
| y | xと同じ次元のベクトルまたは行列、データフレーム。デフォルトでy=x。 |
| use | character型。欠測値がある場合にどのように共分散を計算するか決める。"everything", "all.obs", "complete.obs", "na.or.complete", "pairwise.complete.obs"。 |
| method | character型。どの標本相関係数を計算するか決定する。"pearson"でpearson、"kendall"でKendall、"spearman"でSpearmanの標本相関係数を計算できる。 |
次を実行することで、下の画像のような\(\rho= 0.1, 0.5, 0.8, 0.99\)の多変量正規分布からの乱数をプロットすることができます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | #散布図 par(mfrow = c(2, 2)) par(mar = c(4, 3, 1, 1)) par(oma = c(0, 0, 0, 0)) par(mgp = c(2, 1, 0)) for (rho in rhos) { Sigma <- mkMatrix(rho) samples <- rmvnorm(n, mu, Sigma) x <- samples[, 1] y <- samples[, 2] r <- cor(x, y) plot(x, y, sub = paste0("r=", round(r, 5))) } |
散布図1
相関係数\(\rho\)が高くなるにつれて\(y = x\)の直線の周りに標本が集中することが確認できます。また、標本相関係数\(r\)も相関係数\(\rho\)に近い値をとっており、ちゃんと推定できていることがわかります。
特殊な例
次に、特殊な例をいくつか紹介します。
次を実行すると、下の画像のような標本相関係数\(r\)が\(1\)または\(-1\)となる散布図と、2変数間に関連性があるのにもかかわらず標本相関係数が\(0\)に近くなる散布図を描くことができます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | #r=1, r=-1の例 par(mfrow = c(2, 2)) par(mar = c(4, 3, 1, 1)) par(oma = c(0, 0, 0, 0)) par(mgp = c(2, 1, 0)) signs <- c(1, -1) for (sign in signs) { x <- runif(n, 0, 2) y <- sign * x r <- cor(x, y) plot(x, y, sub = paste0("r=", round(r, 5))) } #r=0でも2変量間の関係がある例 x <- runif(n, -0.5, 0.5) radius <- 0.5 y <- sapply(x, function(z) { (sqrt(- 4 * (z^2 - radius^2))) / 2 }) r <- cor(x, y) plot(x, y, sub = paste0("r=", round(r, 5))) |
散布図2
上の2つの散布図から分かるように、標本相関係数\(r\)が\(1\)または\(-1\)になるのは、\(x_i = y_i,\ i = 1, 2, \ldots, n\)を満たす場合であることが分かります。
一番下の散布図は中心\((0, 0)\)、半径\(0.5\)の半円上に一様乱数を発生させた例です。2変数は円上に分布するという関係があるのにもかかわらず、標本相関係数は\(r = 0.09676\)となっており、2変数間の相関は小さいという解釈ができてしまいます。
このように、ピアソンの積率相関係数は2変数間の線形相関の尺度であるため、2変数間に比例関係以外の関係性がある場合注意が必要です。
Rの標準データセットを用いた例
最後に、Rの標準データセットを用いた相関係数の計算例と散布図のプロット例を紹介します。
irisを用いた例
データセットirisの散布図と標本相関係数の算出の例をまずみていきます。
関数plotを用いることで、簡単に各変数間の12通りの散布図を描くことが可能です。
1 2 3 4 | library(ggplot2) library(GGally) dataset <- iris plot(dataset[, -seq_len(ncol(dataset))[colnames(dataset) == "Species"]]) |
iris散布図
パッケージGGallyを用いると、次のような散布図を作成することも可能です。
1 2 | ggpairs(dataset, columns = seq_len(ncol(dataset) - 1), aes(color = Species, alpha = 0.5), lower = list(continuous = "smooth")) |
iris散布図GGally
関数corを用いて各変数間の12通りの標本相関係数の計算を行うことができますが、少々面倒です。関数corではなくcov2corを用いると便利です。この関数は標本共分散行列を相関行列に変換してくれます（標本共分散行列の対角成分が1になるように変換してくれます）。以下、標本相関係数の計算例です。
1 2 3 4 5 6 7 8 9 10 | > cor(dataset$Sepal.Length, dataset$Sepal.Width) [1] -0.1175698 > cor(dataset$Petal.Length, dataset$Petal.Width) [1] 0.9628654 > cov2cor(var(dataset[, -seq_len(ncol(dataset))[colnames(dataset) == "Species"]])) Sepal.Length Sepal.Width Petal.Length Petal.Width Sepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411 Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259 Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654 Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000 |
anscombeを用いた例
次にanscombeを用いた例を紹介します。
このデータセットは、母集団分布が異なる場合でも標本相関係数すなわち回帰曲線が一致してしまう例をあげています。
単に相関係数を計算して2変数間の関連性について言及するのではなく、母集団分布がどうなっているのかまで確認する必要があることを説明するために作られたデータセットとなります。
1 2 3 4 5 6 7 8 9 | > dataset <- anscombe > head(dataset) x1 x2 x3 x4 y1 y2 y3 y4 1 10 10 10 8 8.04 9.14 7.46 6.58 2 8 8 8 8 6.95 8.14 6.77 5.76 3 13 13 13 8 7.58 8.74 12.74 7.71 4 9 9 9 8 8.81 8.77 7.11 8.84 5 11 11 11 8 8.33 9.26 7.81 8.47 6 14 14 14 8 9.96 8.10 8.84 7.04 |
さっそく、anscombeの標本相関係数の計算をしてみましょう。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | > cor(dataset$x1, dataset$y1) [1] 0.8164205 > cor(dataset$x2, dataset$y1) [1] 0.8164205 > cor(dataset$x3, dataset$y1) [1] 0.8164205 > cor(dataset$x4, dataset$y4) [1] 0.8165214 > cov2cor(var(dataset)) x1 x2 x3 x4 y1 y2 y3 y4 x1 1.0000000 1.0000000 1.0000000 -0.5000000 0.8164205 0.8162365 0.8162867 -0.3140467 x2 1.0000000 1.0000000 1.0000000 -0.5000000 0.8164205 0.8162365 0.8162867 -0.3140467 x3 1.0000000 1.0000000 1.0000000 -0.5000000 0.8164205 0.8162365 0.8162867 -0.3140467 x4 -0.5000000 -0.5000000 -0.5000000 1.0000000 -0.5290927 -0.7184365 -0.3446610 0.8165214 y1 0.8164205 0.8164205 0.8164205 -0.5290927 1.0000000 0.7500054 0.4687167 -0.4891162 y2 0.8162365 0.8162365 0.8162365 -0.7184365 0.7500054 1.0000000 0.5879193 -0.4780949 y3 0.8162867 0.8162867 0.8162867 -0.3446610 0.4687167 0.5879193 1.0000000 -0.1554718 y4 -0.3140467 -0.3140467 -0.3140467 0.8165214 -0.4891162 -0.4780949 -0.1554718 1.0000000 |
いくつかの標本相関係数が一致しているのが分かります。標本相関係数一致した変数の散布図を次のように並べてみます。
anscombe散布図
上記のように2変数間に比例関係が無い場合でも標本相関係数が同じになってしまう例が存在します。全く異なる母集団分布をもつ変数間の標本相関係数を比較する際には注意する必要があることが上のグラフから分かります。比例関係でないのに高い相関を持つという間違った結論を出してしまうことがあります。
まとめ
R言語で相関係数を計算する例を紹介しました。
corを用いることで簡単に2つのデータ間の相関係数（ピアソンの積率相関係数）を算出することが可能です。
また、plotを用いることで散布図を作成することができ、2つ以上の変数の場合でも簡単に複数の散布図を作成することができます。
