Let's say I have an experiment with three within-subject factors, A, B, & C. The data looks like this.

 s  a  b  c
 1  1  1  1
 1  1  1  2
 1  1  2  1
 1  1  2  2
 1  2  1  1
 1  2  1  2
 1  2  2  1
 1  2  2  2

Simple enough. I have 49 subjects. Now, to do this ANOVA in R, I use

m1 <- aov(score ~ a*b*c + Error(subject/(a*b*c)), data)
summary(m1, type=3)
... (clipped) ...
Error: s:b:c
                  Df Sum Sq Mean Sq F value  Pr(>F)   
b:c                1  4.608   4.608   8.121 0.00643 **
Residuals         48 27.236   0.567  

That looks fine, and matches SPSS's repeated measures GLM. All is well.

We can also do a mixed model in R using lme4 and get the exact same results as this, as well as a mixed model done in JMP.

m2 <- lmer(score~a*b*c + (a*b*c|s), data)
library(car); Anova(m2, type=3, test.statistic"F")
... (clipped) ...
                     F Df Df.res    Pr(>F) 
b:c             8.1206  1 48.000  0.006430 ** 

I can do the same thing in SAS and get similar results.

proc mixed data=mixedexample method=reml covtest;
    class a b c s;
    model score = a|b|c; 
    random intercept a|b|c/sub=s;

I can do the ANOVA, by hand in Stata:

anova score a / s|a ///
            b / s#b ///
            c / s#c ///
            a#b / s#a#b ///
            a#c / s#a#c /// 
            b#c / s#b#c /// 

I leave off the full interaction error term so it is the residual. Though, for some reason, the df of the main effect of a is twice the size of the others. But that's not the question I have.

My question is: how do I do the lmer and proc mixed version of the full LMM in Stata? The simple version is

xtmixed score a##b##c || s:, reml

But how do I add the fully crossed error terms to xtmixed in the same way I do in lmer by adding + (a*b*c|s)? The data is balanced with no missing values, so the LMM should be the same as the repeated-measures ANOVA, right? Why can't I do this in Stata?

I barely know the basics of LMMs, but this is a Stata question. I'm just trying to figure out all the different ways of performing these two models. Also, if anyone knows a simpler way of doing the univariate ANOVA in Stata without specifying every single error term by hand? This may not even be the "right" way of doing this procedure, but because I get the same output everywhere else, how do I get Stata to do the same thing as R's lmer, SAS's proc mixed, and JMP?

  • $\begingroup$ Does this do it? xtmixed score a##b##c || s: R.a || s: R.b || s: R.c (I'm really just guessing, based on a very quick skim of the xtmixed manual [Examples 10 and 11 seem relevant]) $\endgroup$ – timbp Mar 6 '13 at 22:33
  • $\begingroup$ It's better but the standard errors are still much larger than SPSS's GLM or R's lmer. :/ b and b*c are big effects in the ANOVA, but aren't coming out as such in any model I try in Stata. SAS, JMP, and lmer seem to have the ability to do this, and it seems that Stata should, but I have yet to find out how. $\endgroup$ – jtth Mar 7 '13 at 4:19

Ok, here's how it's done. Of course, the answer is on UCLA ATS.

This is the "brute force" method that doesn't take into account the structure of your random effects. It works for smaller datasets, but for mine it explodes.

Another way of doing this to get the results from lmer or whatever is to findit desmat and install that (it's like xi3 or other expansion things).


desmat a*b*c, defcon(hel) // or defoncon(var)
xtmixed score _x* || subject: R._x_1 || subject: R._x_2 || subject: R._x_3 || subject: R._x_4 || subject: R._x_5 || subject: R._x_6 || subject: R._x_7, var 

This matches lmer(score ~ a*b*c + (a*b*c|subject), REML=F) exactly.

But it ruins your ability to do contrasts, as the "factors" are no longer factors. The ATS "brute force" way should work for smaller stuff, but it crashes for me after about 15 minutes. You can get close to the lmer by doing something like

xtmixed score a##b##c || subject: R._x_1 || subject: R._x_2 || subject: R._x_3 || subject: R._x_4 || subject: R._x_5 || subject: R._x_6 || subject: R._x_7, var 

That way margins still works, but some stuff is off because the random factors are not dummy coded, but whatever you specify through desmat (or xi3). But you can at least examine interactions.

Now to figure out why standard error is so much larger in mixed models than in ANOVAs. Even in SPSS the so-called "Repeated Measures ANOVA" (which uses the GLM command) gives much larger standard error terms. But that's a separate question.

I hope this helps someone. For now, I'll just write my ANOVAs by hand... sigh. Not even SPSS gets the standard error right, compared to Stata's anova command. Makes a big difference for simple effects.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.