Why does the degre of freedom of one level of a cathegoric variable is different from the other levels in lme?

Question

I am doing mixed models with nested random effects. When I run the summary of the model I see that the number of degrees of freedom of the level "control fishless" of the variable "treatment 1" is 11 while the four other levels' df is 51. I would expect them to be all the same. What could be the cause of this and how can I fix the problem?

Thank you!

My model:

library(nlme)
library(emmeans)

data.score.merge2$treatment1 <- relevel(data.score.merge2$treatment1,"control fish")
mod.full.PC1 <- lme(PC1 ~ season + elevation + treatment1, method = "REML",
                random = ~1|lake/replicate,
                data = data.score.merge2)
summary(mod.full.PC1)
anova(mod.full.PC1)
pairs(emmeans(mod.full.PC1, ~treatment1))

My data (I put them all because it may be because of their distribution or something like that):

structure(list(sample = structure(1:102, .Label = c("Annette 1", 
"Annette 2", "Annette 3", "Annette 4", "Annette 5", "Annette 6", 
"Cobb1", "Cobb2", "Cobb3", "Cobb4", "Cobb5", "Cobb6", "Dog1", 
"Dog2", "Dog3", "Dog4", "Dog5", "Dog6", "Helen1", "Helen2", "Helen3", 
"Helen4", "Helen5", "Helen6", "Herbert 1", "Herbert 2", "Herbert 3", 
"Herbert 4", "Herbert 5", "Herbert 6", "Hidden 10 Months 1", 
"Hidden 10 Months 2", "Hidden 10 Months 3", "Hidden 10 Months 4", 
"Hidden 10 Months 5", "Hidden 10 Months 6", "Hidden 10 Months 7", 
"Hidden 10 Months 8", "Hidden After 1", "Hidden After 2", "Hidden After 3", 
"Hidden After 4", "Hidden After 5", "Hidden After 6", "Hidden After 7", 
"Hidden After 8", "Hidden Before 1", "Hidden Before 2", "Hidden Before 3", 
"Hidden Before 4", "Hidden Before 5", "Hidden Before 6", "Hidden Before 7", 
"Hidden Before 8", "Kootenay Pond 1", "Kootenay Pond 2", "Kootenay Pond 3", 
"Kootenay Pond 4", "Kootenay Pond 5", "Kootenay Pond 6", "Margaret1", 
"Margaret2", "Margaret3", "Margaret4", "Margaret5", "Margaret6", 
"McNair1", "McNair2", "McNair3", "McNair4", "McNair5", "McNair6", 
"Mud1", "Mud2", "Mud3", "Mud4", "Mud5", "Mud6", "Olive1", "Olive2", 
"Olive3", "Olive4", "Olive5", "Olive6", "Ross1", "Ross2", "Ross3", 
"Ross4", "Ross5", "Ross6", "Sentinel 1", "Sentinel 2", "Sentinel 3", 
"Sentinel 4", "Sentinel 5", "Sentinel 6", "Temple1", "Temple2", 
"Temple3", "Temple4", "Temple5", "Temple6"), class = "factor"), 
    replicate = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
    3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
    3L, 1L, 2L, 3L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 
    4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 
    3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
    3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
    3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
    3L), season = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 
    1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
    2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 
    1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
    2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 
    1L, 2L, 2L, 2L), .Label = c("july", "Sept"), class = "factor"), 
    elevation = c(1898L, 1898L, 1898L, 1898L, 1898L, 1898L, 1260L, 
    1260L, 1260L, 1260L, 1260L, 1260L, 1185L, 1185L, 1185L, 1185L, 
    1185L, 1185L, 2400L, 2400L, 2400L, 2400L, 2400L, 2400L, 1506L, 
    1506L, 1506L, 1506L, 1506L, 1506L, 2400L, 2400L, 2400L, 2400L, 
    2400L, 2400L, 2400L, 2400L, 2400L, 2400L, 2400L, 2400L, 2400L, 
    2400L, 2400L, 2400L, 2400L, 2400L, 2400L, 2400L, 2400L, 2400L, 
    2400L, 2400L, 1128L, 1128L, 1128L, 1128L, 1128L, 1128L, 1808L, 
    1808L, 1808L, 1808L, 1808L, 1808L, 1532L, 1532L, 1532L, 1532L, 
    1532L, 1532L, 1600L, 1600L, 1600L, 1600L, 1600L, 1600L, 1470L, 
    1470L, 1470L, 1470L, 1470L, 1470L, 1735L, 1735L, 1735L, 1735L, 
    1735L, 1735L, 2274L, 2274L, 2274L, 2274L, 2274L, 2274L, 2207L, 
    2207L, 2207L, 2207L, 2207L, 2207L), lake = structure(c(1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 
    3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 
    6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
    6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 
    8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 
    10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 
    12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 
    14L, 14L, 14L), .Label = c("Annette", "Cobb", "Dog", "Helen", 
    "Herbert", "Hidden", "Kootenay Pond", "Margaret", "McNair", 
    "Mud", "Olive", "Ross", "Sentinel", "Temple"), class = "factor"), 
    treatment1 = structure(c(5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 5L, 5L, 5L, 5L, 5L, 5L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
    4L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 
    1L, 1L, 1L, 1L), .Label = c("control fish", "10 month after", 
    "2 weeks", "Before", "control fish less"), class = "factor"), 
    PC1 = c(0.44326796769495, 0.391805088607907, 0.366947929408647, 
    0.441091779859967, 0.403451309938093, 0.290713511918726, 
    -0.116054225460659, 0.0837383015520615, -0.138878110788659, 
    -0.0964488183669892, -0.00381046526130373, -0.158529727412508, 
    0.228810116000753, 0.116624496704548, -0.0690976721539107, 
    0.241173083526596, 0.149149975299026, 0.144601131921389, 
    -0.132779778461221, -0.128217883448178, -0.131033607355805, 
    -0.214176904428602, -0.181864972831349, -0.178177705570812, 
    0.172493720206089, 0.208305957239218, 0.20607376930861, 0.238450520739572, 
    0.102790731251281, 0.0160873782197574, 0.302705601248893, 
    0.425953281002366, 0.274172226495253, 0.280315453212515, 
    0.0798037386796905, 0.506141404170783, 0.262641034377766, 
    0.400158058034833, 0.0532097743994839, -0.132559135111449, 
    -0.208628889600582, -0.283005491924375, 0.00488503259743488, 
    0.0532097743994839, -0.200151311512991, -0.19723781894123, 
    -0.2861581962982, -0.192654352006073, -0.260656011895704, 
    -0.222547061478245, -0.198133890113568, -0.201428667334639, 
    -0.303429745739617, -0.25748877830387, 0.438220569263306, 
    0.50013316765322, 0.516180868603105, 0.457174900484235, 0.191265288575061, 
    0.403629612339766, -0.032327404790195, -0.266767528170862, 
    -0.165290371763968, -0.253115189605474, -0.2661674781074, 
    -0.246530476585813, 0.161873213488864, 0.230102615861032, 
    0.247427735270813, 0.49154515223575, 0.332658929465825, 0.375955962214077, 
    -0.164941207044441, -0.14560006272253, -0.120381478479053, 
    -0.229436843745494, -0.135108208826027, -0.240164909215447, 
    0.132065444955584, -0.0264627444494227, -0.00763873712457839, 
    0.109028111315133, 0.0757425066027548, 0.0983200921588997, 
    -0.301499561057097, -0.226818786095393, -0.169404319537732, 
    -0.233550579992025, -0.19970763507962, -0.239831967809433, 
    -0.309404682344026, -0.311166025270994, -0.295989491984986, 
    -0.31253624380418, -0.300847881797013, -0.18667072522063, 
    -0.301875992307235, -0.308143568566698, -0.299702156473456, 
    -0.314275297789266, -0.306728420509431, -0.238861120432659
    )), class = "data.frame", row.names = c(NA, -102L))

