#r #predict #confidence-interval #glmmtmb
#r #прогнозировать #доверительный интервал #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)
Этот код работает нормально, но когда я добавляю в intervals = «confidence» для вычисления доверительных интервалов, кажется, что это не работает. R игнорирует последнюю часть кода, и вычисляются только предсказанные средние значения.
new.dat$prediction <- predict(model_1, new.data = new.dat,
type = "response", re.form = NA, intervals = "confidence")
Почему intervals = «confidence» не работает? Может ли это быть проблемой, связанной с пакетом glmmTMB?
Комментарии:
1. Несколько вещей, которые могут помочь: чтобы просмотреть страницу справки для конкретной версии glmmTMB
predict(), перейдите на?predict.glmmTMB. Вы увидите, что там нетintervalsаргумента. Вы можете сделать приблизительный CI «вручную». В часто задаваемых вопросах Болкера по GLMM здесь приведен пример. Я также использовалggpredict()из пакета ggeffects , чтобы получить приблизительный CI вместе с прогнозами, хотя я признаю, что, как известно, я немного разочаровывался, пытаясь определить точный набор данных, для которого я хочу получить прогнозы.2. @aosmith, на это, возможно, стоит ответить (особенно, если вы разъясните различие между универсальным методом (
?predict), наиболее широко используемым методом (?predict.lm) и методом, используемым операционной системой (?predict.glmmTMB) .) …3. для
?predict.glmmTMBу вас есть опцияse.fit=TRUE, которая значительно упрощает вашу работу по построению доверительных интервалов.
Ответ №1:
Вы можете использовать аргумент se.fit = TRUE , чтобы получить стандартные ошибки прогнозируемых значений, а затем использовать их для вычисления доверительных интервалов.
https://www.rdocumentation.org/packages/glmmTMB/versions/1.0.2.1/topics/predict.glmmTMB
Комментарии:
1. R чувствителен к регистру, это
TRUE, неTrue.2. Спасибо, это было действительно полезно!
Ответ №2:
Есть некоторые пакеты, которые выполняют эту работу за вас, такие как emmeans или ggeffects, или пакет effects (и, вероятно, еще несколько пакетов):
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).
Создано 2020-09-14 пакетом reprex (версия 0.3.0)
Ответ №3:
Я думаю, что другой ответ дает вам обходной путь для получения CI для объектов glmmTMB с использованием se.fit аргумента. Но проблема наличия определенных версий функций для разных типов объектов (определяемых классом объекта) — это то, что вызывало у меня некоторое огорчение в прошлом, поэтому, возможно, стоит подробнее остановиться на этом здесь.
Не вдаваясь в подробности, функции в R, которые являются общими для многих типов объектов, имеют общие версии. Например, если вы ознакомились с документацией для ?predict , вы увидите страницу справки для общей версии функции. Там вы увидите несколько общих утверждений о том, как обычно работает функция, но практически никаких объяснений конкретных аргументов, поскольку доступные аргументы зависят от типа объекта, с которым вы работаете. Описание со страницы справки для generic 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