Dear all,
I am analysing an experiment where the treatments cause heteroscedasticty
in the residuals. Following Pinheiro & Bates (2000 Mixed effects models in
S and S-plus) and Zuur et al (2009 Mixed effects models and extensions in
ecology with R) I try to include separate treatment variances using
varIdent structure. This improves thing according to my residual plots and
a LR test. However, what I don't understand is that when I test for the
effect of my treatment, I don't seem to lose any degrees of freedom for
modelling these separate variances! Does anyone know why this is? Or know
sources where this is explained?

Below a simplified example dataset with 4 observations in each of 5
treatments (the sample size is insufficient to show the differences in
variances so the LR is not significant, but that's not my point). Comparing
both models I see a difference in d.f. among the models for the LR test of
variance structures (6 vs. 10), but when I look at the tests of fixed
effects for both models, both still have 15 residual d.f. even though the
second model needs to estimate 5 variances instead of one.

I am curious about your thoughts!

Best wishes,
Jasper Wubs

> library(nlme)
> m1 <- gls(SM~Treat,weights=NULL,data=data)
> vf1 <- varIdent(form=~1|Treat) #estimate a variance per treatment
> m1v <- gls(SM~Treat,weights=vf1,data=data)
> anova(m1,m1v)

    Model df       AIC       BIC   logLik   Test  L.Ratio p-value
m1      1  6 -85.76002 -81.51172 48.88001
m1v     2 10 -84.85906 -77.77855 52.42953 1 vs 2 7.099036  0.1307

> anova(m1)
Denom. DF: 15
            numDF   F-value p-value
(Intercept)     1 10246.756  <.0001
Treat           4     1.937  0.1564

> anova(m1v)
Denom. DF: 15
            numDF   F-value p-value
(Intercept)     1 27538.274  <.0001
Treat           4     7.733  0.0014

	[[alternative HTML version deleted]]

