class: center, middle, inverse, title-slide # R Club - Day 12 ##
https://opueco.github.io/rclub-slides/2019f-12/slides.html
### Kenji Sato ### 2019/7/24 --- ## 先週やったこと - 重回帰分析 - 回帰係数の検定に関する詳細 - 重み付き線形回帰モデル --- ## 今日やること - 最尤法 - ロジスティック回帰 --- class: section ## 最尤法 --- ## 尤度(ゆうど) コインを10回投げて表が出る回数の確率変数を `\(X\)` とする。表が出る確率を `\(\theta\)` とすれば, `\(X\)` の確率分布は `$${}_{10}C_X \theta^X (1-\theta)^{1-X}$$` 実際にコイン投げをして `\(X = 4\)` であることを観測した。この事象が起こる確率は `$$L(\theta) = {}_{10}C_4 \theta^4 (1-\theta)^{10-4}$$` `\(L(\theta)\)` を尤度という。これは,特定の分布,パラメータを想定したときに,データが観測される確率である。 --- ## 最尤法(さいゆうほう)(1/) `\(\theta\)` を選ぶ方法として, `$$\max L(\theta)$$` となる `\(\theta\)` を選ぶのが最尤法 (Maximum Likelihood, ML) `\(L(\theta)\)` を最大化するかわりに対数尤度 `\(\log L(\theta)\)` を最大化する方が計算が簡単になる場合が多い。 --- ## 最尤法 (2/) 上の例では `$$\log L(\theta) = \log {}_{10}C_4 + 4 \log \theta + 6 \log (1-\theta)$$` なので, `$$\max_\theta \log L(\theta) \Longrightarrow \frac{4}{\theta} - \frac{6}{1-\theta}=0$$` これを `\(\theta\)` について解くと `$$\theta = 0.4$$` --- class: section ## ロジスティック回帰 --- ## 問題設定 線形回帰モデルは,被説明変数の値が連続的だった。 ロジスティック回帰モデルは「2値判別」のために使う。 - 例 - 「試験の合否」を模擬試験・内申書などで予測する - 「病気の有無」を性別・年齢・血圧などで予測する - 「女性の就業状態」を夫の収入,年齢,教育年数,子供の数で予測する --- ## 基本的な考え方 (1/) 観測対象 `\(i=1,\dots,n\)` について - 被説明変数 `\(y_i = 0\)` or `\(1\)` (便宜的に 0-1 とする) と - 説明変数 `\((x_{i1},\dots, x_{ik})\)` を持っている。説明変数のベクトルを `$$x_i = \begin{bmatrix} 1 & x_{i1} & \cdots & x_{ik} \end{bmatrix}$$` --- ## 基本的な考え方 (2/) 係数ベクトル `$$\beta =\begin{bmatrix} \beta_0 \\ \beta_{1} \\ \vdots \\ \beta_{k} \end{bmatrix}$$` をうまく選んで, `$$x_i \beta$$` に関する情報から `\(y_i\)` をできるだけうまく言い当てたい。 --- ## 基本的な考え方 (3/) 0-1 の値を直接推定するのではなく, `$$p_i = \mathrm{Prob}(y_i = 1)$$` と考えて, `\(x_i \beta\)` と `\(p_i\)` を関連付ける方法を探せばよい。そこで, `$$x_i\beta \text{ が大きい} \Longleftrightarrow p_i \text{ が1に近い}$$` `$$x_i\beta \text{ が小さい} \Longleftrightarrow p_i \text{ が0に近い}$$` --- ## `\(x_i \beta\)` と `\(p_i\)` の関連付け (1/) 性質 `$$F' (x) \geq 0$$` `$$\lim_{x \to +\infty} F(x) = 1$$` `$$\lim_{x \to -\infty} F(x) = 0$$` を持つ関数 `\(F\)` を使って `$$p_i = F(x_i \beta)$$` とすればよい。 `\(F\)` (正確には `\(F^{-1}\)`)を**リンク関数**という --- ## `\(x_i \beta\)` と `\(p_i\)` の関連付け (2/) リンク関数 `\(F\)` として標準正規分布の累積分布関数 `\(\Phi\)` を用いるモデルを**プロビット・モデル**という。 `$$\Phi(x) = (2\pi)^{-1/2} \int_{-\infty}^x \exp (-z^2 / 2)dz$$` ```r > plot(pnorm, xlim = c(-8, 8)) ``` <img src="slides_files/figure-html/probit-1.png" height="220" style="display: block; margin: auto;" /> --- ## `\(x_i \beta\)` と `\(p_i\)` の関連付け (3/) `\(F\)` としてロジスティック分布の分布関数を用いるモデルを **ロジット・モデル** (**ロジスティック回帰モデル**)という。 `$$\mathrm{logit}^{-1}(x) = \frac{e^x}{1 + e^x}$$` ```r > curve(exp(x) / (1 + exp(x)), xlim = c(-8, 8)) ``` <img src="slides_files/figure-html/logit-1.png" height="220" style="display: block; margin: auto;" /> --- ## `\(x_i \beta\)` と `\(p_i\)` の関連付け (4/) ロジット `$$p = \frac{e^x}{1 + e^x}$$` を `\(x\)` について解くと, 対数オッズが現れる。 `$$x = \log \left(\frac{p}{1 - p}\right)$$` これを使うと, `\(\beta\)` の線形関数として対数オッズを表現できる。 `$$x_i \beta = \log \left(\frac{p_i}{1 - p_i}\right)$$` --- ## 最尤推定 (1/) `\(y_i\)` のそれぞれは `\(p_i = F(x_i\beta)\)` のベルヌーイ分布であるから,尤度関数は `$$L(\beta; y) = \prod_{i=1}^n p_i^{y_i}(1-p_i)^{1 - y_i}$$` 対数尤度関数は `$$\begin{aligned} \log L(\beta; y) &= \sum_{i=1}^n \left[ y_i \log p_i + (1 - y_i) \log (1 - p_i) \right]\\ &= \sum_{i=1}^n \left[ y_i \log F(x_i\beta) + (1 - y_i) \log (1 - F(x_i\beta)) \right] \end{aligned}$$` --- ## 最尤推定 (2/) `$$\begin{aligned} \frac{\partial}{\partial \beta} \log L(\beta; y) &= \sum_{i=1}^n \left[ \frac{y_i}{F(x_i\beta)} - \frac{1 - y_i}{1 - F(x_i\beta)} \right] F'(x_i\beta) x_i^\top \end{aligned}$$` これがゼロになる `\(\beta\)` を数値計算で探せばよい。 R では `glm()` で実行できる。 --- ## 例: `wooldridge::mroz` (1/) ```r > library(wooldridge) > data(mroz) > > fml <- (inlf ~ nwifeinc + educ + exper + + I(exper^2) + age + kidslt6 + kidsge6) > ml <- glm(fml, data = mroz, family = binomial(link = "logit")) > ml ``` ``` ## ## Call: glm(formula = fml, family = binomial(link = "logit"), data = mroz) ## ## Coefficients: ## (Intercept) nwifeinc educ exper ## 0.425452 -0.021345 0.221170 0.205870 ## I(exper^2) age kidslt6 kidsge6 ## -0.003154 -0.088024 -1.443354 0.060112 ## ## Degrees of Freedom: 752 Total (i.e. Null); 745 Residual ## Null Deviance: 1030 ## Residual Deviance: 803.5 AIC: 819.5 ``` --- ```r > summary(ml) ``` ``` ## ## Call: ## glm(formula = fml, family = binomial(link = "logit"), data = mroz) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.1770 -0.9063 0.4473 0.8561 2.4032 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 0.425452 0.860365 0.495 0.62095 ## nwifeinc -0.021345 0.008421 -2.535 0.01126 * ## educ 0.221170 0.043439 5.091 3.55e-07 *** ## exper 0.205870 0.032057 6.422 1.34e-10 *** ## I(exper^2) -0.003154 0.001016 -3.104 0.00191 ** ## age -0.088024 0.014573 -6.040 1.54e-09 *** ## kidslt6 -1.443354 0.203583 -7.090 1.34e-12 *** ## kidsge6 0.060112 0.074789 0.804 0.42154 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 1029.75 on 752 degrees of freedom ## Residual deviance: 803.53 on 745 degrees of freedom ## AIC: 819.53 ## ## Number of Fisher Scoring iterations: 4 ``` --- ## 例: `wooldridge::mroz` (2/) ```r > mp <- glm(fml, data = mroz, family = binomial(link = "probit")) > mp ``` ``` ## ## Call: glm(formula = fml, family = binomial(link = "probit"), data = mroz) ## ## Coefficients: ## (Intercept) nwifeinc educ exper ## 0.270074 -0.012024 0.130904 0.123347 ## I(exper^2) age kidslt6 kidsge6 ## -0.001887 -0.052852 -0.868325 0.036006 ## ## Degrees of Freedom: 752 Total (i.e. Null); 745 Residual ## Null Deviance: 1030 ## Residual Deviance: 802.6 AIC: 818.6 ``` --- ```r > summary(mp) ``` ``` ## ## Call: ## glm(formula = fml, family = binomial(link = "probit"), data = mroz) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.2156 -0.9151 0.4315 0.8653 2.4553 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 0.2700736 0.5080782 0.532 0.59503 ## nwifeinc -0.0120236 0.0049392 -2.434 0.01492 * ## educ 0.1309040 0.0253987 5.154 2.55e-07 *** ## exper 0.1233472 0.0187587 6.575 4.85e-11 *** ## I(exper^2) -0.0018871 0.0005999 -3.145 0.00166 ** ## age -0.0528524 0.0084624 -6.246 4.22e-10 *** ## kidslt6 -0.8683247 0.1183773 -7.335 2.21e-13 *** ## kidsge6 0.0360056 0.0440303 0.818 0.41350 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 1029.7 on 752 degrees of freedom ## Residual deviance: 802.6 on 745 degrees of freedom ## AIC: 818.6 ## ## Number of Fisher Scoring iterations: 4 ``` --- ## 限界効果 (1/) 係数の解釈には注意が必要。 例えば,`educ` にかかる係数 ≒ 0.13 であるというのを見ただけでは「`educ` が1増えたときに `inlf==1` になる確率の増分」は分からない。 `\(\beta\)` の意味を理解するには確率のモデルに戻る必要がある。 `$$p = F(x \beta)$$` 説明変数 `\(x_{j}\)` が1単位増加したときの確率の変化「限界効果」は `$$\frac{\partial p}{\partial x_{j}} = F'(x\beta) \beta_j$$` --- ## 限界効果 (2/) 限界効果は説明変数の値 `\(x\)` に依存してしまうので,代表値をレポートするのが定石。 - `\(x\)` のすべての観測値について限界効果を計算して,それを平均する - `\(x\)` の平均値を求めて,平均における限界効果を計算する という2つの方法を紹介しよう。 以下ではプロビットの場合を紹介する。ロジットのケースはドリルを参照。 --- ## 限界効果 (3/) #### 限界効果の平均 `$$\mathrm{AME}_j = \frac{1}{n} \sum_{i=1}^n F'(x_{i\bullet}\hat\beta) \hat\beta_j$$` ```r > mean(dnorm(predict(mp))) * coef(mp)[-1] ``` ``` ## nwifeinc educ exper I(exper^2) ## -0.003616176 0.039370095 0.037097345 -0.000567546 ## age kidslt6 kidsge6 ## -0.015895665 -0.261153464 0.010828887 ``` --- ## 限界効果 (4/) #### 平均における限界効果 `$$\mathrm{MEM}_j = F'(\bar{x} \hat{\beta}) \hat{\beta}_j$$` `colMeans()` で各列の平均を計算したあとに, `predict()` にわたして `\(\bar{x}\hat\beta\)` を計算する。 ```r > mp_bar <- as.data.frame.list(colMeans(mroz)) > mean(dnorm(predict(mp, newdata = mp_bar))) * coef(mp)[-1] ``` ``` ## nwifeinc educ exper I(exper^2) ## -0.0045447163 0.0494793222 0.0466229885 -0.0007132772 ## age kidslt6 kidsge6 ## -0.0199772628 -0.3282109539 0.0136094667 ``` --- ## 限界効果 (5/) 注意! 上の計算はモデルに `exper` の 2次項が入っていることを考慮していない。きちんと計算するには **margins** パッケージを使うとよい。 ```r > install.packages("margins") ``` ```r > library(margins) > (ame <- margins(mp)) ``` ``` ## Average marginal effects ``` ``` ## glm(formula = fml, family = binomial(link = "probit"), data = mroz) ``` ``` ## nwifeinc educ exper age kidslt6 kidsge6 ## -0.003616 0.03937 0.02558 -0.0159 -0.2612 0.01083 ``` --- ## 限界効果 (6/) `exper` のAMEは次のように計算している。 ```r ame_exper <- mean(dnorm(predict(mp)) * (coef(mp)["exper"] + 2 * coef(mp)["I(exper^2)"] * mroz$exper)) ame_exper ``` ``` ## [1] 0.02558251 ``` 問題1: なぜこのコードでうまくいくのか考えなさい。 問題2: `exper` のMEMを計算しなさい。 --- ## 問題2 の解答例 ```r > xbar <- as.data.frame.list( + colMeans(get_all_vars(fml, data = mroz)[-1])) > margins(mp, at = xbar) ``` ``` ## Average marginal effects at specified values ``` ``` ## glm(formula = fml, family = binomial(link = "probit"), data = mroz) ``` ``` ## at(nwifeinc) at(educ) at(exper) at(age) at(kidslt6) ## 20.13 12.29 10.63 42.54 0.2377 ## at(kidsge6) nwifeinc educ exper age kidslt6 ## 1.353 -0.004545 0.04948 0.03146 -0.01998 -0.3282 ## kidsge6 ## 0.01361 ``` --- class: section ## モデルの評価 --- ## 判定の精度 (1/) `predict()` と `fitted()` の違いに注意する。 `predict()` で `\(x_i\beta\)` を計算,これに `\(F\)` を適用した `\(F(x_i\beta)\)` は `fitted()` で計算される。 ロジット ```r > invlogit <- function(x) exp(x) / (1 + exp(x)) > sum((invlogit(predict(ml)) - fitted(ml)) ^ 2) ``` ``` ## [1] 0 ``` プロビット ```r > sum((pnorm(predict(mp)) - fitted(mp))^2) ``` ``` ## [1] 0 ``` --- ## 判定の精度 (2/) `\(p = F(x\hat{\beta}) \ge 0.5\)` となるときに `\(y = 1\)` となるように定めてみよう。(ボーダーラインの 0.5 という数字は適当に定める) 判定の正しさを表にする。対角成分が正しく判定できた数。 ```r > tbl <- table(fitted(ml) >= 0.5, mroz$inlf) > tbl ``` ``` ## ## 0 1 ## FALSE 207 81 ## TRUE 118 347 ``` ```r > sum(diag(tbl)) / sum(tbl) ``` ``` ## [1] 0.7357238 ``` --- ## 判定の精度 (3/) 予測性能の評価についてもう少し真面目に考えるならば, - `\(\beta\)` を決めるのに使うデータ(トレーニングデータ) - 性能を評価するデータ(テストデータ) を分けなければいけない。学習データに対する性能がよくても, 他のデータの予測には使えないかもしれない。 R では `sample()`, `caret::createDataPartition()`, `modelr::resample_partition()` などの方法でデータを分割できるので各自試してほしい。 --- ## ROC曲線 (1/) ```r > source("https://git.io/fjP7v") # ROC.R > ROC(fitted(ml), mroz$inlf) ``` <img src="slides_files/figure-html/roc1-1.png" height="400" style="display: block; margin: auto;" /> ``` ## AUC = 0.8014378 th = 0.5623703 ## BER = 0.2595291 OR = 8.199634 ## Actual ## Predicted 0 1 ## FALSE 233 101 ## TRUE 92 327 ``` --- ## ROC曲線 (2/) - `\(\alpha\)` は `\(p \ge \alpha\)` のとき `\(y = 1\)` と判断する閾値 - `\(\alpha\)` を 0から1まで動かして偽陽性率と真陽性率の組み合わせをプロットしたもの。 ```r > tbl # α = 0.5 のとき ``` ``` ## ## 0 1 ## FALSE 207 81 ## TRUE 118 347 ``` ```r > tbl[2, ] / colSums(tbl) ``` ``` ## 0 1 ## 0.3630769 0.8107477 ``` --- ## ROC曲線 (3/) AUC = Area Under Curve は性能の指標の1つ ```r > ROC(fitted(ml), mroz$inlf) ``` <img src="slides_files/figure-html/roc2-1.png" height="200" style="display: block; margin: auto;" /> ``` ## AUC = 0.8014378 th = 0.5623703 ## BER = 0.2595291 OR = 8.199634 ## Actual ## Predicted 0 1 ## FALSE 233 101 ## TRUE 92 327 ``` --- ## AIC (1/) モデルに含める変数をどのように定めるかについて「情報量基準」と呼ばれる概念がいくつか提唱されている。 赤池情報量基準 (AIC) `$$\begin{aligned} \mathrm{AIC} = - 2\log L(\beta; y) + 2(k + 1) \end{aligned}$$` を最小にする変数の組み合わせを選ぶ,というのが1つの基準である。 ```r > y <- mroz$inlf; ftd <- fitted(ml) > aic <- - 2 * sum(y * log(ftd) + + (1 - y) * log(1 - ftd)) + 2 * 8 > aic # AIC(ml) ``` ``` ## [1] 819.5303 ``` --- ## AIC (2/) `step()` を使うと,AIC を減少させる変数の組み合わせを教えてくれる。 ```r > ml2 <- step(ml) > ml3 <- step(ml2) ``` --- ## 練習してみよう ```r devtools::install_github("opueco/R4FunDrill") R4FunDrill::start("Day12") ``` チュートリアルがはじまるよ。