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

[T. Florian Jaeger]

Hi Jan,

I was trying to replicate, but the problem doesn't occur  on my system.
Here's the output from the last few lines of your script

> newdat1[5,]
    age        use        plo        phi
5 -0.36 -0.5430525 -0.7083502 -0.3777548

> eff1$fit[5]
[1] -0.5430525

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

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

As you can see I get the same CIs under both methods.

Here's what I'm running:

platform       x86_64-w64-mingw32
arch           x86_64
os             mingw32
system         x86_64, mingw32
status
major          3
minor          0.2
year           2013
month          09
day            25
svn rev        63987
language       R
version.string R version 3.0.2 (2013-09-25)
nickname       Frisbee Sailing


Effects: 2.3-0
lme4: 1.0-5






On Thu, Dec 19, 2013 at 11:14 PM, Jan Vanhove <jan.vanhove@unifr.ch> wrote:

> 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/
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucsd.edu/pipermail/ling-r-lang-l/attachments/20131222/90b59847/attachment.html