Доверительные интервалы для смешанных моделей с барьерным или нулевым завышением

Я хочу рассчитать CI в смешанных моделях, нулевой раздутый отрицательный бином и модель с препятствиями. Мой код для модели препятствий выглядит так (x1, x2 непрерывный, x3 категориальный):

m1 <- glmmTMB(count~x1+x2+x3+(1|year/class),
          data = bd, zi = ~x2+x3+(1|year/class), family = truncated_nbinom2,
          )

Я использовал confint и получил следующие результаты:

ci <- confint(m1,parm="beta_")
ci
                                          2.5 %       97.5 %     Estimate
cond.(Intercept)                   1.816255e-01  0.448860094  0.285524861
cond.x1                            9.045278e-01  0.972083366  0.937697401
cond.x2                            1.505770e+01 26.817439186 20.094998772
cond.x3high                        1.190972e+00  1.492335046  1.333164894
cond.x3low                         1.028147e+00  1.215828654  1.118056377
cond.x3reg                         1.135515e+00  1.385833853  1.254445909
class:year.cond.Std.Dev.(Intercept)2.256324e+00  2.662976154  2.441845815
year.cond.Std.Dev.(Intercept)      1.051889e+00  1.523719169  1.157153015
zi.(Intercept)                     1.234418e-04  0.001309705  0.000402085
zi.x2                              2.868578e-02  0.166378014  0.069084606
zi.x3high                          8.972025e-01  1.805832900  1.272869874

Правильно ли я вычисляю интервалы? Почему в x3 только одна категория для zi? Если возможно, я также хотел бы знать, возможно ли построить эти КИ.

Спасибо!

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

     class id year count  x1      x2        x3 
956    5 3002 2002      3 15.6    47.9      high
957    5 4004 2002      3 14.3    47.9      low
958    5 6021 2002      3 14.2    47.9      high
959    4 2030 2002      3 10.5    46.3      high
960    4 2031 2002      3 15.3    46.3      high
961    4 2034 2002      3 15.2    46.3      reg

с x1 и x2 непрерывными, x3 трехуровневой категориальной переменной (фактором)

Резюме модели:

summary(m1)
'giveCsparse' has been deprecated; setting 'repr = "T"' for you'giveCsparse' has been deprecated; setting 'repr = "T"' for you'giveCsparse' has been deprecated; setting 'repr = "T"' for you
 Family: truncated_nbinom2  ( log )
Formula:          count ~ x1 + x2 + x3 + (1 | year/class)
Zero inflation:          ~x2 + x3 + (1 | year/class)
Data: bd

     AIC      BIC   logLik deviance df.resid 
 37359.7  37479.7 -18663.8  37327.7    13323 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 class:year(Intercept) 0.79701  0.8928  
 year      (Intercept) 0.02131  0.1460  
Number of obs: 13339, groups:  class:year, 345; year, 15

Zero-inflation model:
 Groups    Name        Variance  Std.Dev. 
 dpto:year (Intercept) 1.024e+02 1.012e+01
 year      (Intercept) 7.842e-07 8.856e-04
Number of obs: 13339, groups:  class:year, 345; year, 15

Overdispersion parameter for truncated_nbinom2 family (): 1.02 

Conditional model:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -1.25343    0.23081  -5.431 5.62e-08 ***
x1             -0.06433    0.01837  -3.501 0.000464 ***
x2              3.00047    0.14724  20.378  < 2e-16 ***
x3high          0.28756    0.05755   4.997 5.82e-07 ***
x3low           0.11159    0.04277   2.609 0.009083 ** 
x3reg           0.22669    0.05082   4.461 8.17e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -7.8188     0.6025 -12.977  < 2e-16 ***
x2              -2.6724     0.4484  -5.959 2.53e-09 ***
x3high           0.2413     0.1784   1.352  0.17635    
x3low           -0.1325     0.1134  -1.169  0.24258    
x3reg           -0.3806     0.1436  -2.651  0.00802 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

CI с broom.mixed

> broom.mixed::tidy(m1, effects="fixed", conf.int=TRUE)
# A tibble: 12 x 9
   effect component term           estimate std.error statistic  p.value conf.low conf.high
   <chr>  <chr>     <chr>             <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
 1 fixed  cond      (Intercept)     -1.25      0.231      -5.43 5.62e- 8  -1.71     -0.801 
 2 fixed  cond      x1              -0.0643    0.0184     -3.50 4.64e- 4  -0.100    -0.0283
 3 fixed  cond      x2               3.00      0.147      20.4  2.60e-92   2.71      3.29  
 4 fixed  cond      x3high           0.288     0.0575      5.00 5.82e- 7   0.175     0.400 
 5 fixed  cond      x3low            0.112     0.0428      2.61 9.08e- 3   0.0278    0.195 
 6 fixed  cond      x3reg            0.227     0.0508      4.46 8.17e- 6   0.127     0.326 
 7 fixed  zi        (Intercept)     -9.88      1.32       -7.49 7.04e-14 -12.5      -7.30  
 8 fixed  zi        x1               0.214     0.120       1.79 7.38e- 2  -0.0206    0.448 
 9 fixed  zi        x2              -2.69      0.449      -6.00 2.01e- 9  -3.57     -1.81  
