【R言語】相関係数の検定・無相関性の検定　関数cor.testの使い方

R言語で相関係数の検定を行う方法を紹介します。

この記事では、t分布に基づく無相関性の検定を実行する関数やその使い方について見ていきます。

フィッシャーのz変換を用いた相関係数の検定、差の検定については次の記事を参照していください。

: 【R言語】フィッシャーのZ変換を用いた相関係数・相関係数の差の検定
R言語でフィッシャーのZ変換を用いた相関係数の検定を行う方法を紹介します。この記事では、正規近似であるフィッシャーのZ変換による相関係数の検定と相関係数の差の検定を解説していきます。相関係数の検定 ...
続きを見る

無相関性の検定だけでなく、フィッシャーのz変換を用いた相関係数の検定も解説します。

この記事で紹介するプログラミングコードは以下からダウンロードできます。

R言語　相関係数の検定

1 ファイル 1.42 KB

ダウンロード

相関係数の検定の詳細については次の記事を参照してください。

: 【統計学】相関係数の検定・無相関性の検定
相関係数の検定を解説する。相関係数の検定の検定統計量や棄却域の導出について解説する。相関係数の分布の導出を行い、検定統計量をどのように構成すればよいかみていく。相関係数については以下の記事を参照 ...
続きを見る

R言語ではじめるプログラミングとデータ分析

Rでデータ解析を始めるならコレ

¥3,300 （2022/06/01 17:20時点 | Amazon調べ）

ポチップ

統計学初学者向けの教材です

¥2,970 （2022/06/01 17:41時点 | Amazon調べ）

Amazon

楽天市場

Yahooショッピング

ポチップ

相関係数の検定

相関係数の検定の関数をいくつか紹介します。

関数cor.test

関数cor.testによって、相関係数の検定（無相関性の検定）を実行することができます。

cor.test(x, y, alternative = c("two.sided", "less", "greater"), method = c("pearson", "kendall", "spearman"), exact = NULL, conf.level = 0.95, continuity = FALSE, ...)

cor.testの引数は以下の通りです。

関数cor.testの引数
x, y	numeric型の観測値のベクトル。xとyの長さは同じでなくてはならない。
alternative	character型。"two.sided"で両側仮説、"greater"で右片側仮説。"less"で左片側仮説を指定する。
method	character型。"pearson"でPearsonの積率相関係数、"kendall"でKendallのτ、"spearman"でSpearmanのρに関する検定を指定できる。
exact	logical型。KendallとSpearmanの検定の際に正確なp値を計算するかどうか。
conf.level	信頼区間の信頼水準。ピアソンの積率相関係数の検定かつサンプルサイズが4以上のときに適用される。
continuity	logical型。trueの場合、KendallとSpearmanの検定の際に連続性補正が適用される。
formula	~ u + vの形から成るformula型のオブジェクト。
data	観測値から成る行列またはデータフレーム。
subset	検定に用いる観測値を指定するためのベクトル。
na.action	データにNAがある場合に実行する関数。デフォルトでna.action。

関数cor.testは以下の無相関性の検定を実行します。任意の相関係数に関する検定ではないことに注意が必要です。この次に紹介する関数を用いることで、任意の相関係数についての検定を実行することが可能です。

相関係数の検定（ピアソンの積率相関係数）

相関係数$\rho$の2変量正規分布から$(x_1, y_1), (x_2, y_2)\ldots, (x_n, y_n)$の2組から成る標本が与えられているとし、ピアソンの標本相関係数を$r$とする。また、対応する確率変数をそれぞれ$(X_1, Y_1), (X_2, Y_2)\ldots, (X_n, Y_n)$、$R$とする。このとき次の仮説検定を考える。\begin{align}&H_0:\ \rho = 0\\ &H_1:\ \rho \neq 0\end{align}この仮説検定の検定統計量は次で与えられる。\begin{align} R \sqrt{\cfrac{n - 2}{1 - R^2}} \sim t_{n-2}. \end{align}また、有意水準$\alpha$の棄却域は次で与えられる。\begin{align}(-\infty, - t_{n - 2, \alpha / 2}) \cup (t_{n - 2, \alpha / 2}, \infty),\end{align}ここに、$t_{n , \alpha}$は自由度$n$のt分布の上側$\alpha$点である。

