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