Rで解析：事前分布で混合モデルの推定を安定させる「blme」パッケージ

「blme」パッケージは、混合モデルの分散や係数に事前分布を置き、推定を安定させるパッケージです。グループの数が少ないときに、グループ間の分散が0に張り付いたり、相関が±1になったりするのを和らげられます。渡すのは、lme4パッケージのlmerやglmerと同じ式とデータフレームで、正規分布の応答も、二項分布やポアソン分布の応答も扱えます。lme4パッケージの上に乗るパッケージなので、lme4で混合モデルを当てはめている人が、書き方をほとんど変えずに試せるのではないかと考えます。

パッケージバージョンは1.0-8。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

＜おすすめのRに関する書籍です＞

R入門 毎週ちょっと、統計解析

少しずつ進められる構成のR入門。まとまった時間が取りにくくても続けやすいようにまとめられている。

www.amazon.co.jp

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("blme")

# パッケージの読み込み
library("blme")
Loading required package: lme4
Loading required package: Matrix

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。例は、宇治茶の茶園8か所で5年間、施肥量を変えて一番茶の収量を記録した架空のデータです。

混合モデルは、茶園ごとの違いを「ランダム効果」として、全体に共通する施肥量の効き目を「固定効果」として同時に推定します。茶園が8か所ほどしかないと、茶園間のばらつきを0と推定してしまうことがあります。blmeパッケージは、そのばらつき（共分散）や固定効果に事前分布を置いて、推定値が端に寄りすぎないようにします。なお、lme4 2.0.0以降の構造化共分散（diag、cs、ar1）にはまだ対応していません。

# 乱数の種を固定
set.seed(2026)

# 宇治茶の茶園8か所
farm <- c("宇治", "宇治田原", "和束", "南山城", "城陽", "京田辺", "木津川", "久御山")

# 茶園と年（1〜5年目）の組み合わせ
tea <- expand.grid(year = 1:5, farm = farm)

# 施肥量（kg/10a）
tea$fert <- round(runif(nrow(tea), 20, 60))

# 茶園ごとの収量の高さと、施肥量の効き方の違い
farm_eff <- rnorm(8, 0, 25)
slope_eff <- rnorm(8, 0, 0.6)

# 一番茶の収量（kg/10a）
tea$yield <- round(400 + farm_eff[as.integer(tea$farm)] +
  (3 + slope_eff[as.integer(tea$farm)]) * tea$fert + rnorm(nrow(tea), 0, 20))

# データの先頭を確認
head(tea)
  year     farm fert yield
1    1     宇治   48   566
2    2     宇治   42   516
3    3     宇治   26   491
4    4     宇治   31   480
5    5     宇治   42   511
6    1 宇治田原   21   496

事前分布付きで線形混合モデルを推定：blmerコマンド

