非ITエンジニアのための Stan と Rstan でベイズな重回帰分析
前回に続いて、Stan と Rstan で重回帰分析を行います。単回帰分析ができれば、重回帰分析はすんなりできます。っていうか、できました。変数が増えただけ、だからです。ほんのちょっぴり、Stan コードを洗練させました(IT系の方には「どこが?」と言われそうだけど)。
Stan で単回帰分析をしたことがない方は、前回のエントリーを見てください。先に見たほうがこのエントリーを理解しやすいと思います。
hikaru1122.hatenadiary.jp
データは再び「アイスクリーム統計学」からお借りしました。ほんとにいい本ですよね。
http://kogolab.chillout.jp/elearn/icecream/chap5/sec1.html
まず分析するデータを上記「アイスクリーム統計学」のページから、Excel にコピペ。タブ区切りのテキスト形式 icecream2.txt で保存して、Rstudio に読み込みます。何も考えなくていいので、楽チンです。でも、きっとRだけで行える方法があるんですよね…。誰が教えていただけないでしょうか?
なお、データ番号、最高気温、最低気温、客数はそれぞれ id、saiko_kion、saitei_kion、kyaku という列名にしています。
> icecream2 id saiko saitei kyaku 1 33 22 382 2 33 26 324 3 34 27 338 4 34 28 317 5 35 28 341 6 35 27 360 7 34 28 339 8 32 25 329 9 28 24 218 10 35 24 402 11 33 26 342 12 28 25 205 13 32 23 368 14 25 22 196 15 28 21 304 16 30 23 294 17 29 23 275 18 32 25 336 19 34 26 384 20 35 27 368
今回の分析で行うのは、次の重回帰式の (ベータ)を推定することです。
推定する は3つあるので、この数を p として、Stan で分析するデータをリスト形式で dat にセットします。Stan で分析するデータはリスト形式はお約束(たぶん)。
dat <- list( N = nrow(icecream2), #サンプルサイズ、20。つまりデータフレーム icecream2 の行数。 saiko_kion = icecream2$saiko, #最高気温。 saitei_kion = icecream2$saitei, #最低気温。 kyaku = icecream2$kyaku, #客数。 p = 3 #推定するベータの数。今回は3。 )
次はStan コード。単回帰分析の拡張版なので、難しくありません。icecream2.stan というファイル名で R の作業ディレクトリに保存します。作業ディレクトリの場所は getwd()
でわかります。
data { int<lower=0> N; #サンプルサイズの数。 real<lower=0> kyaku[N]; #継続したお客さんの数は20回。 real<lower=0> saiko_kion[N]; #計測した最高気温は20個。 real<lower=0> saitei_kion[N]; #計測した最低気温は20個。 int<lower=0> p; #切片および偏回帰係数の合計の数。 } parameters { real beta[p]; #切片beta1、最高気温の偏回帰係数beta2、最低気温の偏回帰係数beta3。 real<lower=0> sigma; #客数が従う正規分布の標準偏差。 } transformed parameters { real mu[N]; #mu は data と parameters セクションで定義した変数を使って定義する。 for (i in 1:N) { mu[i] <- beta[1] + beta[2]*saiko_kion[i] + beta[3]*saitei_kion[i]; #重回帰の式をここで書いて mu に代入。 } } model { kyaku ~ normal(mu, sigma); #客数は正規分布に従う。 }
あとは分析実行。私の環境では35秒ほどかかりました。結果もいっしょに出力します。サンプリング回数などは『基礎からベイズ統計学』の第7章章末問題そのままです。適切なサンプリング回数などはまだよくわからないんですよね…。
fit<-stan(file = scr, model_name = scr, data = dat, chains = cha, warmup = war, iter=ite) print(fit, pars = c("beta"))
結果は次のとおりです。上が Stan によるベイズ推定、下が一般の重回帰分析 lm(kyaku ~ saiko + saitei, data = icecream2)
の結果です。
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat beta[1] -90.57 0.21 41.39 -172.32 -117.40 -90.53 -63.87 -8.07 39158 1 beta[2] 25.96 0.01 1.63 22.72 24.90 25.96 27.01 29.20 33787 1 beta[3] -16.71 0.01 2.21 -21.07 -18.13 -16.71 -15.29 -12.30 33320 1
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -90.640 37.718 -2.403 0.0279 * saiko 25.957 1.486 17.463 2.71e-12 *** saitei -16.703 2.012 -8.301 2.20e-07 ***
だいたい同じになりました。めでたし、めでたし。今回参考にしたページは次のところです。ありがとうございました。