Estimating Multivariate Models with brms

Paul Bürkner

2022-12-14

イントロダクション

今回の vignette では、brms を用いた多変量マルチレベルモデルの指定方法について説明したい。複数の応答変数が含まれ、それぞれが予測変数のセットによって予測されるモデルを多変量と呼ぶ。生物学の例で考えてみよう。Hadfield, Nutall, Osorio, and Owens (2007)は、Eurasian blue tit (https://en.wikipedia.org/wiki/Eurasian_blue_tit)のデータを分析した。彼らはヒナの tarsus 長さと back 色を予測した。雛の半分を別の fosternest に入れ、残りの半分を自分の dam の養育舎にとどめた。これによって、遺伝的要因と環境的要因を分離することができる。さらに、ヒナの hatchdate と sex についての情報もある（後者は94%の動物についてわかっている）。

data("BTdata", package = "MCMCglmm")
head(BTdata)
       tarsus       back  animal     dam fosternest  hatchdate  sex
1 -1.89229718  1.1464212 R187142 R187557      F2102 -0.6874021  Fem
2  1.13610981 -0.7596521 R187154 R187559      F1902 -0.6874021 Male
3  0.98468946  0.1449373 R187341 R187568       A602 -0.4279814 Male
4  0.37900806  0.2555847 R046169 R187518      A1302 -1.4656641 Male
5 -0.07525299 -0.3006992 R046161 R187528      A2602 -1.4656641  Fem
6 -1.13519543  1.5577219 R187409 R187945      C2302  0.3502805  Fem

基本的な多変量解析モデル

まず、比較的単純な多変量正規モデルから始める。

bform1 <- 
  bf(mvbind(tarsus, back) ~ sex + hatchdate + (1|p|fosternest) + (1|q|dam)) +
  set_rescor(TRUE)

fit1 <- brm(bform1, data = BTdata, chains = 2, cores = 2)

モデルコードに見られるように、mvbind の表記で、次のように伝えている。

モデルで見たように、mvbind を使って brms に、tarsus と back の両方が個別の応答変数であることを伝えた。(1|p|fosternest) という用語は、fosternest の切片が変化していることを示す。間に |p| を書くことで、fosternest のすべての変化する効果が相関するようにモデル化されるべきであることを示す。これは、実際には2つのモデル部分、1つは tarsus について、もう1つは back についてなので、理にかなっている。p という指標は任意であり、思いついた他の記号で置き換えることができる（brms のマルチレベル構文の詳細については、help("brmsformula") と vignette("brms_multilevel") を参照）。同様に、(1|q|dam) という用語は、雛の遺伝的母体による相関的な変動効果を示している。あるいは、血統とそれに対応する関連性行列によって遺伝的類似性をモデル化することもできましたが、これはこの vignette の焦点ではない（ vignette("brms_phylogenetics") を参照）。モデルの結果は次のように簡単にまとめることができる。

fit1 <- add_criterion(fit1, "loo")
summary(fit1)
 Family: MV(gaussian, gaussian) 
  Links: mu = identity; sigma = identity
         mu = identity; sigma = identity 
Formula: tarsus ~ sex + hatchdate + (1 | p | fosternest) + (1 | q | dam) 
         back ~ sex + hatchdate + (1 | p | fosternest) + (1 | q | dam) 
   Data: BTdata (Number of observations: 828) 
  Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 2000

Group-Level Effects: 
~dam (Number of levels: 106) 
                                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(tarsus_Intercept)                     0.48      0.05     0.39     0.58 1.00      854
sd(back_Intercept)                       0.25      0.08     0.09     0.39 1.01      271
cor(tarsus_Intercept,back_Intercept)    -0.51      0.23    -0.94    -0.05 1.00      474
                                     Tail_ESS
sd(tarsus_Intercept)                     1463
sd(back_Intercept)                        585
cor(tarsus_Intercept,back_Intercept)      698

~fosternest (Number of levels: 104) 
                                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(tarsus_Intercept)                     0.27      0.05     0.17     0.37 1.00      827
sd(back_Intercept)                       0.35      0.06     0.24     0.46 1.00      584
cor(tarsus_Intercept,back_Intercept)     0.69      0.21     0.21     0.98 1.00      299
                                     Tail_ESS
sd(tarsus_Intercept)                     1085
sd(back_Intercept)                       1148
cor(tarsus_Intercept,back_Intercept)      533

Population-Level Effects: 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
tarsus_Intercept    -0.41      0.07    -0.53    -0.28 1.00     1905     1798
back_Intercept      -0.01      0.06    -0.14     0.12 1.00     2824     1602
tarsus_sexMale       0.77      0.06     0.66     0.88 1.00     4353     1741
tarsus_sexUNK        0.23      0.13    -0.01     0.49 1.00     4032     1494
tarsus_hatchdate    -0.04      0.06    -0.16     0.07 1.00     1587     1231
back_sexMale         0.01      0.07    -0.13     0.14 1.00     4343     1605
back_sexUNK          0.15      0.15    -0.15     0.43 1.00     4015     1594
back_hatchdate      -0.09      0.05    -0.20     0.01 1.00     2142     1792

Family Specific Parameters: 
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_tarsus     0.76      0.02     0.72     0.80 1.00     2913     1659
sigma_back       0.90      0.02     0.86     0.95 1.00     2355     1575

Residual Correlations: 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(tarsus,back)    -0.05      0.04    -0.13     0.02 1.00     2375     1382

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).

