「ある変数Aの値が上がると,他の変数Bの値も上がる」「ある変数Aの値が上がると,他の変数Bの値は下がる」といったように,2つの変数の動きが連動しているような場合,二つの変数の間に相関があるという.そして,その関係性の強さを数値で表したものを相関係数と呼び,相関係数を計算したり,その有意性を評価したりすることを相関分析と呼ぶ.
相関係数は上述しているように,2つの変数の間の関係性の強さを数値で表したものである.数値は-1から1の間の値を取り,1に近いほど強い正の相関があることを示し,-1に近いほど強い負の相関があることを示す.0に近いほど相関が弱い,もしくは無相関であることを示す.
相関係数を算出する式は以下の通りである.
\[ \begin{align} R_{ab} &= \frac{aとbの共分散}{aの標準偏差\times bの標準偏差}\\\\ & = \frac{\frac{1}{n}\sum_{i=1}^{n} (a_i - \bar{a})(b_i - \bar{b})}{\sqrt{\frac{1}{n}\sum_{i=1}^{n} (a_i - \bar{a})^2}\sqrt{\frac{1}{n} \sum_{i=1}^{n} (b_i - \bar{b})^2}}\\ \\ & = \frac{\sum_{i=1}^{n} (a_i - \bar{a})(b_i - \bar{b})}{\sqrt{\sum_{i=1}^{n} (a_i - \bar{a})^2}\sqrt{\sum_{i=1}^{n} (b_i - \bar{b})^2}} \end{align} \]
数式そのものはそこそこ複雑な数式であるが,Rで相関を算出するのは極めて簡単である.以下のような例題を考えよう.
次のデータは,ある高校の生徒たち100人の数学と英語のテストの点数である.このデータを用いて,数学と英語の点数の相関係数を求めよ.
相関係数を算出する関数はcor()である.この関数の引数は2つあり,互いの相関を求めたい変数のベクトルがそれぞれ入る.
data <- read.csv("./practice/example_09_Correl.csv")
head(data)
## math english
## 1 54 50
## 2 58 59
## 3 76 75
## 4 61 59
## 5 61 56
## 6 77 77
summary(data)
## math english
## Min. :37.00 Min. :37.00
## 1st Qu.:55.00 1st Qu.:53.75
## Median :61.00 Median :61.50
## Mean :60.93 Mean :60.46
## 3rd Qu.:67.00 3rd Qu.:67.00
## Max. :82.00 Max. :85.00
nrow(data)
## [1] 100
cor(data$math, data$english)
## [1] 0.879465
なお,出力させる際に桁数を制限したい場合には以下のようにすればよい(->【参考】).
Res<-cor(data$math, data$english)
print(Res, digits=3)
## [1] 0.879
分野によって数値の見方は多少異なるが一般的には以下の基準が用いられる(負の場合には絶対値で考える).
| 相関係数 | 解釈 |
|---|---|
| 0.7~8以上 | 強い相関 |
| 0.5から0.7~8 | 中程度の相関 |
| 0.2~3から0.5 | 弱い相関 |
| 0.2~3未満 | 実質的に無相関 |
ただし,人の心理を扱う場合には,色々な外的影響が想定されることから,基本的に高い相関値となることはあり得ないため,以下のような基準で評価される.
| 相関係数 | 解釈 |
|---|---|
| 0.5以上 | 強い相関 |
| 0.3~0.5 | 中程度の相関 |
| 0.1~0.3 | 弱い相関 |
| 0.1未満 | 実質的に無相関 |
相関係数の有意検定とは「両者の相関係数は本当は0(無相関)である」ということを帰無仮説とする検定である.Pが0.05 (5%) を下回っていれば帰無仮説を棄却して「両者の間に有意な相関がある」と結論付ける.逆にP>=0.05であった場合には「両者は本当は無相関である可能性は否定できない」という結論となる. 相関係数の有意性を確認するにはcor.test()関数を用いる.
先の例題において,数学と英語の点数の相関係数の有意性を確認せよ.
cor.test(data$math, data$english)
##
## Pearson's product-moment correlation
##
## data: data$math and data$english
## t = 18.292, df = 98, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8256767 0.9174075
## sample estimates:
## cor
## 0.879465
この結果から,例題の数学と英語の点数の相関係数は,p値が0.05よりも小さいため,有意な相関があると結論付けられる.
2.2e-16や1.34e+5などの表記は,数字の桁数が大きい場合に用いられる表記であり,2.2e-16は\(2.2 \times
10^{-16}\)(つまり小数第15位まで0がずっと続いている数字)を,1.34e+5は\(1.34 \times 10^{5}\)を意味する.
今回の例のように<2.2e-16となっている場合は,p値が\(2.2 \times 10^{-16}\)
よりも小さいことを意味し,非常に小さな値であることを示している.
注意しないといけないのは,相関係数が有意である,というのは単に「本当の相関は0ではない」という意味でしかない.また,サンプルサイズが100より大きくなると(心理学分野でも500より大きくなると)相関係数が実質的に相関がないという値でも5%水準で有意になり得る.従って,サンプルサイズが大きい場合には,相関係数の有意性の評価に大した意義はなく,数値そのものからの実質的な相関の有無を評価すべきである.
なお,相関係数の有意性はサンプル数だけで決まるので,以下の表を用いて,算出された相関の値が表の値よりも大きな値になっていれば,5%水準,もしくは1%水準で相関係数が有意である,ということができる.
| サンプル数 | 5%有意境界 | 1%有意境界 | サンプル数 | 5%有意境界 | 1%有意境界 | |
|---|---|---|---|---|---|---|
| 3 | 0.997 | 1.000 | 31 | 0.355 | 0.456 | |
| 4 | 0.950 | 0.990 | 32 | 0.349 | 0.449 | |
| 5 | 0.878 | 0.959 | 33 | 0.344 | 0.442 | |
| 6 | 0.811 | 0.917 | 34 | 0.339 | 0.436 | |
| 7 | 0.754 | 0.875 | 35 | 0.334 | 0.430 | |
| 8 | 0.707 | 0.834 | 36 | 0.329 | 0.424 | |
| 9 | 0.666 | 0.798 | 37 | 0.325 | 0.418 | |
| 10 | 0.632 | 0.765 | 38 | 0.320 | 0.413 | |
| 11 | 0.602 | 0.735 | 39 | 0.316 | 0.408 | |
| 12 | 0.576 | 0.708 | 40 | 0.312 | 0.403 | |
| 13 | 0.553 | 0.684 | 45 | 0.294 | 0.380 | |
| 14 | 0.532 | 0.661 | 50 | 0.279 | 0.361 | |
| 15 | 0.514 | 0.641 | 60 | 0.254 | 0.330 | |
| 16 | 0.497 | 0.623 | 70 | 0.235 | 0.306 | |
| 17 | 0.482 | 0.606 | 80 | 0.220 | 0.286 | |
| 18 | 0.468 | 0.590 | 90 | 0.207 | 0.270 | |
| 19 | 0.456 | 0.575 | 100 | 0.197 | 0.256 | |
| 20 | 0.444 | 0.561 | 125 | 0.176 | 0.230 | |
| 21 | 0.433 | 0.549 | 150 | 0.160 | 0.210 | |
| 22 | 0.423 | 0.537 | 175 | 0.148 | 0.194 | |
| 23 | 0.413 | 0.526 | 200 | 0.139 | 0.182 | |
| 24 | 0.404 | 0.515 | 250 | 0.124 | 0.163 | |
| 25 | 0.396 | 0.505 | 300 | 0.113 | 0.149 | |
| 26 | 0.388 | 0.496 | 400 | 0.098 | 0.129 | |
| 27 | 0.381 | 0.487 | 500 | 0.088 | 0.115 | |
| 28 | 0.374 | 0.479 | 750 | 0.072 | 0.094 | |
| 29 | 0.367 | 0.471 | 1000 | 0.062 | 0.081 |
実際に見た目がどの程度の時に,どの程度の相関値となるかを示したのが以下のグラフである.いずれも50サンプルのデータを用いている.
サンプル数50の場合では,人の目には0.7程度でようやく関係があるように見える程度であり,弱い相関とされる0.2では到底相関があるようには見えない.この程度のサンプル数であれば,相関値よりも有意性の方が重要であると言える(サンプル数50の5%有意境界値は先の表より0.279).
3つ以上の変数がある場合,ペアの組み合わせの数だけ相関係数があり得る.Rを用いればそれらを一度に計算することができ,それらは表形式で出力される.出力結果は相関行列と呼ぶ.
次のデータは,リーダーシップに関するアンケート調査から得られたデータ(抜粋)である.このデータを用いて,すべての変数間の相関行列を求めよ.
相関行列を求める関数はcor()である.この関数の引数にはデータフレームを入れる.
data <- read.csv("./practice/example_09_CorMatrix.csv")
head(data)
## 指導.アドバイス 承認.信頼 共感.寄り添い 道具的サポート 日頃からのコミュニ
## 1 3.000000 3.0 3.000000 2.666667 2.666667
## 2 3.000000 3.5 3.000000 3.333333 3.000000
## 3 4.333333 4.0 4.000000 3.000000 4.000000
## 4 2.333333 2.5 3.000000 2.333333 2.000000
## 5 2.333333 2.0 2.000000 2.666667 1.666667
## 6 4.333333 3.0 3.333333 3.333333 3.666667
## 公正 リーダー行動 心理的安全 上司への満足度 離職意思
## 1 2.666667 2.823529 2.571429 3.000000 3
## 2 3.333333 3.176471 3.000000 3.000000 3
## 3 4.000000 3.882353 3.000000 5.000000 1
## 4 3.000000 2.529412 2.714286 2.166667 3
## 5 2.000000 2.117647 2.428571 1.833333 3
## 6 3.333333 3.529412 3.714286 4.000000 2
summary(data)
## 指導.アドバイス 承認.信頼 共感.寄り添い 道具的サポート
## Min. :1.000 Min. :1.00 Min. :1.000 Min. :1.000
## 1st Qu.:2.333 1st Qu.:2.50 1st Qu.:2.333 1st Qu.:2.000
## Median :3.000 Median :3.00 Median :3.000 Median :3.000
## Mean :3.062 Mean :3.01 Mean :2.988 Mean :2.792
## 3rd Qu.:4.000 3rd Qu.:4.00 3rd Qu.:3.667 3rd Qu.:3.667
## Max. :5.000 Max. :5.00 Max. :5.000 Max. :5.000
## 日頃からのコミュニ 公正 リーダー行動 心理的安全
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:2.667 1st Qu.:2.333 1st Qu.:2.529 1st Qu.:2.857
## Median :3.333 Median :3.000 Median :3.059 Median :3.143
## Mean :3.117 Mean :3.000 Mean :2.994 Mean :3.140
## 3rd Qu.:4.000 3rd Qu.:3.667 3rd Qu.:3.647 3rd Qu.:3.429
## Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.000
## 上司への満足度 離職意思
## Min. :1.000 Min. :1.000
## 1st Qu.:2.167 1st Qu.:2.000
## Median :3.000 Median :3.000
## Mean :3.009 Mean :2.825
## 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :5.000 Max. :5.000
nrow(data)
## [1] 200
cor(data)
## 指導.アドバイス 承認.信頼 共感.寄り添い 道具的サポート
## 指導.アドバイス 1.0000000 0.8297383 0.8374576 0.8457553
## 承認.信頼 0.8297383 1.0000000 0.8031177 0.7952222
## 共感.寄り添い 0.8374576 0.8031177 1.0000000 0.7963156
## 道具的サポート 0.8457553 0.7952222 0.7963156 1.0000000
## 日頃からのコミュニ 0.7932627 0.7845982 0.8470742 0.7658898
## 公正 0.8748846 0.8477138 0.8848106 0.8329005
## リーダー行動 0.9365796 0.9020516 0.9350335 0.9099032
## 心理的安全 0.6044248 0.6049806 0.5627454 0.5730466
## 上司への満足度 0.8224389 0.7734356 0.8110860 0.7321629
## 離職意思 -0.3374343 -0.3911060 -0.3782470 -0.2940360
## 日頃からのコミュニ 公正 リーダー行動 心理的安全
## 指導.アドバイス 0.7932627 0.8748846 0.9365796 0.6044248
## 承認.信頼 0.7845982 0.8477138 0.9020516 0.6049806
## 共感.寄り添い 0.8470742 0.8848106 0.9350335 0.5627454
## 道具的サポート 0.7658898 0.8329005 0.9099032 0.5730466
## 日頃からのコミュニ 1.0000000 0.8167715 0.9070783 0.6238759
## 公正 0.8167715 1.0000000 0.9491341 0.6146664
## リーダー行動 0.9070783 0.9491341 1.0000000 0.6459790
## 心理的安全 0.6238759 0.6146664 0.6459790 1.0000000
## 上司への満足度 0.7683507 0.7925804 0.8481132 0.6225417
## 離職意思 -0.2762140 -0.4076889 -0.3724925 -0.4499767
## 上司への満足度 離職意思
## 指導.アドバイス 0.8224389 -0.3374343
## 承認.信頼 0.7734356 -0.3911060
## 共感.寄り添い 0.8110860 -0.3782470
## 道具的サポート 0.7321629 -0.2940360
## 日頃からのコミュニ 0.7683507 -0.2762140
## 公正 0.7925804 -0.4076889
## リーダー行動 0.8481132 -0.3724925
## 心理的安全 0.6225417 -0.4499767
## 上司への満足度 1.0000000 -0.4338144
## 離職意思 -0.4338144 1.0000000
相関行列は対角成分を除いて対称行列となる.また,対角成分は自分を相手とする相関なので必ず1となる.
cor()関数のそのままの出力だと小数以下の表示桁数が多くなってしまうので,以下のようにして表示桁数を制限する.
Res<-cor(data)
print(Res, digits=3)
## 指導.アドバイス 承認.信頼 共感.寄り添い 道具的サポート
## 指導.アドバイス 1.000 0.830 0.837 0.846
## 承認.信頼 0.830 1.000 0.803 0.795
## 共感.寄り添い 0.837 0.803 1.000 0.796
## 道具的サポート 0.846 0.795 0.796 1.000
## 日頃からのコミュニ 0.793 0.785 0.847 0.766
## 公正 0.875 0.848 0.885 0.833
## リーダー行動 0.937 0.902 0.935 0.910
## 心理的安全 0.604 0.605 0.563 0.573
## 上司への満足度 0.822 0.773 0.811 0.732
## 離職意思 -0.337 -0.391 -0.378 -0.294
## 日頃からのコミュニ 公正 リーダー行動 心理的安全
## 指導.アドバイス 0.793 0.875 0.937 0.604
## 承認.信頼 0.785 0.848 0.902 0.605
## 共感.寄り添い 0.847 0.885 0.935 0.563
## 道具的サポート 0.766 0.833 0.910 0.573
## 日頃からのコミュニ 1.000 0.817 0.907 0.624
## 公正 0.817 1.000 0.949 0.615
## リーダー行動 0.907 0.949 1.000 0.646
## 心理的安全 0.624 0.615 0.646 1.000
## 上司への満足度 0.768 0.793 0.848 0.623
## 離職意思 -0.276 -0.408 -0.372 -0.450
## 上司への満足度 離職意思
## 指導.アドバイス 0.822 -0.337
## 承認.信頼 0.773 -0.391
## 共感.寄り添い 0.811 -0.378
## 道具的サポート 0.732 -0.294
## 日頃からのコミュニ 0.768 -0.276
## 公正 0.793 -0.408
## リーダー行動 0.848 -0.372
## 心理的安全 0.623 -0.450
## 上司への満足度 1.000 -0.434
## 離職意思 -0.434 1.000
さらに,データフレームをそのまま入力すると,データフレームに含まれるすべての変数同士の相関が算出されるが,特定の変数同士の相関だけを算出したい場合には,以下のようにする->参考.
Res<-cor(data[c(1,7,8,9,10)])
# 以下のようにしてもよい
# Res<-cor(data[c(1,7:10)]) # 連続している箇所をシーケンス演算子で省略表記した
# あるいは,以下のように直接変数名を書き込んでもよい.
# また下記の通り変数名を指定する際には変数名と変数名の間には改行が入っていても構わない
# Res<-cor(data[c("指導.アドバイス",
# "リーダー行動",
# "心理的安全",
# "上司への満足度",
# "離職意思")])
print(Res, digits=3)
## 指導.アドバイス リーダー行動 心理的安全 上司への満足度 離職意思
## 指導.アドバイス 1.000 0.937 0.604 0.822 -0.337
## リーダー行動 0.937 1.000 0.646 0.848 -0.372
## 心理的安全 0.604 0.646 1.000 0.623 -0.450
## 上司への満足度 0.822 0.848 0.623 1.000 -0.434
## 離職意思 -0.337 -0.372 -0.450 -0.434 1.000
論文などでは,しばしば調査で用いた変数の相関行列が表形式で記述されるが,多くの場合対角成分よりも下のものだけを記述されている.
(出典:Schaubroeck, J., Lam, S. S. K.,
& Peng, A. C. (2011). Cognition-based and affect-based trust as
mediators of leader behavior influences on team performance. Journal of
Applied Psychology, 96(4), 863–871. https://doi.org/10.1037/a0022625 )
なおこの例では,相関行列と合わせて,各尺度のM(平均)とSD(標準偏差)も記述している.論文ではこのような表記の仕方をするのが一般的である.
さらにこの例では,対角行列の1の代わりに(.88)や(.90)といった値が記述されているが,これはクロンバックの\(\alpha\)係数と呼ばれる統計量で,心理学的なアンケート調査(複数の質問項目の回答を集約(合算したり,平均にしたり)して一つの尺度とするもの)でよく用いられる,心理尺度の信頼性(内的整合性,合算することの妥当性)を示す指標である.この値が高いほど,その尺度が信頼性が高いとされる.一般的には0.7以上であれば信頼性が高いとされる.
相関行列を可視化する方法として,ヒートマップが用いられることがある.ヒートマップは,相関係数の値によって色を変えて表示することで,相関係数の強さを視覚的に把握することができる.
以下のコードは,先のCSVファイルから算出した相関行列をヒートマップで表示する例である.\(-1\)を青,\(1\)を赤,\(0\)を白に設定し,間の値はそれに合わせてグラデーションしていくような形で表示させている.
library(ggplot2)
library(reshape2)
## Warning: パッケージ 'reshape2' はバージョン 4.4.2 の R の下で造られました
##
## 次のパッケージを付け加えます: 'reshape2'
## 以下のオブジェクトは 'package:tidyr' からマスクされています:
##
## smiths
data <- read.csv("./practice/example_09_CorMatrix.csv")
cor_matrix <- cor(data)
melted_cor_matrix <- melt(cor_matrix)
ggplot(melted_cor_matrix, aes(Var1, Var2, fill = value)) +
geom_tile() +
scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
今回の例の場合,表示がいずれも1に近いため,赤色に近い色が多くなっているが,一般的にはもっと薄い色が多い.
このようにして可視化することで,相関係数の強さを視覚的に把握することができる.数値的に関連しているからといって,実際にそれらの変数の元となる「事象」が実際に関連しているとは限らない.
アイスの売上高と海辺での人身事故は一般に高い相関係数を示す.実際に2018年~2022年の間での月次の統計値を散布図にしてみよう(参考.なお,統計値はアイスクリームについては総務省統計局 家計調査(月次・356アイスクリーム・シャーベット)から,人身事故の発生件数は海上保安庁の統計年報から取得した.CSVはこちら.
data <- read.csv("./practice/example_09_ice_and_accident.csv")
# データの確認
head(data)
## 年 月 アイスクリームの購入金額_世帯当たり 海浜事故発生件数_全国計
## 1 2018年 1 月 507 75
## 2 2018年 2 月 416 84
## 3 2018年 3 月 607 132
## 4 2018年 4 月 746 125
## 5 2018年 5 月 894 114
## 6 2018年 6 月 1021 117
## 東京の月別平均気温
## 1 4.7
## 2 5.4
## 3 11.5
## 4 17.0
## 5 19.8
## 6 22.0
# 散布図
library(ggplot2)
ggplot(data, aes(x=`アイスクリームの購入金額_世帯当たり`, y=`海浜事故発生件数_全国計`)) +
geom_point() +
labs(title="アイスの売上高と海辺での人身事故の関係", x="アイスの売上高", y="海辺での人身事故の発生件数")
見ての通り,非常に強い相関があることがうかがえる. 実際に相関係数を求めてみよう.
cor(data$`アイスクリームの購入金額_世帯当たり`, data$`海浜事故発生件数_全国計`)
## [1] 0.8582547
cor.test(data$`アイスクリームの購入金額_世帯当たり`, data$`海浜事故発生件数_全国計`)
##
## Pearson's product-moment correlation
##
## data: data$アイスクリームの購入金額_世帯当たり and data$海浜事故発生件数_全国計
## t = 12.736, df = 58, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7727333 0.9131703
## sample estimates:
## cor
## 0.8582547
以下のように高い相関となっている(当然ながら有意な相関値である).しかし,当然ながら,アイスが売れたから海辺での人身事故が発生するわけではないし,海辺での人身事故が起こるとアイスが売れるわけでもない.このような関係を疑似相関と呼ぶ.
考えられるのは,気温が背後要因にある,ということである.気温が高くなるとアイスクリームの購入金額が増えるのは十分に理解できる.同様に気温が高くなると海でのレジャーに出かける人が増えるので,それに比例して事故の発生件数も増えるはずであろう.実際に2018年から2022年の東京の月別平均気温とこれらのデータとの関係を図示すると以下のようになる.
ggplot(data, aes(x=`東京の月別平均気温`, y=`アイスクリームの購入金額_世帯当たり`)) +
geom_point() +
labs(title="東京の月別平均気温とアイスの売上高の関係", x="東京の月別平均気温", y="アイスの売上高")
cor(data$`東京の月別平均気温`, data$`アイスクリームの購入金額_世帯当たり`)
## [1] 0.9114073
ggplot(data, aes(x=`東京の月別平均気温`, y=`海浜事故発生件数_全国計`)) +
geom_point() +
labs(title="東京の月別平均気温と海浜事故発生件数", x="東京の月別平均気温", y="海浜事故発生件数")
cor(data$`東京の月別平均気温`, data$`海浜事故発生件数_全国計`)
## [1] 0.7762446
このように,気温とアイスの売上高,気温と海浜事故発生件数の間にも強い相関があることがわかる.このように,疑似相関は背後にある要因が同じであることが原因で生じることがある.
次の例を見てほしい.2021年度NHK受信料世帯支払率と都道府県別の人口あたりの新型コロナウイルス感染者数の推定値から取得した2022年10月1日時点の累計コロナ罹患者の散布図は以下となり,相関値は-0.85と非常に強い負の相関を示す.
では,NHKの受信料世帯支払率とコロナ罹患者数との間に直接的な関係性があると言えるだろうか?普通に考えて,NHK受信料をキチンと支払う世帯が多くなるとコロナに罹患しなくなる,あるいは,コロナに罹患したらNHK受信料を支払わなくなる,ということは考えにくい.
これに関して考察するとすると, まず,NHK受診料をきちんと支払う世帯が多い地域を見ていると,東北や北陸,山陰などが多い. 実際に支払い率の上位10県は以下の通りである(ちなみに福井は88.5%で11位).
| ランク | 都道府県 | 支払率 |
|---|---|---|
| 1 | 秋田 | 97.9 |
| 2 | 新潟 | 94.9 |
| 3 | 岩手 | 94.6 |
| 4 | 島根 | 94.3 |
| 5 | 山形 | 93.5 |
| 6 | 鳥取 | 92.9 |
| 7 | 青森 | 92.5 |
| 8 | 富山 | 91.9 |
| 9 | 山口 | 91.1 |
| 10 | 岐阜 | 89.4 |
これらの地域に共通しているのは高齢化率が高いことだろう.実際に2021年(令和3年)の高齢化率についてのデータをもとに,高齢化率とNHK受信料支払率の相関を求めてみると,0.68と高い相関があることがわかる.社会問題の一つとして若年層がNHK受信料を支払わない(NHKを視聴しない)ことが問題視されていることから,高齢者ほどNHK受信料を真面目に払っているということは十分に考えられる.
翻って,コロナ感染については,高齢者ほど重症化しやすいと言われていたことから,コロナ対策は高齢者ほど厳格に行っていたであろうし,ワクチンの接種も高齢者を優先して進められていた.こうしたことから,高齢者が多い地域ほどコロナ感染者数が少ないということは十分に考えられる (重症患者に絞れば,逆に高齢者が多い地域ほど重症患者が多い,ということかもしれない).
すなわち,NHK受信料支払率とコロナ感染者数は,互いに「地域の高齢化」という要因と強い関連があると考えられ,この背後要因によって,表面上,NHK受信料支払率とコロナ感染者数の間に高い相関が生じてしまっていたのではないかと考えられる.このように考えると,両者の間の相関(疑似相関)が生じた理屈として納得がいくだろう.
NHK受信料支払率とコロナ感染者数,高齢化率についてのデータこちらの以下の変数の組み合わせについて,散布図と相関係数をそれぞれ求めよ.
このように,相関係数は単なる数値計算で算出されるものなので,理論的にはなんの直接的な関係がないものであっても高い相関が算出されることがある.このような関係のことを疑似相関と呼ぶ. 疑似相関が生じるのは,背後に別の共通要因があり,その共通要因がそれぞれと相関しているためである.
このような場合には,AとBの間の相関係数の高さは共通要因CとA,CとBのそれぞれの相関の高さに由来している.そこで,CからA,Bへの相関の影響を取り除いた,より純粋なAとBの間の相関をもとめたい.そのような場合には以下の数式を用いる.こうして得られた値のことを偏相関と呼ぶ. \(R_{ab}\)はA-Bの相関,\(R_{ca}\)はC-Aの相関,\(R_{cb}\)はC-Bの相関とする.
\[ 偏相関R_{ab|c} = \frac{R_{ab} - R_{ca}\times R_{cb}}{\sqrt{(1-R_{ca}^2)}\sqrt{(1-R_{cb}^2)}} \]
なお,Cの相関の影響を取り除くことを「Cの影響を統制する」と表現する. 相関係数自体には特に因果関係の方向性はないため,上記の図で言えば,A-B間の偏相関のほか,B-C間の偏相関(Aの影響を統制),C-A間の偏相関(Bの影響を統制)という形で,3つの偏相関を考えることができる.同様に,4つの変数を見ている場合には6つの偏相関を,5つの変数を見ている場合には,10の偏相関を考えることができる.このように見ている変数の間でのペアの組み合わせの数だけ偏相関を考えることができる.
偏相関を求めるには,ppcorパッケージのpcor()関数もしくはpcor.test()関数を用いる.pcor()関数は,データフレームを引数に取り,そのデータフレームに含まれる変数間の偏相関行列と各偏相関の有意確率,並びに有意確率を算出するための統計量を返す.
pcor.test()関数は,3つの変数を引数に取り,第3引数に置いた変数の影響を統制した場合の第1引数と第2引数の両変数間の偏相関係数とその有意確率を返す.
いずれも利用にあたっては,「Tools」->「install.packages」からppcorパッケージをインストールしてた上で,library(ppcor)でパッケージを読み込む必要がある.
アイスの売上高と海辺での人身事故のデータを用いて,アイスの売上高と海辺での人身事故の間の偏相関を求めよ.
library(ppcor)
## Warning: パッケージ 'ppcor' はバージョン 4.4.2 の R の下で造られました
## 要求されたパッケージ MASS をロード中です
##
## 次のパッケージを付け加えます: 'MASS'
## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## select
data <- read.csv("./practice/example_09_ice_and_accident.csv")
head(data)
## 年 月 アイスクリームの購入金額_世帯当たり 海浜事故発生件数_全国計
## 1 2018年 1 月 507 75
## 2 2018年 2 月 416 84
## 3 2018年 3 月 607 132
## 4 2018年 4 月 746 125
## 5 2018年 5 月 894 114
## 6 2018年 6 月 1021 117
## 東京の月別平均気温
## 1 4.7
## 2 5.4
## 3 11.5
## 4 17.0
## 5 19.8
## 6 22.0
# 普通の相関行列の算出
Res<-cor(data[c("アイスクリームの購入金額_世帯当たり","海浜事故発生件数_全国計","東京の月別平均気温")])
print(Res, digits=3)
## アイスクリームの購入金額_世帯当たり
## アイスクリームの購入金額_世帯当たり 1.000
## 海浜事故発生件数_全国計 0.858
## 東京の月別平均気温 0.911
## 海浜事故発生件数_全国計 東京の月別平均気温
## アイスクリームの購入金額_世帯当たり 0.858 0.911
## 海浜事故発生件数_全国計 1.000 0.776
## 東京の月別平均気温 0.776 1.000
# 偏相関行列の算出
Res.pcor<-pcor(data[c("アイスクリームの購入金額_世帯当たり","海浜事故発生件数_全国計","東京の月別平均気温")])
print(Res.pcor, digits=3)
## $estimate
## アイスクリームの購入金額_世帯当たり
## アイスクリームの購入金額_世帯当たり 1.000
## 海浜事故発生件数_全国計 0.581
## 東京の月別平均気温 0.758
## 海浜事故発生件数_全国計 東京の月別平均気温
## アイスクリームの購入金額_世帯当たり 0.5812 0.7578
## 海浜事故発生件数_全国計 1.0000 -0.0283
## 東京の月別平均気温 -0.0283 1.0000
##
## $p.value
## アイスクリームの購入金額_世帯当たり
## アイスクリームの購入金額_世帯当たり 0.00e+00
## 海浜事故発生件数_全国計 1.39e-06
## 東京の月別平均気温 3.69e-12
## 海浜事故発生件数_全国計 東京の月別平均気温
## アイスクリームの購入金額_世帯当たり 1.39e-06 3.69e-12
## 海浜事故発生件数_全国計 0.00e+00 8.32e-01
## 東京の月別平均気温 8.32e-01 0.00e+00
##
## $statistic
## アイスクリームの購入金額_世帯当たり
## アイスクリームの購入金額_世帯当たり 0.00
## 海浜事故発生件数_全国計 5.39
## 東京の月別平均気温 8.77
## 海浜事故発生件数_全国計 東京の月別平均気温
## アイスクリームの購入金額_世帯当たり 5.392 8.769
## 海浜事故発生件数_全国計 0.000 -0.214
## 東京の月別平均気温 -0.214 0.000
##
## $n
## [1] 60
##
## $gp
## [1] 1
##
## $method
## [1] "pearson"
# 偏相関の検定
pcor.test(data$`アイスクリームの購入金額_世帯当たり`, data$`海浜事故発生件数_全国計`, data$`東京の月別平均気温`)
## estimate p.value statistic n gp Method
## 1 0.5812047 1.388535e-06 5.39227 60 1 pearson
結果として,元々アイスクリームの売上高と海浜事故発生件数の相関係数は0.712であったが,東京の月別平均気温の影響を統制した場合の偏相関係数は0.542となり,単純な相関係数より偏相関係数が低くなっていることがわかる. ただ,それでもまだ強い相関値が出ているため,その他の共通要因も潜在している可能性がある.その点はさらに別なデータを用いて検証する必要があるだろう.
なお,偏相関行列を見てみると,アイスクリームの売り上げ高で統制した場合の海浜事故発生件数と気温の偏相関係数は -0.0283という値が出ている.この結果から,海浜事故と気温の間には関係がないのではないかという主張をすることが考えられるが,これは誤りである. あくまで因果関係として,気温が高いとアイスクリームの売上高が上がるし,気温が高いと海で遊ぶ人が増えて結果として海浜事故が増える,という関係が論理的に成り立ち得るからこそ,気温で統制した場合のアイスクリームの売上高と海浜事故発生件数の偏相関係数を算出することに意味が出てくる. アイスクリームの売り上げで統制した場合の海浜事故発生件数と気温の偏相関係数を求める,ということは,アイスクリームの売り上げが気温と事故発生件数の共通の背後要因になっているという因果関係を想定していることとなる.そのような因果関係が成立しないのは明らかである.
このように偏相関を検討するときには,あくまで事前にどのような変数間にどのような因果関係があるのかを理論的に考えて仮説を立てておくことが大切である.
先ほどのコロナ感染者数とNHK受信料支払世帯のデータについて,高齢化率を統制した場合の偏相関を求めよ.
疑似相関は,本当は相関がないのに数値上,高い相関の値となる,という現象であるが,逆に,本当は相関があるのに,数値上相関が低い値となるという現象もある.
以下はシミュレーションとして作成したデータであるが,パッと見は相関がないように見えるが,色分けすると正負の相関があり,それが組み合わさったデータになっていることがわかる.
具体的な例としては,日本においては低所得世帯と高所得者世帯とでは,食事の摂取カロリー量にほとんど違いはない(平成 30 年 国民健康・栄養調査結果の概要より),一方で,米国においては低所得世帯と高所得世帯では,食事の摂取カロリー量は低所得世帯の方が高くなりがちであると言われている(子どもの頃の家庭環境と健康格差:肥満の要因分析(研究論文), 米国人の肥満率,10年以内に50%超える恐れ 研究報告(CNNニュースより)).こうしたデータを,「国別である」ということを考慮せずに,まとめて相関分析に掛けた場合,アメリカにおける世帯収入と摂取カロリーの相関関係が,日本の低い相関関係によって希薄化する可能性がある.
このように一見相関がないという結果になっていたとしても,ある属性でデータを切り分けてみると,それぞれの属性で逆方向の相関ががある,ということはたまにある.このため,相関分析を行う際には,手許にある生のデータを実際に眺めてみて「もしかして,この属性で分けたら相関関係が出てくるのでは?」という直感を働かせてみたり,属性ごとで散布図を描いてみるなどして,潜在していた相関関係を探ることが大切である.
このデータは高校生における対人依存欲求(他者と親密な関係を築くことで自身の心の安定を得るとともに,問題に直面したときに他者に助けを求めようとする欲求)と,SNSの利用時間との関係を調べたものである.
※ 本データは次の論文を参考に作成した.
相関係数の差の検定
以下の内容はかなり複雑なものであるとともに,実際にはそう使うものでもないので,興味があれば見るという程度でよい.
相関係数の差が有意かどうかを検定することが可能である.帰無仮説は,2つの相関係数が等しいということである. 差の有意性を検定することによって,二つの要因のいずれがより強くある変数と相関しているかを比較することが可能となる.
相関の差の検定の方法には次の3つの種類がある. いずれもcocorパッケージを用いるので,パッケージのインストールと読み込みが必要である.
互いに異なるサンプルから得られた2つの相関係数の差が有意かどうかを検定するには,cocorパッケージのcocor.indep.groups()関数を用いる.引数は以下の通りである.
ある高校の3年生(男子165名,女子181名)で,数学と英語のテストの点数についての相関を調べたところ,男性の相関係数は0.65,女性の相関係数は0.72であった.この2つの相関係数の差が有意かどうかを検定せよ.
library(cocor)
## Warning: パッケージ 'cocor' はバージョン 4.4.2 の R の下で造られました
# 独立サンプルにおける相関の差の検定
cocor.indep.groups(r1=0.65, r2=0.72, n1=165, n2=181)
##
## Results of a comparison of two correlations based on independent groups
##
## Comparison between r1.jk = 0.65 and r2.hm = 0.72
## Difference: r1.jk - r2.hm = -0.07
## Group sizes: n1 = 165, n2 = 181
## Null hypothesis: r1.jk is equal to r2.hm
## Alternative hypothesis: r1.jk is not equal to r2.hm (two-sided)
## Alpha: 0.05
##
## fisher1925: Fisher's z (1925)
## z = -1.2188, p-value = 0.2229
## Null hypothesis retained
##
## zou2007: Zou's (2007) confidence interval
## 95% confidence interval for r1.jk - r2.hm: -0.1868 0.0423
## Null hypothesis retained (Interval includes 0)
結果の表記は大きく3つのパートに分かれているが,前半は入力した値の情報が記載されているだけである.
検定結果は中段のfisher1925の部分にあり,この中のp-valueが有意確率である.これが有意水準0.05より小さければ,帰無仮説を棄却し,相関係数の差が有意であると結論する(この検定のことをFisherのZ検定もしくはFisherのZ変換を用いた検定と呼ぶ).
下段にあるのは,差の推定値の信頼区間が表記されており,今回の例の場合には,信頼区間は-0.1868から0.0423である.
1つの変数を共有した相関とは,変数J, K, Hという3つの変数があるとき,Jを共通の変数とした(オーバーラップさせた) JとKの相関係数と,JとHの相関係数の差を検定することである.
この場合には,cocorパッケージのcocor.dep.groups.overlap()関数を用いる.引数は以下の通りである.
ある高校の男子生徒165名で,数学と英語と物理のテストの点数についての相関を調べたところ,数学と英語の相関係数は0.65,数学と物理の相関係数は0.72,物理と英語の相関は0.52であった.数学と英語の相関係数と数学と物理の差が有意かどうかを検定せよ.
# 同一サンプルにおける1つの変数が共通している相関の差の検定
cocor.dep.groups.overlap(0.65, 0.72, 0.52, 165)
##
## Results of a comparison of two overlapping correlations based on dependent groups
##
## Comparison between r.jk = 0.65 and r.jh = 0.72
## Difference: r.jk - r.jh = -0.07
## Related correlation: r.kh = 0.52
## Group size: n = 165
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
##
## pearson1898: Pearson and Filon's z (1898)
## z = -1.4127, p-value = 0.1577
## Null hypothesis retained
##
## hotelling1940: Hotelling's t (1940)
## t = -1.4800, df = 162, p-value = 0.1408
## Null hypothesis retained
##
## williams1959: Williams' t (1959)
## t = -1.4156, df = 162, p-value = 0.1588
## Null hypothesis retained
##
## olkin1967: Olkin's z (1967)
## z = -1.4127, p-value = 0.1577
## Null hypothesis retained
##
## dunn1969: Dunn and Clark's z (1969)
## z = -1.4119, p-value = 0.1580
## Null hypothesis retained
##
## hendrickson1970: Hendrickson, Stanley, and Hills' (1970) modification of Williams' t (1959)
## t = -1.4800, df = 162, p-value = 0.1408
## Null hypothesis retained
##
## steiger1980: Steiger's (1980) modification of Dunn and Clark's z (1969) using average correlations
## z = -1.4111, p-value = 0.1582
## Null hypothesis retained
##
## meng1992: Meng, Rosenthal, and Rubin's z (1992)
## z = -1.4105, p-value = 0.1584
## Null hypothesis retained
## 95% confidence interval for r.jk - r.jh: -0.3163 0.0516
## Null hypothesis retained (Interval includes 0)
##
## hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
## z = -1.4099, p-value = 0.1586
## Null hypothesis retained
##
## zou2007: Zou's (2007) confidence interval
## 95% confidence interval for r.jk - r.jh: -0.1712 0.0272
## Null hypothesis retained (Interval includes 0)
さまざまな種類の検定結果が出力されるが,この中でsteiger1980の結果を見るのが一般的である.今回の結果だと,P値は0.1582となっており有意ではない,という結果となっている(この検定のことをStiegerのZ検定と呼ぶ).
加えて,一番下の段
(zou2007)には,差の推定値の信頼区間が表記されている.今回の例の場合には,信頼区間は-
-0.1712から0.0272となっている.
同一サンプルにおいて,互いに異なる変数間の相関の差(J-K間の相関とH-M間の相関)が有意かどうかを検定する場合には,cocorパッケージのcocor.dep.groups.nonoverlap()関数を用いる.引数は以下の通りである.
与えなければならない引数が多く,この順番が入れ替わると結果もかわってしまうので注意が必要である.それぞれの引数の添え字に注意してほしい.
ある高校の女子生徒181名で,数学と物理と英語と国語の点数の相関を調べたところ,次のような表となった.
| 数学 | 物理 | 英語 | 国語 | |
|---|---|---|---|---|
| 数学 | 1 | |||
| 物理 | 0.41 | 1 | ||
| 英語 | 0.53 | 0.62 | 1 | |
| 国語 | 0.32 | 0.43 | 0.72 | 1 |
数学と物理の相関係と,英語と国語の相関の差が有意かどうか検定せよ.
この場合,数学がJ,物理がK,英語がH,国語がMだと考えればよい.
# 同一サンプルにおける互いに異なる変数間の相関の差の検定
cocor.dep.groups.nonoverlap(r.jk=0.42, r.hm=0.72,
r.jh=0.53, r.jm=0.32,
r.kh=0.62, r.km=0.43,
181)
##
## Results of a comparison of two nonoverlapping correlations based on dependent groups
##
## Comparison between r.jk = 0.42 and r.hm = 0.72
## Difference: r.jk - r.hm = -0.3
## Related correlations: r.jh = 0.53, r.jm = 0.32, r.kh = 0.62, r.km = 0.43
## Group size: n = 181
## Null hypothesis: r.jk is equal to r.hm
## Alternative hypothesis: r.jk is not equal to r.hm (two-sided)
## Alpha: 0.05
##
## pearson1898: Pearson and Filon's z (1898)
## z = -4.4920, p-value = 0.0000
## Null hypothesis rejected
##
## dunn1969: Dunn and Clark's z (1969)
## z = -4.6514, p-value = 0.0000
## Null hypothesis rejected
##
## steiger1980: Steiger's (1980) modification of Dunn and Clark's z (1969) using average correlations
## z = -4.6963, p-value = 0.0000
## Null hypothesis rejected
##
## raghunathan1996: Raghunathan, Rosenthal, and Rubin's (1996) modification of Pearson and Filon's z (1898)
## z = -4.6514, p-value = 0.0000
## Null hypothesis rejected
##
## silver2004: Silver, Hittner, and May's (2004) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
## z = -4.6803, p-value = 0.0000
## Null hypothesis rejected
##
## zou2007: Zou's (2007) confidence interval
## 95% confidence interval for r.jk - r.hm: -0.4353 -0.1709
## Null hypothesis rejected (Interval does not include 0)
結果はsteiger1980のP値を読み取る.今回の例だと0.0000となっていおり,0.1%の有意水準でも有意ということになる(必ずしもP=0というわけではなく,小数第4位までは0であった,ということである).
これまでと同様に最下部には差の推定値の信頼区間が表記されている.
平均や分散,比率の比較のときと同様に,3つ以上の相関係数の比較を行う場合には,多重比較の考え方を導入する.すなわち,まずは一対比較を行い,その後,p.adjust()関数を用いて多重比較の補正を行う.
たとえば,J-KとJ-Hの間の相関を比較したい場合に,その背後で別の変数ZがJ-K,J-Hの相関に影響を与えている場合には,Zを統制した偏相関を求め,その差の有意性を検定することができる.その方法は,単にZを統制したJ-K,J-Hの偏相関を求め,その偏相関を上記のcocor.indep.groups()やcocor.dep.groups.overlap(),cocor.dep.groups.nonoverlap()にそれぞれ代入すればよい.
ある組織で従業員満足度調査として,従業員323名(課長級以上の管理職は除く)に対して以下の項目について調査した.
その結果,次のようなデータが得られた.このデータについて以下の課題を行え.