# Large differences between raw (plotted) data and least-squares means from mixed model

I’m analysing data with mixed-models (using the afex package which I believe is based on lme4) from an experiment that had a (continuous-type) score as the dependent variable.

I used a factorial design; the fixed-effects were 3 subject variables (z1, z2, and z3; all continuous and centred), and 2 within-subjects experimental manipulations (one dichotomous x1, the other with three levels x2). My uncorrelated random effects were the participants [ID].

Upon running my mixed model (y ~ x1 + x2 + z1 + z2 + z3 + x1:x2 + x2:z1 + x2:z2 + (1|ID)) I found an interaction between one of the subject variables [z1] and one of the within-subject manipulations [x2]. Fixed effects can be seen below:

                                                       Estimate Std. Error t value
(Intercept)                                           1.718e-02  6.356e-02   0.270
x1b                                                   2.577e-01  3.423e-02   7.531
x2b                                                  -5.242e-02  4.999e-02  -1.048
x2c                                                   9.705e-03  4.999e-02   0.194
z1                                                   -2.254e-01  2.335e-01  -0.965
z2                                                    1.235e-01  2.282e-01   0.541
z3                                                    2.242e-01  3.911e-02   5.734
x1b:x2b                                              -1.253e-01  4.840e-02  -2.588
x1b:x2c                                               2.187e-01  4.840e-02   4.518
x2b:z1                                               -8.068e-01  1.861e-01  -4.335
x2c:z1                                                2.709e-01  1.861e-01   1.456
x2b:z2                                                5.806e-01  1.668e-01   3.480
x2c:z2                                               -1.065e-01  1.668e-01  -0.639


A quick scatterplot for this interaction appears to show that x2b has the flattest slope and x2c has the steepest:

When I use lsmeans and lstrend to try and interpret the interaction, it tells me that level a has a different slope from levels b and c, and that only level a had a significant slope. The lsmeans call:

lsmeans(model.mixed.optim, "x2", by = "z1", at = list("z1" = c(summary(data$z1)[[2]],summary(data$z1)[[5]]))) # compare 1st and 3rd quartiles

I was perplexed as to how this could be the case, until I looked at the least-squares means using lsmip:

I have a few questions:

1. Is lsmeans the correct way to analyse a mixed-model continuous-categorical interaction?

2. Can anyone tell me why the second plot is so different to the first (I know I'm plotting lsmeans vs scatter, but I would have thought they would match up at least a little) -- is this a consequence of specifying participant ID as a random effect?

3. Does anyone have recommended descriptives when reporting the outcome of a mixed model (as I would guess means are not appropriate)?

• What exactly was the model that you fitted that is the basis for the lsmeans results? Second, why are the two values of z1 in the plot -.1854 and .1156? Commented Mar 24, 2015 at 16:56
• I've added the model into the main question. The values of z1 represent the first and third quartiles of the continuous data. Commented Mar 25, 2015 at 14:40

The plot you show cannot have resulted from the model you show in your question. In that model, the slopes for z1 are:

x2 = a :                  -0.2254    (coef of z1)
x2 = b : -.2254 - .8068 = -1.0322    (sum of coefs for z1 and x2b:z1)
x2 = c : -.2254 + .2709 =  0.0455    (sum of coefs for z1 and x2c:z1)


So my answer is that there are two reasons the results are so discrepant

1. Because they are based on a different fitted model. I suggest you go back and re-fit, and find out which results are correct
2. Collinearity among the predictors (that includes both factors with possibly unbalanced levels, and covariates). It is not at all unusual to see the signs of regression coefficients change radically in the presence of other predictors, because the coefficients represent partial effects after taking into account simultaneous linear changes in other predictors.

Firstly you are not using a mixed effect model in your using mixed effect analysis. Secondly you need to decide what analysis you plan to use before caring out the experiment if not your analysis is more akin to data mining. Any advice about the analysis would be wrong. I hope you are a student in which case ask your tutor or supervisor what was planned when you designed the experiment they will be best equipped to give you a good answer. Hope this has been useful.

• Thanks for the clarification re: mixed effect models vs analysis. I had initially conceptualised this analysis as an ANCOVA, but many on this website and elsewhere (and I agree, after having done my own reading) advise that mixed effects analysis would be better-suited. I am specifically looking for advice on how to interpret my interaction and/or the most appropriate way to conduct post-hoc tests based on the output of the mixed effects analysis. Commented Mar 23, 2015 at 18:46
• If this is the case then yes it may be better to use mixed effects has it has its advantages. The problem you have is what to do next there is an obvious and clear interaction - "a" is clearly different to "b" and "c". As mixed effects relatively new in its use there really isn't a "gold standard" for what to do.You could for example create a regression model for "a" and then see if it fits (i'm pretty sure it want) "b", and "c" or you could do something as simple (and relatively easy) and Run a Tukey test (assuming that none where a control) - bit old fashioned I know. Commented Mar 23, 2015 at 20:44