多変量モデルの要約出力は、パラメータが対応する応答変数を接頭辞として持つことを除けば、一変量モデルの出力とよく似ている。ダムでは、足根長さと背中の色は負の相関があるように見えるが、フォスターネストではその逆が真である。これは、この2つの特性に対する遺伝的および環境的要因の影響の差を示している。さらに、出力の底にある小さな残差相関 rescor(tarsus, back) は、足根長さと背の色の間にモデル化されていない依存関係がほとんどないことを示している。この時点では必要ないが、私たちはすでに LOO 情報量基準 fit1 を計算し、保存している。これはモデル比較のために使用する予定である。次に、モデル適合の第一印象を与える、いくつかの事後予測チェックを見てみよう。

pp_check(fit1, resp = "tarsus")

pp_check(fit1, resp = "back")

これはかなり強固に見えるが、tarsus の分布にわずかにモデル化されていない左歪度があることに気づく。 このことについては後でまた説明する。次に、私たちは応答変数のどれくらいの変動が私たちのモデルで説明できるかを調査したいので、 \(R^2\) 係数のベイズ一般化を使用する。

bayes_R2(fit1)
          Estimate  Est.Error      Q2.5     Q97.5
R2tarsus 0.4345953 0.02355886 0.3885075 0.4792487
R2back   0.1976980 0.02713929 0.1436405 0.2510507

明らかに両者の動物特性には説明できない多くのバリエーションがあるが、どうやら背中の色よりも足根管の長さのバリエーションの方が説明できそうだ。

より複雑な多変量解析モデル

さて、tarsus では sex のみを制御し、back では制御せず、hatchdate ではその逆を行いたいとする。これは今の例では特に合理的ではないが、異なる応答変数に異なる数式を指定する方法を説明することができる。mvbind の構文はもう使えないので、もっと冗長なアプローチを使わなければならない。

bf_tarsus <- bf(tarsus ~ sex + (1|p|fosternest) + (1|q|dam))
bf_back <- bf(back ~ hatchdate + (1|p|fosternest) + (1|q|dam))
fit2 <- brm(bf_tarsus + bf_back + set_rescor(TRUE), 
            data = BTdata, chains = 2, cores = 2)

この場合、+ 演算子によって、2つのモデル部分を文字通り追加していることに注意。これは、mvbf(bf_tarsus, bf_back) と書くことと同じである。この構文の詳細については、help("brmsformula") と help("mvbrmsformula") を参照。ここでも、まずモデルを要約する。

fit2 <- add_criterion(fit2, "loo")
summary(fit2)
 Family: MV(gaussian, gaussian) 
  Links: mu = identity; sigma = identity
         mu = identity; sigma = identity 
Formula: tarsus ~ sex + (1 | p | fosternest) + (1 | q | dam) 
         back ~ hatchdate + (1 | p | fosternest) + (1 | q | dam) 
   Data: BTdata (Number of observations: 828) 
  Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 2000

