Включение смещения в модели Пуассона удвоенное количество моделей-кандидатов
Пересмотренный вопрос и код (10.07.2020)
В этот вопрос было внесено несколько редакций по мере продвижения анализа. Я включил самую последнюю информацию здесь, вверху, но для тех, кто застрял на других этапах, см. Информацию ниже.
Резюме: я моделирую отношения между несколькими непрерывными переменными (A, B и C), категориальной переменной (YN) и случайным перехватом (ID). Я включил смещение (T), потому что данные подсчета были собраны во время поведенческих наблюдений разной продолжительности. Также существует термин для нулевой инфляции, основанный на случайном пересечении. Я сравниваю все подмножества этой модели, чтобы определить наиболее подходящие модели.
В настоящее время функция земснаряда выдает ошибку, когда я пытаюсь сравнить все подмножества.
#simulated data
df <- data.frame(
count = sample(size = 300, c(0:6), replace = TRUE, prob = c(0.7, 0.15, 0.05, 0.04, 0.03, 0.02, 0.01)),
A = sample(size = 300, x = 24:65, replace = TRUE),
B = runif(n = 300, min = 7, max = 19),
C = rbeta(n = 300, shape1 = 1, shape2 = 25)*100,
YN = sample(size = 300, c("Y", "N"), replace = TRUE, prob = c(0.2, 0.8)),
T = runif(n = 300, min = 10, max = 20),
ID = sample(size = 300, letters, replace = TRUE)
)
#standardize explanatory variables
dfSTD <- stdize(df,
omit.cols = c("count", "YN", "T", "ID"),
center = TRUE, scale = TRUE)
#specify full model with standardized covariates for analysis
full <- stdizeFit(glmmTMB(count ~ z.A + z.B + z.C + YN + offset(log(T)) + (1|ID),
zi = ~ (1|ID),
family = poisson(),
data = dfSTD,
na.action = "na.fail"),
data = dfSTD,
which = "formula",
evaluate = FALSE)
#compare all possible model subsets
all <- dredge(full,
beta = "sd",
fixed = "cond(offset(log(T)))")
Error in nobs.default(global.model) : no 'nobs' method is available
Я работаю в R версии 3.6.1.
Я моделирую отношения между несколькими непрерывными переменными (A, B и C), категориальной переменной (YN) и случайным перехватом (ID). Я включил смещение (T), потому что данные подсчета были собраны во время поведенческих наблюдений разной продолжительности. Также существует термин для нулевой инфляции, основанный на случайном пересечении.
Я сравниваю все подмножества этой модели, чтобы определить наиболее подходящие модели.
Мой код выглядит следующим образом:
full <- glmmTMB(count ~ A + B + C + YN + offset(log(T)) + (1|ID),
zi = ~ (1|ID),
family = poisson(),
data = df,
na.action = "na.fail")
View(dredge(full, beta = sd))
Но когда смотрю на таблицу, кажется, что модели дублируются, т.е. со смещением и без него. Таблица включает 32 модели, а если не учитывать смещение, то всего 16 моделей. Плюс коэффициенты очень похожи. (показано подмножество этой таблицы)
Я попытался изменить свой код, чтобы смещение всегда было включено, как показано ниже, но это приводит к предупреждениям.
full <- glmmTMB(count ~ A + B + C + YN + (1|ID),
offset = log(T),
zi = ~ (1|ID),
family = poisson(),
data = df,
na.action = "na.fail")
Warning messages:
1: In (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluation
repeated 9 times
РЕДАКТИРОВАТЬ 2020/06/16
Я обновил этот пост с полным обзором модели, когда смещение включено как термин:
Family: poisson ( log )
Formula: count ~ A + B + C + YN + offset(log(T)) +
(1 | ID)
Zero inflation: ~(1 | ID)
Data: df
AIC BIC logLik deviance df.resid
286.4 317.4 -135.2 270.4 344
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
ID (Intercept) 0.1006 0.3172
Number of obs: 352, groups: name, 64
Zero-inflation model:
Groups Name Variance Std.Dev.
ID (Intercept) 2.089e-08 0.0001445
Number of obs: 352, groups: name, 64
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.582134 0.382926 -11.966 < 2e-16 ***
A 0.003517 0.170809 0.021 0.98357
B 0.356596 0.162703 2.192 0.02840 *
C 0.594737 0.095198 6.247 4.17e-10 ***
YN -1.397989 0.538510 -2.596 0.00943 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4744 0.7323 -0.648 0.517
И когда смещение включено как отдельный параметр вне формулы:
Family: poisson ( log )
Formula: count ~ A + B + C + YN + (1 | ID) #NOTE: offset isn't listed here despite being included when I specified the model above
Zero inflation: ~(1 | ID)
Data: df
Offset: log(duration)
AIC BIC logLik deviance df.resid
286.4 317.4 -135.2 270.4 344
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
ID (Intercept) 0.1006 0.3172
Number of obs: 352, groups: name, 64
Zero-inflation model:
Groups Name Variance Std.Dev.
ID (Intercept) 2.089e-08 0.0001445
Number of obs: 352, groups: name, 64
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.582134 0.382926 -11.966 < 2e-16 ***
A 0.003517 0.170809 0.021 0.98357
B 0.356596 0.162703 2.192 0.02840 *
C 0.594737 0.095198 6.247 4.17e-10 ***
YN -1.397989 0.538510 -2.596 0.00943 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4744 0.7323 -0.648 0.517
РЕДАКТИРОВАТЬ 2020/06/24
Пересмотренная функция земснаряда решает проблему с сохранением переменной смещения во всех моделях.
View(dredge(full,
beta = "sd"), #note sd was not in quotes above!
fixed = "cond(offset(log(time)))") #retains offset in all models
Однако это привело к новой проблеме:
Warning message:
In dredge(full, beta = "sd", fixed = "cond(offset(log(time)))") :
do not know how to standardize coefficients of 'glmmTMB', argument 'beta' ignored
@Kamil Barton предлагает ниже, что мне нужно использовать strike и stdizeFit для устранения этого предупреждающего сообщения. Но пока я не могу заставить этот код работать без другого сообщения об ошибке. Кроме того, я обеспокоен тем, что стандартизация моих переменных ответа (как предлагается ниже) - неправильный путь, и, поскольку стандартизация временной переменной приводит к отрицательным значениям, это больше не может регистрироваться в модели (дает значения NA).
Тем не менее, вот моя попытка реализовать предложенный Камилом подход с использованием почти такого же подхода, что и в файле справки:
dfSTD <- stdize(select(df, counts, A, B, C, YN, T, ID),
center = TRUE, scale = TRUE)
full <- stdizeFit(glmmTMB(z.counts ~ z.A + z.B + z.C + z.YN +
offset(log(z.T)) + (1|ID),
zi = ~ (1|ID),
family = poisson(),
data = dfSTD,
na.action = "na.fail"),
which = "formula", evaluate = FALSE)
Error in stdizeFit(glmmTMB(z.counts ~ z.A + z.B + z.C + z.YN + :
argument "data" is missing, with no default
Мы будем очень признательны за любые мысли по устранению этих последних ошибок, будь то подход, предложенный Камилом, или то, что я сделал выше.
РЕДАКТИРОВАТЬ 2020/07/10
Мне удалось решить указанную выше проблему с помощью следующего фрагмента кода:
dfSTD <- stdize(df, omit.cols = c("counts", "YN", "T", "ID"),
center = TRUE, scale = TRUE)
full <- stdizeFit(glmmTMB(counts ~ z.A + z.B + z.C + YN +
offset(log(T)) + (1|ID),
zi = ~ (1|ID),
family = poisson(),
data = dfSTD,
na.action = "na.fail"),
data = dfSTD, which = "formula", evaluate = FALSE)
Однако теперь, когда я использую функцию земснаряда, у меня возникает другая ошибка.
allEVCounts <- dredge(full,
beta = "sd",
fixed = "cond(offset(log(T)))")
Error in nobs.default(global.model) : no 'nobs' method is available
Любая помощь в определении того, как и где обратиться к методу "nobs", была бы оценена.
2 ответа
dredgeобрабатывает член смещения так же, как фиксированные члены модели, когда он включен в формулу. Чтобы сохранить смещение во всех моделях, добавьте аргументfixed = "<offset term name>", который в данном случае offset(log(T)):
dredge(full, beta = "sd", fixed = "cond(offset(log(T)))")
Обратите внимание, что в вашем коде beta = sd (без кавычек) рассматривается как beta = "none"потому что вы передаете функцию sd не символьная строка "sd" и это недопустимое значение.
Предупреждения об "оценке функции NA/NaN" исходят от glmmTMB, некоторые подмодели могут не сходиться со смещением. Задаватьoptions(warn = 1) (немедленные предупреждающие сообщения) и используйте dredge с участием trace = TRUE чтобы узнать какие.
Решение от Kamil Bartoń : в
stdizeFit, установленный
evaluate=TRUE , в противном случае он производит вызов, а не соответствующий объект модели, что не является тем, что вы (и
dredge) ожидал.