[R] lme and Assay data: Test for block effect when block is systematic - anova/summary goes wrong
Søren Højsgaard
Soren.Hojsgaard at agrsci.dk
Tue Feb 7 01:02:14 CET 2006
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
More information about the R-help
mailing list