Russ Lenth · Accepted Answer · 2020-04-28T22:09:06.427

Note the relationship between treatment1 and lake:

> with(data.score.merge2, table(lake, treatment1))
               treatment1
lake            control fish 10 month after 2 weeks Before control fish less
  Annette                  0              0       0      0                 6
  Cobb                     6              0       0      0                 0
  Dog                      6              0       0      0                 0
  Helen                    6              0       0      0                 0
  Herbert                  0              0       0      0                 6
  Hidden                   0              8       8      8                 0
  Kootenay Pond            0              0       0      0                 6
  Margaret                 6              0       0      0                 0
  McNair                   6              0       0      0                 0
  Mud                      6              0       0      0                 0
  Olive                    6              0       0      0                 0
  Ross                     6              0       0      0                 0
  Sentinel                 0              0       0      0                 6
  Temple                   6              0       0      0                 0

Three of the treatments appear only with Hidden lake, and no other treatment is observed on more than one lake. So comparisons among levels 10 month after, 2 weeks, and Before are within-lake comparisons, and all other comparisons are between-lake. The results I get are:

> pairs(emmeans(mod.full.PC1, ~treatment1))
 contrast                           estimate     SE df t.ratio p.value
 control fish - 10 month after       -0.6384 0.2210 84 -2.888  0.0385 
 control fish - 2 weeks              -0.2206 0.2209 84 -0.999  0.8552 
 control fish - Before               -0.0816 0.2210 84 -0.369  0.9960 
 control fish - control fish less    -0.2442 0.1142 11 -2.138  0.2710 
 10 month after - 2 weeks             0.4178 0.0470 84  8.885  <.0001 
 10 month after - Before              0.5568 0.0438 84 12.701  <.0001 
 10 month after - control fish less   0.3942 0.2314 11  1.704  0.4698 
 2 weeks - Before                     0.1390 0.0470 84  2.957  0.0321 
 2 weeks - control fish less         -0.0235 0.2313 11 -0.102  1.0000 
 Before - control fish less          -0.1626 0.2314 11 -0.703  0.9516 

