[R] model.tables

[Berwin Turlach]

>>>>> "JM" == John Maindonald <john.maindonald at anu.edu.au> writes:

  JM> At 08:02 30/11/99 +0000, Prof Brian D Ripley wrote:
  >> On Tue, 30 Nov 1999, spoon <spoon at hilbert.maths.utas.edu.au> wrote:
  >> 
  >>> Hi,
  >>> Is this a bug or do I just not understand model.tables?
  >>> [...]
  >>> Or am I just completely misinterpreting something basic?
  >> 
  >> Basically, yes. This is an incompletely replicated design, and d and e
  >> occur on different litters.
  >> 
  >> The results are identical to the S-PLUS original. I think you are
  >> probably looking for what dummy.coef gives you.
  >> [...]

  JM> So what is it that model.tables() gives?
Good question, in R (0.65.1) the documentation of model.tables says:

   WARNING:
        The implementation is incomplete, and only the simpler
        cases have been tested thoroughly.

  JM> S-PLUS says they are estimates.  Of what?  These are not the
  JM> marginal means of the fitted values, ignoring other factors.  Do
  JM> they mean anything at all?  The issue of unbalance does not
  JM> arise here.
The documentation in R also states that

    Details:

        For `type = "effects"' give tables of the coefficients
        for each term, optionally with standard errors.

        For `type = "means"' give tables of the mean response
        for each combinations of levels of the factors in a
        term.

Hence, my understanding would be that if `type="means"' is specified
the marginal means of the observations ignoring other factors is
given.  But this doesn't seem to be the case, at least not with an
unbalanced design.  You can use either your dreamed up example or the
data that spoon has posted.

> bdes_data.frame(blk=rep(1:3,rep(2,3)),
+    trt=c(1,2,2,3,3,1), y=c(.4,.6,1,1.4,2,1.2))
> bdes$blk_factor(bdes$blk,labels=LETTERS[1:3])
> bdes$trt_factor(bdes$trt,labels=letters[1:3])
> tapply(bdes$y,bdes$trt,mean)
  a   b   c 
0.8 0.8 1.7 
> bdes.aov_aov(y~blk+trt,data=bdes,projections=T)
> model.tables(bdes.aov,"mean")
Tables of means
Grand mean
    
1.1 

 blk 
  A   B   C 
0.5 1.2 1.6 

 trt 
   a    b    c 
0.85 1.05 1.40

However, model.tables seems to be sensitive regarding the way the
factors have been specified.  If I change the order in the model, I
get the results that I expect.  But note that the results for the
factor blk change. 

> bdes1.aov_aov(y~trt+blk,data=bdes,projections=T)
> model.tables(bdes1.aov,"mean")
Tables of means
Grand mean
    
1.1 

 trt 
  a   b   c 
0.8 0.8 1.7 

 blk 
   A    B    C 
0.80 1.05 1.45 


Note, if I take the data that is on the R help page for model.tables
(which is balanced) then I get the results that I expect:

> N <- c(0,1,0,1,1,1,0,0,0,1,1,0,1,1,0,0,1,0,1,0,1,1,0,0)
> P <- c(1,1,0,0,0,1,0,1,1,1,0,0,0,1,0,1,1,0,0,1,0,1,1,0)
> K <- c(1,0,0,1,0,1,1,0,0,1,0,1,0,1,1,0,0,0,1,1,1,0,1,0)
> yield <- c(49.5,62.8,46.8,57.0,59.8,58.5,55.5,56.0,62.8,55.8,69.5,
+ 55.0, 62.0,48.8,45.5,44.2,52.0,51.5,49.8,48.8,57.2,59.0,53.2,56.0)
> npk <- data.frame(block=gl(6,4), N=factor(N), P=factor(P),
+ K=factor(K), yield=yield)
> npk.aov <- aov(yield ~ block + N*P*K, npk)
> model.tables(npk.aov, "means")
Tables of means
Grand mean
       
54.875 

 block 
    1     2     3     4     5     6 
54.03 57.45 60.77 50.12 50.52 56.35 

 N 
    0     1 
52.07 57.68 

 P 
    0     1 
55.47 54.28 

 K 
    0     1 
56.87 52.88 

 N:P 
   P
N   0     1    
  0 51.72 52.42
  1 59.22 56.15

 N:K 
   K
N   0     1    
  0 52.88 51.25
  1 60.85 54.52

 P:K 
   K
P   0     1    
  0 57.60 53.33
  1 56.13 52.43
> mean(npk[,"yield"])
[1] 54.875
> tapply(npk$yield,npk$block,"mean")
     1      2      3      4      5      6 
54.025 57.450 60.775 50.125 50.525 56.350 
> tapply(npk$yield,npk$N,"mean")
       0        1 
52.06667 57.68333 
> tapply(npk$yield,npk$P,"mean")
       0        1 
55.46667 54.28333 
> tapply(npk$yield,npk$K,"mean")
       0        1 
56.86667 52.88333 
> tapply(npk$yield,list(npk$N,npk$P),"mean")
         0        1
0 51.71667 52.41667
1 59.21667 56.15000
> tapply(npk$yield,list(npk$N,npk$K),"mean")
         0        1
0 52.88333 51.25000
1 60.85000 54.51667
> tapply(npk$yield,list(npk$P,npk$K),"mean")
         0        1
0 57.60000 53.33333
1 56.13333 52.43333


Cheers,

        Berwin
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._