Effects package provides a very fast and convenient way for plotting linear mixed effect model results obtained through lme4 package. The effect function calculates confidence intervals (CIs) very quickly, but how trustworthy are these confidence intervals?

For example:



fm1  <- lmer(strength ~ batch + (1 | cask), Pastes)
effs <- as.data.frame(effect(c("batch"), fm1))
ggplot(effs, aes(x = batch, y = fit, ymin = lower, ymax = upper)) + 
  geom_rect(xmax = Inf, xmin = -Inf, ymin = effs[effs$batch == "A", "lower"],
        ymax = effs[effs$batch == "A", "upper"], alpha = 0.5, fill = "grey") +
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

enter image description here

According to CIs calculated using effects package, batch "E" does not overlap with batch "A".

If I try the same using confint.merMod function and the default method:

a <- fixef(fm1)
b <- confint(fm1)
# Computing profile confidence intervals ...
# There were 26 warnings (use warnings() to see them)

b <- data.frame(b)
b <- b[-1:-2,]

b1 <- b[[1]]
b2 <- b[[2]]

dt <- data.frame(fit   = c(a[1],  a[1] + a[2:length(a)]), 
                 lower = c(b1[1],  b1[1] + b1[2:length(b1)]), 
                 upper = c(b2[1],  b2[1] + b2[2:length(b2)]) )
dt$batch <- LETTERS[1:nrow(dt)]

ggplot(dt, aes(x = batch, y = fit, ymin = lower, ymax = upper)) +
  geom_rect(xmax = Inf, xmin = -Inf, ymin = dt[dt$batch == "A", "lower"], 
        ymax = dt[dt$batch == "A", "upper"], alpha = 0.5, fill = "grey") + 
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

enter image description here

I see that all of the CIs overlap. I also get warnings indicating that the function failed to calculate trustworthy CIs. This example, and my actual dataset, makes me to suspect that effects package takes shortcuts in CI calculation that might not entirely be approved by statisticians. How trustworthy are the CIs returned by effect function from effects package for lmer objects?

What have I tried: Looking into the source code, I noticed that effect function relies on Effect.merMod function, which in turn directs to Effect.mer function, which looks like this:

function (focal.predictors, mod, ...) 
    result <- Effect(focal.predictors, mer.to.glm(mod), ...)
    result$formula <- as.formula(formula(mod))
<environment: namespace:effects>

mer.to.glm function seems to calculate Variance-Covariate Matrix from the lmerobject:


function (mod) 
mod2$vcov <- as.matrix(vcov(mod))

This, in turn, is probably used in Effect.default function to calculate CIs (I might have misunderstood this part):

     z <- qnorm(1 - (1 - confidence.level)/2)
        V <- vcov.(mod)
        eff.vcov <- mod.matrix %*% V %*% t(mod.matrix)
        rownames(eff.vcov) <- colnames(eff.vcov) <- NULL
        var <- diag(eff.vcov)
        result$vcov <- eff.vcov
        result$se <- sqrt(var)
        result$lower <- effect - z * result$se
        result$upper <- effect + z * result$se

I do not know enough about LMMs to judge whether this is a right approach, but considering the discussion around confidence interval calculation for LMMs, this approach appears suspiciously simple.

  • 1
    When you have long lines of code, I would greatly appreciate it if you break them into multiple lines so we don't have to scroll to see it all. – rvl Oct 2 '14 at 17:30
  • 1
    @rvl The code should be more readable now. – Mikko Oct 3 '14 at 7:48
up vote 48 down vote accepted

All of the results are essentially the same (for this particular example). Some theoretical differences are:

  • as @rvl points out, your reconstruction of CIs without taking account of covariance among parameters is just wrong (sorry)
  • confidence intervals for parameters can be based on Wald confidence intervals (assuming a quadratic log-likelihood surface): lsmeans, effects, confint(.,method="Wald"); except for lsmeans, these methods ignore finite-size effects ("degrees of freedom"), but in this case it barely makes any difference (df=40 is practically indistinguishable from infinite df)
  • ... or on profile confidence intervals (the default method; ignores finite-size effects but allows for non-quadratic surfaces)
  • ... or on parametric bootstrapping (the gold standard -- assumes the model is correct [responses are Normal, random effects are Normally distributed, data are conditionally independent, etc.], but otherwise makes few assumptions)

I think all of these approaches are reasonable (some are more approximate than others), but in this case it barely makes any difference which one you use. If you're concerned, try out several contrasting methods on your data, or on simulated data that resemble your own, and see what happens ...

