RとStanではじめるベイズ統計モデリングによるデータ分析入門-Part4:GLMM-
UPDATE: 2024-01-17 12:47:10.246772
はじめに
このノートは「ベイズ統計」に関する何らかの内容をまとめ、ベイズ統計への理解を深めていくために作成している。今回は「RとStanではじめるベイズ統計モデリングによるデータ分析入門」を写経していく。基本的には気になった部分を写経しながら、ところどころ自分用の補足をメモすることで、自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。
4.1.1 本章の目的と概要
階層ベイズモデルの具体例として過分散が生じているデータに対する一般化線形混合モデル(GLMM)の推定を行う。以降で必要なパッケージやデータを読み込んでおく。
library(tidyverse)
library(rstan)
library(brms)
library(bayesplot)
library(patchwork)
options(max.print = 999999)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
PATH <- 'https://raw.githubusercontent.com/logics-of-blue/book-r-stan-bayesian-model-intro/master/book-data/'
fish_num_climate_2 <- read.csv(paste0(PATH,'4-1-1-fish-num-2.csv'))
fish_num_climate_3 <- read.csv(paste0(PATH,'4-2-1-fish-num-3.csv'))
fish_num_climate_4 <- read.csv(paste0(PATH,'4-3-1-fish-num-4.csv'))4.1.3 分析の準備
湖で1時間釣りをしたときの釣果数、天気、気温のサンプルデータを利用する。サンプルサイズは100。このデータの特徴は、釣り人、釣り道具、釣りをしている時の湖の様子は全く統一されていないため、観測されている天気や気温以外の「計測されていない要因」が理由で釣果数が変化することを想定しないといけない。
fish_num_climate_2$id <- as.factor(fish_num_climate_2$id)
head(fish_num_climate_2, n = 10)## fish_num weather temperature id
## 1 0 cloudy 5.0 1
## 2 1 cloudy 24.2 2
## 3 6 cloudy 11.5 3
## 4 0 cloudy 9.8 4
## 5 1 cloudy 18.1 5
## 6 1 cloudy 18.1 6
## 7 0 cloudy 3.7 7
## 8 5 cloudy 8.8 8
## 9 9 cloudy 17.3 9
## 10 3 cloudy 18.9 104.1.4 通常のポアソン回帰モデルを適用した結果
観測できていない要因はさておき、とりあえず天気と気温を利用して釣果数をポアソン回帰モデルで推定する。
glm_pois_brms <- brm(
formula = fish_num ~ weather + temperature, # modelの構造を指定
family = poisson(), # ポアソン分布を使う
data = fish_num_climate_2, # データ
seed = 1, # 乱数の種
prior = c(set_prior("", class = "Intercept"))# 無情報事前分布にする
)推定結果はこちら。
glm_pois_brms## Family: poisson
## Links: mu = log
## Formula: fish_num ~ weather + temperature
## Data: fish_num_climate_2 (Number of observations: 100)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.30 0.16 -0.02 0.59 1.00 2701 2353
## weathersunny -0.73 0.13 -0.99 -0.48 1.00 2347 2711
## temperature 0.06 0.01 0.04 0.08 1.00 2781 2499
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).予測区間を可視化する。広めの99％予測区間の外側にも多くのデータが存在しており、\(\lambda\)で平均と分散を決めるポアソン分布では対応できていない。想定された分散よりも大きくデータがばらついている。
set.seed(1)
eff_glm_pre <- conditional_effects(
glm_pois_brms,
method = "predict",
effects = "temperature:weather",
probs = c(0.005, 0.995))
plot(eff_glm_pre, points = T)
4.1.5 過分散対処のためのGLMMの構造
モデルの改良を試みる。1時間釣りをした時の釣果数は、気温、天気によって一般化線形モデルで表現する。釣果数はポアソン分布に従い、\(\lambda\)が気温と湿度で変化すると想定していた。ただ、これではうまく説明できなかった。
\[ \begin{eqnarray} log(\lambda_{i}) &=& \beta_{0} + \beta_{1}x_{i1} + \beta_{2}x_{i2} \\ y_{i} &\sim& Poisson(\lambda_{i}) \\ \end{eqnarray} \]
調査のたびに調査した人、道具、技術、湖の状態は変化するが観測できていない。そのため、これらの変化を表現できるようにモデルにランダム効果を取り込む。つまり、調査ごとに変化するランダムな影響\(r_{i}\)を組み込み。ランダムな影響\(r_{i}\)は平均0、分散\(\sigma^{2}_{r}\)の正規分布に従うとする。\(\sigma^{2}_{r}\)はランダム効果\(r_{i}\)というパラメタの分散を表すパラメタであるため、ハイパーパラメタとも呼ばれる。ランダム効果\(r_{i}\)は100人分存在するため、100個の事後分布が得られるが、好き勝手な分布にならないように、この分散パラメタの大きさを限定することで、ランダム効果\(r_{i}\)の大きさを限定する。
\[ \begin{eqnarray} r_{i} &\sim& Normal(0, \sigma^{2}_{r}) \\ log(\lambda_{i}) &=& \beta_{0} + \beta_{1}x_{i1} + \beta_{2}x_{i2} + r_{i} \\ y_{i} &\sim& Poisson(\lambda_{i}) \\ \end{eqnarray} \]
\(\beta_{0},\beta_{1},\beta_{2}\)というパラメタは固定効果と呼ばれ、新たに追加したランダム効果とは異なる。ランダム効果は何らかの確率分布に従って生成される。固定効果とランダム効果が含まれるモデルを混合モデル、一般家線形混合モデルと呼ぶ。このように、上位の層の確率変数の実現値が、下位の層の確率分布のパラメタとなるモデルを階層ベイズモデルとも呼ぶ。
4.1.10 brmsによるGLMMの推定
(1|id)という表記を追加することで、ランダム効果を追加できる。縦棒の左の1は切片を意味する。縦棒の右がグループ名を意味する。つまり、idごとで切片にランダム効果を入れることを意味する。
glmm_pois_brms <- brm(
formula = fish_num ~ weather + temperature + (1|id), # ランダム効果
family = poisson(), # ポアソン分布を使う
data = fish_num_climate_2, # データ
seed = 1, # 乱数の種
prior = c(set_prior("", class = "Intercept"),
set_prior("", class = "sd")) # 無情報事前分布にする
)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 14.0.0 (clang-1400.0.29.202)’
## using SDK: ‘MacOSX13.1.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1Group-Level Effectsはランダム効果に関わるパラメタで、Population-Level Effectsは固定効果を表す。
# 結果の表示
# plot(glmm_pois_brms)
# stancode(glmm_pois_brms)
glmm_pois_brms## Family: poisson
## Links: mu = log
## Formula: fish_num ~ weather + temperature + (1 | id)
## Data: fish_num_climate_2 (Number of observations: 100)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Group-Level Effects:
## ~id (Number of levels: 100)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.10 0.15 0.83 1.44 1.00 1073 1645
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.46 0.33 -1.14 0.15 1.01 1319 2427
## weathersunny -0.72 0.29 -1.31 -0.15 1.00 963 1760
## temperature 0.08 0.02 0.04 0.11 1.01 862 1632
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).パッケージを使わずに自ら推定する場合は下記を参照。
################################
# Stan
################################
# # ダミー変数を作る
# formula_pois <- formula(fish_num ~ weather + temperature)
# design_mat <- model.matrix(formula_pois, fish_num_climate_2)
# sunny_dummy <- as.numeric(design_mat[, "weathersunny"])
#
# # データの作成
# data_list_1 <- list(
# N = nrow(fish_num_climate_2),
# fish_num = fish_num_climate_2$fish_num,
# temp = fish_num_climate_2$temperature,
# sunny = sunny_dummy
# )
# # 結果の表示
# data_list_1
#
# # MCMCの実行
# glmm_pois_stan <- stan(
# file = "4-1-1-glmm-pois.stan",
# data = data_list_1,
# seed = 1
# )
#
# # 収束の確認
# mcmc_rhat(rhat(glmm_pois_stan))
#
# # 参考：トレースプロットなど
# mcmc_sample <- rstan::extract(glmm_pois_stan, permuted = FALSE)
# mcmc_combo(
# mcmc_sample,
# pars = c("Intercept", "b_sunny", "b_temp", "sigma_r", "lp__"))
#
#
# # 結果の表示
# print(glmm_pois_stan,
# pars = c("Intercept", "b_sunny", "b_temp", "sigma_r"),
# probs = c(0.025, 0.5, 0.975))
################################
# Model
################################
// data {
// int N; // サンプルサイズ
// int fish_num[N]; // 釣獲尾数
// vector[N] sunny; // 晴れダミー変数
// vector[N] temp; // 気温データ
// }
//
// parameters {
// real Intercept; // 切片
// real b_temp; // 係数(気温)
// real b_sunny; // 係数(晴れの影響)
// vector[N] r; // ランダム効果
// real<lower=0> sigma_r; // ランダム効果の標準偏差
// }
//
// transformed parameters{
// vector[N] lambda = Intercept + b_sunny*sunny + b_temp*temp + r;
// }
//
// model {
// r ~ normal(0, sigma_r);
// fish_num ~ poisson_log(lambda);
// }4.2.1 本章の目的と概要
先程はすべてのデータに対してランダム効果を含めるモデルを推定したが、ここではランダム効果をグループ単位で含めるモデルを推定する。つまり、同一グループ内では、同じランダム効果を加わえるモデルを推定する。ランダム効果の加え方が異なるだけで、ランダム切片モデルと同じである。
4.2.2 分析の準備
湖で1時間釣りをしたときの釣果数、天気、気温のサンプルデータを利用する。サンプルサイズは100。先程のデータとの違いは、調査員が10人(A-J)おり、10回分の自分が担当した釣果数に関しては調査員が識別できるようになっている。調査員ごとに釣りの技術が異なるので、このようなバイアスを排除したいという意図がある。
fish_num_climate_3 %>%
group_by(human) %>%
summarise(
cnt = n(),
m_fish_num = mean(fish_num)
)## # A tibble: 10 × 3
## human cnt m_fish_num
## <chr> <int> <dbl>
## 1 A 10 5.5
## 2 B 10 1.9
## 3 C 10 4.6
## 4 D 10 0.9
## 5 E 10 1.7
## 6 F 10 0.8
## 7 G 10 3.3
## 8 H 10 1.9
## 9 I 10 2.3
## 10 J 10 1.94.2.3 ランダム切片モデルの構造
\(k\)はA-Jの10人分で、ランダム効果は平均0、分散\(\sigma^{2}_{r}\)に従うと仮定する。ランダム効果としては\(r_{A}=+1, r_{B}=-2,...,r_{J}=+3\)のような形で、釣果数に影響を与える。
\[ \begin{eqnarray} r_{k} &\sim& Normal(0, \sigma^{2}_{r}) \\ log(\lambda_{i}) &=& \beta_{0} + \beta_{1}x_{i1} + \beta_{2}x_{i2} + r_{k} \\ y_{i} &\sim& Poisson(\lambda_{i}) \\ \end{eqnarray} \]
4.2.4 ランダム効果の使い所
ランダム効果を固定効果として取り込むことで、切片を調整することはもちろん可能ではある。分析の目的は、天気や気温と釣果数の関係を調べることであり、各調査員の釣り技術を分析することではない。固定効果として調べることも可能だけれども、その必要はないので、ランダム効果としている。
4.2.5 brmsによるランダム切片モデルの推定
(1|human)と指定すればOK。考え方を変えると(1|id)でもidが100グループあると考えれば良い。
glmm_pois_brms_human <- brm(
formula = fish_num ~ weather + temperature + (1|human),
family = poisson(),
data = fish_num_climate_3,
seed = 1,
prior = c(set_prior("", class = "Intercept"),
set_prior("", class = "sd"))
)推定結果はこちら。ランダム効果のばらつきの大きさ\(\sigma_{r}\)は0.65となっている。
# plot(glmm_pois_brms_human)
# stanplot(glmm_pois_brms_human, type = "rhat")
glmm_pois_brms_human## Family: poisson
## Links: mu = log
## Formula: fish_num ~ weather + temperature + (1 | human)
## Data: fish_num_climate_3 (Number of observations: 100)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Group-Level Effects:
## ~human (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.65 0.22 0.35 1.20 1.00 820 1100
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.87 0.30 -1.50 -0.31 1.00 1520 1529
## weathersunny -0.52 0.13 -0.78 -0.26 1.00 3570 2933
## temperature 0.10 0.01 0.08 0.12 1.00 4290 3179
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).各々の調査者の影響の大きさを調べたい場合は、ranef()関数を利用する。ランダム効果としては\(r_{A}=0.76,
r_{B}=0.08,...,r_{J}=-0.15\)のような形で、釣果数に影響を与えている。\(r_{F}=-0.73\)であるため、Fさんは釣りが苦手であることがわかる。
ranef(glmm_pois_brms_human)## $human
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## A 0.76154307 0.2606235 0.27306232 1.30732897
## B 0.08633475 0.2993386 -0.48967083 0.69408085
## C 0.69400206 0.2668357 0.20209372 1.25840029
## D -0.61530646 0.3367841 -1.31656524 0.03750431
## E -0.11743562 0.2990496 -0.72568979 0.45232235
## F -0.73877968 0.3476499 -1.45848467 -0.06341159
## G 0.45411672 0.2738553 -0.04625609 1.02393953
## H -0.31098554 0.2951315 -0.88107846 0.26629226
## I 0.02825932 0.2853584 -0.50965747 0.60159890
## J -0.15326425 0.2954807 -0.73720697 0.466154074.2.6 回帰曲線の図示
ランダム切片効果モデルを可視化する。
eff_glmm_human1 <- conditional_effects(
glmm_pois_brms_human,
effects = "temperature",
re_formula = NULL)
plot(eff_glmm_human1, points = TRUE)
effects = "temperature:weather"と指定することで交互作用を表現でき、conditionsを指定することでグループごとに可視化できる。グループごとに可視化することで調査員の能力を反映したグラフを確認できる。
conditions <- data.frame(human = LETTERS[1:10])
eff_glmm_human2 <- conditional_effects(
glmm_pois_brms_human,
effects = "temperature:weather",
re_formula = NULL,
conditions = conditions)
plot(eff_glmm_human2, points = TRUE)
4.3.1 本章の目的と概要
最後はランダム係数モデルを推定する。係数にランダム効果を入れることで、説明変数の固定効果の強さが増減するモデルを表現できる。
4.3.2 分析の準備
釣果数と気温の関係を記録しているデータを利用する。調査員Jに関しては記録が4しかない。
fish_num_climate_4 %>%
group_by(human) %>%
summarise(
cnt = n(),
m_fish_num = mean(fish_num)
)## # A tibble: 10 × 3
## human cnt m_fish_num
## <chr> <int> <dbl>
## 1 A 10 3.7
## 2 B 10 9.7
## 3 C 10 10.7
## 4 D 10 6.2
## 5 E 10 18.9
## 6 F 10 9.5
## 7 G 10 8.2
## 8 H 10 8.8
## 9 I 10 2.8
## 10 J 4 10.24.3.3 交互作用を用いたモデル化
釣りをした人と気温の違いによって、釣果数が変化するモデルであれば、交互作用(カテゴリ×数量)を想定したモデルでも良い。ただ、交互作用モデルは調査員間の関係を考慮できない。
glm_pois_brms_interaction <- brm(
formula = fish_num ~ temperature * human,
family = poisson(),
data = fish_num_climate_4,
seed = 1,
prior = c(set_prior("", class = "Intercept"))
)可視化するとわかりやすい。番号10は調査員Jのことであり、4つしか記録がない。そのため、この4つのデータからのみ推定を行うため、傾きがマイナスになっている。
# glm_pois_brms_interaction
# stanplot(glm_pois_brms_interaction, type = "rhat")
eff_1 <- marginal_effects(glm_pois_brms_interaction,
effects = "temperature",
conditions = conditions)
plot(eff_1, points = TRUE)
ただ、この結果は分析担当者の感覚としては違和感がある場合があったとする。データが多ければプラスの傾きで推定することができたはずだとする。ランダム効果を用いることで、全体と似たような傾向を表現するように調整できる。これををランダム効果の縮約という。これは、ランダム効果が特定の確率分布に従って生成されるという仮定を立てているためである。これはランダム係数効果に限った話ではない。
4.3.5 ランダム係数モデルの構造
ランダム係数効果は、ランダム切片効果と共に利用されることが多いため、ここでは双方を含めたモデルを推定する。
\[ \begin{eqnarray} r_{k} &\sim& Normal(0, \sigma^{2}_{r}) \\ \tau_{k} &\sim& Normal(0, \sigma^{2}_{\tau}) \\ log(\lambda_{i}) &=& \beta_{0} + (\beta_{1} + \tau_{k})x_{i1} + r_{k} \\ y_{i} &\sim& Poisson(\lambda_{i}) \\ \end{eqnarray} \]
4.3.6 brmsによるランダム係数モデルの推定
ランダム係数とランダム切片効果モデルを推定する際はtemperature||humanと指定する。縦棒の左が「ランダム効果を与えられる変数」で、縦棒の右がグループの変数を表す。
縦棒を減らして、temperature|humanとも指定できる。この場合、ランダム切片とランダム係数に相関を認めているモデルであり、先程のモデルはランダム切片とランダム係数に相関がないと考えていることになる。釣りが上手い=ランダム切片がプラス、気温があがる=ランダム係数もプラスとなるような関係が想定できるなら、こちらのモデルを利用するべき。
glmm_pois_brms_keisu <- brm(
formula = fish_num ~ temperature + (temperature||human),
family = poisson(),
data = fish_num_climate_4,
seed = 1,
iter = 6000,
warmup = 5000,
control = list(adapt_delta = 0.97, max_treedepth = 15)
)4.3.7 回帰曲線の図示
可視化すると番号10の調査員Jの傾きはプラスになっている。これはランダム効果で縮約された結果である。
# glmm_pois_brms_keisu
# plot(glmm_pois_brms_keisu)
# prior_summary(glmm_pois_brms_keisu)
# stanplot(glmm_pois_brms_keisu, type = "rhat")
eff_2 <- marginal_effects(glmm_pois_brms_keisu,
re_formula = NULL,
effects = "temperature",
conditions = conditions)
plot(eff_2, points = TRUE)