F-test and degrees of freedom based on Satterthwaite approximation

An approximate F-test based on the Satterthwaite approach.

SATmodcomp(
  largeModel,
  smallModel,
  betaH = 0,
  details = 0,
  eps = sqrt(.Machine$double.eps)
)

# S3 method for class 'lmerMod'
SATmodcomp(
  largeModel,
  smallModel,
  betaH = 0,
  details = 0,
  eps = sqrt(.Machine$double.eps)
)

Arguments

largeModel

An lmer model

smallModel

An lmer model or a restriction matrix

betaH

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.

details

If larger than 0 some timing details are printed.

eps

A small number.

Details

Notice: It cannot be guaranteed that the results agree with other implementations of the Satterthwaite approach!

References

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/

See also

getKR, lmer, vcovAdj, PBmodcomp, KRmodcomp

Author

Søren Højsgaard, sorenh@math.aau.dk

Examples

(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
SATmodcomp(fm1, "Days")
#> large : Reaction ~ Days + (Days | Subject)
#> <environment: 0x55a80f476368>
#> small : Reaction ~ (Days | Subject)
#> <environment: 0x55a80f476368>
#>      statistic    ndf ddf   p.value    
#> [1,]    45.853  1.000  17 3.264e-06 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
SATmodcomp(fm1, ~.-Days)
#> large : Reaction ~ Days + (Days | Subject)
#> <environment: 0x55a80f476368>
#> small : Reaction ~ (Days | Subject)
#> <environment: 0x55a80f476368>
#>      statistic    ndf ddf   p.value    
#> [1,]    45.853  1.000  17 3.264e-06 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
L1 <- cbind(0, 1) 
SATmodcomp(fm1, L1)
#> Error in attr(x, "formula") %||% {    form <- switch(mode(x), `NULL` = structure(list(), class = "formula"),         character = eval(str2expression(x)), call = eval(x),         stop("invalid formula"))    environment(form) <- env    form}: invalid formula
SATmodcomp(fm1, fm0)
#> large : Reaction ~ Days + (Days | Subject)
#> <environment: 0x55a80f476368>
#> small : Reaction ~ (Days | Subject)
#> <environment: 0x55a80f476368>
#>      statistic    ndf ddf   p.value    
#> [1,]    45.853  1.000  17 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
SATmodcomp(fm2, "(Days+I(Days^2))")
#> large : Reaction ~ Days + I(Days^2) + (Days | Subject)
#> <environment: 0x55a80f476368>
#> small : Reaction ~ (Days | Subject)
#> <environment: 0x55a80f476368>
#>      statistic    ndf    ddf   p.value    
#> [1,]    23.754  2.000 51.563 4.876e-08 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
SATmodcomp(fm2, ~. - Days - I(Days^2))
#> large : Reaction ~ Days + I(Days^2) + (Days | Subject)
#> <environment: 0x55a80f476368>
#> small : Reaction ~ (Days | Subject)
#> <environment: 0x55a80f476368>
#>      statistic    ndf    ddf   p.value    
#> [1,]    23.754  2.000 51.563 4.876e-08 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
L2 <- rbind(c(0, 1, 0), c(0, 0, 1))
SATmodcomp(fm2, L2)
#> Error in attr(x, "formula") %||% {    form <- switch(mode(x), `NULL` = structure(list(), class = "formula"),         character = eval(str2expression(x)), call = eval(x),         stop("invalid formula"))    environment(form) <- env    form}: invalid formula
SATmodcomp(fm2, fm0)
#> large : Reaction ~ Days + I(Days^2) + (Days | Subject)
#> <environment: 0x55a80f476368>
#> small : Reaction ~ (Days | Subject)
#> <environment: 0x55a80f476368>
#>      statistic    ndf    ddf   p.value    
#> [1,]    23.754  2.000 51.563 4.876e-08 ***
#> ---
#> 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
SATmodcomp(fm2, "I(Days^2)")
#> large : Reaction ~ Days + I(Days^2) + (Days | Subject)
#> <environment: 0x55a80f476368>
#> small : Reaction ~ Days + (Days | Subject)
#> <environment: 0x55a80f476368>
#>      statistic    ndf ddf p.value
#> [1,]    1.6558 1.0000 143  0.2003
SATmodcomp(fm2, ~. - I(Days^2))
#> large : Reaction ~ Days + I(Days^2) + (Days | Subject)
#> <environment: 0x55a80f476368>
#> small : Reaction ~ Days + (Days | Subject)
#> <environment: 0x55a80f476368>
#>      statistic    ndf ddf p.value
#> [1,]    1.6558 1.0000 143  0.2003
L3 <- rbind(c(0, 0, 1))
SATmodcomp(fm2, L3)
#> Error in attr(x, "formula") %||% {    form <- switch(mode(x), `NULL` = structure(list(), class = "formula"),         character = eval(str2expression(x)), call = eval(x),         stop("invalid formula"))    environment(form) <- env    form}: invalid formula
SATmodcomp(fm2, fm1)
#> large : Reaction ~ Days + I(Days^2) + (Days | Subject)
#> <environment: 0x55a80f476368>
#> small : Reaction ~ Days + (Days | Subject)
#> <environment: 0x55a80f476368>
#>      statistic    ndf ddf p.value
#> [1,]    1.6558 1.0000 143  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