Как подобрать доверительные интервалы с помощью функции прогнозирования для glmmTMB
Я запускаю смешанную модель с использованием пакета glmmTMB и использую функцию прогнозирования для вычисления прогнозируемых средних значений, используя следующий код:
запустить модель
model_1 <- glmmTMB(Step.rate ~ Treatment*Week +
(1|Treatment.Group/Lamb.ID) + (1|Plot),
data = data.df, family = nbinom1)
создать новый фрейм данных
new.dat <- data.frame(Treatment = data.df$Treatment,
Week = data.df$Week, Plot = data.df$Plot,
Treatment.Group = data.df$Treatment.Group,
Lamb.ID = data.df$Lamb.ID)
прогнозировать средние значения
new.dat$prediction <- predict(model_1, new.data = new.dat,
type = "response", re.form = NA)
Этот код работает нормально, но когда я добавляю интервалы = "доверие" для вычисления доверительных интервалов, он, похоже, не работает. R игнорирует последнюю часть кода, и вычисляются только предсказанные средние.
new.dat$prediction <- predict(model_1, new.data = new.dat,
type = "response", re.form = NA, intervals = "confidence")
Почему интервалы = "уверенность" не работают? Может ли это быть проблема связанная с пакетом glmmTMB?
2 ответа
Вы можете использовать аргумент se.fit = TRUE чтобы получить стандартные ошибки предсказанных значений, а затем использовать их для вычисления доверительных интервалов.
https://www.rdocumentation.org/packages/glmmTMB/versions/1.0.2.1/topics/predict.glmmTMB
Есть несколько пакетов, которые делают работу за вас, например emmeans или ggeffects, или пакет эффектов (и, вероятно, еще несколько пакетов):
library(ggeffects)
library(glmmTMB)
library(emmeans)
data("Salamanders")
m <- glmmTMB(count ~ spp * mined + sample + (1 | site), family = nbinom1, data = Salamanders)
emmeans(m, c("spp", "mined"), type = "response")
#> spp mined response SE df lower.CL upper.CL
#> GP yes 0.0368 0.0373 627 0.00504 0.269
#> PR yes 0.1099 0.0661 627 0.03368 0.358
#> DM yes 0.3842 0.1397 627 0.18808 0.785
#> EC-A yes 0.1099 0.0660 627 0.03377 0.357
#> EC-L yes 0.3238 0.1222 627 0.15437 0.679
#> DES-L yes 0.4910 0.1641 627 0.25468 0.947
#> DF yes 0.5561 0.1764 627 0.29822 1.037
#> GP no 2.2686 0.4577 627 1.52646 3.372
#> PR no 0.4582 0.1515 627 0.23940 0.877
#> DM no 2.4201 0.4835 627 1.63472 3.583
#> EC-A no 0.8931 0.2373 627 0.53005 1.505
#> EC-L no 3.2017 0.6084 627 2.20451 4.650
#> DES-L no 3.4921 0.6517 627 2.42061 5.038
#> DF no 1.8495 0.3948 627 1.21623 2.813
#>
#> Confidence level used: 0.95
#> Intervals are back-transformed from the log scale
ggpredict(m, c("spp", "mined"))
#>
#> # Predicted counts of count
#> # x = spp
#>
#> # mined = yes
#>
#> x | Predicted | SE | 95% CI
#> --------------------------------------
#> GP | 0.04 | 1.01 | [0.01, 0.27]
#> PR | 0.11 | 0.60 | [0.03, 0.36]
#> DM | 0.38 | 0.36 | [0.19, 0.78]
#> EC-A | 0.11 | 0.60 | [0.03, 0.36]
#> EC-L | 0.32 | 0.38 | [0.15, 0.68]
#> DF | 0.56 | 0.32 | [0.30, 1.04]
#>
#> # mined = no
#>
#> x | Predicted | SE | 95% CI
#> --------------------------------------
#> GP | 2.27 | 0.20 | [1.53, 3.37]
#> PR | 0.46 | 0.33 | [0.24, 0.88]
#> DM | 2.42 | 0.20 | [1.64, 3.58]
#> EC-A | 0.89 | 0.27 | [0.53, 1.50]
#> EC-L | 3.20 | 0.19 | [2.21, 4.65]
#> DF | 1.85 | 0.21 | [1.22, 2.81]
#>
#> Adjusted for:
#> * sample = 2.50
#> * site = NA (population-level)
#> Standard errors are on the link-scale (untransformed).
Создано 14 сентября 2020 г. пакетом REPEX (v0.3.0)
Я думаю, что другой ответ дает вам обходной путь для получения CI для объектов glmmTMB с использованием se.fitаргумент. Но проблема наличия определенных версий функций для разных типов объектов (определяемых классом объекта) - это то, что вызывало у меня некоторое горе в прошлом, поэтому, возможно, стоит подробнее остановиться здесь.
Не вдаваясь в подробности, можно сказать, что функции в R, общие для многих типов объектов, имеют общие версии. Если вы дошли до документации по?predict, например, вы увидите страницу справки для общей версии функции. Там вы увидите несколько общих утверждений о том, как обычно работает функция, но мало или совсем не объясните конкретные аргументы, поскольку доступные аргументы зависят от типа объекта, с которым вы работаете. Описание со страницы справки для общегоpredict():
Прогноз - это общая функция для прогнозов на основе результатов различных функций подгонки модели. Функция вызывает определенные методы, которые зависят от класса первого аргумента.
Конкретные функции подгонки модели могут иметь определенные версииpredict()которые работают с результирующим объектом модели. Например, есть конкретныйpredict() для моделей, подходящих для lm(). Объекты, возвращенные изlm()относятся к классу lm. Вы можете увидеть документацию по версии функции для объектов lm по адресу?predict.lm. Именно эта функция содержитintervalsаргумент для расчета доверительных интервалов и интервалов прогноза. Хотя многие из нас начинают с объектов lm и поэтому учатсяintervals, оказывается, многие (большинство?) других predict() у функций нет этой опции.
Ключ для перехода на страницу справки по конкретным predict()Функция, которую вы используете, - это знать класс объекта, возвращаемого используемой функцией подгонки модели. Например, модели подходятglmmTMB()относятся к классу glmmTMB, поэтому вы можете перейти к?predict.glmmTMB. Модели подходятlme4::lmer()относятся к классу merMod, поэтому вы можете перейти к?predict.merMod. Если вы не знаете класс, возвращаемый функцией подгонки модели, похоже, вы часто можете найти информацию в документации в разделе " Значение ". Это верно дляlm() а также lmer(), как минимум.
Наконец, если вам нужно знать, имеет ли определенный класс объекта конкретную версию связанной с ним общей функции, вы можете посмотреть методы, доступные для этого класса, с помощью methods()функция. Пример для lm:
methods(class = "lm")
[1] add1 alias anova case.names coerce
[6] confint cooks.distance deviance dfbeta dfbetas
[11] drop1 dummy.coef effects extractAIC family
[16] formula hatvalues influence initialize kappa
[21] labels logLik model.frame model.matrix nobs
[26] plot predict print proj qqnorm
[31] qr residuals rstandard rstudent show
[36] simulate slotsFromS3 summary variable.names vcov