kr_modcomp.RdAn approximate F-test based on the Kenward-Roger approach.
KRmodcomp(largeModel, smallModel, betaH = 0, details = 0, ...)
# S3 method for class 'lmerMod'
KRmodcomp(largeModel, smallModel, betaH = 0, details = 0, ...)
# S3 method for class 'gls'
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. If smallModel is a model object then betaH=0.
If larger than 0 some timing details are printed.
Additional arguments, currently not used.
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 linear combination (a weighted sum)
of known matrices can be compared.
The smallModel is the model to be tested against the largeModel.
The largeModel is a model fitted with lmer(). A technical
detail: The model must be fitted with REML=TRUE. If the model is
fitted 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.
The smallModel can be one of several things:
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.
a restriction matrix L specifying the hypothesis
$$L \beta = L \beta_H$$
where L is a k x p matrix (there are k restrictions and p is
the number of fixed effect parameters (the length of
fixef(largeModel)) and beta_H is a p column vector.
A formula or a text string specifying what is to be removed from the larger model to form the smaller model.
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!
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.
(fm0 <- lmer(Reaction ~ (Days|Subject), sleepstudy))
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ (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
(fm1 <- 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
(fm2 <- lmer(Reaction ~ Days + I(Days^2) + (Days|Subject), sleepstudy))
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ Days + I(Days^2) + (Days | Subject)
#> Data: sleepstudy
#> REML criterion at convergence: 1742.816
#> Random effects:
#> Groups Name Std.Dev. Corr
#> Subject (Intercept) 24.761
#> Days 5.925 0.06
#> Residual 25.534
#> Number of obs: 180, groups: Subject, 18
#> Fixed Effects:
#> (Intercept) Days I(Days^2)
#> 255.449 7.434 0.337
## Test for no effect of Days in fm1, i.e. test fm0 under fm1
KRmodcomp(fm1, "Days")
#> Reaction ~ Days + (Days | Subject)
#> Reaction ~ (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 45.853 1.000 17.000 3.264e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
KRmodcomp(fm1, ~.-Days)
#> Reaction ~ Days + (Days | Subject)
#> Reaction ~ (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 45.853 1.000 17.000 3.264e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
L1 <- cbind(0, 1)
KRmodcomp(fm1, L1)
#> Reaction ~ Days + (Days | Subject)
#> Reaction ~ -1 + .X1 + (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 45.853 1.000 17.000 3.264e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
KRmodcomp(fm1, fm0)
#> Reaction ~ Days + (Days | Subject)
#> Reaction ~ (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 45.853 1.000 17.000 3.264e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(fm1, fm0)
#> refitting model(s) with ML (instead of REML)
#> Data: sleepstudy
#> Models:
#> fm0: Reaction ~ (Days | Subject)
#> fm1: Reaction ~ Days + (Days | Subject)
#> npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
#> fm0 5 1785.5 1801.4 -887.74 1775.5
#> fm1 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
## Test for no effect of Days and Days-squared in fm2, i.e. test fm0 under fm2
KRmodcomp(fm2, "(Days+I(Days^2))")
#> Reaction ~ Days + I(Days^2) + (Days | Subject)
#> Reaction ~ (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 23.364 2.000 39.817 1.938e-07 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
KRmodcomp(fm2, ~. - Days - I(Days^2))
#> Reaction ~ Days + I(Days^2) + (Days | Subject)
#> Reaction ~ (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 23.364 2.000 39.817 1.938e-07 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
L2 <- rbind(c(0, 1, 0), c(0, 0, 1))
KRmodcomp(fm2, L2)
#> Reaction ~ Days + I(Days^2) + (Days | Subject)
#> Reaction ~ -1 + .X1 + (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 23.364 2.000 39.817 1.938e-07 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
KRmodcomp(fm2, fm0)
#> Reaction ~ Days + I(Days^2) + (Days | Subject)
#> Reaction ~ (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 23.364 2.000 39.817 1.938e-07 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(fm2, fm0)
#> refitting model(s) with ML (instead of REML)
#> Data: sleepstudy
#> Models:
#> fm0: Reaction ~ (Days | Subject)
#> fm2: Reaction ~ Days + I(Days^2) + (Days | Subject)
#> npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
#> fm0 5 1785.5 1801.4 -887.74 1775.5
#> fm2 7 1764.3 1786.6 -875.14 1750.3 25.194 2 3.382e-06 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## Test for no effect of Days-squared in fm2, i.e. test fm1 under fm2
KRmodcomp(fm2, "I(Days^2)")
#> Reaction ~ Days + I(Days^2) + (Days | Subject)
#> Reaction ~ Days + (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 1.6558 1.0000 143.0000 0.2003
KRmodcomp(fm2, ~. - I(Days^2))
#> Reaction ~ Days + I(Days^2) + (Days | Subject)
#> Reaction ~ Days + (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 1.6558 1.0000 143.0000 0.2003
L3 <- rbind(c(0, 0, 1))
KRmodcomp(fm2, L3)
#> Reaction ~ Days + I(Days^2) + (Days | Subject)
#> Reaction ~ -1 + .X1 + .X2 + (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 1.6558 1.0000 143.0000 0.2003
KRmodcomp(fm2, fm1)
#> Reaction ~ Days + I(Days^2) + (Days | Subject)
#> Reaction ~ Days + (Days | Subject)
#> stat df ddf p.value
#> KR_Ftest 1.6558 1.0000 143.0000 0.2003
anova(fm2, fm1)
#> refitting model(s) with ML (instead of REML)
#> Data: sleepstudy
#> Models:
#> fm1: Reaction ~ Days + (Days | Subject)
#> fm2: Reaction ~ Days + I(Days^2) + (Days | Subject)
#> npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
#> fm1 6 1763.9 1783.1 -875.97 1751.9
#> fm2 7 1764.3 1786.6 -875.14 1750.3 1.6577 1 0.1979