GLMMをやるパッケージを二つほど。

@kosugitti
2015/11/28

GLMMをやるパッケージ

GLMM=様々な分布＋変量効果の線形モデル。 lme4パッケージもいいけどMCMCですよね，ってことになると

glmmstan (@simizu706)
rstanarm
brms

package {rstanarm}

Applied Regression Modeling via RStanの略で，Rstanのラッパーパッケージですね。使えるのは次のようなもの。

stan_lm and stan_aov
stan_glm
stan_glmer and stan_lmer
stan_polr

入れてみよう

私の環境は次の通り。

OS X El Capitan
R version 3.2.2
Rstudio Version 0.99.486

Xcodeやrstan2.8.1, devtoolsが問題なく入っていればいけると思います。

インストール

local=FALSEとしておかないといけないみたい。

library(devtools)
install_github("stan-dev/rstanarm", local = FALSE)

stan_aov

stan_aov(yield ~ block + N*P*K, data = npk, contrasts = "contr.poly",
         prior = R2(0.5), seed = 12345)

分散分析的結果の出し方

anova

stan_glm

familyはglmerに従うみたい。

counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
fit2 <- stan_glm(counts ~ outcome + treatment, family = poisson(link="log"),
                 prior = normal(0, 2.5), prior_intercept = normal(0, 10), 
                 cores = 1)
plot(fit2, ci_level = 0.95, outer_level = 0.99, show_density = TRUE)

出力例

arm

package{brms}

Bayesian Regression Models using Stanの頭文字でbrmsパッケージ。 rstanが入っている必要がありますが，入れ方そのものは

install.packages("brms")

でオーケー。

使い方もご他聞にもれずlme4パッケージと同じ書式で式を書くだけです。

デモってみましょう。

こちらのサイトに全部書いてあります。

temp <- c(11.9,14.2,15.2,16.4,17.2,18.1,18.5,19.4,22.1,22.6,23.4,25.1)
units <- c(185L,215L,332L,325L,408L,421L,406L,412L,522L,445L,544L,614L)
log_units <- log(units)
n <- length(units)
market.size <- rep(800, n)

まずは回帰

library(brms) #ライブラリ呼び出し
rstan_options(auto_write = TRUE) #呪文
options(mc.cores = parallel::detectCores()) #呪文

# Linear Gaussian model
lin.mod <- brm(units ~ temp, family="gaussian")
lin.mod

からの対数変換した正規モデル

# Log-transformed Linear Gaussian model
log.lin.mod <- brm(log_units ~ temp, family="gaussian")
plot(log.lin.mod)

ポアソンとかね。

# Poisson model
pois.mod <- brm(units ~ temp, family="poisson")
plot(pois.mod)

家族がいっぱい！

使える族が多いのが利点 family

FamilyとLink(1/3)

Family	Link	用途
gaussian	idenitity	線形回帰
student	idenitity	ロバストな線形回帰
cauchy	idenitity	ロバストな線形回帰
poisson	logit	カウントデータの回帰
negbinomial	logit	カウントデータの回帰
binomial	logit	ロジスティック回帰
bernoulli	logit	ロジスティック回帰
categorical	logit	多項ロジスティック回帰

FamilyとLink(2/3)

Family	Link	用途
cumulative	logit	順序回帰
cratio	logit	順序回帰(contiuation ratio)
sratio	logit	順序回帰(stopping ratio)
acat	logit	順序回帰(adjacent category)
gamma	log	生存分析(ガンマ分布)
weibull	log	生存分析(ワイブル分布)
exponential	log	生存分析(対数分布)
inverse.gaussian	log	生存分析(逆正規分布)

FamilyとLink(3/3)

Family	Link	用途
hurdle_poisson	log	ゼロ切断データ分析(ポアソン分布)
hurdle_negbinomial	log	ゼロ切断データ分析(負の二項分布)
hurdle_gamma	log	ゼロ切断データ分析(ガンマ分布)
zero_inflated_poisson	log	ゼロ過剰データ分析(ポアソン分布)
zero_inflated_negbinoial	log	ゼロ過剰データ分析(負の二項分布)

オプションとか

事前分布はset_priorで指定できます。

set_prior("normal(0,10)", class = "b") # 固定効果に正規分布を

#固定効果のthreat変数に正規分布を
set_prior("normal(1,2)", class = "b", coef = "treat") 

# 変量効果に半コーシーを設定
set_prior("cauchy(0,2)", class = "sd", 
                     group = "subject", coef = "Intercept")

使える分布が多くて幸せ

例として最初に挙がっているのはポアソン回帰によるてんかん患者の発作について，固定効果にt分布，変量効果に正規分布（ただし標準偏差は半コーシー）の事前分布をおいたモデルをやっています。

fit_e <- brm(count ~ log_Age_c + log_Base4_c * Trt_c
             + (1|patient) + (1|visit), 
             data = epilepsy,
             family = "poisson", 
             # 事前分布をまとめて設定
             prior = c(set_prior("student_t(5,0,10)", class = "b"),
                       set_prior("cauchy(0,2)", class = "sd")))

結果出力例

result

結果はshinystanに丸投げだ！

launch_shiny(fit_e)

Stan周りの設定

library(brms)
# Linear Gaussian model
lin.mod <- brm(units ~ temp, family="gaussian",n.chains=4,n.iter=500,
               n.warmup = 50, n.thin=1)

GLMMをやるパッケージを二つほど。