kr-modcomp.RdAn approximate F-test based on the Kenward-Roger approach.
KRmodcomp(largeModel, smallModel, betaH = 0, details = 0)
# S3 method for lmerMod
KRmodcomp(largeModel, smallModel, betaH = 0, details = 0)An lmer model
An lmer model or a restriction matrix
A number or a vector of the beta of the hypothesis, e.g. L
beta=L betaH. betaH=0 if smallModel is a model object and not a restriction matrix.
If larger than 0 some timing details are printed.
The model object must be fitted with restricted maximum
likelihood (i.e. with REML=TRUE). If the object is fitted with
maximum likelihood (i.e. with REML=FALSE) then the model is
refitted with REML=TRUE before the p-values are calculated. Put
differently, the user needs not worry about this issue.
An F test is calculated according to the approach of Kenward and Roger
(1997). The function works for linear mixed models fitted with the
lmer function of the lme4 package. Only models where the
covariance structure is a sum of known matrices can be compared.
The largeModel may be a model fitted with lmer either using
REML=TRUE or REML=FALSE. The smallModel can be a model
fitted with lmer. It must have the same covariance structure as
largeModel. Furthermore, its linear space of expectation must be a
subspace of the space for largeModel. The model smallModel
can also be a restriction matrix L specifying the hypothesis \(L
\beta = L \beta_H\), where \(L\) is a \(k \times p\) matrix and
\(\beta\) is a \(p\) column vector the same length as
fixef(largeModel).
The \(\beta_H\) is a \(p\) column vector.
Notice: if you want to test a hypothesis \(L \beta = c\) with a \(k\) vector \(c\), a suitable \(\beta_H\) is obtained via \(\beta_H=L c\) where \(L_n\) is a g-inverse of \(L\).
Notice: It cannot be guaranteed that the results agree with other implementations of the Kenward-Roger approach!
This functionality is not thoroughly tested and should be used with care. Please do report bugs etc.
Ulrich Halekoh, Søren Højsgaard (2014)., A Kenward-Roger Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed Models - The R Package pbkrtest., Journal of Statistical Software, 58(10), 1-30., https://www.jstatsoft.org/v59/i09/
Kenward, M. G. and Roger, J. H. (1997), Small Sample Inference for Fixed Effects from Restricted Maximum Likelihood, Biometrics 53: 983-997.
(fmLarge <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ Days + (Days | Subject)
#> Data: sleepstudy
#> REML criterion at convergence: 1743.628
#> Random effects:
#> Groups Name Std.Dev. Corr
#> Subject (Intercept) 24.741
#> Days 5.922 0.07
#> Residual 25.592
#> Number of obs: 180, groups: Subject, 18
#> Fixed Effects:
#> (Intercept) Days
#> 251.41 10.47
## removing Days
(fmSmall <- lmer(Reaction ~ 1 + (Days|Subject), sleepstudy))
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ 1 + (Days | Subject)
#> Data: sleepstudy
#> REML criterion at convergence: 1769.845
#> Random effects:
#> Groups Name Std.Dev. Corr
#> Subject (Intercept) 25.53
#> Days 11.93 -0.18
#> Residual 25.59
#> Number of obs: 180, groups: Subject, 18
#> Fixed Effects:
#> (Intercept)
#> 257.8
anova(fmLarge, fmSmall)
#> refitting model(s) with ML (instead of REML)
#> Data: sleepstudy
#> Models:
#> fmSmall: Reaction ~ 1 + (Days | Subject)
#> fmLarge: Reaction ~ Days + (Days | Subject)
#> npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
#> fmSmall 5 1785.5 1801.4 -887.74 1775.5
#> fmLarge 6 1763.9 1783.1 -875.97 1751.9 23.537 1 1.226e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
KRmodcomp(fmLarge, fmSmall)
#> large : Reaction ~ Days + (Days | Subject)
#> small : Reaction ~ 1 + (Days | Subject)
#> stat ndf ddf F.scaling p.value
#> Ftest 45.853 1.000 17.000 1 3.264e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## The same test using a restriction matrix
L <- cbind(0, 1)
KRmodcomp(fmLarge, L)
#> large : Reaction ~ Days + (Days | Subject)
#> L =
#> [,1] [,2]
#> [1,] 0 1
#> stat ndf ddf F.scaling p.value
#> Ftest 45.853 1.000 17.000 1 3.264e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## Same example, but with independent intercept and slope effects:
m.large <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), data = sleepstudy)
m.small <- lmer(Reaction ~ 1 + (1|Subject) + (0+Days|Subject), data = sleepstudy)
anova(m.large, m.small)
#> refitting model(s) with ML (instead of REML)
#> Data: sleepstudy
#> Models:
#> m.small: Reaction ~ 1 + (1 | Subject) + (0 + Days | Subject)
#> m.large: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
#> npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
#> m.small 4 1783.6 1796.4 -887.8 1775.6
#> m.large 5 1762.0 1778.0 -876.0 1752.0 23.6 1 1.186e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
KRmodcomp(m.large, m.small)
#> large : Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
#> small : Reaction ~ 1 + (1 | Subject) + (0 + Days | Subject)
#> stat ndf ddf F.scaling p.value
#> Ftest 45.046 1.000 18.188 1 2.571e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1