Group-Level Effects: 
~dam (Number of levels: 106) 
                                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(tarsus_Intercept)                     0.48      0.05     0.39     0.58 1.00      795
sd(back_Intercept)                       0.25      0.08     0.10     0.39 1.00      283
cor(tarsus_Intercept,back_Intercept)    -0.49      0.22    -0.91    -0.07 1.00      368
                                     Tail_ESS
sd(tarsus_Intercept)                     1297
sd(back_Intercept)                        441
cor(tarsus_Intercept,back_Intercept)      636

~fosternest (Number of levels: 104) 
                                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(tarsus_Intercept)                     0.27      0.05     0.17     0.37 1.01      561
sd(back_Intercept)                       0.35      0.06     0.23     0.46 1.00      436
cor(tarsus_Intercept,back_Intercept)     0.68      0.21     0.21     0.97 1.00      172
                                     Tail_ESS
sd(tarsus_Intercept)                      991
sd(back_Intercept)                        752
cor(tarsus_Intercept,back_Intercept)      631

Population-Level Effects: 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
tarsus_Intercept    -0.41      0.07    -0.54    -0.28 1.00      842     1068
back_Intercept      -0.00      0.06    -0.12     0.11 1.00     1175     1420
tarsus_sexMale       0.77      0.06     0.65     0.89 1.00     2113     1485
tarsus_sexUNK        0.23      0.13    -0.02     0.49 1.00     2107     1335
back_hatchdate      -0.09      0.05    -0.19     0.01 1.00     1591     1452

Family Specific Parameters: 
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_tarsus     0.76      0.02     0.71     0.80 1.00     1708     1447
sigma_back       0.90      0.03     0.85     0.95 1.00     1794     1366

Residual Correlations: 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(tarsus,back)    -0.05      0.04    -0.13     0.02 1.00     2096     1235

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).

ここでは、初期モデルから特定の効果を除外することで、モデルの適合度がどのように変化するかを見てみよう。

loo(fit1, fit2)
Output of model 'fit1':

Computed from 2000 by 828 log-likelihood matrix

         Estimate   SE
elpd_loo  -2126.3 33.7
p_loo       175.9  7.5
looic      4252.5 67.5
------
Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     810   97.8%   243       
 (0.5, 0.7]   (ok)        15    1.8%   112       
   (0.7, 1]   (bad)        3    0.4%   19        
   (1, Inf)   (very bad)   0    0.0%   <NA>      
See help('pareto-k-diagnostic') for details.

Output of model 'fit2':

Computed from 2000 by 828 log-likelihood matrix

         Estimate   SE
elpd_loo  -2124.5 33.7
p_loo       174.3  7.4
looic      4249.0 67.5
------
Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     814   98.3%   88        
 (0.5, 0.7]   (ok)        13    1.6%   129       
   (0.7, 1]   (bad)        1    0.1%   34        
   (1, Inf)   (very bad)   0    0.0%   <NA>      
See help('pareto-k-diagnostic') for details.

Model comparisons:
     elpd_diff se_diff
fit2  0.0       0.0   
fit1 -1.8       1.3

どうやら、モデルの適合度には特筆すべき差はないようである。したがって、sex と hatchdate を両方の応答変数についてモデル化する必要は本当にないが、それらを含めても害はない（だから、私はそれらを含めるだけだろう）。

brms の多変量解析シンタックスの能力を垣間見るために、同時に様々な方向にモデルを変えてみよう。tarsus のわずかな左歪度を思い出してみよう。今度は、gaussian の代わりに skew_normal の系列を使用してモデル化する。多変量正規（またはstudent-t）モデルではないので、残差相関の推定はもはや不可能である。このことは、set_rescor 関数を用いて明示す。さらに、hatchdate の非線形スプラインをフィッティングして、back と hatchdate の関係が本当に以前仮定したように線形かどうかを調べる。その上で、tarsus の残差分散をオスとメスで別々にモデル化する。

bf_tarsus <- bf(tarsus ~ sex + (1|p|fosternest) + (1|q|dam)) +
  lf(sigma ~ 0 + sex) + skew_normal()
