【Rでベイズ統計モデリング#14】ランダム係数モデル

記事の目的

GLMMであるランダム係数モデルのベイズ推定を、RとStanを使用して実装していきます。データの作成から実装するので、コピペで再現することが可能です。

0 前準備

0.1 今回のモデル

0.2 ワーキングディレクトリの設定

以下の画像のようにワーキングディレクトリを設定します。設定したディレクトリに、RファイルとStanファイルを保存します。

1 ライブラリ

# 1 ライブラリ
library(dplyr)
library(ggplot2)
library(rstan)
library(bayesplot)

set.seed(1)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

2 データ

2.1 コード

# 2 データ
気温 <- rnorm(100, 20,5) %>% round(1) 
休日 <- rbinom(100, 1, 2/7)
店舗 <- runif(100, 1, 4) %>% round()
店舗ランダム効果 <- rnorm(4, 0, 1)
data　<- data.frame(気温, 休日, 店舗) %>% 
  mutate(r=ifelse(店舗==1, 店舗ランダム効果[1], ifelse(店舗==2, 店舗ランダム効果[2],
                                               ifelse(店舗==3, 店舗ランダム効果[3],店舗ランダム効果[4] ))))%>%
  mutate(lambda=exp(-2+0.2*気温+(0.5+r)*休日))
data$売り上げ個数 <- rpois(100, data$lambda)
data %>% select(気温, 休日, 店舗, 売り上げ個数) %>% head()

plot <- ggplot() +
  geom_point(data=data,aes(x=気温, y=売り上げ個数, color=factor(休日))) + 
  theme_classic(base_family = "HiraKakuPro-W3") +
  theme(text=element_text(size=25))+
  labs(x="気温", y="売り上げ個数", title="データ") +
  scale_color_manual("資格",values=c("red","blue")) + 
  facet_wrap(.~ 店舗)
plot

2.2 結果

11行目の結果	20行目の結果

3 Stanの利用

3.1 Stanファイル

data {
  int N;
  int y[N];
  vector[N] x1;
  vector[N] x2;
  
  int K;
  int<lower=1, upper=K> x3[N];
}

parameters {
  vector[3] b;
  vector[K] r;
  
  real<lower=0> sigma;
}

transformed parameters{
  vector[N] lambda;
  lambda = exp(b[1] + b[2]*x1 + (b[3] + r[x3]).*x2);
}

model{
  r ~ normal(0, sigma);
  y ~ poisson(lambda);
}

3.2 Stanを利用するRのコード

# 3 stanの使用
data_list <- list(
  N = nrow(data), K = 4, y = data$売り上げ個数,
  x1 = data$気温, x2 = data$休日, x3 = data$店舗
)

mcmc_result <- stan(
  file="13ランダム係数モデル.stan",
  data=data_list,
  seed=10,
  iter = 2000, warmup = 200, chains = 4, thin=1
)

4 分析結果

4.1 コード

## 4.1 推定結果
print(mcmc_result, probs = c(0.025, 0.5, 0.975), pars=c("b","sigma","r"))

## 4.2 収束の確認
mcmc_sample <- rstan::extract(mcmc_result, permuted=FALSE)
mcmc_combo(mcmc_sample, pars=c("b[1]","b[2]","b[3]","sigma"))