(PS: I wouldn't put too much weight on the fact that the confidence intervals of A and E don't overlap. You'd have to do a proper pairwise comparison procedure to make reliable inferences about the differences between this particular pair of estimates ...)

95% CIs:

enter image description here

Comparison code:

fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
c0 <- confint(fm2,method="Wald")
c1 <- confint(fm2)
c2 <- confint(fm2,method="boot")
c3 <- with(effect("batch",fm2),cbind(lower,upper))
c4 <- with(summary(lsmeans(fm2,spec="batch")),cbind(lower.CL,upper.CL))
tmpf <- function(method,val) {
library(ggplot2); theme_set(theme_bw())
allCI <- rbind(tmpf("lme4_wald",c0),

  • 2
    I accept this answer as it is right to the point and gives a nice comparison between different methods. However, check out rlv's excellent answer for more information. – Mikko Oct 3 '14 at 17:22
  • Thank you for excellent and very helpful answer. Do I understand correctly that one cannot use CIs to compare groups/batches, but it is possible to compare effects. Say that I had two treatments, several individuals and several measurements within individuals. I would use individuals as random effect as each of them would contain x measurements. Then I wanted to know whether these two treatments resulted to a different response. Could I use effects package and CI overlap in this case? – Mikko Oct 3 '14 at 17:41
  • 4
    This is a more general question which is relevant to any standard model-based approach. Might be worth a separate question. (1) In general the way one answers questions about differences between treatments is to set up the model so that the difference between the focal treatments is a contrast (i.e., an estimated parameter) in the model, and then to calculate a p-value or check whether the confidence intervals at a particular alpha-level include zero. (continued) – Ben Bolker Oct 3 '14 at 18:57
  • 4
    (2) overlapping CIs is at best a conservative, and approximate, criterion for differences between parameters (there are several published papers on this topic). (3) There's a separate/orthogonal issue with pairwise comparisons, which is that one has to control appropriately for multiplicity and non-independence of the comparisons (this can be done, e.g. by the methods in the multcomp package, but it requires at least a little bit of care) – Ben Bolker Oct 3 '14 at 18:58
  • 1
    For what? You might want to ask a new question. – Ben Bolker May 9 '17 at 5:31

It looks like what you have done in the second method is to have computed confidence intervals for the regression coefficients, then transformed those to obtain CIs for the predictions. This ignores the covariances between the regression coefficients.

Try fitting the model without an intercept, so that the batch effects will actually be the predictions, and confint will return the intervals you need.

Addendum 1

I did exactly what I suggested above:

> fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
> confint(fm2)
Computing profile confidence intervals ...
           2.5 %    97.5 %
.sig01  0.000000  1.637468
.sigma  2.086385  3.007380
batchA 60.234772 64.298581
batchB 57.268105 61.331915
batchC 60.018105 64.081915
batchD 57.668105 61.731915
batchE 53.868105 57.931915
batchF 59.001439 63.065248
batchG 57.868105 61.931915
batchH 61.084772 65.148581
batchI 56.651439 60.715248
batchJ 56.551439 60.615248

These intervals seem to jibe with the results from effects.

Addendum 2

Another alternative is the lsmeans package. It obtains degrees of freedom and an adjusted covariance matrix from the pbkrtest package.

> library("lsmeans")
> lsmeans(fm1, "batch")
Loading required namespace: pbkrtest
 batch   lsmean       SE    df lower.CL upper.CL
 A     62.26667 1.125709 40.45 59.99232 64.54101
 B     59.30000 1.125709 40.45 57.02565 61.57435
 C     62.05000 1.125709 40.45 59.77565 64.32435
 D     59.70000 1.125709 40.45 57.42565 61.97435
 E     55.90000 1.125709 40.45 53.62565 58.17435
 F     61.03333 1.125709 40.45 58.75899 63.30768
 G     59.90000 1.125709 40.45 57.62565 62.17435
 H     63.11667 1.125709 40.45 60.84232 65.39101
 I     58.68333 1.125709 40.45 56.40899 60.95768
 J     58.58333 1.125709 40.45 56.30899 60.85768

Confidence level used: 0.95 

These are even more in line with the effect results: the standard errors are identical, but effect uses different d.f. The confint results in Addendum 1 are even narrower than the asymptotic ones based on using $\pm1.96\times\mbox{se}$. So now I think those are not very trustworthy.

Results from effect and lsmeans are similar, but with an unbalanced multi-factor situation, lsmeans by default averages over unused factors with equal weights, whereas effect weights by the observed frequencies (available as an option in lsmeans).

  • Thank you for this solution. Intervals are now more similar, although not exactly the same. Your answer still does not answer the question whether CIs from effects package can be trusted for lmer objects. I am considering using the results in a publication and want to be sure that CIs are calculated using an approved method for LMMs. – Mikko Oct 3 '14 at 7:58
  • Would you please tell : in your Addendum 1 the first two parameters .sig01 and .sigma produces by confint , are those confidence interval for variance ? or confidence interval of standard deviation ? – ABC Jul 30 '15 at 8:38
  • They are CIs for whatever parameters are labeled that way in the model. You should look at the documentation for lmer for a definitive answer. However, people usually use notations like sigma to refer to standard deviations, and sigma.square or sigma^2 to refer to variances. – rvl Jul 30 '15 at 13:31
  • Is it better to use lmertest, lsmeans or mertools? – skan May 9 '17 at 0:55

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.