Solved – Specify contrasts for lme with interactions

contrastslme4-nlme

I have used lme to generate a mixed effect model of the response of cells to a certain stimulus. The stimulus is applied 3 times in a row (coded by the Exposure factor that can be 1, 2, or 3) and the Genotype factor has 3 levels (wt, A, and B)

My model looks like:

model <- lme(AUC ~ Exposure + Sex + Genotype + 
            Exposure:Sex + Exposure:Genotype + Basal, 
         random = ~ 1 | Date / Experiment / Cell, 
         data = mydata)

anova returns significant interactions between the exposure number and both sex and genotype

                  numDF denDF   F-value p-value
(Intercept)           1   168   55.1581  <.0001
Exposure              2   168   10.1711  0.0001
Sex                   1    22    4.1581  0.0536
Genotype              2    22    3.7246  0.0404
Basal                 1    41 1142.2234  <.0001
Exposure:Sex          2   168    5.5801  0.0045
Exposure:Genotype     4   168    2.0960  0.0835

For the fixed effects, summary returns:

                             Value Std.Error  DF  t-value p-value
(Intercept)               204417.8 109088.33 168  1.87387  0.0627
Exposure2                -192542.9  58653.05 168 -3.28274  0.0013
Exposure3                -332725.9  58653.05 168 -5.67278  0.0000
SexM                     -232599.7 104210.57  22 -2.23202  0.0361
GenotypeA                -184772.3 125723.74  22 -1.46967  0.1558
GenotypeB                 -40073.2 128715.16  22 -0.31133  0.7585
Basal                       1650.4     48.83  41 33.79680  0.0000
Exposure2:SexM            135000.3  58133.66 168  2.32224  0.0214
Exposure3:SexM            203637.4  58133.66 168  3.50292  0.0006
Exposure2:GenotypeA       106377.5  71472.21 168  1.48838  0.1385
Exposure3:GenotypeA       159548.1  71472.21 168  2.23231  0.0269 
Exposure2:GenotypeB       101246.7  70397.58 168  1.43821  0.1522 
Exposure3:GenotypeB       191111.2  70397.58 168  2.71474  0.0073

Am I correct in interpreting the interaction terms of this output as:

There is a statistically significant difference in response of M vs F both during exposure 2 and 3.

There is a statistically significant difference between genotype A and the reference genotype (wt) for exposure 3 (p=0.02) but not for exposure 2 (p=0.13)?
And similarly for genotype B?

What if I wanted to compare exposure 2 to exposure 3?

I read about using estimable in the gmodels package, but I struggle to understand its syntax.

Best Answer

If you look at the summary of your fixed effects portion of the model, you can label each row as follows:

                                    Value Std.Error  DF  t-value p-value
beta0  (Intercept)               204417.8 109088.33 168  1.87387  0.0627
beta1  Exposure2                -192542.9  58653.05 168 -3.28274  0.0013
beta2  Exposure3                -332725.9  58653.05 168 -5.67278  0.0000
beta3  SexM                     -232599.7 104210.57  22 -2.23202  0.0361
beta4  GenotypeA                -184772.3 125723.74  22 -1.46967  0.1558
beta5  GenotypeB                 -40073.2 128715.16  22 -0.31133  0.7585
beta6  Basal                       1650.4     48.83  41 33.79680  0.0000
beta7  Exposure2:SexM            135000.3  58133.66 168  2.32224  0.0214
beta8  Exposure3:SexM            203637.4  58133.66 168  3.50292  0.0006
beta9  Exposure2:GenotypeA       106377.5  71472.21 168  1.48838  0.1385
beta10 Exposure3:GenotypeA       159548.1  71472.21 168  2.23231  0.0269 
beta11 Exposure2:GenotypeB       101246.7  70397.58 168  1.43821  0.1522 
beta12 Exposure3:GenotypeB       191111.2  70397.58 168  2.71474  0.0073

So the fixed effects portion of your model looks like this:

beta0 + beta1*Exposure2 + beta2*Exposure3 + beta3*SexM +     
beta4*GenotypeA + beta5*GenotypeB + beta6*Basal + 
beta7*Exposure2*SexM + beta8*Exposure3*SexM + 
beta9*Exposure2*GenotypeA + beta10*Exposure3*GenotypeA + 
beta11*Exposure2*GenotypeB + beta12*Exposure3*GenotypeB

where the symbol * denotes multiplication and all the beta's are true fixed effects (i.e., unknown but estimable from the data via the values listed in the Value column of your fixed effects summary).

If you are interested in describing the effect of Sex, for instance, all you have to do is to find all the terms in the model which include the dummy variable SexM (which is equal to 1 for Males and 0 for Females) and group them together. The coefficient of SexM obtained after this grouping denotes the effect of Sex:

(beta3 + beta7*Exp2 + beta8*Exp3)*SexM        (1)

From the above, we can see that the effect of Sex depends on the value of Exp. Recall that Exp2 and Exp3 are dummy variables defined as follows: Exp2 = 1 if Exp = 2 and 0 otherwise; Exp3 = 1 if Exp = 3 and 0 otherwise. We can exploit this to spell out the effect of SexM for:

Exp = 1 (that is, for Exp2 = 0 and Exp3 = 0);

Exp = 2 (that is, for Exp2 = 1 and Exp3 = 0);

Exp = 3 (that is, for Exp2 = 0 and Exp3 = 1).

When Exp = 1, substituting Exp2 = 0 and Exp3 = 0 in expression (1) above yields that the coefficient of SexM is beta3. This is the effect of Sex on your outcome variable when Exp = 1, all else being equal.

When Exp = 2, substituting Exp2 = 1 and Exp3 = 0 in expression (1) above yields that the coefficient of SexM is beta3 + beta7. This is the effect of Sex on your outcome variable when Exp = 2, all else being equal.

When Exp = 3, substituting Exp2 = 0 and Exp3 = 1 in expression (1) above yields that the coefficient of SexM is beta3 + beta8. This is the effect of Sex on your outcome variable when Exp = 3, all else being equal.

So if you wanted to set up linear combinations of parameters which encapsulate the effect of SexM for each value of Exp, all you need to do is to specify these combinations so they reflect the parameters you are interested in: beta3, beta3 + beta7 and beta3 + beta8. (Again, the betas are true fixed effects, not estimated fixed effects.)

Since the fixed effects portion of your model includes the parameters beta0 through beta12, you are going to set up each combination via a row vector which includes 13 components. The first component of this vector corresponds to beta0, the second component corresponds to beta1, ..., the last component corresponds to beta12.

In R, beta3 is a combination of all the fixed effects model parameters with weights given by the components of the row vector c1:

c1 <- rep(0, 13)
names(c1) <- paste0("beta",0:12)
c1[names(c1)=="beta3"] <- 1
c1

In other words, beta3 = 0*beta0 + 0*beta1 + 0*beta2 + 1*beta3 + 0*beta4 + 0*beta5 + 0*beta6 + 0*beta7 + 0*beta8 + 0*beta9 + 0*beta10 + 0*beta11 + 0*beta12.

The linear combination of model parameters beta0 through beta12 that will encapsulate the parameter beta3 + beta7 is given by:

c2 <- rep(0, 13)
names(c2) <- paste0("beta",0:12)
c2[names(c2)=="beta3"] <- 1
c2[names(c2)=="beta7"] <- 1
c2

In other words, beta3 + beta7 = 0*beta0 + 0*beta1 + 0*beta2 + 1*beta3 + 0*beta4 + 0*beta5 + 0*beta6 + 1*beta7 + 0*beta8 + 0*beta9 + 0*beta10 + 0*beta11 + 0*beta12.

The linear combination of model parameters beta0 through beta12 that will encapsulate the parameter beta3 + beta8 is given by:

c3 <- rep(0, 13)
names(c3) <- paste0("beta",0:12)
c3[names(c3)=="beta3"] <- 1
c3[names(c3)=="beta8"] <- 1
c3

such that beta3 + beta8 = 0*beta0 + 0*beta1 + 0*beta2 + 1*beta3 + 0*beta4 + 0*beta5 + 0*beta6 + 0*beta7 + 1*beta8 + 0*beta9 + 0*beta10 + 0*beta11 + 0*beta12.

If you now want to simultaneously test hypotheses such as:

H0: beta3 = 0 versus Ha: beta3 != 0 

H0: beta3 + beta7 = 0 versus Ha: beta3 + beta7 != 0 

H0: beta3 + beta8 = 0 versus Ha: beta3 + beta8 != 0

you can achieve this by using the multcomp package in R. (Here, != 0 stands for "not equal to zero".) These hypotheses will enable you to determine whether Sex has an effect for any of the specific levels of Exp (i.e., if males differ from females with respect to the average response, all else being the same). This package can also be used to compute simultaneous confidence intervals for beta3, beta3 + beta7 and beta3 + beta8.

All you need to do to use multcomp in conjunction with your lme model is something like this:

library(multcomp)

c <- rbind(c1, c2, c3)

g <- glht(model, linfct = c)

s <- summary(g, test=adjusted("holm"))

s

ci <- confint(summary(g, test=adjusted("holm")))

ci

If you don't want to adjust p-values or confidence levels for multiplicity, just use adjusted("none") in the above.

Of course, there are R packages which will do all of the above for you with a minimal number of commands. But it helps to know how to test your own simple effects (e.g., beta3, beta3 + beta7, beta3 + beta8) when interested in probing a significant interaction.

In the approach I presented here, contrasts are specified via linear combinations of all model parameters. There are other ways to specify contrasts, which I will leave to others on this forum to address.

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.

Best Answer

Related Solutions

Solved – Specifying contrasts for lme in 3×3 repeated-measures design

Treatment contrasts:

Repeated contrasts (sliding differences):

R – How to Set Custom Contrasts with lmer in R

Custom contrasts

Helmert contrast with contr.helmert

Related Question

Helmert contrast with `contr.helmert`