[R-sig-ME] output of MCMCglmm
John Maindonald
john.maindonald at anu.edu.au
Wed Mar 16 23:23:16 CET 2011
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
More information about the R-sig-mixed-models
mailing list