R model.tables() incorrect means – possible bug?

Question

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.

head(data):

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

Code:

x <- aov(zodpovedny ~ pohlavie*skupina,data=data)
model.tables(x,"means")

Problem illustrated:

Answer 1

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.

Answer 2

@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:

library(reshape2)
set.seed(5)
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.)