# When to use Kenward-Roger, not standard glmmTMB

DF <- data.frame(ID = c(5,6,8,10,10,11,12,14,14,15,15,16,17,17,18,23,25,26),
time = factor(c(1,24,24,1,24,1,1,6,24,6,96,6,6,96,96,24,24,24)),
CB = c(78.66,123.87,121.93,109.32,110.29,179.81,181.07,102.37,141.4,50.13,92.36,82.64,106.15,194.5,128.6,95.33,152.88,124.16),
M = c(103,124.1,97.3,102.1,106.8,80.6,89.3,107.2,96.5,131.3,101.9,122.5,110.4,95.2,94.4,114.2,112.6,97.5),
CC = c(120.1,148.1,126.3,116.7,133.4,159.1,139.1,111.9,138.2,82.3,131.7,107.1,116.3,162,111.8,107,160.6,129.3))

I have repeated measures data of 14 individuals distributed unevenly over 4 time points, and I'm analyzing the effect of time on CC. I have fit the data with both glmmTMB and lmer, the latter package has more complex tricks that make results more conservative, but I'm puzzled by the output and I'm wondering if I need it at all.

library(lme4)
library(glmmTMB)
library(emmeans)
model.glmmTMB <- glmmTMB(CC ~ time + (1|ID), data=DF, REML=T)
emm.glmmTMB <- emmeans(model.glmmTMB, specs=~time)
ph.glmmTMB <- contrast(emm.glmmTMB, "consec", simple=list("time"), adjust="bonferroni")
emmplot <- emmip(emm.glmmTMB, ~time, CI=T)

Judging from emmplot, I see some significant effects that I would expect based on the EMM. Even though my N isn't very small, the data is very sparsely distributed over my time points, so I tried another approach with lme4's Kenward-Roger df correction.

model.lmer <- lmer(CC ~ time + (1|ID), data=DF, REML=T)
emm.lmer <- emmeans(model.lmer, specs=~time)
ph.lmer <- contrast(emm.lmer, "consec", simple=list("time"), adjust="bonferroni")

Looking at ph.lmer, it's clear the correction has wiped any effect previously off the map. We can also see why: it completely decimated the df. However this doesn't seem to happen when regressing against CB or M, and it also doesn't happen when I log transform CC. Due to my limited knowledge of LMMs, I'm unable to figure out which strategy I have to apply now.

I have two questions stemming from this issue:

1. I'm puzzled by the complete lack of effects in the face of clear differences in EMMs after applying KR, and the fact that df has been reduced to ~1 makes me wonder if the algorithm tripped up somewhere. But perhaps I'm hurting my head over this for nothing. I read that KR was developed for situations where n<<p, which not really the case here, but I want to be sure: is the Kenward-Roger correction appropriate/necessary in my situation?
2. If the reduced df is correct, how is it possible that df isn't reduced as much when regressing M, or CB (whose sd is much higher)?

If this is your actual full dataset and you think it's important to account for subject effects, you're in real trouble, and whether or not to use K-R is the least of your problems. For most subjects, you have only one observation,so the subject effects seriously confound the time effects, and it's a fool's errand to try to sort them out.

Your only reasonable option with this dataset is to fit a regular linear model. with time as the sole predictor. Get consulting help to design your next study so that it has a chance of providing useful information.

There is a saying about trying to make a silk purse out of a sow's ear; it applies here.

### Follow-up

I decided to see more specifically how four models compare -- (1) a simple linear fit, (2) the mixed model considered here, (3) the same model as (2), only with fixed subject effects, and (4) the first model with added indicator variables for the four subjects having repeat measurements.

DF <- data.frame(ID = c(5,6,8,10,10,11,12,14,14,15,15,16,17,17,18,23,25,26),
time = factor(c(1,24,24,1,24,1,1,6,24,6,96,6,6,96,96,24,24,24)),
CB = c(78.66,123.87,121.93,109.32,110.29,179.81,181.07,102.37,141.4,50.13,92.36,82.64,106.15,194.5,128.6,95.33,152.88,124.16),
M = c(103,124.1,97.3,102.1,106.8,80.6,89.3,107.2,96.5,131.3,101.9,122.5,110.4,95.2,94.4,114.2,112.6,97.5),
CC = c(120.1,148.1,126.3,116.7,133.4,159.1,139.1,111.9,138.2,82.3,131.7,107.1,116.3,162,111.8,107,160.6,129.3))

DF = transform(DF,
i10 = 0 + (ID == 10),
i14 = 0 + (ID == 14),
i15 = 0 + (ID == 15),
i17 = 0 + (ID == 17))

require(lme4)
require(emmeans)

DF.lm = lm(CC ~ time, data = DF)
DF.lmer = lmer(CC ~ time + (1|ID), data = DF)
DF.fix = lm(CC ~ time + factor(ID), data = DF)
DF.lmx = lm(CC ~ time + i10 + i14 + i15 + i17, data = DF)

emmeans(DF.lm, consec ~ time)
## $$emmeans ## time emmean SE df lower.CL upper.CL ## 1 134 9.31 14 113.8 154 ## 6 104 9.31 14 84.4 124 ## 24 135 7.04 14 119.6 150 ## 96 135 10.75 14 112.1 158 ## ## Confidence level used: 0.95 ## ##$$contrasts
##  contrast estimate   SE df t.ratio p.value
##  6 - 1     -29.350 13.2 14 -2.230  0.1083
##  24 - 6     30.300 11.7 14  2.597  0.0556
##  96 - 24     0.467 12.8 14  0.036  1.0000
##
## P value adjustment: mvt method for 3 tests
emmip(DF.lm, ~ time)

emmeans(DF.lmer, consec ~ time)
## $$emmeans ## time emmean SE df lower.CL upper.CL ## 1 119 6.41 13.4 105 133 ## 6 108 6.30 13.8 94 121 ## 24 135 5.99 13.7 122 148 ## 96 155 6.58 12.5 140 169 ## ## Degrees-of-freedom method: kenward-roger ## Confidence level used: 0.95 ## ##$$contrasts
##  contrast estimate   SE   df t.ratio p.value
##  6 - 1       -11.7 4.54 1.15 -2.583  0.3368
##  24 - 6       27.3 3.16 1.11  8.661  0.1061
##  96 - 24      19.7 3.88 1.12  5.085  0.1788
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: mvt method for 3 tests
emmip(DF.lmer, ~ time)

emmeans(DF.fix, consec ~ time)
## $$emmeans ## time emmean SE df lower.CL upper.CL ## 1 118 2.20 1 90.1 146 ## 6 108 1.96 1 83.5 133 ## 24 135 1.19 1 119.6 150 ## 96 156 2.51 1 124.1 188 ## ## Results are averaged over the levels of: ID ## Confidence level used: 0.95 ## ##$$contrasts
##  contrast estimate   SE df t.ratio p.value
##  6 - 1        -9.6 3.70  1 -2.595  0.3385
##  24 - 6       26.3 2.62  1 10.052  0.0924
##  96 - 24      21.2 3.20  1  6.632  0.1392
##
## Results are averaged over the levels of: ID
## P value adjustment: mvt method for 3 tests
emmip(DF.fix, ~ time)

emmeans(DF.lmx, consec ~ time)
## $$emmeans ## time emmean SE df lower.CL upper.CL ## 1 146 21.4 10 98.7 194 ## 6 107 15.3 10 72.6 141 ## 24 145 19.7 10 100.9 189 ## 96 138 16.9 10 100.4 176 ## ## Results are averaged over the levels of: i10, i14, i15, i17 ## Confidence level used: 0.95 ## ##$$contrasts
##  contrast estimate   SE df t.ratio p.value
##  6 - 1      -39.64 16.3 10 -2.429  0.0829
##  24 - 6      38.14 14.2 10  2.684  0.0545
##  96 - 24     -6.84 16.0 10 -0.428  0.9418
##
## Results are averaged over the levels of: i10, i14, i15, i17
## P value adjustment: mvt method for 3 tests
emmip(DF.lmx, ~ time)

Created on 2021-02-09 by the reprex package (v0.3.0)

Note that models (1) and (4) have similar EMMs with the first, 3rd, and 4th times about equal and the 2nd one lower; but the SEs for the EMMs in model (4) are much higher. Models (2) and (3) have similar predictions, with the one the last time way higher. As mentioned in the OP, the df for the comparisons in model (2) are near 1 -- really indicating that we can't do those comparisons. MIn model (3), the df are exactly 1. The SEs of the means are way low but with the small d.f., the CIs are actually wider than for model (2). In model (4), we account for the subject variations only where we can assess them -- which "costs" us 4 d.f. for error. We still have 10 df left. The tests of comparisons are similar to those of model (1), and the SEs of the comparisons are not too much higher.

To be honest, I wouldn't attach a whole lot of meaning to these interpretations, but it seems to me that adding indicators only for subjects having repeat measurements is one way to account for subject variations in a situation like this.

• That makes sense. These data are the results of additional experiments which were very prone to error (with no chance for a do-over) and a lot of them failed as a result. I thought regular linear modeling was out of the question because some data points are correlated and so I hoped LMM would be a viable option to account for that. Thanks for your reply! Feb 8 at 8:06