Thanks a lot for your reply. This helps a lot!
Just to confirm, using lm this model will give me the mean yield value for each cell in the two way array. Now if I want to obtain the mean
of group means (like a SS type III approach) using LME (since I have random effects in the model) how can I parametrize this?
I could definitivelly use xtabs in a two-way case but in my case I have 2 other (continuous) covariates that are potential confounders in
the model so I need to keep them to obtain the corrected means.
I added a continuous variable (NewVar) to the dataset Newxmp11.07 and obtaned a model with the covariate (fm4) and another without it (fm3)


Newxmp11.07<-fix(xmp11.07)
str(Newxmp11.07)
`data.frame':   36 obs. of  4 variables:
$ Yield  : num  10.5 9.2 7.9 12.8 11.2 13.3 12.1 12.6 14 10.8 ...
$ Variety: Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
$ Density: Factor w/ 4 levels "1","2","3","4": 1 1 1 2 2 2 3 3 3 4 ...
$ NewVar : num  10 9 7 12 11 11 12 11 15 16 ...

fm3<-gls(Yield~-1+Variety + Density, xmp11.07)
fm4<-gls(Yield~-1+Variety + Density+NewVar, Newxmp11.07)
fm3
Coefficients:
Variety1  Variety2  Variety3  Density2  Density3  Density4
8.922222  9.797222 15.713889  2.911111  4.300000  2.433333

Degrees of freedom: 36 total; 30 residual
Residual standard error: 1.239243

fm4
Coefficients:
  Variety1    Variety2    Variety3    Density2    Density3    Density4
8.75757265  9.70589316 15.43347009  2.88152564  4.32186752  2.40117521
    NewVar
0.01157692

Degrees of freedom: 36 total; 29 residual
Residual standard error: 1.252991

fm4 gives me the mean of the group means for all the varieties but apparently it gives me the treatment contrasts for the densities. If I change the order of the factors in the model specification I get

coef(fm5<-gls(Yield~-1+ Density+Variety+NewVar, Newxmp11.07))
  Density1    Density2    Density3    Density4    Variety2    Variety3
8.75757265 11.63909829 13.07944017 11.15874786  0.94832051  6.67589744
    NewVar
0.01157692

This, just like fm4 will include the original intercept value in Density 1 which is not the actual density 1 mean. What am I missing? I am sorry if these questions are very basic but I want to make sure that I understand what I am doing. I guess that this is the price that I am paying for having used in the past packages like SAS where you just ask for lsmeans and the software will give you a "black box" answer!

Best regards

Francisco


From: Douglas Bates <[EMAIL PROTECTED]>
To: "Francisco Vergara" <[EMAIL PROTECTED]>
CC: [EMAIL PROTECTED]
Subject: Re: [R] Example of cell means model
Date: 15 Oct 2003 08:40:33 -0500

This is an example from chapter 11 of the 6th edition of Devore's
engineering statistics text.  It happens to be a balanced data set in
two factors but the calculations will also work for unbalanced data.

I create a factor called 'cell' from the text representation of the
Variety level and the Density level using '/' as the separator
character. The coefficients for the linear model 'Yield ~ cell - 1'
are exactly the cell means.

> library(Devore6)
> data(xmp11.07)
> str(xmp11.07)
`data.frame': 36 obs. of 3 variables:
$ Yield : num 10.5 9.2 7.9 12.8 11.2 13.3 12.1 12.6 14 10.8 ...
$ Variety: Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
$ Density: Factor w/ 4 levels "1","2","3","4": 1 1 1 2 2 2 3 3 3 4 ...
> xmp11.07$cell = with(xmp11.07, factor(paste(Variety, Density, sep = '/')))
> xtabs(~ Variety + Density, xmp11.07)
Density
Variety 1 2 3 4
1 3 3 3 3
2 3 3 3 3
3 3 3 3 3
> means = xtabs(Yield ~ Variety+Density, xmp11.07)/xtabs(~ Variety + Density, xmp11.07)
> means
Density
Variety 1 2 3 4
1 9.200000 12.433333 12.900000 10.800000
2 8.933333 12.633333 14.500000 12.766667
3 16.300000 18.100000 19.933333 18.166667
> coef(fm1 <- lm(Yield ~ cell - 1, xmp11.07))
cell1/1 cell1/2 cell1/3 cell1/4 cell2/1 cell2/2 cell2/3 cell2/4
9.200000 12.433333 12.900000 10.800000 8.933333 12.633333 14.500000 12.766667
cell3/1 cell3/2 cell3/3 cell3/4
16.300000 18.100000 19.933333 18.166667


The residual sum of squares for this model is the same as that for the
model with crossed terms

> deviance(fm1)
[1] 38.04
> deviance(fm2 <- lm(Yield ~ Variety * Density, xmp11.07))
[1] 38.04

but the coefficients are quite different because they represent a
different parameterization.

> coef(fm2)
      (Intercept)          Variety2          Variety3          Density2
       9.20000000       -0.26666667        7.10000000        3.23333333
         Density3          Density4 Variety2:Density2 Variety3:Density2
       3.70000000        1.60000000        0.46666667       -1.43333333
Variety2:Density3 Variety3:Density3 Variety2:Density4 Variety3:Density4
       1.86666667       -0.06666667        2.23333333        0.26666667

I hope this answers your question.  Sorry for the delay in
responding to you.

"Francisco Vergara" <[EMAIL PROTECTED]> writes:

> Many thanks for your reply.  An example would be very helpful to make
> sure that I understand correctly what you described.

_________________________________________________________________
Surf and talk on the phone at the same time with broadband Internet access. Get high-speed for as low as $29.95/month (depending on the local service providers in your area). https://broadband.msn.com


______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help

Reply via email to