10 fixed  zi        x3high           0.232     0.178       1.30 1.93e- 1  -0.117     0.582 
11 fixed  zi        x3low           -0.135     0.113      -1.19 2.36e- 1  -0.357     0.0878
12 fixed  zi        x4reg           -0.382     0.144      -2.66 7.74e- 3  -0.664    -0.101 

person MKie45    schedule 12.04.2021    source источник
comment
(1) Да, насколько я могу судить. (2) трудно сказать без минимально воспроизводимого примера, я согласен, это кажется странным. (3) вы можете использовать broom.mixed::tidy(m1, effects="fixed", conf.int=TRUE) для получения данных в формате, удобном для ggplot2 (или даже dotwhisker::dwplot(m1, effects="fixed"))   -  person Ben Bolker    schedule 12.04.2021
comment
Я добавил образец данных и описание переменных. Проблема в том, что x3 имеет три категории, а confint сообщает только об одной. В любом случае, я попробую broom.mixed.   -  person MKie45    schedule 12.04.2021
comment
Я посмотрю, смогу ли я продвинуться дальше, но я подозреваю, что нет. Мне/нам нужен воспроизводимый пример (см. ссылку минимально воспроизводимый пример, а также stackoverflow.com/ questions/5963269/ для определения "воспроизводимого" и совета по построению MRE. (Скорее всего, проще всего будет дать ссылку на ваш полный набор данных...) . Можете ли вы также показать результаты summary(m1), а также любые предупреждающие сообщения, которые могли быть напечатаны... ?   -  person Ben Bolker    schedule 12.04.2021
comment
(1) Я не понимаю, как вы получаете коэффициенты x3 для высоких, низких и регулярных, если это трехуровневый фактор: что такое levels(bd$x3)? (2) в качестве перекрестной проверки, каковы результаты broom.mixed::tidy(bd, effects="fixed", conf.int=TRUE)?   -  person Ben Bolker    schedule 13.04.2021
comment
Я включил краткое описание модели. У меня была ошибка, x3 - четырехуровневая категориальная переменная. Я попробовал другую модель с категориальными и непрерывными переменными и с ziformula, и она также не сообщает о полных переменных в zi. В любом случае, я попытаюсь создать MRE для трехуровневой многоуровневой модели, а затем обновлю пост.   -  person MKie45    schedule 13.04.2021
comment
Я сделал один для вас.   -  person Ben Bolker    schedule 13.04.2021


Ответы (1)


tl;dr насколько я могу судить, это ошибка в confint.glmmTMB (и, возможно, во внутренней функции glmmTMB:::getParms). А пока broom.mixed::tidy(m1, effects="fixed") должен делать то, что вы хотите. (Сейчас в разрабатываемой версии на GitHub находится исправление, должно ли оно когда-нибудь появиться в CRAN? Скоро...)

Воспроизводимый пример:

настроить данные

set.seed(101)
n <- 1e3
bd <- data.frame(
    year=factor(sample(2002:2018, size=n, replace=TRUE)),
    class=factor(sample(1:20, size=n, replace=TRUE)),
    x1 = rnorm(n),
    x2 = rnorm(n),
    x3 = factor(sample(c("low","reg","high"), size=n, replace=TRUE),
                levels=c("low","reg","high")),
    count = rnbinom(n, mu = 3, size=1))

поместиться

library(glmmTMB)
m1 <- glmmTMB(count~x1+x2+x3+(1|year/class),
          data = bd, zi = ~x2+x3+(1|year/class), family = truncated_nbinom2,
          )

доверительные интервалы

confint(m1, "beta_")  ## wrong/ incomplete
broom.mixed::tidy(m1, effects="fixed", conf.int=TRUE)  ## correct

Вы можете подумать о том, какие доверительные интервалы вам нужны:

  • CI Wald (по умолчанию) намного быстрее вычисляются и, как правило, подходят, если (1) ваш набор данных велик и (2) вы не оцениваете какие-либо параметры в логарифмической/логарифмической шкале, которые находятся рядом с границами.
  • CI профиля вероятности более точны, но намного медленнее
person Ben Bolker    schedule 12.04.2021
comment
Я попробовал broom.mixed, и CI выглядит нормально. Извините, что не смог вовремя сделать MRE самостоятельно! Но спасибо большое, теперь мне стало понятнее. Я собираюсь загрузить в пост аккуратный вывод для перекрестной проверки. - person MKie45; 13.04.2021