bf_back <- bf(back ~ s(hatchdate) + (1|p|fosternest) + (1|q|dam)) +
  gaussian()

fit3 <- brm(
  bf_tarsus + bf_back + set_rescor(FALSE),
  data = BTdata, chains = 2, cores = 2,
  control = list(adapt_delta = 0.95)
)

ここでも、モデルを要約し、いくつかの事後予測的なチェックを見る。

fit3 <- add_criterion(fit3, "loo")
summary(fit3)
 Family: MV(skew_normal, gaussian) 
  Links: mu = identity; sigma = log; alpha = identity
         mu = identity; sigma = identity 
Formula: tarsus ~ sex + (1 | p | fosternest) + (1 | q | dam) 
         sigma ~ 0 + sex
         back ~ s(hatchdate) + (1 | p | fosternest) + (1 | q | dam) 
   Data: BTdata (Number of observations: 828) 
  Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 2000

Smooth Terms: 
                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sds(back_shatchdate_1)     2.05      1.07     0.34     4.70 1.01      412      214

Group-Level Effects: 
~dam (Number of levels: 106) 
                                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(tarsus_Intercept)                     0.47      0.05     0.38     0.58 1.00      653
sd(back_Intercept)                       0.24      0.07     0.11     0.37 1.02      238
cor(tarsus_Intercept,back_Intercept)    -0.51      0.22    -0.91    -0.05 1.01      307
                                     Tail_ESS
sd(tarsus_Intercept)                      803
sd(back_Intercept)                        536
cor(tarsus_Intercept,back_Intercept)      365

~fosternest (Number of levels: 104) 
                                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(tarsus_Intercept)                     0.27      0.06     0.16     0.37 1.00      411
sd(back_Intercept)                       0.31      0.06     0.20     0.42 1.00      431
cor(tarsus_Intercept,back_Intercept)     0.65      0.23     0.17     0.98 1.01      156
                                     Tail_ESS
sd(tarsus_Intercept)                      916
sd(back_Intercept)                        918
cor(tarsus_Intercept,back_Intercept)      377

Population-Level Effects: 
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
tarsus_Intercept        -0.41      0.07    -0.54    -0.28 1.00      665      641
back_Intercept           0.00      0.05    -0.10     0.10 1.00      897      988
tarsus_sexMale           0.77      0.06     0.66     0.88 1.00     2124     1550
tarsus_sexUNK            0.22      0.12    -0.02     0.46 1.00     2043     1321
sigma_tarsus_sexFem     -0.30      0.04    -0.39    -0.21 1.00     1629     1604
sigma_tarsus_sexMale    -0.24      0.04    -0.32    -0.17 1.00     1756     1567
sigma_tarsus_sexUNK     -0.39      0.13    -0.63    -0.15 1.00     1804     1636
back_shatchdate_1       -0.34      3.42    -6.33     7.39 1.00      702      769

Family Specific Parameters: 
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_back       0.90      0.02     0.85     0.95 1.00     1844     1580
alpha_tarsus    -1.21      0.45    -1.93     0.10 1.00     1077      547

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).

tarsus の（対数）残差標準偏差は、雌雄の区別がつかないヒナでは、雄や雌のヒナに比べ、やや大きいことがわかる。さらに、tarsus の負の alpha (歪度)パラメータから、残差は確かに少し左巻きであることがわかる。最後に

conditional_effects(fit3, "hatchdate", resp = "back")

を見ると、hatchdate on back、孵化日の経過とともに波状に変化しているような非線形な関係があることがわかる。

多変量解析モデルには、この vignette では説明しないが、さらに多くのモデリング・オプションがある。例えば、自己相関構造、ガウス過程、明示的な非線形予測変数（例えば、help("brmsformula") または vignette("brms_multilevel") を参照）などがある。実際、一変量モデルの柔軟性は、多変量モデルでもほぼ全て保持されている。

参考文献

Hadfield JD, Nutall A, Osorio D, Owens IPF (2007). Testing the phenotypic gambit: phenotypic, genetic and environmental correlations of colour. Journal of Evolutionary Biology, 20(2), 549-557.