Pearsonの相関係数だけでなく、引数methodを"Spearman"または"Kendall"に指定することで、以下のSpearmanやKendallの相関係数に基づく検定を行うことができます。母集団分布が正規性を持たないとき、ノンパラメトリック手法であるSpearmanやKendallの相関係数の検定を行います。

実行例

上で紹介した関数の実行例をみていきます。

関数cor.test

次のデータセットirisを用いてcor.testを解説していきます。

library(ggplot2)
dataset <- iris

1 2	library(ggplot2) dataset <- iris

irisのSepal.LengthとPetal.Lengthの標本相関係数と散布図は以下の通りです。

> r <- cor(dataset$Sepal.Length, dataset$Petal.Length)
> r
[1] 0.8717538

> r <- cor(dataset$Sepal.Length, dataset$Petal.Length)

> r

[1] 0.8717538

ggplot(dataset, aes(x = Sepal.Length, y = Petal.Width)) + 
  geom_point(size = 2) +
  stat_ellipse()

ggplot(dataset, aes(x = Sepal.Length, y = Petal.Width)) +

geom_point(size = 2) +

stat_ellipse()

次のように、標本相関係数の関数corと同じ書き方で相関係数の検定を行うことができます。cor.testの第1引数と第2引数に1つ目の標本と2つ目の標本を指定するだけです。

testResult <- cor.test(dataset$Sepal.Length, dataset$Petal.Length)
cor.test(~ Sepal.Length + Petal.Length, data = dataset) #別の書き方

1 2	testResult <- cor.test(dataset$Sepal.Length, dataset$Petal.Length) cor.test(~ Sepal.Length + Petal.Length, data = dataset) #別の書き方

2行目のようにformulaを使った書き方もあります。好きな方を適宜使ってください。

testResultを参照すると次の検定結果がコンソールに表示されます。

> testResult

	Pearson's product-moment correlation

data:  dataset$Sepal.Length and dataset$Petal.Length
t = 21.646, df = 148, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8270363 0.9055080
sample estimates:
      cor 
0.8717538

> testResult

Pearson's product-moment correlation

data: dataset$Sepal.Length and dataset$Petal.Length

t = 21.646, df = 148, p-value < 2.2e-16

alternative hypothesis: true correlation is not equal to 0

95 percent confidence interval:

0.8270363 0.9055080

sample estimates:

cor

0.8717538

ポイント

「p-value < 2.2e-16」とあるように有意水準0.05で相関係数は0であるとはいえないことが分かりました。すなわち、Sepal.LengthとPetal.Lengthには相関があることが分かりました。

t検定などで紹介したように、p値や信頼区間などの検定結果はtestResultの後ろに$を付けることで抽出することができます。以下のように、検定統計量の値、自由度、p値などの値を取得することができます。

stat <- testResult$statistic       #検定統計量
df <- testResult$parameter         #自由度
pValue <- testResult$p.value       #p値
estimate <- testResult$estimate    #相関係数の推定値
nullValue <- testResult$null.value #帰無仮説の下での相関係数の値
ci <- testResult$conf.int          #相関係数の信頼区間

stat <- testResult$statistic #検定統計量

df <- testResult$parameter #自由度

pValue <- testResult$p.value #p値

estimate <- testResult$estimate #相関係数の推定値

nullValue <- testResult$null.value #帰無仮説の下での相関係数の値

ci <- testResult$conf.int #相関係数の信頼区間

検定結果をcsvファイルに保存したいときは、上の検定結果をデータフレームに格納し、write.csvでデータフレームを出力するよいでしょう。