# 茶園ごとに切片と施肥量の傾きが違うモデル（既定のウィシャート分布の事前分布）
fit1 <- blmer(yield ~ fert + (1 + fert | farm), data = tea)
fit1
Cov prior  : farm ~ wishart(df = 4.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev  : 6.2894

Linear mixed model fit by REML ['blmerMod']
Formula: yield ~ fert + (1 + fert | farm)
   Data: tea
REML criterion at convergence: 360.4929
Random effects:
 Groups   Name        Std.Dev. Corr
 farm     (Intercept) 45.7466
          fert         0.8954  -0.60
 Residual             16.3080
Number of obs: 40, groups:  farm, 8
Fixed Effects:
(Intercept)         fert
    416.831        2.969

# 施肥量の効き目（固定効果）
fixef(fit1)
(Intercept)        fert
  416.83118     2.96863

# 茶園間のばらつき（ランダム効果の標準偏差と相関）
VarCorr(fit1)
 Groups   Name        Std.Dev. Corr
 farm     (Intercept) 45.7466
          fert         0.8954  -0.603
 Residual             16.3080

# 切片だけが茶園ごとに違うモデル。茶園間の標準偏差にガンマ分布の事前分布を置く
fit2 <- blmer(yield ~ fert + (1 | farm), data = tea, cov.prior = gamma)
fit2
Cov prior  : farm ~ gamma(shape = 2.5, rate = 0, posterior.scale = sd, common.scale = TRUE)
Prior dev  : -2.1034

Linear mixed model fit by REML ['blmerMod']
Formula: yield ~ fert + (1 | farm)
   Data: tea
REML criterion at convergence: 360.4945
Random effects:
 Groups   Name        Std.Dev.
 farm     (Intercept) 36.11
 Residual             17.91
Number of obs: 40, groups:  farm, 8
Fixed Effects:
(Intercept)         fert
    419.632        2.883

＜おすすめのRに関する書籍です＞

データサイエンティストのための最新知識と実践 Rではじめよう! [モダン]なデータ分析

モダンな書き方でデータ分析を進めるための一冊。今どきのRの流儀を押さえたい人に。

www.amazon.co.jp

事前分布付きで一般化線形混合モデルを推定：bglmerコマンド

例は、九条ねぎの畑10か所で、殺菌剤を使った区と使わない区を3つずつ設け、各区20株のうちさび病が出た株を数えた架空のデータです。病気が出たかどうかの割合なので、二項分布のモデルを当てはめます。

事前分布の種類やパラメータの決め方は、パッケージのヘルプのbmerDist-classに一覧があります。

# 乱数の種を固定
set.seed(2026)

# 九条ねぎの畑10か所
field <- c("九条", "上鳥羽", "久世", "淀", "向日", "長岡京", "大山崎", "八幡", "久御山", "城陽")

# 畑ごとに6区。殺菌剤を使わない区は0、使った区は1
negi <- data.frame(
  field = rep(field, each = 6),
  fungicide = rep(c(0, 1), times = 30)
)

# 畑ごとの病気の出やすさと、殺菌剤の効き目から発病の確率を作る
p <- plogis(-0.5 + rep(rnorm(10, 0, 0.3), each = 6) - 1.5 * negi$fungicide)

# 各区20株のうち、さび病が出た株の数
negi$n <- 20
negi$sick <- rbinom(60, negi$n, p)

# 比べるため、lme4のglmerで当てはめる
fit_g0 <- glmer(cbind(sick, n - sick) ~ fungicide + (1 | field),
  data = negi, family = binomial
)

# 畑間の標準偏差が0近くに寄る
VarCorr(fit_g0)
 Groups Name        Std.Dev.
 field  (Intercept) 0.05154

# bglmerで当てはめる（既定のウィシャート分布の事前分布）
fit_g1 <- bglmer(cbind(sick, n - sick) ~ fungicide + (1 | field),
  data = negi, family = binomial
)

# 畑間の標準偏差が0から離れる
VarCorr(fit_g1)
 Groups Name        Std.Dev.
 field  (Intercept) 0.19065

# 固定効果にt分布の事前分布を置く
fit_g2 <- bglmer(cbind(sick, n - sick) ~ fungicide + (1 | field),
  data = negi, family = binomial, fixef.prior = t
)
fit_g2
Cov prior  : field ~ wishart(df = 3.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Fixef prior: t(df = 3, mean = c(0, 0), diag(scale) = c(100, 6.25), common.scale = FALSE)
Prior dev  : 15.833

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [bglmerMod]
 Family: binomial  ( logit )
Formula: cbind(sick, n - sick) ~ fungicide + (1 | field)
   Data: negi
      AIC       BIC    logLik -2*log(L)  df.resid
 229.9236  236.2067 -111.9618  223.9236        57
Random effects:
 Groups Name        Std.Dev.
 field  (Intercept) 0.1904
Number of obs: 60, groups:  field, 10
Fixed Effects:
(Intercept)    fungicide
    -0.6117      -1.7269

＜おすすめのRに関する書籍です＞

StanとRでベイズ統計モデリング Wonderful R 2 | 松浦健太郎, 石田基広

Stanによるベイズ統計モデリングの実践書。現実のデータ解析を念頭に、モデリングの考え方から丁寧に解説する。

この記事が誰かの役に立ちますように。