Извлечение апостериорных мод и достоверных интервалов из вывода glmmTMB

Question

Я обычно работаю с lme4 пакет, но glmmTMB Пакет становится все более подходящим для работы с очень сложными данными (подумайте о чрезмерной дисперсии и / или нулевой инфляции).

Есть ли способ извлечь апостериорные моды и достоверные интервалы из glmmTMBмодели, аналогично тому, как это делается дляlme4модели (пример здесь).

Детали:

Я работаю с данными подсчета (доступны здесь), которые не завышены и диспергированы и имеют случайные эффекты. Пакет, который лучше всего подходит для работы с такими данными, - этоglmmTMB(подробности здесь). (Обратите внимание на два выброса:euc0==78 а также np_other_grass==20).

Данные выглядят так:

euc0 ea_grass ep_grass np_grass np_other_grass month year precip season   prop_id quad
 3      5.7      0.0     16.7            4.0     7 2006    526 Winter    Barlow    1
 0      6.7      0.0     28.3            0.0     7 2006    525 Winter    Barlow    2
 0      2.3      0.0      3.3            0.0     7 2006    524 Winter    Barlow    3
 0      1.7      0.0     13.3            0.0     7 2006    845 Winter    Blaber    4
 0      5.7      0.0     45.0            0.0     7 2006    817 Winter    Blaber    5
 0     11.7      1.7     46.7            0.0     7 2006    607 Winter    DClark    3

В glmmTMB модель:

model<-glmmTMB(euc0 ~ ea_grass + ep_grass + np_grass + np_other_grass + (1|prop_id), data = euc, family= nbinom2) #nbimom2 lets var increases quadratically
summary(model)
confint(model) #this gives the confidence intervals

Как я обычно извлекаю апостериорную моду и достоверные интервалы для lmer/glmer модель:

#extracting model estimates and credible intervals
sm.model <-arm::sim(model, n.sim=1000)
smfixef.model = sm.model@fixef
smfixef.model =coda::as.mcmc(smfixef.model)
MCMCglmm::posterior.mode(smfixef.model)  #mode of the distribution
coda::HPDinterval(smfixef.model)  #credible intervals

#among-brood variance
bid<-sm.model@ranef$prop_id[,,1]
bvar<-as.vector(apply(bid, 1, var)) #between brood variance posterior distribution
bvar<-coda::as.mcmc(bvar)
MCMCglmm::posterior.mode(bvar) #mode of the distribution
coda::HPDinterval(bvar) #credible intervals

Answer 1

Большая часть ответа:

1. Получить многомерную Нормальную выборку параметров условной модели довольно просто (я думаю, что это то, чтоarm::sim() делается.

library(MASS)
pp <- fixef(model)$cond
vv <- vcov(model)$cond
samp <- MASS::mvrnorm(1000, mu=pp, Sigma=vv)

(затем используйте оставшуюся часть описанного выше метода).

2. Я немного скептически отношусь к тому, что ваш второй пример делает то, что вы хотите. Дисперсия условных режимов не обязательно является хорошей оценкой дисперсии между группами (например, см. Здесь). Кроме того, я нервничаю из-за полузависимого байесовского подхода (например, почему нет априорных значений? Зачем смотреть на апостериорную моду, которая редко является значимым значением в байесовском контексте?), Хотя я сам иногда использую аналогичные подходы!) Однако использовать результаты glmmTMB для правильного анализа цепей Маркова методом Монте-Карло не так уж сложно:

library(tmbstan)
library(rstan)
library(coda)
library(emdbook) ## for lump.mcmc.list(), or use runjags::combine.mcmc()

t2 <- system.time(m2 <- tmbstan(model$obj))
m3 <- rstan::As.mcmc.list(m2)
lattice::xyplot(m3,layout=c(5,6))
m4 <- emdbook::lump.mcmc.list(m3)
coda::HPDinterval(m4)

