Dear all, I am analyzing 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 seperate 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 loose any degrees of freedom for modelling these seperate variances! Does anyone know why this is?
Below a simple example dataset with 4 observations in each of 5 treatments (the sample size is insufficient to show the differences in variance 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. eventhough 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 -- View this message in context: http://r-sig-ecology.471788.n2.nabble.com/heteroscedasticity-and-d-f-tp7578507.html Sent from the r-sig-ecology mailing list archive at Nabble.com. _______________________________________________ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology