lme and Assay data: Test for block effect when block is systematic - anova/summary goes wrong

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

lme and Assay data: Test for block effect when block is systematic - anova/summary goes wrong

Søren Højsgaard
Consider the Assay data where block, sample within block and dilut within block is random.
This model can be fitted with (where I define Assay2 to get an ordinary data frame rather
than a grouped data object):
 
Assay2 <- as.data.frame(Assay)
fm2<-lme(logDens~sample*dilut, data=Assay2,
  random=list(Block = pdBlocked(list(pdIdent(~1), pdIdent(~sample-1),pdIdent(~dilut-1))) ))
 
Now, block has only 2 levels so I prefer to treat it as fixed:
 
fm3<-lme(logDens~Block+sample*dilut, data=Assay2,
  random=list(Block = pdBlocked(list(pdIdent(~sample-1),pdIdent(~dilut-1))) ))
 
This works fine until I try a summary() or anova(), e.g.
 
> anova(fm3)
             numDF denDF  F-value p-value
(Intercept)      1    29 824.2612  <.0001
Block              1     0   1.5320     NaN
sample            5    29  11.2133  <.0001
dilut               4    29 420.7901  <.0001
sample:dilut    20    29   1.6069  0.1193
Warning messages:
1: NaNs produced in: pf(q, df1, df2, lower.tail, log.p)
2: NAs introduced by coercion
 
The output from SAS (when I request residual denominator dfs)
     Type 1 Tests of Fixed Effects
                   Num     Den
Effect              DF      DF    F Value    Pr > F
Block                 1      29       1.53    0.2257
sample               5      29      11.21    <.0001
dilut                  4      29     420.79    <.0001
sample*dilut      20      29       1.61    0.1193
 
The test for block effect is usually not of interest. Moreover, I don't want to enroll into the "degree of freedom police force", but 0 dfs is certainly wrong. Maybe it would be better if lme simply used the residual df's (when no better choice is available??)
 
Best regards
Søren

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html