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.
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.