The following example from the JSS paper on pkbrtest shows an LRT meant to test for the effect of the whole-plot treatment (harvest):

data("beets", package = "pbkrtest")
sug <- lmer(sugpct ~ block + sow + harvest + (1 | block:harvest),
            data = beets, REML = FALSE)
sug_no.harv <- update(sug, . ~ . - harvest)
sug_no.sow <- update(sug, . ~ . - sow)

anova(sug, sug_no.harv)

At the end of Section 3.2, the authors write that, in addition to inadequacies stemming from the finite sample approximation, "the test for no effect of harvesting time is misleading because the hierarchical structure of the data has not been appropriately accounted for."

Really? But the two error terms are in there!

In this balanced data situation, a another analysis is:

beets$bh <- with(beets, interaction(block, harvest))
summary(aov(sugpct ~ block + sow + harvest + Error(bh), data = beets))

Which they claim accounts for the hierarchical structure. The p-value in the second situation is much larger, which they attribute to this.

How can the likelihood ratio test be a truly general tool if it can't test hypotheses in the presence of multiple error terms? Aside from poor finite sample size properties, is the LRT shown above logically okay?

  • 1
    $\begingroup$ I think this is the df issue. If you run anova(sug, ddf='Kenward-Roger') after library('lmerTest'), then you will get exactly the same p-value as in RM-ANOVA. anova(sug, ddf='Satterthwaite') gives much lower p-value, close to the one from anova(sug, sug_no.harv). $\endgroup$
    – amoeba
    Mar 4, 2018 at 16:56
  • $\begingroup$ Just realized that the paper you referenced is exactly on Kenward-Roger approximation. I guess I should read it first :-) $\endgroup$
    – amoeba
    Mar 4, 2018 at 17:00
  • 1
    $\begingroup$ Note that KR procedure is not only changing the estimate of df, it also changes the statistic itself (at least as far as I understand). $\endgroup$
    – amoeba
    Mar 4, 2018 at 17:45
  • 1
    $\begingroup$ There are no differences btw Kenward-Roger and Satterthwaite in this case. KR assumes REML=TRUE and using this in the lmer() call makes anova(sug, ddf="Ken"), anova(sug, ddf="Sat") and the aov call all agree exactly. $\endgroup$ May 8, 2018 at 18:34
  • 1
    $\begingroup$ KR is in general adjusting the value of the F-statistic, but in a "simple" balanced case as this one the F-statistic is not adjusted (and it is the same as that used with the Satterthwaite method and the same as the one that comes out of equating observed and expected mean squares as shown in the paper) $\endgroup$ May 8, 2018 at 18:37

2 Answers 2


The way I think about this there are several 'layers' of asymptotic adjustments going on here.

First note that the likelihood ratio test and the Wald test are asymptotically equivalent (though the LRT can be seen as a small adjustment to the Wald test to adjust for non-quadratic behavior of the log-likelihood function.) The LR and equivalent Wald tests are:

> sug <- lmer(sugpct ~ block + sow + harvest + (1 | block:harvest),
+             data = beets)
> sug_no.harv <- update(sug, . ~ . - harvest)
> anova(sug, sug_no.harv, refit=TRUE) # LR
refitting model(s) with ML (instead of REML)
Data: beets
sug_no.harv: sugpct ~ block + sow + (1 | block:harvest)
sug: sugpct ~ block + sow + harvest + (1 | block:harvest)
            Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
sug_no.harv  9 -69.084 -56.473 43.542  -87.084                             
sug         10 -79.998 -65.986 49.999  -99.998 12.914      1  0.0003261 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Corresponding Wald Statistic:

> z <- coef(summary(sug))["harvestharv2", "t value"]
> z^2 # Wald chi-square statistic comparable to the LR statistic
[1] 15.21053
> pchisq(z^2, df=1, lower.tail = FALSE)
[1] 9.616582e-05
> 2 * pnorm(abs(z), lower.tail = FALSE) # equivalent
[1] 9.616582e-05

Now we can make the first finite-sample adjustment and use 20 df as we would in the corresponding fixed-effects model:

> 2 * pt(abs(z), df=20, lower.tail = FALSE)
[1] 0.0008886302

or we can take the full step and adjust for the hierarchical structure of the data and use 2 ddf:

> 2 * pt(abs(z), df=2, lower.tail = FALSE)
[1] 0.05989785

It is the latter adjustment from 20 to 2 - not the one from infinity to 20 that really matters and I think that is the one the pbkrtest authors were thinking of here. With ~200 whole plots there is no need for adjustments anymore and LR, asymptotic Wald and exact Wald (=F-test) all agree.

PS: Note that the REML LR test does not make sense because it depends on the coding of the fixed-effects in the model (and if, say, sum-to-zero contrasts are used a different test is obtained), so the LR test has to use ML -- not REML. On the other hand in the Wald test, the REML adjustment to the variance parameter estimation is important and non-ignorable given the size of this dataset: it ensures the correct variance of the harvest contrast in the example.

  • $\begingroup$ I like the simplicity of "just take the T-statistic and change the degrees of freedom," and it reminds me of textbook ANOVA cases where you created ratios of sums of squares based on what made sense. But I know pbkrtest is doing a lot more than that and I'm wondering if this is "cutting corners" in any way. It feels like it takes us from pure ML to a hybrid, which isn't necessarily bad but I still want to learn more. $\endgroup$
    – Ben Ogorek
    May 10, 2018 at 12:48

Will award best answer for theoretical justification, but a quick simulation makes me believe that the authors of pbrktest misspoke in Section 3.2. Indeed, a whole-plot treatment likelihood ratio test (LRT) is quite bad for small whole plot sample sizes, but performance improves as the number of whole plots increases. Below are whole plot p-value histograms from a Monte Carlo simulation where there is no whole plot treatment effect and we would expect a uniform distribution.

Five whole plots Two hundred whole plots

You try different whole plot and split plot sample sizes below. Of course, in an experimental design setting, the number of whole plots may indeed be small, making the methods in the pbkrtest package especially useful.


N <- 5 # whole plots. Use multiple of 2 for balanced data set
J <- 6  # subplots within w.p.s. Use multiple of 2 for balanced data set

B <- 5000 # Monte Carlo reps
p_values <- c()
for (rep in 1:B) {
  if (rep %% 100 == 0) cat("On Monte Carlo rep ", rep, "\n")
  u_i <- 5.5 * rnorm(N)
  e_i <- 2.9 * rnorm(N * J)

  sp_df <- expand.grid(subunit = 1:J, unit = 1:N)
  sp_df$wp_trt <- sp_df$unit <= N / 2
  sp_df$sp_trt <- sp_df$subunit <= J / 2

  sp_df$y <- 30 + u_i[sp_df$unit] + 1.5 * sp_df$sp_trt + e_i

  lmer_full_wp <- lmer(y ~ factor(wp_trt) + (1 | unit), data = sp_df,
                       REML = FALSE)
  lmer_red_wp <- lmer(y ~ (1 | unit), data = sp_df, REML = FALSE)

  my_anova <- anova(lmer_full_wp, lmer_red_wp)
  p_values <- c(p_values, my_anova["Pr(>Chisq)"][2, 1])

# With no whole plot treatment effect, the p-value distr should be uniform
hist(p_values, main = paste("Distribution of p-values for", N,
                            "whole plots", "and", J, "subplots / w.p."))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.