[R] A model-building strategy in mixed-effects modelling

[Stats Wolf]

Dear all,

Consider a completely randomized block design (let's use data(Oats)
irrespoctive of the split-plot design it was arranged in). Look:

library(nlme)
fit <- lme(yield ~ nitro, Oats, random = ~1|Block, method="ML")
fit2 <- lm(yield ~ nitro + Block, Oats)
anova(fit, fit2)
gives this:
     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit      1  4 624.3245 633.4312 -308.1623
fit2     2  8 611.9309 630.1442 -297.9654 1 vs 2 20.39366   4e-04

Clearly, considering block a random term is worse than considering it
a fixed term. Let's see if blocking should be included in the model at
all:

fit3 <- lm(yield ~nitro, Oats)
anova(fit2,fit3)

which gives a very small P value in favor of fit2, which suggests the
block term should be included. So, I go for the second model, with
block considered fixed.

Is this indeed how I should generally proceed when choosing the
optimum model for a situation that calls for mixed effects? Of course,
the example above is overly simplistic, yet such situations can occur
-- from a complex model with a couple of random terms one can finally
get to a simple fixed-effects model. Please comment.

Thanks in advance,

Stats Wolf