Results are averaged over the levels of: season 
Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 5 estimates

[I get 84 d.f. where you say 51 -- that is a mystery.] I don't see, either, why there are not 11 d.f. for comparisons with control fish. These results use the containment method.

I tried this also with the approximate Satterthwaite method:

> pairs(emmeans(mod.full.PC1, ~treatment1, mode = "appx"))
 contrast                           estimate     SE    df t.ratio p.value
 control fish - 10 month after       -0.6384 0.2210  9.88 -2.888  0.0938 
 control fish - 2 weeks              -0.2206 0.2209  9.87 -0.999  0.8502 
 control fish - Before               -0.0816 0.2210  9.88 -0.369  0.9954 
 control fish - control fish less    -0.2442 0.1142  9.98 -2.138  0.2769 
 10 month after - 2 weeks             0.4178 0.0470 84.79  8.885  <.0001 
 10 month after - Before              0.5568 0.0438 84.79 12.701  <.0001 
 10 month after - control fish less   0.3942 0.2314  9.89  1.704  0.4735 
 2 weeks - Before                     0.1390 0.0470 84.79  2.957  0.0320 
 2 weeks - control fish less         -0.0235 0.2313  9.88 -0.102  1.0000 
 Before - control fish less          -0.1626 0.2314  9.89 -0.703  0.9512 

Results are averaged over the levels of: season 
Degrees-of-freedom method: appx-satterthwaite 
P value adjustment: tukey method for comparing a family of 5 estimates

These d.f. make a lot more sense, being distinctly higher only for the three treatments on Hidden Lake. All I can guess is that there is come bug in the containment method for d.f. connected with the fact that control fish is the reference level. I will investigate that.

Addendum

More on containment... If I change the contrasts (i.e., how the dummy variables are constructed), that can affect the d.f. from the containment method:

> mfp = update(mod.full.PC1, contrasts = list(treatment1 = "contr.sum"))
> pairs(emmeans(mfp, ~treatment1, mode = "cont"))
 contrast                           estimate     SE df t.ratio p.value
 control fish - 10 month after       -0.6384 0.2210 11 -2.888  0.0874 
 control fish - 2 weeks              -0.2206 0.2209 11 -0.999  0.8507 
 control fish - Before               -0.0816 0.2210 11 -0.369  0.9954 
 control fish - control fish less    -0.2442 0.1142 11 -2.138  0.2710 
 10 month after - 2 weeks             0.4178 0.0470 84  8.885  <.0001 
 10 month after - Before              0.5568 0.0438 84 12.701  <.0001 
 10 month after - control fish less   0.3942 0.2314 11  1.704  0.4698 
 2 weeks - Before                     0.1390 0.0470 84  2.957  0.0321 
 2 weeks - control fish less         -0.0235 0.2313 11 -0.102  1.0000 
 Before - control fish less          -0.1626 0.2314 11 -0.703  0.9516 

Results are averaged over the levels of: season 
Degrees-of-freedom method: containment 
P value adjustment: tukey method for comparing a family of 5 estimates

Note that only 3 comparisons have the larger d.f., more consistent with my intuition and the Satterthwaite results. My guess right now is that this has to do with the reference level for a factor. With the default "contr.treatment" contrasts, the reference level doesn't involve any predictors associated with the factor in question, and this messes up the d.f. calculations for comparisons with the reference level. I would speculate that this bug can't be fixed, and the workaround is to use other than "contr.treatment" contrasts when containment d.f. is used. But I will investigate further.

I just realised that I deleated the nested random effect by mistake. My model was actualy: mod.full.PC1 <- lme(PC1 ~ season + elevation + treatment1, method = "REML", random = ~1|lake/replicate, data = data.score.merge2) When I take off the "replicate" I also get 85. Sorry for the confusion. — Julien Beaulieu, Apr 28 '20 at 13:51

Why does the degre of freedom of one level of a cathegoric variable is different from the other levels in lme?

1 Answers1

Addendum