Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I was doing ANOVA in SPSS and then in R and to my huge surprise, the results from the latest version of R were incorrect.

When I use the function model.tables(x,"means") to get descriptive statistics, the independent variable means by the second dependent are slightly incorrect (e.g. 129 instead of 130.27).

My question is what could cause the problem? I am a novice to R but using the same data, SPSS gets the result correctly, so something is obviously wrong.


  skupina pohlavie zodpovedny
1       1        1        152
2       1        1        118
3       2        2         88
4       2        1        140


x <- aov(zodpovedny ~ pohlavie*skupina,data=data)

Problem illustrated:

This is unfortunate.

share|improve this question
Can you post your complete data? – mnel Dec 13 '12 at 23:32
@mnel The data do not belong to me and I need to ask for a permission first, so maybe eventually. All zodpovedny values are indeed integers from an interval [0,200] and when I use the mean() function independently, R returns the correct value of 130.2692 instead of 129. Additionally, when I switch the formula from pohlavie*skupina to skupina*pohlavie, values for pohlavie are correct while values for skupina are not. Interestingly, as I have made many trials and restarted R a few times during the process, it even returned 141.88 instead of the former and correct 141.9. – Harold Cavendish Dec 14 '12 at 0:03
There isn't much difference between 141.88 and 141.9 apart from the number of decimal places you are reporting. – mnel Dec 14 '12 at 0:45
Why do you assume that R and SPSS are using "the same procedure"? If you can't post your original data, you could at least post a reproducible example with some made-up data. It's hard to help you if you don't. – mark999 Dec 14 '12 at 0:53
No. You can't. What you can do however is work out what R is doing by looking at the source code. – mnel Dec 14 '12 at 0:56
up vote 5 down vote accepted

As you point out, the individual cell means match, but where you see the problem is in the marginal means. There are multiple ways to calculate the marginal means. Suppose that the data has information on sex (male/female) and age (old/young) and we want to calculate the margin for sex. One approach is to ignore the age variable and just take the mean of all the males and the mean of all the females. Another approach is to find the mean of males by averaging the mean of old males and the mean of young males (add the 2 means and divide by 2). In a balanced design those 2 methods will give the same answer (can be shown with simple algebra), but in the unbalanced case they will usually give different answers because the weight that each data point contributes to the overall mean is different. With model based means you can get different weightings from the 2 I mentioned (I used them for examples as simple ways to understand). I expect in your case that R and SPSS are likely using different approaches.

share|improve this answer

@mnel is correct in that because of the unbalanced design, the order of the terms matter in the output of model.tables.

ADDED: In the help file for aov, we read that it "is designed for balanced designs, and the results can be hard to interpret without balance." So if you want simple descriptive statistics, better to ask for them directly.

Now, it would be better if you had posted a full data set yourself, even if you had to make one an alternate one that showed the same problem. But you got lucky and a curious reader wanted to know what was going on, so I did that for you. Here's a sample data set:

d <- expand.grid(a=factor(LETTERS[1:2]), b=factor(letters[1:2]))
d <- d[rep(1:4, c(15,9,11,10)),]
d$y <- round(rnorm(nrow(d), mean=10, sd=2),1)

And we see that the order of the terms in the model matters (output truncated):

> model.tables(aov(y ~ a*b, data=d), "means")
 a      A      B
    10.43  9.921
 b      a      b
    9.843  10.64

> model.tables(aov(y ~ b*a, data=d), "means")
 b       a      b
     9.867  10.61
 a       A      B
     10.46  9.877

The first term in the model agrees with the actual mean and the other is different.

> tapply(d$y, d$a, mean)
        A         B 
10.426923  9.921053 
> tapply(d$y, d$b, mean)
        a         b 
 9.866667 10.609524 

Note that I say different, not wrong. It's telling you something correct about the model. I'm not sure what, actually, but I'm curious enough that I may look into the code for model.tables to see what. (Or maybe not, it's getting late.)

share|improve this answer
Sorry about not providing the data and thank you for making it for me. It would be much appreciated if you found out why is it different and explained it, along with how it can not be wrong when SPSS returns correct means using the same model. – Harold Cavendish Dec 14 '12 at 9:30
What is the correct method to do ANOVA for unbalanced data in R, then? – Harold Cavendish Dec 14 '12 at 17:58
That comment isn't just about R, it's true about ANOVA in general. aov and lm are the right tools to use, you just have to know what the output from them means. But that's true in any software. The most common error is perhaps not knowing what kind of tests are being performed; to use SAS terms, Type I, II, or III. The summary of an aov fit will be Type I; to get other types, a common recommendation is Anova from the car package. – Aaron Dec 14 '12 at 18:05

Watch out: the model.tables() function only works with balanced designs. If you want to have the marginal means for unbalanced design you should use the popMeans() function. Imagine you have the following model:

Check.Model <- aov(dependent ~ factor1 + factor2,

If you would want the marginal means for the levels of factor1 (i.e., averaged over the levels of factor2) in an unbalanced design, you should use the popMeans() function from the doBy package:

popMeans(Check.Model, eff=c("factor1"))
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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