R でロジスティック回帰 - glm, MCMCpack
前回 に続き、今回も書籍 「 データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学) 」のサンプルを使って GLM とベイズ統計を試してみます。
題材は、6章「GLM の応用範囲をひろげる -ロジスティック回帰など-」のサンプルデータを使ったロジスティック回帰です。
- (1) GLM によるロジスティック回帰 （glm 関数）
- (2) MCMCpack を使ったベイズ統計によるロジスティック回帰（MCMCmetrop1R 関数）
ここで、MCMCpack のロジスティック回帰に MCMClogit 関数を使えると思ったのですが、使い方がよく分からなかったので MCMCmetrop1R 関数のみ利用しました。
サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20131229/
はじめに
MCMCpack パッケージを install.packages("MCMCpack") で R の実行環境へインストールしておきます。
ロジスティック回帰を試すためのデータは書籍のサポート web サイトから取得できます。
データ data4a.csv
N,y,x,f 8,1,9.76,C 8,6,10.48,C 8,5,10.83,C ・・・
データ内容は下記のようになっており、個体 i それぞれにおいて 「 個の観察種子のうち生きていて発芽能力があるものは
個、死んだ種子は
個」 となっています。
| 項目 | 内容 |
|---|---|
| N | 観察種子数 |
| y | 生存種子数 |
| x | 植物の体サイズ |
| f | 施肥処理 （C: 肥料なし, T: 肥料あり） |
体サイズ x と肥料による施肥処理 f が種子の生存する確率（ある個体 i から得られた 1個の種子が生存している確率）にどのように影響しているかをロジスティック回帰で解析します。
(1) GLM によるロジスティック回帰 （glm 関数）
まず、glm 関数を使ったロジスティック回帰を行って予測値の曲線を描画します。
特徴は以下のようになります。
- family に binomial （二項分布）を指定 （デフォルトのリンク関数が logit となる）
- family が binomial で応答変数 （
cbind(y, N - y)の箇所） が 2値（0, 1）以外なら cbind する必要あり （2値なら y だけでよい）
下記では、生存する確率 の関数が
、
線形予測子 で解析する事になります。
logisticGlm.R
d <- read.csv('data4a.csv') d.res <- glm(cbind(y, N - y) ~ x + f, data = d, family = binomial) summary(d.res) png("logisticGlm.png") # C: 肥料なしを 赤、T: 肥料ありを 青 で描画 plot(d$x, d$y, col = c("red", "blue")[d$f]) xx <- seq(min(d$x), max(d$x), length = 1000) ft <- factor("T", levels = levels(d$f)) fc <- factor("C", levels = levels(d$f)) # 下記でも可 #ft <- factor("T", levels = c("C", "T")) #fc <- factor("C", levels = c("C", "T")) # T: 肥料ありの予測値を算出 qq.t <- predict(d.res, data.frame(x = xx, f = ft), type="response") # C: 肥料なしの予測値を算出 qq.c <- predict(d.res, data.frame(x = xx, f = fc), type="response") # T: 肥料ありの予測値の曲線を 緑 で描画 lines(xx, max(d$N) * qq.t, col = "green") # C: 肥料なしの予測値の曲線を 黄 で描画 lines(xx, max(d$N) * qq.c, col = "yellow") dev.off()
predict 関数を使えば glm 結果の回帰モデルから新しいデータによる予測値を求めることができますので、"T:肥料あり" と "C:肥料なし" のそれぞれの予測値を求めて曲線を描画しています。
predict の結果は、生存する確率 の予測値ですので、縦軸が生存種子数
のグラフへ描画するには観察種子数 N （ここでは 8）を乗算する事になります。
ここで、"T:肥料あり" と "C:肥料なし" は factor 関数を使ってファクタ（因子）として作成しています。
d$f と同じ水準を持たせる必要がありますので levels 指定しています。
また、predict には glm で使った説明変数と同じ名前（上記の x と f）を使ったデータを渡す点に注意が必要です。
実行
実行すると下記のような結果になりました。
> R CMD BATCH logisticGlm.R
実行結果 logisticGlm.Rout
・・・
> d.res <- glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)
>
> summary(d.res)
Call:
glm(formula = cbind(y, N - y) ~ x + f, family = binomial, data = d)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.2584 -0.7948 0.1753 0.8757 3.1589
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -19.5361 1.4138 -13.82 <2e-16 ***
x 1.9524 0.1389 14.06 <2e-16 ***
fT 2.0215 0.2313 8.74 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 499.23 on 99 degrees of freedom
Residual deviance: 123.03 on 97 degrees of freedom
AIC: 272.21
Number of Fisher Scoring iterations: 5
・・・
この結果より、線形予測子 のパラメータ推定値
が
-19.5361, 1.9524, 2.0215 となります。
(2) MCMCpack を使ったベイズ統計によるロジスティック回帰（MCMCmetrop1R 関数）
次は、MCMCpack の MCMCmetrop1R 関数を使ったロジスティック回帰です。
今回のモデルの対数尤度は で、
、
なので、
これらを関数化 （下記の func） して MCMCmetrop1R 関数へ渡しています。
ちなみに、 です。（combination）
前回は、項目毎のデータ（d$x や d$y）を func 関数へ渡すようにしていましたが、今回はデータ d を丸ごと渡すようにしました。 MCMCmetrop1R(・・・, data = d, ・・・)
ファクタ（因子）のデータ （d$f） を直接計算に使えないようなので as.numeric で数値化して使っています。
なお、 の部分は
dbinom と log 関数を使うとシンプルになります。
パラメータ
の初期値は
theta.init で c(0, 0, 0) と指定しました。
logisticMcmcMetrop.R
library(MCMCpack) d <- read.csv('data4a.csv') func <- function(beta, data) { z <- beta[1] + beta[2] * data$x + beta[3] * as.numeric(data$f) q <- 1.0 / (1.0 + exp(-z)) sum(log(choose(data$N, data$y)) + data$y * log(q) + (data$N - data$y) * log(1 - q)) # 下記でも可 # sum(log(dbinom(data$y, data$N, q))) } d.res <- MCMCmetrop1R(func, theta.init = c(0, 0, 0), data = d, logfun = TRUE) summary(d.res) png("logisticMcmcMetrop.png") plot(d$x, d$y, col = c("red", "blue")[d$f]) xx <- seq(min(d$x), max(d$x), length = 1000) ft <- factor("T", levels = levels(d$f)) fc <- factor("C", levels = levels(d$f)) zt <- mean(d.res[,1]) + mean(d.res[,2]) * xx + mean(d.res[,3]) * as.numeric(ft) zc <- mean(d.res[,1]) + mean(d.res[,2]) * xx + mean(d.res[,3]) * as.numeric(fc) lines(xx, max(d$N) * 1.0 / (1.0 + exp(-zt)), col="green") lines(xx, max(d$N) * 1.0 / (1.0 + exp(-zc)), col="yellow") dev.off()
MCMCmetrop1R 結果の各パラメータの mean の値を使って、
肥料の有無（T と C）の生存確率の予測値をそれぞれ算出して、
観察種子数 N （ここでは 8）を乗算し曲線を描画しています。
実行
実行すると下記のような結果になります。
> R CMD BATCH logisticMcmcMetrop.R
実行結果 logisticMcmcMetrop.Rout
・・・
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
The Metropolis acceptance rate was 0.44590
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
>
> summary(d.res)
Iterations = 501:20500
Thinning interval = 1
Number of chains = 1
Sample size per chain = 20000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
[1,] -21.694 1.5583 0.0110185 0.037186
[2,] 1.966 0.1386 0.0009798 0.003322
[3,] 2.027 0.2268 0.0016035 0.005383
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
var1 -24.836 -22.733 -21.667 -20.596 -18.766
var2 1.703 1.870 1.963 2.059 2.248
var3 1.594 1.873 2.022 2.172 2.494
・・・
Mean の値が glm の結果と似たような値になりました。 グラフにすると違いが分からないくらいです。
の分布
の分布
の分布
