R help - Oct 2005 - The mathematics inside lme()

Hello all!

Consider a dataset with a grouping structure, Group (factor)
Several treatments, Treat (factor)
Some sort of yield, Yield (numeric)
Something, possibly important, measured for each group; GroupCov (numeric)

To look for fixed effects from Treat on Yield, a first attempt could be:

m1 <- lm(Yield ~ Treat)

which gives, in a symmetric situation, the same estimated fixed effects as:

m2 <- lme(Yield ~ Treat,
          random =~1| Group)

but m2 is a much better model with safer significances.

Now I want to evaluate GroupCov as a covariate to Treat. I can then start
with either m1 or m2 as base, but what is most correct when GroupCov has
only one value for each Group?

m3 <- lm(Yield ~ Treat + GroupCov + Treat:GroupCov)

gives the same fixed effects as

m4 <- lme(Yield ~ Treat + GroupCov + Treat:GroupCov,
          random =~1| Group)

but this time the prob.values for GroupCov is much stronger in m3 than in
m4. Needless to say, anova(m3,m4) tells that m4 is a better *model* than
m3. But is it better for my purpose? Trying an old-fashioned style model
with only fixed effects? (Don??t shout at me, it is only a dirty test of
the system):

m5 <- lm(Yield ~ Group + Treat + GroupCov + Treat:GroupCov)

is accepted by lm() but here GroupCov is silently removed from the
analysis by lm(). I accept this removal, but I get even more suspicious
that something fishy is going on in m4. My gut-feeling is that m3 is the
right starting point but I have got a general recommendation always to
start with m4-type calls when evaluating numeric covariates. I need one
(or two) "doctor second opinion" on this.

How is the mathematics inside lme() working? Is some part of the variation
I want to catch as an effect from GroupCov already removed by the random
call, or why do I get better significances in the pure fixed call? Could
these sigificances be some sort of artefact?

Cheers
/CG

PS.
I sent basically this question, but in a more special case and with
another header, to the list two days ago. Nobody was interested, hopefully
this is more tasty ;-)
DS.


-- 
CG Pettersson, MSci, PhD Stud.
Swedish University of Agricultural Sciences (SLU)
Dep. of Crop Production Ekology. Box 7043.
SE-750 07 Uppsala, Sweden
cg.pettersson at evp.slu.se

> Now I want to evaluate GroupCov as a covariate to Treat. I 
> can then start with either m1 or m2 as base, but what is most
> correct when GroupCov has only one value for each Group?
> 
> m3 <- lm(Yield ~ Treat + GroupCov + Treat:GroupCov)
> 
> gives the same fixed effects as
> 
> m4 <- lme(Yield ~ Treat + GroupCov + Treat:GroupCov,
>           random =~1| Group)
> 
> but this time the prob.values for GroupCov is much stronger 
> in m3 than in m4. Needless to say, anova(m3,m4) tells that m4
> is a better *model* than m3. But is it better for my purpose?
Well, I do not actually know what your purpose is ...
... but in my oppinion the second model is much better (and I am tempted to
say the first one is wrong).

The crucial point here is that there is only one value of GroupCov in each
Group. Thus the number of replications that provide degrees of freedom for
the effect of GroupCov is the number of groups. m4 adjusts for this fact,
has a lower df for the GroupCov and thus a lower p-value.

In m3, you model as if all your observations are independent for all
variables. This is actually the case for none but may become more visible
for GroupCov because this variable is constant for all units within group
(and thus this value is certainly not independent).

Cheers, Lorenz
- 
Lorenz Gygax, Dr. sc. nat.
Centre for proper housing of ruminants and pigs
Swiss Federal Veterinary Office
agroscope FAT T??nikon, CH-8356 Ettenhausen / Switzerland