分散共分散行列と相関係数行列を用いることで、様々な定量的手法が可能となる。
実際の手法については「中級レベル」で解説するとして、要するに、
これら行列の意味を理解することは極めて重要である。
そもそも、この二つの行列は、どのようなものであったか？
たしか、分散共分散行列は、各変数の分散が対角成分に並び、
それ以外の成分には「共分散」が入っているような「対称行列」であった。
一方、相関係数行列は、対角成分に「1」が並び、
それ以外の部分には「相関係数」が入っているような行列であった。
相関係数行列は、変数間の関係を一目で確認できるので、色々と重宝するのだった。
ふむふむ。そのような話を確かにした。何か「問題」でもあるのか？
実は、2変数の相関係数の話では気づかなかったのであるが、
多変数の相関係数行列にしてみると、ある種の疑惑が生まれる。
本当に、二つの変数の関係「のみ」を表しているのか？
実は、他の変数の影響を受けているのではないか？
という疑惑である。
例えば、3つの変数、「
」、「
」、「
」の3つの変数が存在した場合、
この3つの変数の相関係数行列は以下のようになるはず。
ふむ。何も不審な点は見当たらないように見えるが、
ひょっとすると、「」と「」の関係の中には、
「」の影響が「潜在的」に存在するかもしれない。
2変数の間の相関係数の話では、意識する必要は無いが、
たしかに、多変数となる相関係数行列では、そういった懸念が生じる。
では、そのような別の変数の影響を取り除くにはどうすれば良いか？
この問題が、本日のメインテーマ。
とりあえず、いつものように、いつもの如く、チュートリアルデータを準備する。
今回のデータは、総務省統計局のデータ（e-stat）のデータを用いて、
チュートリアル用の「擬似」データを作成したもの。
今回のデータは、居住世帯の有無と一ヶ月辺りの平均家賃の関係を見るためのデータ。
家賃のデータは、19階級に区分し家賃を集計したものであったので、
各階級の平均を一ヶ月の家賃とし、世帯数を掛けたもので都道府県別の平均家賃を算出。
かなり「大雑把な加工」をしているので、正式な分析には用いてはならない。
あくまで、手法を学ぶための「擬似データ」であるので、
ここから実際の解釈を導いても、「深い意味」はない。
一応、注意事項を整理したところで、いつもどおりにデータを読み込む。
# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",", header=TRUE)
# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",", header=TRUE)
今回のチュートリアルのために用意したデータは以下の通り。
OCU,UNO,PRI
Hokkaido,2340300,390100,40687
Aomori-ken,493500,87400,111566
Iwate-ken,470700,78800,17002
Miyagi-ken,869700,144200,16084
Akita-ken,380300,57100,41254
Yamagata-ken,383000,49700,10523
Fukushima-ken,699700,108600,11879
Ibaraki-ken,1036200,187600,25207
Tochigi-ken,708700,131300,36003
Gumma-ken,725300,130400,27494
Saitama-ken,2688000,341100,22045
Chiba-ken,2344500,373100,94198
Tokyo-to,5939900,840500,60863
Kanagawa-ken,3612200,455600,246996
Niigata-ken,810700,119000,113336
Toyama-ken,368800,55500,24592
Ishikawa-ken,421600,76400,10485
Fukui-ken,259700,49000,15203
Yamanashi-ken,314600,83600,6005
Nagano-ken,758300,188000,10698
Gifu-ken,712600,123100,25608
Shizuoka-ken,1359400,238500,22521
Aichi-ken,2764400,368400,54336
Mie-ken,680900,110100,122595
Shiga-ken,491300,76300,17786
Kyoto-fu,1086800,183400,14849
Osaka-fu,3685100,660900,39764
Hyogo-ken,2169400,351300,150179
Nara-ken,502500,90100,68843
Wakayama-ken,382100,85700,12267
Tottori-ken,208600,38600,10868
Shimane-ken,249900,45900,6686
Okayama-ken,734700,131900,7459
Hiroshima-ken,1147600,208700,27600
Yamaguchi-ken,584100,107600,48448
Tokushima-ken,297000,58600,21579
Kagawa-ken,372700,73700,10305
Ehime-ken,574000,107100,12877
Kochi-ken,312800,64900,22953
Fukuoka-ken,2034000,340800,11205
Saga-ken,286100,36800,103313
Nagasaki-ken,539200,91900,9725
Kumamoto-ken,663800,105700,20633
Oita-ken,467200,79300,29546
Miyazaki-ken,443800,65900,22566
Kagoshima-ken,718200,133100,17694
Okinawa-ken,504400,62100,29919
このデータにおいて、
cov(X) # 分散共分散行列
cor(X) # 相関係数行列
以下が、その結果。
> cov(X) # 分散共分散行列
OCU UNO PRI
OCU 1.297147e+12 185407240994 21146986138
UNO 1.854072e+11 27363809917 2614407571
PRI 2.114699e+10 2614407571 2151559501
> cor(X) # 相関係数行列
OCU UNO PRI
OCU 1.0000000 0.9841103 0.4002922
UNO 0.9841103 1.0000000 0.3407284
PRI 0.4002922 0.3407284 1.0000000
さて、相関係数行列を見てみると、
居住世帯数（OCU）と非居住世帯数（UNO）の相関が非常に高く、
平均家賃（PRI）と他の変数の関係は強く無さそうに見える。
相関の有無は、ひとまずは置いておいて、
平均家賃に注目してやると、居住世帯数が多いと賃料が高く、
同様に、非居住世帯数が多くて「も」平均賃料が高くなる傾向がある。
うん？非居住世帯が多くなるほど、月の平均賃料が高くなるというのは、
直感的に合わないような気がする。
居住・非居住の世帯数の関係が影響しているかも。
問題は、2変数間の関係に、別の変数の影響が存在している可能性があること。
では、この「潜在的」な影響は、どのようにしたら見えてくるのか？
ここで、相関係数とは別の方法で2変数の関係を考えてみる。
別の方法とは、すなわち、単回帰と呼ばれる方法である。
何のことか解らない人は、もう一度、やり直し。
たしか、単回帰は、2つの変数の関係について、
一方の変数で、もう一方の変数を「説明」するという考え方だった。
この方法を用いることで、他の変数の潜在的な影響を知ることができるかもしれない。
OCUで「説明できない部分」のPRIの値と、
OCUで「説明できない部分」のUNOの値との、
関係を「比較」すれば良いのではないか？とにかくやってみる。
手を動かした方が解りやすいので、まずは、普通に回帰分析を行なってみる。
Rでは、単回帰を行うためのlm()関数があるので、今回は、この関数を使ってみる。
# 平均家賃を固定し、単回帰を行う
X.lm.3_1 <- lm(X[,3]~X[,1]) # PRIをOCUで説明する
X.lm.2_1 <- lm(X[,2]~X[,1]) # UNOをOCUで説明する
本当は、PRIとUNOを比較したいが、OCUの影響が気になる。
したがって、PRIとUNOをそれぞれ「y軸」に固定し、OCUからの説明を考える。
# グラフを並べて描くためのオマジナイ
par(mfrow=c(1,2))
# PRIをy軸とし、OCUをx軸とした散布図
plot(X[,1], X[,3], xlab="OCU", ylab="PRI", ylim=c(0, max(X[,3])))
# 回帰直線を重ねて描く
par(new=TRUE)
plot(X[,1], X.lm.3_1$fitted.value, type="l", ylim=c(0, max(X[,3])), xaxt="n", yaxt="n", xlab="", ylab="")
# UNOをy軸とし、OCUをx軸とした散布図
plot(X[,1], X[,2], xlab="OCU", ylab="UNO", ylim=c(0, max(X[,2])))
# 回帰直線を重ねて描く
par(new=TRUE)
plot(X[,1], X.lm.2_1$fitted.value, type="l", ylim=c(0, max(X[,2])), xaxt="n", yaxt="n", xlab="", ylab="")
次に、各々の「残差の状況」を確認してみる。
# グラフを並べて描くためのオマジナイ
par(mfrow=c(1,2))
# PRIとOCUの残差の検討（標準化残差）
plot(X[,1], scale(X.lm.3_1$residuals), ylim=c(-4,4), xlab="OCU", ylab="PRI")
abline(h=0)
abline(h=c(-2,2), lty=2, col="grey")
# UNOとOCUの残差の検討（標準化残差）
plot(X[,1], scale(X.lm.2_1$residuals), ylim=c(-4,4), xlab="OCU", ylab="UNO")
abline(h=0)
abline(h=c(-2,2), lty=2, col="grey")
この図では、残差の状況を見やすくするために残差を標準化している。
これは「標準化残差」と呼ばれる。「±2」の線は、おおよその「外れ値」の基準線。
詳しい話は、正規分布の話で述べる。
さてさて、たしか、理屈上は、説明できない部分が「残差」と呼ばれるものだった。
つまり、OCUの影響を取り除くには、「残差同士の相関」を見れば良い。
# 残差同士の相関係数
cor(X.lm.3_1$residuals,X.lm.2_1$residuals)
以下は実行結果。
> # 残差同士の相関係数
> cor(X.lm.3_1$residuals,X.lm.2_1$residuals)
[1] -0.3269775
なるほど。OCUの影響を抜いた状況で相関係数を見てみると、
月の平均賃料と非居住世帯数の関係は、「負の関係」となった。
つまり、非世帯数が増加するほど、平均賃料は低くなる傾向がある。
この結果は、影響を抜かない相関係数よりも直感に即しているように思える。
すなわち、人が住んでいない住宅が多いということは、
需要に対する過剰供給が背後に存在し得るので賃料が低くなる。
試しに、残差同士の関係を散布図で可視化すると以下のようになる。
# とりあえず、変数を定義し直す
pri.res <- X.lm.3_1$residuals
uno.res <- X.lm.2_1$residuals
# 残差同士の単回帰（偏回帰）
res.lm <- lm(uno.res ~ pri.res)
# 残差同士の関係を散布図で表す
plot(pri.res, uno.res, ylim=c(min(uno.res),max(uno.res)), xlab="", ylab="")
# 一応、回帰直線も追加
par(new=TRUE)
plot(pri.res, res.lm$fitted.value, type="l", xlab="", ylab="", xaxt="n", yaxt="n", ylim=c(min(uno.res),max(uno.res)), col="red")
# 最後に図を整える
pcor <- round(cor(pri.res, uno.res),3)
plm <- round(res.lm$coefficient[2],3)
leg1 <- paste("Partial Correlation", pcor, sep=": ")
leg2 <- paste("Partial Regression Coefficient", plm, sep=": ")
legend("topright", legend=c(leg1, leg2))
title(xlab="Residuals for building rental price", ylab="Residuals for unoccupied buildings")
abline(h=0, v=0)
実際の手法については「中級レベル」で解説するとして、要するに、
これら行列の意味を理解することは極めて重要である。
そもそも、この二つの行列は、どのようなものであったか？
たしか、分散共分散行列は、各変数の分散が対角成分に並び、
それ以外の成分には「共分散」が入っているような「対称行列」であった。
一方、相関係数行列は、対角成分に「1」が並び、
それ以外の部分には「相関係数」が入っているような行列であった。
相関係数行列は、変数間の関係を一目で確認できるので、色々と重宝するのだった。
ふむふむ。そのような話を確かにした。何か「問題」でもあるのか？
実は、2変数の相関係数の話では気づかなかったのであるが、
多変数の相関係数行列にしてみると、ある種の疑惑が生まれる。
本当に、二つの変数の関係「のみ」を表しているのか？
実は、他の変数の影響を受けているのではないか？
という疑惑である。
例えば、3つの変数、「
」、「」、「」の3つの変数が存在した場合、この3つの変数の相関係数行列は以下のようになるはず。
ふむ。何も不審な点は見当たらないように見えるが、
ひょっとすると、「
」と「」の関係の中には、「
」の影響が「潜在的」に存在するかもしれない。2変数の間の相関係数の話では、意識する必要は無いが、
たしかに、多変数となる相関係数行列では、そういった懸念が生じる。
では、そのような別の変数の影響を取り除くにはどうすれば良いか？
この問題が、本日のメインテーマ。
とりあえず、いつものように、いつもの如く、チュートリアルデータを準備する。
今回のデータは、総務省統計局のデータ（e-stat）のデータを用いて、
チュートリアル用の「擬似」データを作成したもの。
今回のデータは、居住世帯の有無と一ヶ月辺りの平均家賃の関係を見るためのデータ。
家賃のデータは、19階級に区分し家賃を集計したものであったので、
各階級の平均を一ヶ月の家賃とし、世帯数を掛けたもので都道府県別の平均家賃を算出。
かなり「大雑把な加工」をしているので、正式な分析には用いてはならない。
あくまで、手法を学ぶための「擬似データ」であるので、
ここから実際の解釈を導いても、「深い意味」はない。
一応、注意事項を整理したところで、いつもどおりにデータを読み込む。
# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",", header=TRUE)
# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",", header=TRUE)
今回のチュートリアルのために用意したデータは以下の通り。
OCU,UNO,PRI
Hokkaido,2340300,390100,40687
Aomori-ken,493500,87400,111566
Iwate-ken,470700,78800,17002
Miyagi-ken,869700,144200,16084
Akita-ken,380300,57100,41254
Yamagata-ken,383000,49700,10523
Fukushima-ken,699700,108600,11879
Ibaraki-ken,1036200,187600,25207
Tochigi-ken,708700,131300,36003
Gumma-ken,725300,130400,27494
Saitama-ken,2688000,341100,22045
Chiba-ken,2344500,373100,94198
Tokyo-to,5939900,840500,60863
Kanagawa-ken,3612200,455600,246996
Niigata-ken,810700,119000,113336
Toyama-ken,368800,55500,24592
Ishikawa-ken,421600,76400,10485
Fukui-ken,259700,49000,15203
Yamanashi-ken,314600,83600,6005
Nagano-ken,758300,188000,10698
Gifu-ken,712600,123100,25608
Shizuoka-ken,1359400,238500,22521
Aichi-ken,2764400,368400,54336
Mie-ken,680900,110100,122595
Shiga-ken,491300,76300,17786
Kyoto-fu,1086800,183400,14849
Osaka-fu,3685100,660900,39764
Hyogo-ken,2169400,351300,150179
Nara-ken,502500,90100,68843
Wakayama-ken,382100,85700,12267
Tottori-ken,208600,38600,10868
Shimane-ken,249900,45900,6686
Okayama-ken,734700,131900,7459
Hiroshima-ken,1147600,208700,27600
Yamaguchi-ken,584100,107600,48448
Tokushima-ken,297000,58600,21579
Kagawa-ken,372700,73700,10305
Ehime-ken,574000,107100,12877
Kochi-ken,312800,64900,22953
Fukuoka-ken,2034000,340800,11205
Saga-ken,286100,36800,103313
Nagasaki-ken,539200,91900,9725
Kumamoto-ken,663800,105700,20633
Oita-ken,467200,79300,29546
Miyazaki-ken,443800,65900,22566
Kagoshima-ken,718200,133100,17694
Okinawa-ken,504400,62100,29919
このデータにおいて、
- OCUは居住世帯、
- UNOは非居住世帯、
- PRIは一ヶ月の平均家賃。
cov(X) # 分散共分散行列
cor(X) # 相関係数行列
以下が、その結果。
> cov(X) # 分散共分散行列
OCU UNO PRI
OCU 1.297147e+12 185407240994 21146986138
UNO 1.854072e+11 27363809917 2614407571
PRI 2.114699e+10 2614407571 2151559501
> cor(X) # 相関係数行列
OCU UNO PRI
OCU 1.0000000 0.9841103 0.4002922
UNO 0.9841103 1.0000000 0.3407284
PRI 0.4002922 0.3407284 1.0000000
居住世帯数（OCU）と非居住世帯数（UNO）の相関が非常に高く、
平均家賃（PRI）と他の変数の関係は強く無さそうに見える。
相関の有無は、ひとまずは置いておいて、
平均家賃に注目してやると、居住世帯数が多いと賃料が高く、
同様に、非居住世帯数が多くて「も」平均賃料が高くなる傾向がある。
うん？非居住世帯が多くなるほど、月の平均賃料が高くなるというのは、
直感的に合わないような気がする。
居住・非居住の世帯数の関係が影響しているかも。
問題は、2変数間の関係に、別の変数の影響が存在している可能性があること。
では、この「潜在的」な影響は、どのようにしたら見えてくるのか？
ここで、相関係数とは別の方法で2変数の関係を考えてみる。
別の方法とは、すなわち、単回帰と呼ばれる方法である。
何のことか解らない人は、もう一度、やり直し。
たしか、単回帰は、2つの変数の関係について、
一方の変数で、もう一方の変数を「説明」するという考え方だった。
この方法を用いることで、他の変数の潜在的な影響を知ることができるかもしれない。
OCUで「説明できない部分」のPRIの値と、
OCUで「説明できない部分」のUNOの値との、
関係を「比較」すれば良いのではないか？とにかくやってみる。
手を動かした方が解りやすいので、まずは、普通に回帰分析を行なってみる。
Rでは、単回帰を行うためのlm()関数があるので、今回は、この関数を使ってみる。
# 平均家賃を固定し、単回帰を行う
X.lm.3_1 <- lm(X[,3]~X[,1]) # PRIをOCUで説明する
X.lm.2_1 <- lm(X[,2]~X[,1]) # UNOをOCUで説明する
本当は、PRIとUNOを比較したいが、OCUの影響が気になる。
したがって、PRIとUNOをそれぞれ「y軸」に固定し、OCUからの説明を考える。
par(mfrow=c(1,2))
# PRIをy軸とし、OCUをx軸とした散布図
plot(X[,1], X[,3], xlab="OCU", ylab="PRI", ylim=c(0, max(X[,3])))
# 回帰直線を重ねて描く
par(new=TRUE)
plot(X[,1], X.lm.3_1$fitted.value, type="l", ylim=c(0, max(X[,3])), xaxt="n", yaxt="n", xlab="", ylab="")
# UNOをy軸とし、OCUをx軸とした散布図
plot(X[,1], X[,2], xlab="OCU", ylab="UNO", ylim=c(0, max(X[,2])))
# 回帰直線を重ねて描く
par(new=TRUE)
plot(X[,1], X.lm.2_1$fitted.value, type="l", ylim=c(0, max(X[,2])), xaxt="n", yaxt="n", xlab="", ylab="")
par(mfrow=c(1,2))
# PRIとOCUの残差の検討（標準化残差）
plot(X[,1], scale(X.lm.3_1$residuals), ylim=c(-4,4), xlab="OCU", ylab="PRI")
abline(h=0)
abline(h=c(-2,2), lty=2, col="grey")
# UNOとOCUの残差の検討（標準化残差）
plot(X[,1], scale(X.lm.2_1$residuals), ylim=c(-4,4), xlab="OCU", ylab="UNO")
abline(h=0)
abline(h=c(-2,2), lty=2, col="grey")
この図では、残差の状況を見やすくするために残差を標準化している。
これは「標準化残差」と呼ばれる。「±2」の線は、おおよその「外れ値」の基準線。
詳しい話は、正規分布の話で述べる。
さてさて、たしか、理屈上は、説明できない部分が「残差」と呼ばれるものだった。
つまり、OCUの影響を取り除くには、「残差同士の相関」を見れば良い。
# 残差同士の相関係数
cor(X.lm.3_1$residuals,X.lm.2_1$residuals)
以下は実行結果。
> # 残差同士の相関係数
> cor(X.lm.3_1$residuals,X.lm.2_1$residuals)
[1] -0.3269775
月の平均賃料と非居住世帯数の関係は、「負の関係」となった。
つまり、非世帯数が増加するほど、平均賃料は低くなる傾向がある。
この結果は、影響を抜かない相関係数よりも直感に即しているように思える。
すなわち、人が住んでいない住宅が多いということは、
需要に対する過剰供給が背後に存在し得るので賃料が低くなる。
試しに、残差同士の関係を散布図で可視化すると以下のようになる。
pri.res <- X.lm.3_1$residuals
uno.res <- X.lm.2_1$residuals
# 残差同士の単回帰（偏回帰）
res.lm <- lm(uno.res ~ pri.res)
# 残差同士の関係を散布図で表す
plot(pri.res, uno.res, ylim=c(min(uno.res),max(uno.res)), xlab="", ylab="")
# 一応、回帰直線も追加
par(new=TRUE)
plot(pri.res, res.lm$fitted.value, type="l", xlab="", ylab="", xaxt="n", yaxt="n", ylim=c(min(uno.res),max(uno.res)), col="red")
# 最後に図を整える
pcor <- round(cor(pri.res, uno.res),3)
plm <- round(res.lm$coefficient[2],3)
leg1 <- paste("Partial Correlation", pcor, sep=": ")
leg2 <- paste("Partial Regression Coefficient", plm, sep=": ")
legend("topright", legend=c(leg1, leg2))
title(xlab="Residuals for building rental price", ylab="Residuals for unoccupied buildings")
abline(h=0, v=0)
さて、このように、ある2つの変数の関係を観察するために、
別の変数の影響を取り除いた相関係数のことを、
「偏相関係数（Partial Correlation）」と呼ぶ。
また、ここでは、残差同士の単回帰を行い、その回帰直線も描いているが、
残差同士の単回帰の係数のことを「偏回帰係数（Partial Regression Coefficients）」と呼ぶ。
ちなみに、残差同士の関係なので切片は「0」となる。
さて、今回は、単回帰の残差の関係から偏相関を導いたが、
一般的に、偏相関係数は元の相関係数から、以下のように計算できる。
なお、
は、の影響を抜いた と との偏相関を意味し、
、
、
はそれぞれ、元の相関係数を意味している。
「
」は、ギリシア文字で「ロー」と発音する。「ピー」では無い。
さて、要するに、知りたい2つの変数の相関係数から、
潜在的に影響しているかもしれない別の1つの変数の影響を取り除き、
結果の値が、「±1」の間に収まるように分母の部分で基準化している。
今回のように、3変数の場合は、直感的に理解しやすいのであるが、
4変数以上になると計算が非常に厄介になってくる。行列を使えば、
式自体はコンパクトになるが...余裕の無い人は以下の式の話は飛ばしても良い。
まず、元の相関係数行列における、変数「i」と「j」との相関係数を「
」 とし、
その相関係数行列の逆行列における「i」と「j」の関係を 「
」する。
ちなみに、逆行列の添字を右下ではなく「右上」に置く。ちょっと、注意。
このとき、4変数以上の偏相関係数行列は、次のように求める。
逆行列の対応する要素を、2つの対角要素の平方根で割って基準化し、
さらに、「符号を反転」させている。逆行列を使った数学の魔法!!
この方法で、偏相関係数を計算することができるのである。
まぁ、実際には、手計算で偏相関係数を求める人は少ないであろう。
Rでは、この偏相関係数を簡単に計算できるため、通常はこれを用いるのであるが、
残念ながら、偏相関行列を計算してくれるパッケージは最初からは入っていないので、
「corpcor」パッケージをインストールする。
# corpcor パッケージのインストール
install.packages("corpcor", dep=TRUE)
上手くインストールできたら、実際に計算してみる。
# ライブラリを読み込む
library(corpcor)
# 偏相関係数行列を計算する。
X.pcor <- cor2pcor(cor(X))
# 変数名が付かないので付け直す。
rownames(X.pcor) <- rownames(cor(X))
colnames(X.pcor) <- colnames(cor(X))
以下が、この実行結果。
> X.pcor
OCU UNO PRI
OCU 1.0000000 0.9839439 0.3892439
UNO 0.9839439 1.0000000 -0.3269775
PRI 0.3892439 -0.3269775 1.0000000
PRIとUNOの関係を見てみると、確かに、残差同士の相関係数と一致している。
相関係数行列は、変数間の関係を観察する上で非常に重要であるが、
さらに、細かい状況を把握するためには、
偏相関係数行列で見た方が、状況が良く解る場合もある。
別の変数の影響を取り除いた相関係数のことを、
「偏相関係数（Partial Correlation）」と呼ぶ。
また、ここでは、残差同士の単回帰を行い、その回帰直線も描いているが、
残差同士の単回帰の係数のことを「偏回帰係数（Partial Regression Coefficients）」と呼ぶ。
ちなみに、残差同士の関係なので切片は「0」となる。
さて、今回は、単回帰の残差の関係から偏相関を導いたが、
一般的に、偏相関係数は元の相関係数から、以下のように計算できる。
なお、
は、の影響を抜いた と との偏相関を意味し、、、 はそれぞれ、元の相関係数を意味している。「
」は、ギリシア文字で「ロー」と発音する。「ピー」では無い。さて、要するに、知りたい2つの変数の相関係数から、
潜在的に影響しているかもしれない別の1つの変数の影響を取り除き、
結果の値が、「±1」の間に収まるように分母の部分で基準化している。
今回のように、3変数の場合は、直感的に理解しやすいのであるが、
4変数以上になると計算が非常に厄介になってくる。行列を使えば、
式自体はコンパクトになるが...余裕の無い人は以下の式の話は飛ばしても良い。
まず、元の相関係数行列における、変数「i」と「j」との相関係数を「
」 とし、その相関係数行列の逆行列における「i」と「j」の関係を 「
」する。ちなみに、逆行列の添字を右下ではなく「右上」に置く。ちょっと、注意。
このとき、4変数以上の偏相関係数行列は、次のように求める。
逆行列の対応する要素を、2つの対角要素の平方根で割って基準化し、
さらに、「符号を反転」させている。逆行列を使った数学の魔法!!
この方法で、偏相関係数を計算することができるのである。
まぁ、実際には、手計算で偏相関係数を求める人は少ないであろう。
Rでは、この偏相関係数を簡単に計算できるため、通常はこれを用いるのであるが、
残念ながら、偏相関行列を計算してくれるパッケージは最初からは入っていないので、
「corpcor」パッケージをインストールする。
# corpcor パッケージのインストール
install.packages("corpcor", dep=TRUE)
上手くインストールできたら、実際に計算してみる。
# ライブラリを読み込む
library(corpcor)
# 偏相関係数行列を計算する。
X.pcor <- cor2pcor(cor(X))
# 変数名が付かないので付け直す。
rownames(X.pcor) <- rownames(cor(X))
colnames(X.pcor) <- colnames(cor(X))
以下が、この実行結果。
> X.pcor
OCU UNO PRI
OCU 1.0000000 0.9839439 0.3892439
UNO 0.9839439 1.0000000 -0.3269775
PRI 0.3892439 -0.3269775 1.0000000
相関係数行列は、変数間の関係を観察する上で非常に重要であるが、
さらに、細かい状況を把握するためには、
偏相関係数行列で見た方が、状況が良く解る場合もある。
分散共分散行列と相関係数行列はを→分散共分散行列と相関係数行列を だと思います。
返信削除あっ、そうですね。直します。
削除以下のxlab="OCU"は正しくは"UNO"ではないでしょうか？
返信削除# UNOとOCUの残差の検討（標準化残差）
plot(X[,1], scale(X.lm.2_1$residuals), ylim=c(-4,4), xlab="OCU", ylab="PRI")
abline(h=0)
abline(h=c(-2,2), lty=2, col="grey")
間違えました。ylab="PRI"が正しくはylab="UNO"ではないでしょうか？
返信削除毎度、精査してくださりありがとうございます。ご指摘の通りです。
削除とりあえず、コードの方は直しておきました。図は、後ほど改めて。