Bug 17236 - Degrees of freedom in lme() -- incorrect handling of intercept
Summary: Degrees of freedom in lme() -- incorrect handling of intercept
Status: UNCONFIRMED
Alias: None
Product: R
Classification: Unclassified
Component: Models (show other bugs)
Version: R 3.3.*
Hardware: All All
: P5 normal
Assignee: R-core
URL:
Depends on:
Blocks:
 
Reported: 2017-03-16 18:03 UTC by Russ Lenth
Modified: 2017-03-17 02:01 UTC (History)
0 users

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Russ Lenth 2017-03-16 18:03:52 UTC
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.
Comment 1 Russ Lenth 2017-03-17 02:01:23 UTC
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.