この間までは、相関係数は意外と扱いが難しいよーという話をしました。
相関係数の目安と解釈と注意点~舐めてかかると痛い目を見る話~ - バナナでもわかる話
最小二乗法による相関係数の導出~相関係数の意外な解釈法 - バナナでもわかる話
今回は、因果に関する相関分析を具体例を使ってやってみようと思うのですが、
共立出版の「マーケティング・データ分析の基礎」
の第8章の話を軸に、本の中で考慮不足な点を補いながら解説していこうかと思います。
データはここ（https://www.kyoritsu-pub.co.jp/bookdetail/9784320123663）から取ってきています。
データの説明
まず、データですが
x=read.csv("regression_data.csv",header=TRUE,comment.char="#") head(x) city sake TaJan 1 新潟市 15259 2.8 2 金沢市 13717 3.8 3 福島市 13245 1.6 4 松江市 12716 4.3 5 秋田市 12684 0.1 6 長野市 12511 -0.6
sakeが清酒の消費量、TaJanが平均気温ですね。清酒の消費量はここから取ってきています。
統計局ホームページ/家計調査（二人以上の世帯） 品目別都道府県庁所在市及び政令指定都市ランキング（平成27年（2015年）〜29年（2017年）平均）
プロットするとこんな感じ
attach(x) library(maptools) plot(TaJan,sake,xlab="1月の平均気温",ylab="清酒の消費量ml") pointLabel(TaJan,sake,labels=city,cex=0.7)
本書の主張
で、ここでこの教科書だと、プロットしてみた感じだと平均気温が低いほど清酒の消費量が多そうだということで、相関係数を出していて
下の結果から分かるように相関係数は-0.68、相関係数が0かどうかの仮説検定を行ったら、帰無仮説も棄却され、めでたく相関はありそうです！
> cor(TaJan,sake) [1] -0.6845565 > cor.test(TaJan,sake) Pearson's product-moment correlation data: TaJan and sake t = -6.2996, df = 45, p-value = 1.12e-07 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.8120826 -0.4946238 sample estimates: cor -0.6845565
という話でした。
ここまでが教科書での相関分析の話。後はこれに回帰分析を適用しています。
本当か？
しかし、この分析は納得のいく結果でしょうか。
前回前々回の記事を思い出していただくと、因果を考える場合の相関係数は少なくとも次の点に注意する必要がありました。
・線形の仮定があること
・外れ値の影響は気を付ける必要があること
プロットを見る限り、線形の仮定はまあ大よそ良さそうですが、右下の那覇市、少し気になりませんか？
正直那覇市を取り除くか否かで結構相関係数の値が変わりそうだし、平均気温15度近辺は那覇市のデータしかないわけなので、那覇市に過適合していそうにも見えます。
このケースの外れ値の判別はクックの距離という指標があって、距離が大きいほど外れ値の指標になります。
Rでは次のように求めることが出来ます。
> cooks.distance(lm(TaJan~sake)) 1 2 3 4 5 1.816426e-01 9.817568e-02 4.112882e-03 5.848597e-02 8.984344e-03 6 7 8 9 10 2.426579e-02 8.071903e-03 5.929246e-03 2.338768e-02 5.472206e-03 11 12 13 14 15 6.706070e-02 9.477868e-07 3.262687e-02 4.933207e-02 8.683578e-05 16 17 18 19 20 1.090427e-02 6.473217e-04 4.274496e-03 1.120386e-02 2.563119e-03 21 22 23 24 25 9.809439e-05 3.050151e-03 1.242908e-01 8.352329e-04 2.789172e-03 26 27 28 29 30 1.107317e-02 1.336510e-03 9.814185e-04 5.739481e-03 1.088626e-04 31 32 33 34 35 1.232378e-02 4.506737e-04 9.236139e-04 2.251216e-03 2.463987e-03 36 37 38 39 40 1.031285e-04 1.529464e-06 7.187894e-04 4.674627e-03 1.492860e-03 41 42 43 44 45 2.624656e-02 7.554864e-04 1.273579e-04 8.424634e-03 6.364789e-03 46 47 3.593068e-02 8.781210e-01
47番のクックの距離が0.88を取っているのが那覇市です。他の市よりも大きな値を取っていることが分かります。
つまり、クックの距離を基準にすれば、那覇市は明らかに外れ値なわけです。相関係数は最小二乗法により求められるので、外れ値には弱い。
取り除かねばなりません。
ただ、基準では外れ値と判別されるわけですが、何故外れ値なのでしょうか。この辺をはっきりさせずに「外れ値だから～～～」と取り除くのは少々危険です。
お酒の知識と日本の文化圏の知識が不足しているので断言は出来ませんが、日本国の中でも沖縄は中々独特な文化圏を持っていますよね。成り立ちもかなり異なりますし、その辺の影響から、他の市とは違う動きをしていると考えられます。また歴史が異なるという話で行くと、北海道も気になりますが、確かにプロットを見てみると札幌も変な位置にありますよね。
この辺を考えると、少なくとも那覇市は外れ値と判断するのは妥当かなと思います。
※本来のしっかりした分析であれば、この辺の文化に関するもっと詳細な検討を加えてから外れ値とする必要があります。
さて、外した結果を見てみましょう。
xx=x[-47,] > cor(xx$TaJan,xx$sake) [1] -0.6168623 > > > cor.test(xx$TaJan,xx$sake) Pearson's product-moment correlation data: xx$TaJan and xx$sake t = -5.1988, df = 44, p-value = 4.988e-06 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.7693838 -0.3977998 sample estimates: cor -0.6168623
相関係数は-0.617まで減りました。
回帰分析
教科書では、この後相関だけでは、影響度がわからないので回帰分析もしましょうということで回帰分析を行っています。
やり方は下のような感じです。
res=lm(sake~TaJan) plot(TaJan,sake,xlab="1月の平均気温",ylab="清酒の消費量ml") abline(res$coef[1],res$coef[2])
でも、よく見てください。右下の那覇市に過適合した結果の直線に見えませんか？那覇市の場所に合わせてピッタリ直線が引かれていますよね。
ちなみに回帰分析の結果はこんな感じ
> summary(res) Call: lm(formula = sake ~ TaJan) Residuals: Min 1Q Median 3Q Max -5058.4 -1308.1 -202.1 896.1 5760.7 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 11228.6 546.7 20.54 < 2e-16 *** TaJan -618.0 98.1 -6.30 1.12e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2199 on 45 degrees of freedom Multiple R-squared: 0.4686, Adjusted R-squared: 0.4568 F-statistic: 39.68 on 1 and 45 DF, p-value: 1.12e-07
単回帰のわりに決定係数も0.47で意外と良さそうに見えますが、これって正直、私は那覇市にひっぱられてうまく直線が引けた結果にしか見えないわけです。
そこで那覇市を取って回帰分析をやり直してみました。
> res2=lm(xx$sake~xx$TaJan) > summary(lm(xx$sake~xx$TaJan)) Call: lm(formula = xx$sake ~ xx$TaJan) Residuals: Min 1Q Median 3Q Max -5102.7 -1301.9 -215.3 885.3 5755.5 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 11251.0 605.6 18.580 < 2e-16 *** xx$TaJan -624.1 120.0 -5.199 4.99e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2223 on 44 degrees of freedom Multiple R-squared: 0.3805, Adjusted R-squared: 0.3664 F-statistic: 27.03 on 1 and 44 DF, p-value: 4.988e-06
決定係数は0.38まで落ちました。一個変数を取り除いただけなのに、この落ち方は流石に大きすぎます。
つまり、先ほどの決定係数の高さは結局那覇市のデータに過適合した結果だと考えられますね。
相関係数を見る場合でも回帰分析を行う場合でも同じですが、外れ値に引っ張られている影響だとか、
決定係数を見るのであれば説明変数が多すぎることによる影響は厳しく考慮すべきだと思っています。
そして、今回の結果を解釈するのであれば、確かに緩く相関関係はありそうだが、沖縄や北海道の外れ値からも分かる通り、地域差がありそうなので、地域差を考慮したモデル(空間計量等)を用いた方がより、納得のいく結果が得られそうだと言えます。
このように、簡単なモデルで色々試してから、じゃあ次にこういう結果が見えそうだから、もう少し複雑なモデルを使ってみようというサイクルが個人的には大切だと思っていて、外れ値も考慮せず相関係数が出た！相関がある！終わり！はマズイのかなと思っております。
