Lesson4.Regression!
Kosugitti
2017年11月10日
毎回やる儀式
- インターネットに接続しておいてください（このサイトが見えていればOK）。
- Rstudioを起動し，この授業のプロジェクトを開いてください。(File>Open Project)
- File>New File>R scriptで新しいコードファイルを作ってください。
- 野球のデータセットがbaseballオブジェクトに入っているものとします。
baseball <- read.csv("ball2017.csv",na.strings="*",fileEncoding = "UTF-8")回帰分析はすごいよ
今日のモデルのRQ例
- 野球選手の成績と年収はどう関係しているのだろうか？
- 野球選手の成績と年収は，基本的に関係するんだろうけど，チームごとの違いが結構大きいよね？（階層線形モデル）
- 心理療法の研究で，抑うつに対して相対的に効果があったのは認知行動療法と薬物療法，どっち？（実験計画，メタ分析）
- 収入やクレジットカードの履歴で，クレジットリスクの低い人と高い人を区別できるだろうか？（判別分析）
- ある入院患者が心筋梗塞と診断されました。血圧の上下，心拍数，脈拍数，動脈圧などのデータがあったとして，医療スタッフがこのデータから患者の生存可能性を推定できないだろうか？（生存分析）
- 心理学者が無意味綴りを覚えさせる認知実験を考える。被験者は三つのグループに分けられ，刺激提示が異なる方法で行われる。何十回と刺激を受けて，どれだけ記憶に残っているかが検討されるが，群の分け方（刺激の与え方）で結果は違うだろうか？この時の刺激は同じものもあり違うものもある。刺激の「無意味さの程度」は整えてあるつもりだけど・・・（階層線形モデル，あるいは反復測定分散分析）
- ダイエットなう。この三ヶ月ほど体重の増減を毎朝記録しているけど，明日の体重はどうなるかな？（時系列分析，自己相関の高い回帰分析）
- 風が吹けば桶屋が儲かるらしい。というのもつまり、大風で土ぼこりが立つ→土ぼこりが目に入って、盲人が増える→盲人は三味線を買う（当時の盲人が就ける職に由来）→三味線に使う猫皮が必要になり、ネコが殺される→ネコが減ればネズミが増える→ネズミは桶をかじる→桶の需要が増え桶屋が儲かる，らしい。本当かね？（パス解析）
え，こんなにも色々なことができるの？と思われたかもしれません。ちょっと大げさに表現していますが，これらはいずれも今日の「回帰分析」を応用したRQです。RQのうしろ，カッコ内の分析名がより正しいモデル名なのですが，いずれも「回帰分析」の応用系なのです。
今日のモデルのRQ例解説
- ヒットを打つひとは年俸が高い？ホームランを打つ人は年俸が高い？こうした関係は，グラフを書いてみたり，相関係数を見たりすると検証できますね。しかしそこに「モデル」の観点を入れると，「ヒット一本うつと年俸がいくら違うか？」という話ができるようになります！
- 基本的には連続変数と連続変数の関係にモデル＝数式を当てはめることで，この目的を達成します。
- データがネストしていたり，カテゴリカル変数（連続変数でない）であったりすると，モデルは様々に展開していきます（上記RQのカッコ内）。
- モデルを繰り返し適用することで，経路（パス）を考えることもできます。風が吹けば桶屋が儲かる論の検証ですね。
- 前回のMDSの結果と合わせて考えることで，次元 を説明する変数はなにか，というのをデータから考えることもできます。
回帰分析ってなんなの？
回帰分析とは，ある目的となる変数を説明する関数を当てはめることです。一般に，\[y=ax+b\] という線形関数を当てはめます。
身長と体重を例にして説明しましょう。 この二つの変数の散布図を書いてみます。
library(ggplot2)
g <- ggplot(baseball,aes(x=height,y=weight)) + geom_point()
plot(g)
右肩あがりの関係がありそうです。線を入れてみます。
g <- g + geom_smooth(method="lm",se=FALSE)
plot(g)
この線の引き方が，線形回帰をするということです。回帰モデルをデータに当てはめる，とも言います。 lmは，liner model(線形モデル)という意味です。
この線形モデルを数式で表現するのが次の方法です。
reg <- lm(weight~height,data=baseball)
reg##
## Call:
## lm(formula = weight ~ height, data = baseball)
##
## Coefficients:
## (Intercept) height
## -160.906 1.377すなわち，体重＝1.377* 身長 -160.906 という関係にあるということがわかります。言い換えると，身長が1cm高いと体重が1.37kg大きいということです(このデータの中では）。また。Intercept(切片)は身長が0だった時の(!)値です。
もちろん，すべてのデータがこの直線上に乗っているわけではありません。むしろほとんど外れていますが，この線が「当てはまっている」といえるのはどうしてでしょうか？ここでは，最小二乗基準 によってこの数字を決めています。最小二乗法による回帰分析，ということもあります。 推定値の決め方は他にもある，ということを頭の片隅に置いておいてください。
さて，最小二乗基準で推定値が得られはしましたが，誤差も多いので，これをどう評価すればよいでしょうか。 ちなみに，推定値や誤差は結果のオブジェクトに入っています。
reg$residuals #誤差。残差ともいう。## 1 2 3 4 5 6
## 3.4372215 15.8287594 13.0456836 9.3214092 -4.0701287 14.9592721
## 7 8 9 10 11 12
## -11.0701287 6.1908966 14.1908966 4.0750843 15.5530338 6.1908966
## 13 14 15 16 17 18
## -7.0701287 6.2920085 4.1908966 11.9445717 2.0603840 10.0750843
## 19 20 21 22 23 24
## -7.5480781 -3.9249157 -1.9249157 -1.9249157 7.2055970 -2.1712406
## 25 26 27 28 29 30
## -2.5333778 9.7129472 -10.5480781 -0.6638904 2.2055970 -12.9249157
## 31 32 33 34 35 36
## 0.0897847 -6.3017532 -12.1859409 -10.9249157 -8.5627785 -7.4322659
## 37 38 39 40 41 42
## 5.2055970 -1.1712406 1.0750843 2.6835464 -1.9249157 -13.5480781
## 43 44 45 46 47 48
## -9.5480781 7.2055970 -13.0554283 -11.9249157 0.3214092 -4.9249157
## 49 50 51 52 53 54
## 9.2055970 -0.1712406 -7.8238038 -7.4175655 -4.7944030 -13.7944030
## 55 56 57 58 59 60
## 4.6982468 2.9592721 21.9592721 1.5824345 12.9445717 3.0750843
## 61 62 63 64 65 66
## -8.4175655 -1.9249157 -8.9249157 1.3361096 -3.4175655 10.4666222
## 67 68 69 70 71 72
## -3.5480781 3.1908966 -5.0407279 -3.8091034 3.3214092 5.2055970
## 73 74 75 76 77 78
## -14.0407279 -4.9396160 4.3361096 -2.0554283 1.6982468 -1.0407279reg$fitted.values #推定値## 1 2 3 4 5 6 7
## 96.56278 84.17124 108.95432 89.67859 102.07013 80.04073 102.07013
## 8 9 10 11 12 13 14
## 93.80910 93.80910 86.92492 103.44697 93.80910 102.07013 111.70799
## 15 16 17 18 19 20 21
## 93.80910 91.05543 97.93962 86.92492 85.54808 86.92492 86.92492
## 22 23 24 25 26 27 28
## 86.92492 82.79440 84.17124 74.53338 77.28705 85.54808 78.66389
## 29 30 31 32 33 34 35
## 82.79440 86.92492 75.91022 88.30175 95.18594 86.92492 96.56278
## 36 37 38 39 40 41 42
## 92.43227 82.79440 84.17124 86.92492 99.31645 86.92492 85.54808
## 43 44 45 46 47 48 49
## 85.54808 82.79440 91.05543 86.92492 89.67859 86.92492 82.79440
## 50 51 52 53 54 55 56
## 84.17124 104.82380 81.41757 82.79440 82.79440 88.30175 80.04073
## 57 58 59 60 61 62 63
## 80.04073 81.41757 91.05543 86.92492 81.41757 86.92492 86.92492
## 64 65 66 67 68 69 70
## 78.66389 81.41757 74.53338 85.54808 93.80910 80.04073 93.80910
## 71 72 73 74 75 76 77
## 89.67859 82.79440 80.04073 97.93962 78.66389 91.05543 88.30175
## 78
## 80.04073この推定値と実際のデータとの相関係数が当てはまりの良さ(適合度)の一つの指標です。相関係数が高ければ高いほど，良い予測ができていると言えます。 相関係数は-1から1までの値を取りますが，適合度指標としては相関係数の二乗を使います(R-squared)。 この指標も実は結果のオブジェクトに既に入っています。
summary(reg)##
## Call:
## lm(formula = weight ~ height, data = baseball)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.0407 -5.9865 -0.4176 5.0788 21.9593
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -160.9058 30.0502 -5.355 8.81e-07 ***
## height 1.3768 0.1662 8.285 3.13e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.247 on 76 degrees of freedom
## Multiple R-squared: 0.4746, Adjusted R-squared: 0.4677
## F-statistic: 68.65 on 1 and 76 DF, p-value: 3.128e-120.46ほどですね。これは実際のデータとしては，とても良い当てはまりだと言えます。
うんちく
- summary関数の出力において，推定値がどこにあるか確認してください。
- アスタリスクマークが出ていますが，これはこの回帰係数が「統計的に有意である」ことを示しています。
- R-squaredは\(R^2\)とも書き，説明される分散の大きさを表しています。
- F統計量が出ています。これは回帰分析が分散分析と数学的には同じであることの表れです！
重回帰分析
説明する変数が増えたら重回帰分析(Multiple Regression Analysis)と名前が変わります。 モデルの書き方はそれほど変わりません。
reg2 <- lm(HR~height+weight,data=baseball)
reg2##
## Call:
## lm(formula = HR ~ height + weight, data = baseball)
##
## Coefficients:
## (Intercept) height weight
## -16.1248 -0.1603 0.6568summary(reg2)##
## Call:
## lm(formula = HR ~ height + weight, data = baseball)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.614 -5.449 -2.470 5.274 19.062
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -16.1248 34.0112 -0.474 0.637
## height -0.1603 0.2211 -0.725 0.471
## weight 0.6568 0.1106 5.937 8.42e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.954 on 75 degrees of freedom
## Multiple R-squared: 0.4309, Adjusted R-squared: 0.4158
## F-statistic: 28.4 on 2 and 75 DF, p-value: 6.584e-10このことから，(同じ身長の人であれば)体重が0.65kg増えるとホームランが一本増えることになります。身長はHRとは関係なかったようですね！
説明変数が複数になると，解釈はちょっと注が必要です。従属変数に対して個別に影響しているのではなく，他の説明変数を統制した時の影響力の大きさである，と解釈する必要があることに注意してください。このことについては特に豊田秀樹編著「もう一つの重回帰分析」東京図書で言及されていますので，副読本として紹介しておきます。
さて，説明変数に年収を入れて見ましょう。
reg2.2 <- lm(HR~height+weight+pay,data=baseball)
reg2.2##
## Call:
## lm(formula = HR ~ height + weight + pay, data = baseball)
##
## Coefficients:
## (Intercept) height weight pay
## -6.5900043 -0.1950482 0.5697236 0.0002792summary(reg2.2)##
## Call:
## lm(formula = HR ~ height + weight + pay, data = baseball)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.979 -4.920 -1.905 4.408 16.498
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.590e+00 3.189e+01 -0.207 0.83686
## height -1.950e-01 2.068e-01 -0.943 0.34862
## weight 5.697e-01 1.064e-01 5.356 9.21e-07 ***
## pay 2.792e-04 8.078e-05 3.457 0.00091 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.43 on 74 degrees of freedom
## Multiple R-squared: 0.51, Adjusted R-squared: 0.4902
## F-statistic: 25.68 on 3 and 74 DF, p-value: 1.721e-11年収(pay)も統計的に有意な説明変数のようですが，係数が0.00002792(2.792e-04)ととても小さな値になっています。この小さい値は，影響が小さいことを意味するのではなく，payが1万円(元のデータの単位が「万円」)上がると0.00002792本ホームランが増える，ということを意味しているのですが，ちょっとわかりにくいですね。
また，0.00002792と体重の係数0.5697ではどちらが大きいか？という比較も意味をなしません。単位が違うからです。 そこで，単位にとらわれずに重要度を比較できるようにするために，データの標準化を考えます。
baseball2 <- subset(baseball,select=c("HR","weight","height","pay")) #データの一部抜き出し
baseball2z <- scale(baseball2) #標準化
baseball2z <- data.frame(baseball2z) #データフレーム型にしておく
summary(baseball2z)## HR weight height pay
## Min. :-1.2173 Min. :-1.9418 Min. :-1.7228 Min. :-1.0212
## 1st Qu.:-0.9290 1st Qu.:-0.7696 1st Qu.:-0.6619 1st Qu.:-0.7809
## Median :-0.3524 Median :-0.1724 Median :-0.1315 Median :-0.4248
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.6807 3rd Qu.: 0.6238 3rd Qu.: 0.5316 3rd Qu.: 0.5878
## Max. : 2.1463 Max. : 3.0125 Max. : 3.0512 Max. : 3.0469こうした上で，あらためて。
reg2.3 <- lm(HR~height+weight+pay,data=baseball2z)
reg2.3##
## Call:
## lm(formula = HR ~ height + weight + pay, data = baseball2z)
##
## Coefficients:
## (Intercept) height weight pay
## -2.268e-16 -1.060e-01 6.189e-01 3.014e-01summary(reg2.3)##
## Call:
## lm(formula = HR ~ height + weight + pay, data = baseball2z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6317 -0.4729 -0.1831 0.4236 1.5855
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.268e-16 8.085e-02 0.000 1.00000
## height -1.060e-01 1.124e-01 -0.943 0.34862
## weight 6.189e-01 1.155e-01 5.356 9.21e-07 ***
## pay 3.014e-01 8.721e-02 3.457 0.00091 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.714 on 74 degrees of freedom
## Multiple R-squared: 0.51, Adjusted R-squared: 0.4902
## F-statistic: 25.68 on 3 and 74 DF, p-value: 1.721e-11こうすると単位が整っていますから，大きさを比較することができます。体重の係数は0.6189，年収の係数は0.3014ですから，体重の方が説明する力（影響力）は大きいと考えることができます。
うんちく
- 有意性検定，\(R^2\)は標準化したもの(reg2.3)としていないもの(reg2.2)とで変わりがないことを確認してください。
- 標準化したデータの係数(reg3)のことを標準化係数，していない生データの係数(reg2)のことを非標準化係数と言います。
- 説明変数が一つの回帰分析(reg1)の係数は回帰係数regression coefficientと言いますが，説明変数が二つ以上ある重回帰分析の係数は偏回帰係数partial regression coefficientと言います。
- partialとは，変数の影響を統計的に統制する作業のことを言います。
MDSと回帰分析の面白い関係
前回はMDSで選手をプロットしたりしたのでした。ちょっと再現して見ましょう。
身長と体重だけのサブセットを作って選手の散布図(!)を書いて見ます。
subdat <- subset(baseball,select=c("height","weight")) #二変数だけ抽出
subdat.dist <- dist(subdat) # 距離行列を計算
player.map <- cmdscale(subdat.dist,2) # MDSで座標を得る
plot(player.map,type="n")
text(player.map,levels(baseball$name))
たくさんの選手がいて，野球に詳しくない人には特徴がわかりにくいかもしれません。 しかし，MDSによって一人一人に座標＝連続的な数字が割り当てられたのですから，その座標値を予測する変数があれば，次元の意味がわかるかもしれません！
baseball$dim1 <- player.map[,1]
reg4 <- lm(dim1~HR, data=baseball)
summary(reg4)##
## Call:
## lm(formula = dim1 ~ HR, data = baseball)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.573 -3.629 1.138 5.721 19.802
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.3876 1.6662 5.634 2.84e-07 ***
## HR -0.7411 0.1019 -7.273 2.68e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.305 on 76 degrees of freedom
## Multiple R-squared: 0.4103, Adjusted R-squared: 0.4026
## F-statistic: 52.89 on 1 and 76 DF, p-value: 2.682e-10これはホームランの本数を独立変数にして，座標を従属変数にした時の回帰係数です。-0.7なのでホームランが増えるほど左のほうにプロットされるようです。ホームラン次元なのかもしれませんねえ。
本日の課題
- 野球のデータを使って，本塁打数(HR)を打率(AB)で単回帰分析してください。回帰係数と自由度調節済み\(R^2\)がいくらだったか報告すること。
- 次の図を描いてください。
- （今後のために）次の図を描いてください。すなわち，チームごとに線を引くのです。
- 野球のデータを使って，本塁打数(HR)を従属変数，打率(AB)，三振の数(K)を独立変数とした重回帰分析を行なってください。その際，標準化されてない回帰係数と，標準化された回帰係数，自由度調節済み\(R^2\)を報告すること。
- 野球のデータを使って，年収を一番よく説明する変数は何か，重回帰分析で求めてください。