using profile and boot method within confint option, with glmer model

Question

I am using glmer with a logit link for a gaussian error model.

When I try obtaining the confidence intervals, using either profile or the boot method with the confint option, I obtain an error for use of profile likelihood, and with bootstrapping:

> 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)

I have looked at online suggestions on how to overcome the profile issue, by installing the lme4 development tool and I have also eliminated all NAs from the dataset. However, in both instances, I still receive the same two error messages.

Is anybody able to offer any help as to whether this is a lme4 issue, or whether it's more fundamental.

Here is code to produce the first 20 observations and the model format:

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)

The response is a proportion (0,1), as is the covariate. I have used the logit link to keep the expected values of the response bound between (0,1), rather than being unbound with the normal, though this data fits a LMM reasonably well.

I will also analyse the difference between the two proportions, for which I expect some differences (and fitted values) to violate the (0,1) boundary - so I will model this with either an identity link with gaussian error, or logit link on the scaled difference (coercing this to be (0,1)).

Also, given that there may only be 1-5 records per site, I will naturally consider running linear regression or GLM (proportions), and treat site as a fixed effect in the modelling process, if the estimate of the random effects variance is zero (or very small).

Answer 1

For your example it works with bootstrap:

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

Take into consideration that under your proposed model, although your estimation will be always between 0 and 1, it is expected to observe values lower than 0 and greater than 1.

Answer 2

You've identified a bug with the current version of lme4, partially fixed on Github now (it works with method="boot"). (devtools::install_github("lme4/lme4") should work to install it, if you have development tools installed.)

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

Profiling/profile confidence intervals still don't work, but with a more meaningful error:

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

Bootstrapping does work. There are lot of warnings, and a lot of refitting attempts failed, but I'm hoping that's because of the small size of the data set provided.

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

I want to suggest a couple of other options. You can use glmmADMB or glmmTMB, and these platforms also allow you to use a Beta distribution to model proportions. I would consider modeling the difference between proportion types by "melting" the data (so that there is a prop column and a prop_type column) and including prop_type as a predictor (possibly with an individual-level random effect identifying which proportions were measured on the same individual) ...

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"))

It sounds like your data might be more appropriate for a Beta model?

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"))

Compare results (fancier than it needs to be)

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())