resultTable <- data.frame(stat = stat, d.f. = df, p.value = pValue,
                          c.i.lower = ci[1], c.i.upper = ci[2],
                          estimate = estimate, null.value = nullValue, row.names = NULL)
write.csv(resultTable, "相関係数の検定.csv", row.names = FALSE))

resultTable <- data.frame(stat = stat, d.f. = df, p.value = pValue,

c.i.lower = ci[1], c.i.upper = ci[2],

estimate = estimate, null.value = nullValue, row.names = NULL)

write.csv(resultTable, "相関係数の検定.csv", row.names = FALSE))

上を実行すると次の画像の検定結果のcsvファイルが出力されます。

また、cor.testの引数methodを"spearman"や"kendall"にすることでノンパラメトリックな相関係数の検定を行うことができます。

データセットUSJudgeRatingsについて、CONTとINTGのヒストグラムは次のようになります。

library(gridExtra)
dataset <- USJudgeRatings

normality_CONT <- shapiro.test(dataset$CONT)
normality_INTG <-shapiro.test(dataset$INTG)
hist_CONT <- ggplot(dataset, aes(x = CONT)) +
               geom_histogram(alpha = 0.5) +
             　annotate("text", x = 6.5, y = 5, label = paste0("normality test p-value =", round(normality_CONT$p.value, 5)))
hist_INTG <- ggplot(dataset, aes(x = INTG)) +
               geom_histogram(position = "identity", alpha = 0.5) + 
               annotate("text", x = 6.5, y = 4, label = paste0("normality test p-value =", round(normality_INTG$p.value, 5)))
grid.arrange(hist_CONT, hist_INTG)

library(gridExtra)

dataset <- USJudgeRatings

normality_CONT <- shapiro.test(dataset$CONT)

normality_INTG <-shapiro.test(dataset$INTG)

hist_CONT <- ggplot(dataset, aes(x = CONT)) +

geom_histogram(alpha = 0.5) +

　annotate("text", x = 6.5, y = 5, label = paste0("normality test p-value =", round(normality_CONT$p.value, 5)))

hist_INTG <- ggplot(dataset, aes(x = INTG)) +

geom_histogram(position = "identity", alpha = 0.5) +

annotate("text", x = 6.5, y = 4, label = paste0("normality test p-value =", round(normality_INTG$p.value, 5)))

grid.arrange(hist_CONT, hist_INTG)

上の図のようにデータが正規性を持つかどうか微妙なときは、SpearmanやKendallの相関係数に基づく検定を行う必要があります。

以下のようにPearsonの相関係数のときと同様に、SpearmanとKendallの標本相関係数の算出および相関係数の検定を行うことが可能です。

testSpearman <- cor.test(dataset$CONT, dataset$INTG, method = "spearman")
testKendall <- cor.test(dataset$CONT, dataset$INTG, method = "kendall")

1 2	testSpearman <- cor.test(dataset$CONT, dataset$INTG, method = "spearman") testKendall <- cor.test(dataset$CONT, dataset$INTG, method = "kendall")

SpearmanとKendallの相関係数の検定の結果は以下のようになります。

> stat <- testSpearman$statistic
> pValue <- testSpearman$p.value
> estimate <- testSpearman$estimate
> nullValue <- testSpearman$null.value
> resultTable_Spearman <- data.frame(stat = stat, p.value = pValue,
+                                    estimate = estimate, null.value = nullValue)
> resultTable_Spearman
      stat   p.value   estimate null.value
S 15581.27 0.2576129 -0.1764773          0
>
> stat <- testKendall$statistic
> pValue <- testKendall$p.value
> estimate <- testKendall$estimate
> nullValue <- testKendall$null.value
> resultTable_Kendall <- data.frame(stat = stat, p.value = pValue,
+ estimate = estimate, null.value = nullValue)
> resultTable_Kendall
stat p.value estimate null.value
z -1.103591 0.2697706 -0.120344 0