[R-sig-ME] output of MCMCglmm

[John Maindonald]

The best way to understand this is to write down the formula that
gives the fitted values.

Remember that fitted values are constructed by adding together
parameter estimates.  If the mean is incorporated into the estimates
for the first factor, it is not available for incorporation into the estimates
for any later model term.

If the model is balanced, you can get means for all levels of all factors
from model.tables()  If the model is not balanced, what you get
corresponds to a sequential breakdown a/c model terms, and can be
very misleading.

Observe:
================================================================

> model.tables(aov(aftertaste ~ panelist+product, data=appletaste),  type="means")
Tables of means
Grand mean

61.05 

panelist 
panelist
   a      b      c      d      e      f      g      h      i      j      k 
98.33  57.33 103.33  73.33  72.00  63.67  38.00  55.33  86.33  91.33  21.00 
   l      m      n      o      p      q      r      s      t 
39.00  67.67  58.33  37.00  66.00  50.67  68.00  53.33  21.00 

product 
product
298   493   649   937 
70.29 88.58 58.41 26.92 
> model.tables(aov(aftertaste ~ product+panelist, data=appletaste), type="means")
Tables of means
Grand mean

61.05 

product 
product
298   493   649   937 
68.33 94.07 51.80 30.00 

panelist 
panelist
   a      b      c      d      e      f      g      h      i      j      k 
95.25  54.25 100.25  70.25  68.92  66.09  40.43  57.76  88.76  93.76  32.01 
   l      m      n      o      p      q      r      s      t 
50.01  78.67  69.34  48.01  55.65  40.32  57.65  42.98  10.65 

================================================================

Not only are the product means different, but the estimates of differences
of product effects are different.  This is a balanced incomplete design.  
Here, the SEDs are the same for all pairwise differences of treatments.
Model designs can (and do) get a lot worse, wrt issues of getting sensible 
model summaries.

The output you should believe is what you get from:

summary.lm(aov(aftertaste ~ product+panelist, data=appletaste))
or
summary.lm(aov(aftertaste ~ 0+product+panelist, data=appletaste))
or 
summary.lm(aov(aftertaste ~ panelist+product, data=appletaste))
or . . .  
[You can also change the model contrasts used; see help(C)]

It is a bit disturbing that model.tables() is now allowed to input without
warning output that, relative to what the user might reasonably expect,
is complete nonsense.  

NB that termplot() does give means that are adjusted for all other model 
terms, but as it stands only for main effect models.  This is what it would 
be sensible for a proper model.tables() function to give, albeit what to
do in the presence of interactions has to be thought through carefully.

John Maindonald             email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.
http://www.maths.anu.edu.au/~johnm

On 16/03/2011, at 11:58 PM, Szymek Drobniak wrote:

> Hi all,
> I've one small question: after removing the Intercept from a model
> with >2 fixed effects I'm getting in the output all levels of the
> first factor, but only N-1 levels of the remaining factors. Do I
> understand correctly that they're presented as differences with
> reference levels, which would be the first (alphabetically) levels in
> each fixed effect? Is it possible to obtain the output that will list
> the means (and not differences) for all levels of all fixed factors?
> Cheers,
> szymek
> 
> Institute of Environmental Sciences, Jagiellonian University
> ul. Gronostajowa 7, 30-387 Kraków, POLAND
> tel.: +48 12 664 52 19 fax: +48 12 664 69 12
> www.eko.uj.edu.pl/drobniak
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models