Generating similar estimates of interactions in afex, lsmeans, and lme4 packages

Question

I would like to know if there is a way get the same estimates of an interaction effect in afex & lsmeans packages as in lmer. The toy data below is for two groups with different intercepts and slopes.

set.seed(1234)
A0 <- rnorm(4,2,1)
B0 <- rnorm(4,2+3,1)
A1 <- rnorm(4,6,1)
B1 <- rnorm(4,6+2,1)
A2 <- rnorm(4,10,1)
B2 <- rnorm(4,10+1,1)
A3 <- rnorm(4,14,1)
B3 <- rnorm(4,14+0,1)
score <- c(A0,B0,A1,B1,A2,B2,A3,B3)
id <- factor(rep(1:8,times = 4, length = 32))
time <- factor(rep(0:3, each = 8, length = 32))
timeNum <- as.numeric(rep(0:3, each = 8, length = 32))
group <- factor(rep(c("A","B"), times =2, each = 4, length = 32))
df <- data.frame(id, group, time, timeNum, score)
df  

And here is the plot

(ggplot(df, aes(x = time, y = score, group = group)) + 
    stat_summary(fun.y = "mean", geom = "line", aes(linetype = group)) +
    stat_summary(fun.y = "mean", geom = "point", aes(shape = group), size = 3) +
    coord_cartesian(ylim = c(0,18)))

When I run a standard lmer on the data looking for an estimate of the difference in change in score over time between groups.

summary(modelLMER <- lmer(score ~ group * timeNum + (timeNum|id), df))

I get an estimate for the group*time interaction of -1.07, which means that the increase in score for a one-unit increase in time is ~1 point less in group B than group A. This estimate matches the preset differences I built into the dataset.

What I would like to know is how to do a similar thing in the afex and lsmeans packages.

library(afex)
library(lsmeans)

First I generated the afex model object

modelLM <- aov_ez(id="id", dv="score", data=df, between="group", within="time", 
                 type=3, return="lm")

Then passed that into the lsmeans function

lsMeansLM <- lsmeans(modelLM, ~rep.meas:group) 

My goal is to generate an accurate estimate of the group*time interaction in afex and lsmeans. To do so requires specifying custom contrast matrices based on the split specified in the lsmeans function above.

groupMain = list(c(-1,-1,-1,-1,1,1,1,1)) # group main effect
linTrend = list(c(-3,-1,1,3,-3,-1,1,3)) # linear trend
linXGroup = mapply("*", groupMain, linTrend) # group x linear trend interaction

Then I made a master list

contrasts <- list(groupMain=groupMain, linTrend=linTrend, linXGroup=linXGroup)

Which I passed into the contrast function in lsmeans.

contrast(lsMeansLM, contrasts)

The F and p values in the output match those for the automatic tests for linear trend and for the group difference in linear trend generated from a mixed ANCOVA in SPSS. However the mixed ANCOVA does not generate an estimate.

The estimate of the effect using the procedure above, instead of being approx. -1, like in the lmer (and matching the difference I built into the data) is approx. -10, which is wildly inaccurate.

I assume it has something to do with how I am coding the contrast coefficients. I know if I normalise the coefficients of the groupMain matrix by dividing all coefficients by four that yields an accurate estimate of the main effect of group averaged across all timepoints. But I have no idea how to get an accurate estimate either of linear trend averaged across groups (linTrend), or an accurate estimate of the difference in linear trend across groups (linXGroup).

I am not sure if this question is more suitable for here or Cross Validated. I figured here first because it seems to be software related, but I know there are probably deeper issues involved. Any help would be much appreciated.

Answer 1

The issue here is that timeNum is a numeric predictor. Therefore, the interaction is a comparison of slopes. Note this:

> lstrends(modelLMER, ~group, var = "timeNum")
 group timeNum.trend       SE  df lower.CL upper.CL
 A          4.047168 0.229166 6.2 3.490738 4.603598
 B          2.977761 0.229166 6.2 2.421331 3.534191

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
> pairs(.Last.value)
 contrast estimate        SE  df t.ratio p.value
 A - B    1.069407 0.3240897 6.2     3.3  0.0157

There's your 1.07 - the opposite sign because the comparison is in the other direction.

I will further explain that the lsmeans result you describe in the question is a comparison of the two group means, not an interaction contrast. lsmeans uses a reference grid:

> ref.grid(modelLMER)
'ref.grid' object with variables:
    group = A, B
    timeNum = 1.5

and as you can see, timeNum is being held fixed at its mean of 1.5. The LS means are predictions for each group at timeNum = 1.5 -- often called the adjusted means; and the difference is thus the difference between those two adjusted means.

Regarding the discrepancy claimed in obtaining your linear contrast of about 10.7: The linear contrast coefficients c(-3,-1,1,3) give you a multiple of the slope of the line. To get the slope, you need to divide by sum(c(-3,-1,1,3)^2) -- and also multiply by 2, because the contrast coefficients increment by 2.

Answer 2

Thanks to the invaluable help of @rvl I was able to solve this. Here is the code.

In order to generate the correct contrast matrices we first need to normalise them

(mainMat <- c(-1,-1,-1,-1,1,1,1,1)) # main effects matrix
(trendMat <- c(-3,-1,1,3,-3,-1,1,3) # linear trend contrast coefficients 
(nTimePoints <- 4) # number of timePoints
(mainNorm <- 1/nTimePoints)
(nGroups <- 2) # number of between-Ss groups
(trendIncrem <- 2) # the incremental increase of each new trend contrast coefficient    
(trendNorm <- trendIncrem/(sum(trendMat^2))) # normalising the trend coefficients

Now we create several contrast matrices in the form of lists. These are normalised using the objects we created above

(groupMain = list(mainMat*mainNorm)) # normalised group main effect
(linTrend = list(trendMat*trendNorm)) # normalised linear trend
(linXGroup = list((mainMat*trendMat)*(nGroups*trendNorm))) # group x linear trend interaction

Now pass those lists of matrices into a master list

contrasts <- list(groupMain=groupMain, linTrend=linTrend, linXGroup=linXGroup)

And pass that master list into the contrasts function in lsmeans

contrast(lsMeansLM, contrasts)

This is the output

 contrast                                               estimate        SE df t.ratio p.value
 c(-0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25)  1.927788 0.2230903  6   8.641  0.0001
 c(-0.15, -0.05, 0.05, 0.15, -0.15, -0.05, 0.05, 0.15)  3.512465 0.1609290  6  21.826  <.0001
 c(0.3, 0.1, -0.1, -0.3, -0.3, -0.1, 0.1, 0.3)         -1.069407 0.3218581  6  -3.323  0.0160

How do we check if these are accurate estimates?

Note first that the estimate of the group*time interaction is now approximately the same value as is returned by

summary(modelLMER)

The 'main effect' trend (for want of a better descriptor), which is the rate of change in score across the four time points averaged across both levels of group, is 3.51. If we change the coding of the group factor to simple coding via

contrasts(df$group) <- c(-.5,.5)

and run summary(modelLMER) again, the time estimate will now be 3.51.

Finally for the main effect of group, that is, the difference in score between groups averaged across all time points. We can run

pairs(lsmeans(modelLM,"group"))

And this will be -1.92. Thank you @rvl. A great answer. Using afex and lsmeans we have now forced a mixed ANCOVA that treats the repeated measures variable as categorical to give us estimates of group differences in trend and main effects that match those returned by a mixed-effects model where the repeated measures variable is continuous, and with p- and F-values that match those of SPSS.