Расчет доверительных интервалов для прогноза glmmTMB с нулевой инфляцией - PullRequest
0 голосов
/ 01 января 2019

Я пытаюсь вычислить доверительные интервалы прогнозов, сделанных на основе модели glmmTMB с нулевой инфляцией.

Я просмотрел несколько вопросов, опубликованных на github, и оригинальную статью, описывающую glmmTMB.Тем не менее, похоже, что glmmTMB немного изменился, и я очень озадачен тем, какой метод использовать.

Вот что я рассмотрел до сих пор: - https://github.com/glmmTMB/glmmTMB/issues/324 - https://github.com/glmmTMB/glmmTMB/issues/378 - https://www.biorxiv.org/content/biorxiv/suppl/2017/05/01/132753.DC1/132753-2.pdf

В настоящее время я просто беру прогноз и вычисляю его по шкале ответов (поскольку это должно включать и термин ZI), преобразовывая обратно в логарифмическую шкалу, вычисляя доверительный интервал 0,95и преобразование обратно в шкалу ответов.(см. код ниже).

library(glmmTMB)
data("Salamanders")


fit = glmmTMB(count~spp + cover + mined +(1|site), 
        ziformula=~spp + mined,
        dispformula = ~DOY,
        data = Salamanders, 
        family=nbinom2)


newdata <- expand.grid(
  site = "new",
  spp = unique(Salamanders$spp),
  mined = factor(c("yes", "no"), levels =  levels(Salamanders$mined)),
  cover = seq(from = min(Salamanders$cover),
              to = max(Salamanders$cover),
              length = 25
              ),
  DOY = mean(Salamanders$DOY)
  )

preds <- predict(fit, newdata, se=T, allow.new.levels = T, type='response')
newdata$pred=preds$fit
newdata$se = preds$se.fit
newdata$ulimit = exp(log(newdata$pred)+(qnorm(0.975)*log(newdata$se)))
newdata$llimit = exp(log(newdata$pred)-(qnorm(0.975)*log(newdata$se)))

library(ggplot2)

ggplot(data=newdata, aes(x=cover, y = pred))+
  geom_ribbon(aes(ymin=llimit, ymax=ulimit, fill=mined), alpha = 0.25)+
  geom_line(aes(color = mined), size=1)+
  facet_wrap(~spp)

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

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...