[R-lang] Differences between manual computation of confidence bands for logistic ME regression and 'effects' package?

[Jan Vanhove]

Dear all,

I've been trying out John Fox's 'effects' package for plotting (partial) 
fixed effects of logistic mixed-effects models. The effects package also 
allows the plotting of confidence bands, but by my reckoning, these seem 
to be rather narrow compared to the confidence bands that I computed 
manually following the advice on http://glmm.wikidot.com/faq. Could 
someone please explain to me where the difference comes from?

What follows is an example using a model with only one random effect 
term. My own data are modelled with crossed random intercepts and 
produce a much wider discrepancy between the two approaches. I'd be 
happy to provide the data and code for that example, too, if that would 
be helpful.

Thanks,
Jan

### Example 1 ###

# Analysing the 'Contraception' dataset in the mlmRev package.

# Load packages and data

library(lme4) # version 1.0-5

library(mlmRev)

data(Contraception)

# Fit ME model

mod1 <- glmer(use ~ age + (1|district)

               , data = Contraception

               , family = binomial)

# Compute 95% confidence band by hand

# Create new dataset

newdat1 <- expand.grid(age=seq(min(Contraception$age), max(Contraception$age), length.out=11)

                        , use = 0)

# Calculate effect of age

mm1 <- model.matrix(terms(mod1), newdat1)

newdat1$use <- mm1 %*% fixef(mod1)

#

pvar1 <- diag(mm1 %*% tcrossprod(vcov(mod1), mm1))

newdat1 <- data.frame(

   newdat1

   , plo = newdat1$use - qnorm(.975)*sqrt(pvar1)

   , phi = newdat1$use + qnorm(.975)*sqrt(pvar1)

)

# Plot effect + 95% CI

plot(newdat1$age

      , plogis(newdat1$use) # transform log-odds into probabilities

      , type="l"

      , lwd=3

      , ylim=c(0.25, 0.5))

points(newdat1$age

        , plogis(newdat1$plo)

        , type="l"

        , lty=2)

points(newdat1$age

        , plogis(newdat1$phi)

        , type="l"

        , lty=2)

# Using effects package

library(effects)

# Compute effects + 95% CB

eff1 <- effect('age'

                , mod1

                , xlevels=list(age=seq(min(Contraception$age), max(Contraception$age), length.out=11))

                , confidence.level=0.95)

# Plot

plot(eff1

      , ylim=c(0.25,0.50)

      , rescale.axis=FALSE)

# The difference between the two plots is subtle, but digging into the eff1 object, we find that the effect() function produces somewhat smaller CIs.

# Manual computation:

newdat1[5,]

#     age        use        plo        phi

# 5 -0.36 -0.5430525 -0.7083502 -0.3777548

# effect():

eff1$fit[5]

# -0.5430525 # i.e. same predicted value

eff1$fit[5] + qnorm(0.975) * eff1$se[5] # -0.4516169

eff1$fit[5] - qnorm(0.975) * eff1$se[5] # -0.6344881, i.e. narrower confidence band

-- 
Jan Vanhove
Mehrsprachigkeitsforschung und Fremdsprachendidaktik
Universität Freiburg, Schweiz
http://homeweb.unifr.ch/VanhoveJ/Pub/

Conference "Negotiating methodological challenges in linguistic research", February 2014, http://irg2014.org/