с использованием профиля и метода загрузки в пределах опции confint, с моделью мерцания

Question

Я использую glmer с логической ссылкой для модели ошибок gaussian.

Когда я пытаюсь получить доверительные интервалы, используя либо профиль или метод загрузки с опцией confint, я получаю сообщение об ошибке для использования профилей вероятности, и с самозагрузкой:

> Profile: Computing profile confidence intervals ... Error in 
> names(opt) <- profnames(fm, signames) : 'names' attribute [2] must 
> be the same length as the vector [1] 
> 
> Boot: Error in if (const(t, min(1e-08, mean(t, na.rm = TRUE)/1e+06))) 
> { : missing value where TRUE/FALSE needed In addition: Warning 
> message: In bootMer(object, FUN = FUN, nsim = nsim, ...) : some 
> bootstrap runs failed (9999/10000) 

Я посмотрел на сайте предложения по устранению проблемы с профилем, путем установки инструмента разработки lme4, и я также устранил все НС из набора данных. Однако в обоих случаях я все равно получаю те же два сообщения об ошибках.

Может ли кто-нибудь предложить какую-либо помощь в отношении того, является ли это проблемой lme4 или является ли это более фундаментальным.

Вот код, чтобы произвести первые 20 наблюдения и формат модели:

df2 <- data.frame(
prop1 = c(0.46, 0.471, 0.458, 0.764, 0.742, 0.746, 0.569, 0.45, 0.491, 0.467, 0.464, 
     0.556, 0.584, 0.478, 0.456, 0.46, 0.493, 0.704, 0.693, 0.651), 
prop2 = c(0.458, 0.438, 0.453, 0.731, 0.738, 0.722, 0.613, 0.498, 0.452, 0.451, 0.447, 
     0.48, 0.576, 0.484, 0.473, 0.467, 0.467, 0.722, 0.707, 0.709), 
site = c(1,1,2,3,3,3,4,4,4,4,4,5,5,5,6,6,7,8,8,8) 
) 
df2$site <- factor(df2$site) 

model <- glmer(prop2 ~ prop1 + (1|site), 
      data=df2, family=gaussian(link="logit")) 
summary(model) 

Реакция является доля (0,1), как это ковариата. Я использовал ссылку logit, чтобы поддерживать ожидаемые значения ответа, связанные между (0,1), а не быть несвязанным с нормальным, хотя эти данные хорошо подходят для LMM.

Я также проанализирую разницу между двумя пропорциями, для которых я ожидаю, что некоторые отличия (и установленные значения) нарушат границу (0,1), поэтому я буду моделировать это как с линией связи с гауссовской ошибкой, или logit на масштабированной разности (принудительное выполнение этого параметра (0,1)).

Кроме того, учитывая, что на сайт может быть только 1-5 записей, я, естественно, рассмотрю возможность использования линейной регрессии или GLM (пропорции) и обрабатывать сайт как фиксированный эффект в процессе моделирования, если оценка случайного дисперсия эффектов равна нулю (или очень мала).

Answer 1

Вы обнаружили ошибку с текущей версией lme4, частично зафиксированной на Github (она работает с method="boot"). (devtools::install_github("lme4/lme4") должны работать, чтобы установить его, если у вас установлены инструменты разработки.)

library(lme4) 
fit_glmer <- glmer(prop2 ~ prop1 + (1|site), 
     data=df2, family=gaussian(link="logit")) 

Профилирование/профиль доверительные интервалы до сих пор не работают, но с более значимой ошибкой:

try(profile(fit_glmer)) 
## Error in profile.merMod(fit_glmer) : 
## can't (yet) profile GLMMs with non-fixed scale parameters 

Бутстрапирование делает Работа. Есть много предупреждений, и много попыток повторной установки не удалось, но я надеюсь, это из-за небольшого размера предоставленного набора данных.

bci <- suppressWarnings(confint(fit_glmer,method="boot",nsim=200)) 

Я хочу предложить несколько других вариантов. Вы можете использовать glmmADMB или glmmTMB, и эти платформы также позволяют использовать бета-распределение для пропорций модели.Я хотел бы рассмотреть вопрос о моделировании разницы между типами пропорций путем «плавления» данных (так что есть столбец prop и столбец prop_type) и включающий prop_type в качестве предиктора (возможно, с помощью случайного эффекта на индивидуальном уровне, определяющего, какие пропорции были измерены на то же лицо) ...

library(glmmADMB) 
fit_ADMB <- glmmadmb(prop2 ~ prop1 + (1|site), 
     data=df2, family="gaussian",link="logit") 
## warning message about 'sd.est' not defined -- only a problem 
## for computing standard errors of predictions, I think? 

library(glmmTMB) 
fit_TMB <- glmmTMB(prop2 ~ prop1 + (1|site), 
     data=df2, family=list(family="gaussian",link="logit")) 

Похоже, что ваши данные могут быть более подходящими для бета-модели?

fit_ADMB_beta <- glmmadmb(prop2 ~ prop1 + (1|site), 
     data=df2, family="beta",link="logit") 
fit_TMB_beta <- glmmTMB(prop2 ~ prop1 + (1|site), 
     data=df2, 
     family=list(family="beta",link="logit")) 

Сравните результаты (любитель, чем это должно быть)

library(broom) 
library(plyr) 
library(dplyr) 
library(dwplot) 

tab <- ldply(lme4:::namedList(fit_TMB_beta, 
         fit_ADMB_beta, 
         fit_ADMB, 
         fit_TMB, 
         fit_glmer), 
     tidy,.id="model") %>% 
    filter(term != "(Intercept)") 
dwplot(tab) + 
    facet_wrap(~term,scale="free", 
         ncol=1)+theme_set(theme_bw()) 

Answer 2

Для примера он работает с начальной загрузкой:

confint(model, method = "boot") 
#     2.5 %  97.5 % 
# .sig01  12.02914066 44.71708844 
# .sigma  0.03356588 0.07344978 
# (Intercept) -5.26207985 1.28669024 
# prop1  1.01574201 6.99804555 

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