MichaelChirico
commented
4 years ago Note how the degrees of freedom work with the aov() function. I show this because its calculations are pretty closely akin with those based on groupings such as are done in lme():
data(Oats, package = "nlme")Oats.aov <- aov(yield ~ Variety * factor(nitro) + + Error(Block/Variety), data = Oats)
Oats.aov1 <- aov(yield ~ -1 + Variety * factor(nitro) + + Error(Block/Variety), data = Oats)
summary(Oats.aov)Error: Block Df Sum Sq Mean Sq F value Pr(>F) Residuals 5 15875 3175
Error: Block:Variety Df Sum Sq Mean Sq F value Pr(>F) Variety 2 1786 893.2 1.485 0.272 Residuals 10 6013 601.3
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
factor(nitro) 3 20020 6673 37.686 2.46e-12
Variety:factor(nitro) 6 322 54 0.303 0.932
Residuals 45 7969 177
In the model with intercept, the Block stratum, with 5 d.f., has no model terms listed - though I maintain the intercept is in that stratum. The remaining d.f. are consistent with those obtained by anova(Oats.lme).
summary(Oats.aov1)Error: Block
Df Sum Sq Mean Sq F value Pr(>F)
Variety 1 778336 778336 245.1 1.93e-05
Residuals 5 15875 3175
Error: Block:Variety Df Sum Sq Mean Sq F value Pr(>F) Variety 2 1786 893.2 1.485 0.272 Residuals 10 6013 601.3
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
factor(nitro) 3 20021 6674 37.686 2.46e-12
Variety:factor(nitro) 6 322 54 0.303 0.932
Residuals 45 7969 177
In the model without intercept, the strata have the same respective degrees of freedom. The sums of squares and d.f. in the second and third strata are identical to those for the model with intercept, but in the top stratum, Variety appears with 1 d.f. This extra degrees of freedom and sum of squares is in fact attributable to the intercept, as can be verified by explicitly specifying the intercept:
Oats$one = 1
Oats.aovone <- aov(yield ~ -1 + one + Variety * factor(nitro) +
Error(Block/Variety), data = Oats)
summary(Oats.aovone)Error: Block
Df Sum Sq Mean Sq F value Pr(>F)
one 1 778336 778336 245.1 1.93e-05
Residuals 5 15875 3175
... etc. ...
I believe that it is this sum of squares that is being ignored by lme() when calculating its degrees of freedom with the no-intercept model. And, it further supports the fact that the error stratum for the intercept is Blocks, not residual.
METADATA
- Comment author - Russ Lenth
- Timestamp - 2017-03-17 02:01:23 UTC
github-actions[bot]
I believe that lme() handles the intercept incorrectly in calculating degrees of freedom. For example, consider these two models with and without intercept:
Oats.lme <- lme(yield∼ Variety * factor(nitro), random =∼ 1 | Block/Variety, data = Oats)
Oats.lme1 <- lme(yield∼ Variety * factor(nitro) - 1, random =∼ 1 | Block/Variety, data = Oats)
Compare their anova tables:
R> anova(Oats.lme) numDF denDF F-value p-value (Intercept) 1 45 245.14299 <.0001 Variety 2 10 1.48534 0.2724 factor(nitro) 3 45 37.68562 <.0001 Variety:factor(nitro) 6 45 0.30282 0.9322
R> anova(Oats.lme1) numDF denDF F-value p-value Variety 3 10 82.70456 <.0001 factor(nitro) 3 46 37.68562 <.0001 Variety:factor(nitro) 6 46 0.30282 0.9323
The d.f. for nitro-related terms differ, even though both models have the same fitted values and 12 fixed parameters.
Also, I believe the intercept d.f. are incorrect in the first model. As I understand it, these degrees of freedom are obtained using a "containment" method, whereby, due to the grouping structure, the d.f. for a term can be no greater than the d.f. for any term that contains it. In this case, the d.f. for Variety and factor(nitro) are both <= the d.f. for their interaction. Since the intercept is contained in all terms, its d.f. should not be more than 10. In fact, I believe that the intercept in this particular example should have 5 d.f., because there are 6 blocks at the coarsest level of grouping, hence 6 - 1 = 5 d.f.
I can support these assertions further by comparing with the equivalent model (with intercept) fitted using lme4::lmer() with "contr.sum" contrasts, and applying pbkrtest::get_ddf_Lb() to obtain d.f. for each single coefficient. I obtain 5 df for the intercept, 10 df for each Variety coefficient, 45 df for the other coefficients. I can provide those details if you think it will be helpful. (While the K-R method is not the same as containment, it should yield the same results in this example since it is a balanced design.)
These issues can be repaired, I believe, by appropriate changes to the function nlme:::getFixDF(). I can't really follow too precisely what is going on in there, but that function does work through the grouping structure to obtain its results, but somehow handles the intercept incorrectly.
METADATA