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

contrastslme4-nlmemultiple-comparisonsrrepeated measures

I would like to test the effects of two categorical variables, each with 3 levels, on some continuous data. Three groups of participants completed the same task for 3 types of stimuli, so it's a 3×3 repeated measures design. I'm using lme. The fixed-effects part of the model is the interaction between participant group ("play") and stimuli ("ins"), and the random-effects part is subject grouped within stimulus category.

contrasts(iqr$ins) <- contr.sum(3)
contrasts(iqr$play) <- contr.sum(3)
iqr.lme <- lme(iqr ~ play*ins, random = ~1|subject/ins, data = iqr)

I'm using sum contrasts and would like to see comparisons between each level of the "play" and "ins" variables. In other words, I would like to use a contrast matrix for each variable that looks like the one below. If I understand correctly, -1 is the reference level and 1 is the level compared to it. Specifying contr.sum only results in contrasts for the first two columns being given, however.

          [,1]  [,2]  [,3]
clarinet    1     0     1
piano       0     1    -1
violin     -1    -1     0

Here is the fixed effects part of the summary output:

Fixed effects: iqr ~ play + ins + play * ins 
                 Value  Std.Error DF   t-value p-value
(Intercept)  2.9768519 0.11599299 66 25.664066  0.0000
play1       -0.0185185 0.16403886 33 -0.112891  0.9108
play2        0.2037037 0.16403886 33  1.241801  0.2231
ins1         0.1342593 0.09005357 66  1.490882  0.1408
ins2        -0.3240741 0.09005357 66 -3.598681  0.0006
play1:ins1  -0.3009259 0.12735498 66 -2.362891  0.0211
play2:ins1   0.0601852 0.12735498 66  0.472578  0.6381
play1:ins2   0.1574074 0.12735498 66  1.235974  0.2208
play2:ins2  -0.1481481 0.12735498 66 -1.163269  0.2489

If I understand correctly, ins1 is the contrast for clarinet vs. violin (using the contrast matrix above), and ins2 is the contrast for piano vs. violin. The problem is that I would also like to see a contrast for clarinet vs. piano.

My solution was to run a second model using a contrast matrix that I specified myself. I thought this would give me the missing contrast and replicate the clarinet vs. violin values from the original output.

cst1 <- cbind(c(-1,1,0), c(-1,0,1))
contrasts(iqr$ins) <- cst1
contrasts(iqr$play) <- cst1
iqr.lme <- lme(iqr ~ play*ins, random = ~1|subject/ins, data = iqr)
contrasts(iqr$ins)

          [,1]  [,2]
clarinet   -1    -1
piano       1     0
violin      0     1

Here, ins1 should correspond to clarinet vs. piano, right? But the values are the same as the ins2 contrast in my original output, indicating that the contrast matrix I specified was blatantly ignored and ins1 instead corresponds to piano vs. violin. The values for play2 and ins2 are new, and I could just assume that they indicate the contrast that was missing from the original output, but I would really rather understand what is going on than base my results on dubious assumptions!

Fixed effects: iqr ~ play * ins 
                 Value  Std.Error DF   t-value p-value
(Intercept)  2.9768519 0.11599299 66 25.664066  0.0000
play1        0.2037037 0.16403886 33  1.241801  0.2231
play2       -0.1851852 0.16403886 33 -1.128910  0.2671
ins1        -0.3240741 0.09005357 66 -3.598681  0.0006
ins2         0.1898148 0.09005357 66  2.107799  0.0389
play1:ins1  -0.1481481 0.12735498 66 -1.163269  0.2489
play2:ins1  -0.0092593 0.12735498 66 -0.072704  0.9423
play1:ins2   0.0879630 0.12735498 66  0.690691  0.4922
play2:ins2  -0.2314815 0.12735498 66 -1.817608  0.0737

Am I misunderstanding how contrast matrices are read or specified? Or is there actually a way to get all the contrasts I want listed at once? I am fairly new to R and very new to nlme, so I could very well be missing something obvious. I hope this is enough information for someone to be able to help. Thanks a lot!

Best Answer

If you have a factor with $n$ levels, you can test $n - 1$ effects (at maximum) associated with this factor in one model.

Have a look at the sum contrasts. The function contr.sum(3) returns the following matrix:

  [,1] [,2]
1    1    0
2    0    1
3   -1   -1

The columns denote the contrasts. For a factor with $3$ levels, you obtain $2$ contrasts. Your interpretation of the contrasts is not correct. The first one tests the first level of the factor against levels 2 and 3. The second contrasts tests the second levels against levels 1 and 3.

You manually specified a contrast matrix with cbind(c(-1,1,0), c(-1,0,1)):

     [,1] [,2]
[1,]   -1   -1
[2,]    1    0
[3,]    0    1

This matrix also allows two tests. You test level 2 against levels 1 and 3. This is identical to the second test in the first matrix. Hence, you obtain the identical coefficient in both models. The second test of this contrast matrix is level 3 against levels 1 and 2.

If you want to compare levels against each other, there are two possibilities:

Treatment contrasts:
```
contr.treatment(3)

  2 3
1 0 0
2 1 0
3 0 1
```
Here, you test the first level against (a) the second and (b) the third one.
Repeated contrasts (sliding differences):
```
library(MASS)
contr.sdif(3)

         2-1        3-2
1 -0.6666667 -0.3333333
2  0.3333333 -0.3333333
3  0.3333333  0.6666667
```
Here, you test (a) the first level against the second one and (b) the second level against the third one. This matches the column names.

Please note that no contrast matrix will allow you to have more than $2$ contrasts for a factor with $3$ levels.

If you want to manually create an appropriate contrast matrix based on specific tests, you can use the ginv function from the MASS package:

cst1 <- cbind(c(-1,1,0), c(-1,0,1)) # your contrast matrix

library(MASS)
t(ginv(cst1))

       [,1]       [,2]
[1,] -0.3333333 -0.3333333
[2,]  0.6666667 -0.3333333
[3,] -0.3333333  0.6666667

The resulting matrix allows the tests you are looking for, i.e., level 2 against level 1 and level 3 against level 1.

For more details on how you can create a contrast matrix based on manually defined comparisons, have a look at this answer.

Related Solutions

Paired t-test – Special Case of Linear Mixed-Effect Modeling

The equivalence of the models can be observed by calculating the correlation between two observations from the same individual, as follows:

As in your notation, let $Y_{ij} = \mu + \alpha_i + \beta_j + \epsilon_{ij}$, where $\beta_j \sim N(0, \sigma_p^2)$ and $\epsilon_{ij} \sim N(0, \sigma^2)$. Then $Cov(y_{ik}, y_{jk}) = Cov(\mu + \alpha_i + \beta_k + \epsilon_{ik}, \mu + \alpha_j + \beta_k + \epsilon_{jk}) = Cov(\beta_k, \beta_k) = \sigma_p^2$, because all other terms are independent or fixed, and $Var(y_{ik}) = Var(y_{jk}) = \sigma_p^2 + \sigma^2$, so the correlation is $\sigma_p^2/(\sigma_p^2 + \sigma^2)$.

Note that the models however are not quite equivalent as the random effect model forces the correlation to be positive. The CS model and the t-test/anova model do not.

EDIT: There are two other differences as well. First, the CS and random effect models assume normality for the random effect, but the t-test/anova model does not. Secondly, the CS and random effect models are fit using maximum likelihood, while the anova is fit using mean squares; when everything is balanced they will agree, but not necessarily in more complex situations. Finally, I'd be wary of using F/df/p values from the various fits as measures of how much the models agree; see Doug Bates's famous screed on df's for more details. (END EDIT)

The problem with your R code is that you're not specifying the correlation structure properly. You need to use gls with the corCompSymm correlation structure.

Generate data so that there is a subject effect:

set.seed(5)
x <- rnorm(10)
x1<-x+rnorm(10)
x2<-x+1 + rnorm(10)
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), 
                    rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Then here's how you'd fit the random effects and the compound symmetry models.

library(nlme)
fm1 <- lme(y ~ x, random=~1|subj, data=myDat)
fm2 <- gls(y ~ x, correlation=corCompSymm(form=~1|subj), data=myDat)

The standard errors from the random effects model are:

m1.varp <- 0.5453527^2
m1.vare <- 1.084408^2

And the correlation and residual variance from the CS model is:

m2.rho <- 0.2018595
m2.var <- 1.213816^2

And they're equal to what is expected:

> m1.varp/(m1.varp+m1.vare)
[1] 0.2018594
> sqrt(m1.varp + m1.vare)
[1] 1.213816

Other correlation structures are usually not fit with random effects but simply by specifying the desired structure; one common exception is the AR(1) + random effect model, which has a random effect and AR(1) correlation between observations on the same random effect.

EDIT2: When I fit the three options, I get exactly the same results except that gls doesn't try to guess the df for the term of interest.

> summary(fm1)
...
Fixed effects: y ~ x 
                 Value Std.Error DF   t-value p-value
(Intercept) -0.5611156 0.3838423  9 -1.461839  0.1778
xx2          2.0772757 0.4849618  9  4.283380  0.0020

> summary(fm2)
...
                 Value Std.Error   t-value p-value
(Intercept) -0.5611156 0.3838423 -1.461839  0.1610
xx2          2.0772757 0.4849618  4.283380  0.0004

> m1 <- lm(y~ x + subj, data=myDat)
> summary(m1)
...
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.3154     0.8042  -0.392  0.70403   
xx2           2.0773     0.4850   4.283  0.00204 **

(The intercept is different here because with the default coding, it's not the mean of all subjects but instead the mean of the first subject.)

It's also of interest to note that the newer lme4 package gives the same results but doesn't even try to compute a p-value.

> mm1 <- lmer(y ~ x + (1|subj), data=myDat)
> summary(mm1)
...
            Estimate Std. Error t value
(Intercept)  -0.5611     0.3838  -1.462
xx2           2.0773     0.4850   4.283

Model Selection in Longitudinal Data – Testing the Need for Random-Effects Terms in Longitudinal Data Analysis

The likelihood ratio test is slightly incorrect (in general, conservative) for testing the significance of a random effect, because the null value ($\sigma^2=0$) is at the boundary of the feasible space, but in this case there is overwhelmingly strong evidence against the null hypothesis. The model with random effects of individual is 15713-6772=8941 log-likelihood units better; twice the log-likelihood value is $\chi^2$ distributed, so the direct p-value calculation would give you ...

pchisq(2*8941,df=1,lower.tail=FALSE,log.p=TRUE)/log(10)
## -3885.251

... a p-value of approximately $10^{-3885}$.

You should really consider a random-slope model (random = ~time|id) as well.

Update: relative to the random-intercept model, the random-slopes model is again much better. The improvement is now 935 log-likelihood units, which doing the equivalent calculation as above corresponds to a rejection of the null hypothesis (among-individual variation in slope is equal to zero) with a p-value of "only" $10^{-408}$.

Best Answer

Treatment contrasts:

Repeated contrasts (sliding differences):

Related Solutions

Paired t-test – Special Case of Linear Mixed-Effect Modeling

Model Selection in Longitudinal Data – Testing the Need for Random-Effects Terms in Longitudinal Data Analysis

Related Question