Может быть полезно знать, что theta столбец m4 представляет собой журнал стандартного отклонения между группами...

(Видеть vignette("mcmc", package="glmmTMB") для получения дополнительной информации...)

Answer 2

Я думаю, что Бен уже ответил на ваш вопрос, поэтому мой ответ не добавляет особого смысла в обсуждение... Может быть, только одно, как вы писали в своих комментариях, что вас интересуют различия внутри и между группами. Вы можете получить эту информацию черезparameters::random_parameters()(если я правильно понял, что вы искали). См. Пример ниже, который сначала генерирует смоделированные выборки из многомерной нормы (как в примере с Беном), а затем дает вам сводку дисперсий случайного эффекта...

library(readr)
library(glmmTMB)
library(parameters)
library(bayestestR)
library(insight)

euc_data <- read_csv("D:/Downloads/euc_data.csv")
model <-
  glmmTMB(
    euc0 ~ ea_grass + ep_grass + np_grass + np_other_grass + (1 | prop_id),
    data = euc_data,
    family = nbinom2
  ) #nbimom2 lets var increases quadratically


# generate samples
samples <- parameters::simulate_model(model)
#> Model has no zero-inflation component. Simulating from conditional parameters.


# describe samples
bayestestR::describe_posterior(samples)
#> # Description of Posterior Distributions
#> 
#> Parameter      | Median |           89% CI |    pd |        89% ROPE | % in ROPE
#> --------------------------------------------------------------------------------
#> (Intercept)    | -1.072 | [-2.183, -0.057] | 0.944 | [-0.100, 0.100] |     1.122
#> ea_grass       | -0.001 | [-0.033,  0.029] | 0.525 | [-0.100, 0.100] |   100.000
#> ep_grass       | -0.050 | [-0.130,  0.038] | 0.839 | [-0.100, 0.100] |    85.297
#> np_grass       | -0.020 | [-0.054,  0.012] | 0.836 | [-0.100, 0.100] |   100.000
#> np_other_grass | -0.002 | [-0.362,  0.320] | 0.501 | [-0.100, 0.100] |    38.945


# or directly get summary of sample description
sp <- parameters::simulate_parameters(model, ci = .95, ci_method = "hdi", test = c("pd", "p_map"))
sp
#> Model has no zero-inflation component. Simulating from conditional parameters.
#> # Description of Posterior Distributions
#> 
#> Parameter      | Coefficient | p_MAP |    pd |              CI
#> --------------------------------------------------------------
#> (Intercept)    |      -1.037 | 0.281 | 0.933 | [-2.305, 0.282]
#> ea_grass       |      -0.001 | 0.973 | 0.511 | [-0.042, 0.037]
#> ep_grass       |      -0.054 | 0.553 | 0.842 | [-0.160, 0.047]
#> np_grass       |      -0.019 | 0.621 | 0.802 | [-0.057, 0.023]
#> np_other_grass |       0.019 | 0.999 | 0.540 | [-0.386, 0.450]

plot(sp) + see::theme_modern()
#> Model has no zero-inflation component. Simulating from conditional parameters.

# random effect variances
parameters::random_parameters(model)
#> # Random Effects
#> 
#> Within-Group Variance         2.92 (1.71)
#> Between-Group Variance
#>   Random Intercept (prop_id)   2.1 (1.45)
#> N (groups per factor)
#>   prop_id                       18
#> Observations                   346

insight::get_variance(model)
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may be unreliable.
#> $var.fixed
#> [1] 0.3056285
#> 
#> $var.random
#> [1] 2.104233
#> 
#> $var.residual
#> [1] 2.91602
#> 
#> $var.distribution
#> [1] 2.91602
#> 
#> $var.dispersion
#> [1] 0
#> 
#> $var.intercept
#>  prop_id 
#> 2.104233

^{Создано 2020-05-26 пакетом REPEX (v0.3.0)}