Lesson4.Regression!
Kosugitti
2016年11月04日
毎回やる儀式
- インターネットに接続しておいてください（このサイトが見えていればOK）。
- Rstudioを起動し，この授業のプロジェクトを開いてください。(File>Open Project)
- File>New File>R scriptで新しいコードファイルを作ってください。
- 野球のデータセットがbaseballオブジェクトに入っているものとします。
baseball <- read.csv("baseball2016.csv",na.strings="*")回帰分析はすごいよ
今日のモデルの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")
plot(g)
この線の引き方が，線形回帰をするということです。回帰モデルをデータに当てはめる，とも言います。 lmは，liner model(線形モデル)という意味です。
この線形モデルを数式で表現するのが次の方法です。
reg <- lm(weight~height,data=baseball)
reg##
## Call:
## lm(formula = weight ~ height, data = baseball)
##
## Coefficients:
## (Intercept) height
## -125.518 1.172すなわち，体重＝1.172 * 身長 -125.5118 という関係にあるということがわかります。言い換えると，身長が1cm高いと体重が1.17kg大きいということです(このデータの中では）。また。Intercept(切片)は身長が0だった時の(!)値です。
もちろん，すべてのデータがこの直線上に乗っているわけではありません。むしろほとんど外れていますが，この線が「当てはまっている」といえるのはどうしてでしょうか？ここでは，最小二乗基準 によってこの数字を決めています。最小二乗法による回帰分析，ということもあります。 推定値の決め方は他にもある，ということを頭の片隅に置いておいてください。
さて，最小二乗基準で推定値が得られはしましたが，誤差も多いので，これをどう評価すればよいでしょうか。 ちなみに，推定値や誤差は結果のオブジェクトに入っています。
reg$residuals #誤差。残差ともいう。## 1 2 3 4 5
## 5.52516074 7.55718176 17.80512571 3.90118876 2.80512571
## 6 7 8 9 10
## 6.04117125 9.18115373 15.41719927 4.04117125 8.69716424
## 11 12 13 14 15
## 16.69716424 -1.33485678 11.46111870 8.69716424 5.35315724
## 16 17 18 19 20
## 6.69716424 14.04117125 11.55718176 -5.09881124 -2.44281824
## 21 22 23 24 25
## 7.90118876 -0.44281824 -0.44281824 8.07319227 -1.09881124
## 26 27 28 29 30
## -5.89478671 3.04117125 5.76120628 -0.41079722 3.07319227
## 31 32 33 34 35
## -11.44281824 -1.89478671 3.55718176 -9.47483926 -9.44281824
## 36 37 38 39 40
## -5.64684276 -3.95882875 6.07319227 3.07319227 1.38517825
## 41 42 43 44 45
## 0.00915023 -0.44281824 -0.44281824 -12.27081474 -5.78682525
## 46 47 48 49 50
## -8.27081474 8.07319227 -10.95882875 -5.30283576 5.21317475
## 51 52 53 54 55
## 10.07319227 0.90118876 -8.67886378 -6.75480423 -3.92680773
## 56 57 58 59 60
## -13.92680773 6.38517825 3.24519577 22.41719927 2.24519577
## 61 62 63 64 65
## 11.04117125 4.55718176 -7.75480423 3.55718176 -0.44281824
## 66 67 68 69 70
## -6.44281824 -2.58280073 -2.75480423 -4.23879372 -2.27081474
## 71 72 73 74 75
## 5.69716424 -4.58280073 1.69716424 10.55718176 4.21317475
## 76 77 78 79 80
## 6.07319227 -8.58280073 -13.58280073 -1.81884627 4.58920278
## 81 82 83 84 85
## 0.04117125 -0.58280073 1.32113622 -3.30283576 -7.67886378
## 86 87 88 89 90
## 0.18115373 11.66514322 -19.02287079 2.52516074 -4.81884627
## 91 92 93 94 95
## 9.55718176 -1.53888130 -2.16285327 -14.39889881 -4.19487429
## 96 97 98 99 100
## 5.66514322 12.18115373 14.46111870 12.66514322 -0.30283576
## 101 102 103 104 105
## -1.92680773 -8.44281824 9.86916775 3.21317475 -7.09881124
## 106 107 108 109 110
## -4.44281824 -8.44281824 1.69716424 -2.27081474 3.35315724
## 111 112 113 114 115
## -2.13083225 1.07319227 -14.78682525 -0.95882875 5.10521329
## 116 117 118 119 120
## -4.41079722 8.41719927 -0.47483926 1.55718176 -6.61482175
## 121 122 123 124 125
## -5.44281824 2.79322730 -12.47483926 -16.61482175 -0.44281824
## 126 127 128 129 130
## -1.09881124 1.86916775 -10.95882875 1.79322730 -3.44281824
## 131 132 133 134 135
## -6.61482175 -4.13083225 -16.36687780 2.52516074 -12.44281824
## 136 137 138 139 140
## -7.47483926 -8.47483926 1.07319227 -5.92680773 0.07319227
## 141 142 143
## -6.30283576 -4.44281824 10.86916775reg$fitted.values #推定値## 1 2 3 4 5 6 7
## 92.47484 85.44282 104.19487 83.09881 104.19487 88.95883 94.81885
## 8 9 10 11 12 13 14
## 79.58280 88.95883 91.30284 91.30284 98.33486 106.53888 91.30284
## 15 16 17 18 19 20 21
## 93.64684 91.30284 88.95883 85.44282 83.09881 85.44282 83.09881
## 22 23 24 25 26 27 28
## 85.44282 85.44282 81.92681 83.09881 74.89479 88.95883 77.23879
## 29 30 31 32 33 34 35
## 78.41080 81.92681 85.44282 74.89479 85.44282 92.47484 85.44282
## 36 37 38 39 40 41 42
## 93.64684 88.95883 81.92681 81.92681 86.61482 95.99085 85.44282
## 43 44 45 46 47 48 49
## 85.44282 84.27081 87.78683 84.27081 81.92681 88.95883 91.30284
## 50 51 52 53 54 55 56
## 87.78683 81.92681 83.09881 100.67886 80.75480 81.92681 81.92681
## 57 58 59 60 61 62 63
## 86.61482 80.75480 79.58280 80.75480 88.95883 85.44282 80.75480
## 64 65 66 67 68 69 70
## 85.44282 85.44282 85.44282 79.58280 80.75480 77.23879 84.27081
## 71 72 73 74 75 76 77
## 91.30284 79.58280 91.30284 85.44282 87.78683 81.92681 79.58280
## 78 79 80 81 82 83 84
## 79.58280 94.81885 78.41080 88.95883 79.58280 100.67886 91.30284
## 85 86 87 88 89 90 91
## 100.67886 94.81885 98.33486 103.02287 92.47484 94.81885 85.44282
## 92 93 94 95 96 97 98
## 106.53888 97.16285 112.39890 104.19487 98.33486 94.81885 106.53888
## 99 100 101 102 103 104 105
## 98.33486 91.30284 81.92681 85.44282 90.13083 87.78683 83.09881
## 106 107 108 109 110 111 112
## 85.44282 85.44282 91.30284 84.27081 93.64684 90.13083 81.92681
## 113 114 115 116 117 118 119
## 87.78683 88.95883 74.89479 78.41080 79.58280 92.47484 85.44282
## 120 121 122 123 124 125 126
## 86.61482 85.44282 70.20677 92.47484 86.61482 85.44282 83.09881
## 127 128 129 130 131 132 133
## 90.13083 88.95883 70.20677 85.44282 86.61482 90.13083 105.36688
## 134 135 136 137 138 139 140
## 92.47484 85.44282 92.47484 92.47484 81.92681 81.92681 81.92681
## 141 142 143
## 91.30284 85.44282 90.13083この推定値と実際のデータとの相関係数が当てはまりの良さ(適合度)の一つの指標です。相関係数が高ければ高いほど，良い予測ができていると言えます。 相関係数は-1から1までの値を取りますが，適合度指標としては相関係数の二乗を使います(R-squared)。 この指標も実は結果のオブジェクトに既に入っています。
summary(reg)##
## Call:
## lm(formula = weight ~ height, data = baseball)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.0229 -5.2008 -0.4428 5.1592 22.4172
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -125.5178 18.3500 -6.84 2.22e-10 ***
## height 1.1720 0.1008 11.63 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.704 on 141 degrees of freedom
## Multiple R-squared: 0.4896, Adjusted R-squared: 0.486
## F-statistic: 135.3 on 1 and 141 DF, p-value: < 2.2e-160.5ほどですね。これは実際のデータとしては，とても良い当てはまりだと言えます。
うんちく
- summary関数の出力において，推定値がどこにあるか確認してください。
- アスタリスクマークが出ていますが，これはこの回帰係数が「統計的に有意である」ことを示しています。
- R-squaredは\(R^2\)とも書き，説明される分散の大きさを表しています。
- F統計量が出ています。これは回帰分析が分散分析と数学的には同じであることの表れです！
重回帰分析
説明する変数が増えたら重回帰分析(Multiple Regression Analysis)と名前が変わります。 モデルの書き方はそれほど変わりません。
reg2 <- lm(weight~age+height,data=baseball)
reg2##
## Call:
## lm(formula = weight ~ age + height, data = baseball)
##
## Coefficients:
## (Intercept) age height
## -134.7002 0.3667 1.1607summary(reg2)##
## Call:
## lm(formula = weight ~ age + height, data = baseball)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.0067 -5.3351 0.1374 4.4695 21.4743
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -134.70024 18.41031 -7.317 1.8e-11 ***
## age 0.36669 0.14887 2.463 0.015 *
## height 1.16072 0.09912 11.711 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.57 on 140 degrees of freedom
## Multiple R-squared: 0.5108, Adjusted R-squared: 0.5038
## F-statistic: 73.09 on 2 and 140 DF, p-value: < 2.2e-16年齢が1年上だと体重が0.3kg増えるようで，身長が1kg高いと1.16cm高いようです。一応どちらも統計的にも有意であるようですが，どちらが重要なファクターなのでしょう？単位が違うから比べられないと思いませんか？
そこで，データの標準化を考えます。
baseball2 <- subset(baseball,select=c("weight","age","height")) #データの一部抜き出し
baseball2z <- scale(baseball2) #標準化
baseball2z <- data.frame(baseball2z) #データフレーム型にしておく
summary(baseball2z)## weight age height
## Min. :-2.02513 Min. :-2.49175 Min. :-2.3347
## 1st Qu.:-0.72233 1st Qu.:-0.61885 1st Qu.:-0.7761
## Median :-0.07093 Median : 0.08349 Median :-0.3085
## Mean : 0.00000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.67353 3rd Qu.: 0.66878 3rd Qu.: 0.4709
## Max. : 3.18607 Max. : 2.42463 Max. : 3.2765こうした上で，あらためて。
reg3 <- lm(weight~age+height,baseball2z)
reg3##
## Call:
## lm(formula = weight ~ age + height, data = baseball2z)
##
## Coefficients:
## (Intercept) age height
## 1.003e-15 1.458e-01 6.930e-01summary(reg3)##
## Call:
## lm(formula = weight ~ age + height, data = baseball2z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.76871 -0.49647 0.01279 0.41592 1.99833
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.003e-15 5.890e-02 0.000 1.000
## age 1.458e-01 5.918e-02 2.463 0.015 *
## height 6.930e-01 5.918e-02 11.711 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7044 on 140 degrees of freedom
## Multiple R-squared: 0.5108, Adjusted R-squared: 0.5038
## F-statistic: 73.09 on 2 and 140 DF, p-value: < 2.2e-16こうすると単位が整っていますから，大きさを比較することができます。年齢の係数は0.145,身長の係数は0.693なので，身長の方が体重との関係が強いことがわかります。
うんちく
- 有意性検定，\(R^2\)は標準化したもの(reg3)としていないもの(reg2)とで変わりがないことを確認してください。
- 標準化したデータの係数(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
## -36.471 -5.534 1.157 7.620 26.004
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3402 1.3894 3.843 0.000183 ***
## HR -0.5241 0.1030 -5.087 1.14e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.88 on 141 degrees of freedom
## Multiple R-squared: 0.1551, Adjusted R-squared: 0.1491
## F-statistic: 25.87 on 1 and 141 DF, p-value: 1.142e-06これはホームランの本数を独立変数にして，座標を従属変数にした時の回帰係数です。-0.5なのでホームランが増えるほど左のほうにプロットされるようです。もっとも\(R^2\)が小さいので，それほど説明力はありませんが・・・。
とまれ，このような使い方もできる，というのは知っておいてください。
本日の課題
- 野球のデータにおける，連続変数を使って回帰分析，重回帰分析を実践して見てください。
- 前回のMDSの座標と，野球データの変数を使って回帰分析をして見てください。
発展課題
- とあるサイト(http://yakkun.com/xy/status_list.htm) からとってきたデータをmoodleにアップしました。このデータの一部を使って書いた図が下になります。
pokemon <- read.csv("pokemon.csv",header=T)
dat <- subset(pokemon[1:30,],select=c("HP","Kougeki","Bougyo","Tokko","Tokubou","Subayasa"))
result.map <- cmdscale(dist(dat),2)
#パッケージのインストール
#install.packages("ggiraph")
library(ggplot2)
plot.dat <- transform(result.map)
plot.dat$name <- pokemon[1:30,]$name
g2 <- ggplot(plot.dat,aes(x=X1,y=X2,color=name))+geom_point()
g2
この二軸をどのように解釈すれば良いでしょうか。エビデンスとともに（統計的な根拠を元に）論じてください。