Solved – Confidence-intervals for conditions tested with a mixed-effects model

confidence intervalmixed modelr

Sorry for a possibly ignorant question.

I have fit a mixed-effects model using the lmer function from the lme4 package, and the main fixed effect (a factor with three levels) in the model was significant according to a run with pvals.fnc (from the languageR package).

To illustrate the effect in an appealing way, I would like to plot three bars with the "baseline condition" (the intercept including the effect level 1), the "experimental condition 1" (the intercept + the effect level 2), the "experimental condition 2" (intercept + effect level 3), and the corresponding confidence intervals for these conditions.

However, how do I do that? The function plotLMER.fnc used to do that, but has stopped working for me (lme4_0.999375-32; languageR_1.0).
Model: mylmer <- lmer(outcome ~ (1|subject) + (1|item) + Factor, data)

So I have MCMC output in the form of:

$fixed
              Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)
(Intercept)     0.4728   0.4718     0.2250     0.7368 0.0010   0.0008
Factor2        -0.0420  -0.0420    -0.1732     0.0931 0.5414   0.5710
Factor3        -0.1643  -0.1631    -0.3153    -0.0112 0.0328   0.0508

My idea was to construct the conditions from this data by simply using the intercept as the baseline condition, and construct the experimental conditions by adding the effect of each factor level to the intercept. The new (intercept + effect combined) CI would be constructed by turning the HPD interval into a standard deviation and then use the square root of the sum of squared standard deviations (like here)

Is this approach appropriate? If I do it, the CIs become become fairly large and the effect no longer looks significant (conflicting with the pvals.fnc output).

Is there a better way (code example would be great)?

If not, is there a "least worst" solution that could satisfy readers who strongly request the standard bar plot with CIs?

Problems in the back of my head, which I am too ignorant to properly formulate:
CI not equal to HPD.
Transforming CI to SD is inappropriate, as CI may be asymmetric.
Problems with determining CI at all in mixed models (like here).

Best Answer

The DRAFT r-sig-mixed-models FAQ details (in the "Predictions and/or confidence (or prediction) intervals on predictions" section) how to obtain predictions and confidence intervals for cells in the design of a mixed effects model. The ezPredict() function in the ez package wraps the code for the lme4 case (well, obtaining predictions and variances, leaving the user to decide their own CI).

Related Solutions

Solved – Relationship between paired t test and simple mixed model

I think the problem is the way the paired t test is computed. Try this:

t.test(all_by_sub$var[all_by_sub$condition==1], all_by_sub$var[all_by_sub$condition==0], paired=TRUE)

This gives:

data:  all_by_sub$var[all_by_sub$condition == 1] and all_by_sub$var[all_by_sub$condition == 0]
t = 2.0529, df = 37, p-value = 0.0472
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.0009400428 0.1435297005
sample estimates:
mean of the differences 
             0.07223487

Obviously, t.test computes x minus y in the paired t test. That's why the sign of the estimate and the t value was reversed. Beyond that, both the estimate (0.072 vs. 0.082) and the t value (2.05 vs. 2.19) of the mixed model are very close to the results of the t test:

            Estimate Std. Error      df t value Pr(>|t|)
(Intercept)    0.118      0.028 106.772   4.159    0.000
condition      0.082      0.037 462.992   2.192    0.029

R – How to Set Custom Contrasts with lmer in R

For the following steps, we need the data frame in the long format. The data frame dat contains the dependent variable result, the categorical predictor cond (levels: a, b, and c), and the random factor s.

library(tidyr)
dat <- gather(temp, cond, result, a, b, c)

In the following, I will illustrate two approaches to create a contrast matrix corresponding to the conditions you want to compare:

$a - \frac{b+c}{2}$
$b - c$

Custom contrasts

The matrix mat corresponds to the level differences.

mat <- rbind(c(1, -0.5, -0.5),     # a vs. (b + c) / 2
             c(0, 1, -1))          # b vs. c

To create the actual contrast matrix, we compute the generalized inverse with ginv (from MASS).

library(MASS)
cMat <- ginv(mat)
#            [,1]          [,2]
# [1,]  0.6666667 -7.130169e-17
# [2,] -0.3333333  5.000000e-01
# [3,] -0.3333333 -5.000000e-01

This contrast matrix cMat can be used in lmer.

library(lme4)
res <- lmer(result ~ cond + (1|s), data = dat, 
            contrasts = list(cond = cMat))
coef(summary(res))    
#              Estimate Std. Error    t value
# (Intercept) -2.948115  0.0946025 -31.163182
# cond1        1.351517  0.2006822   6.734612
# cond2        1.153918  0.2317279   4.979625

As you can see, the fixed-effect estimates correspond to the differences specified above. Furthermore, the intercept represents the overall mean.

Helmert contrast with `contr.helmert`

You can also use the built-in contr.helmert function to create the contrast matrix.

cHelmert <- contr.helmert(3)
#   [,1] [,2]
# 1   -1   -1
# 2    1   -1
# 3    0    2

However, the order does not correspond to the one you specified in the question. Hence, we have to reverse the order of columns and rows. The first column corresponds to b vs. a and the second one corresponds to c vs. the mean of b and a.

cHelmert2 <- cHelmert[c(3:1), 2:1]
#   [,1] [,2]
# 3    2    0
# 2   -1    1
# 1   -1   -1

Compare the contrast matrix cHelmert2 to cMat. You will notice that the columns are scaled versions of the other matrix.

The result of lmer is:

library(lme4)
res2 <- lmer(result ~ cond + (1|s), data = dat, 
             contrasts = list(cond = cHelmert2))
coef(summary(res2))    
#               Estimate Std. Error    t value
# (Intercept) -2.9481150 0.09460250 -31.163182
# cond1        0.4505056 0.06689407   6.734612
# cond2        0.5769590 0.11586393   4.979625

This contrast matrix allows fo the same comparisons as the custom contrast matrix. However, since the values in the matrix are different, the fixed-effects coefficients are different too. Not suprisingly, the $t$-values are the same.