Print and Summary Method Utilities for Mixed Effects

The print, summary methods (including the print for the summary() result) in lme4 are modular, using about ten small utility functions. Other packages, building on lme4 can use the same utilities for ease of programming and consistency of output.

Notably see the Examples.

llikAIC() extracts the log likelihood, AIC, and related statics from a Fitted LMM.

formatVC() “format()”s the VarCorr matrix of the random effects – for print()ing and show()ing; it is also the “workhorse” of .prt.VC(), and returns a character matrix.

.prt.*() all use cat and print to produce output.

llikAIC(object, devianceFUN = devCrit, chkREML = TRUE,
        devcomp = object@devcomp)

methTitle(dims)

.prt.methTit(mtit, class)
.prt.family (famL)
.prt.resids (resids, digits, title = "Scaled residuals:", ...)
.prt.call   (call, long = TRUE)
.prt.aictab (aictab, digits = 1)
.prt.grps   (ngrps, nobs)
.prt.warn   (optinfo, summary = FALSE, ...)

.prt.VC (varcor, digits, comp = "Std.Dev.", corr = any(comp == "Std.Dev."),
         formatter = format, ...)
formatVC(varcor, digits = max(3, getOption("digits") - 2),
         comp = "Std.Dev.", corr = any(comp == "Std.Dev."),
         formatter = format,
         useScale = attr(varcor, "useSc"), ...)

Arguments

object: a LMM model fit
devianceFUN: the function to be used for computing the deviance; should not be changed for lme4 created objects.
chkREML: optional logical indicating if object maybe a REML fit.

devcomp: for lme4 always the equivalent of object@devcomp; here a list

dims: for lme4 always the equivalent of object@devcomp$dims, a named vector or list with components "GLMM", "NLMM", "REML", and "nAGQ" of which the first two are logical scalars, and the latter two typically are FALSE or numeric.

mtit: the result of methTitle(object)
class: typically class(object).

famL: a list with components family and link, each a character string; note that standard R family objects can be used directly, as well.

resids: numeric vector of model residuals.
digits: non-negative integer of (significant) digits to print minimally.
title: character string.
...: optional arguments passed on, e.g., to residuals().

call: the call of the model fit; e.g., available via (generic) function getCall().
long: logical indicating if the output may be long, e.g., printing the control part of the call if there is one.

aictab: typically the AICtab component of the result of llikAIC().

varcor: typically the result of VarCorr().
comp: optional character vector of length 1 or 2, containing "Std.Dev." and/or "Variance", indicating the columns to use.
corr: logical indicating if correlations or covariances should be used for vector random effects.
formatter: a function used for formatting the numbers.

ngrps: integer (vector), typically the result of ngrps(object).
nobs: integer; the number of observations, e.g., the result of nobs.

optinfo: typically object @ optinfo, the optimization infos, including warnings if there were.
summary: logical

useScale: (logical) whether the parent model estimates a scale parameter.

Value

llikAIC() returns a list with components

logLik: which is logLik(object), and
AICtab: a “table” of AIC, BIC, logLik, deviance and df.residual() values.

Examples

## Create a few "lme4 standard" models ------------------------------
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
fmM <- update(fm1, REML=FALSE) # -> Maximum Likelihood
fmQ <- update(fm1, . ~ Days + (Days | Subject))

gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
             data = cbpp, family = binomial)
gmA <- update(gm1, nAGQ = 5)


(lA1 <- llikAIC(fm1))
#> $logLik
#> 'log Lik.' -871.8141 (df=6)
#> 
#> $AICtab
#>     REML 
#> 1743.628 
#> 
(lAM <- llikAIC(fmM))
#> $logLik
#> 'log Lik.' -875.9697 (df=6)
#> 
#> $AICtab
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#> 1763.9393 1783.0971 -875.9697 1751.9393  174.0000 
#> 
(lAg <- llikAIC(gmA))
#> $logLik
#> 'log Lik.' -50.00568 (df=5)
#> 
#> $AICtab
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#> 110.01137 120.13813 -50.00568 100.01137  51.00000 
#> 

(m1 <- methTitle(fm1 @ devcomp $ dims))
#> [1] "Linear mixed model fit by REML"
(mM <- methTitle(fmM @ devcomp $ dims))
#> [1] "Linear mixed model fit by maximum likelihood "
(mG <- methTitle(gm1 @ devcomp $ dims))
#> [1] "Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)"
(mA <- methTitle(gmA @ devcomp $ dims))
#> [1] "Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 5)"

.prt.methTit(m1, class(fm1))
#> Linear mixed model fit by REML ['lmerMod']
.prt.methTit(mA, class(gmA))
#> Generalized linear mixed model fit by maximum likelihood (Adaptive
#>   Gauss-Hermite Quadrature, nAGQ = 5) [glmerMod]

.prt.family(gaussian())
#>  Family: gaussian  ( identity )
.prt.family(binomial())
#>  Family: binomial  ( logit )
.prt.family( poisson())
#>  Family: poisson  ( log )

.prt.resids(residuals(fm1), digits = 4)
#> Scaled residuals: 
#>      Min       1Q   Median       3Q      Max 
#> -101.179  -11.859    0.592   11.859  132.547 
#> 
.prt.resids(residuals(fmM), digits = 2)
#> Scaled residuals: 
#>    Min     1Q Median     3Q    Max 
#> -100.9  -11.9    0.7   11.9  132.5 
#> 

.prt.call(getCall(fm1))
#> Formula: Reaction ~ Days + (Days | Subject)
#>    Data: sleepstudy
.prt.call(getCall(gm1))
#> Formula: cbind(incidence, size - incidence) ~ period + (1 | herd)
#>    Data: cbpp

.prt.aictab ( lA1 $ AICtab ) # REML
#> REML criterion at convergence: 1743.6
.prt.aictab ( lAM $ AICtab ) # ML --> AIC, BIC, ...
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#>    1763.9    1783.1    -876.0    1751.9       174 

V1 <- VarCorr(fm1)
m <- formatVC(V1)
stopifnot(is.matrix(m), is.character(m), ncol(m) == 4)
print(m, quote = FALSE) ## prints all but the first line of .prt.VC() below:
#>  Groups   Name        Std.Dev. Corr 
#>  Subject  (Intercept) 24.7407       
#>           Days         5.9221  0.066
#>  Residual             25.5918       
.prt.VC( V1, digits = 4)
#> Random effects:
#>  Groups   Name        Std.Dev. Corr
#>  Subject  (Intercept) 24.741       
#>           Days         5.922   0.07
#>  Residual             25.592       
## Random effects:
##  Groups   Name        Std.Dev. Corr
##  Subject  (Intercept) 24.740
##           Days         5.922   0.07
##  Residual             25.592
p1 <- capture.output(V1)
p2 <- capture.output( print(m, quote=FALSE) )
pX <- capture.output( .prt.VC(V1, digits = max(3, getOption("digits")-2)) )
stopifnot(identical(p1, p2),
          identical(p1, pX[-1])) # [-1] : dropping 1st line

(Vq <- VarCorr(fmQ)) # default print()
#>  Groups   Name        Std.Dev. Corr 
#>  Subject  (Intercept) 24.7407       
#>           Days         5.9221  0.066
#>  Residual             25.5918       
print(Vq, comp = c("Std.Dev.", "Variance"))
#>  Groups   Name        Std.Dev. Variance Corr 
#>  Subject  (Intercept) 24.7407  612.100       
#>           Days         5.9221   35.072  0.066
#>  Residual             25.5918  654.940       
print(Vq, comp = c("Std.Dev.", "Variance"), corr=FALSE)
#>  Groups   Name        Std.Dev. Variance Cov    
#>  Subject  (Intercept) 24.7407  612.100         
#>           Days         5.9221   35.072    9.604
#>  Residual             25.5918  654.940         
print(Vq, comp = "Variance")
#>  Groups   Name        Variance Cov    
#>  Subject  (Intercept) 612.100         
#>           Days         35.072    9.604
#>  Residual             654.940         

.prt.grps(ngrps = ngrps(fm1),
          nobs  = nobs (fm1))
#> Number of obs: 180, groups:  Subject, 18
## --> Number of obs: 180, groups:  Subject, 18

.prt.warn(fm1 @ optinfo) # nothing .. had no warnings
.prt.warn(fmQ @ optinfo) # (ditto)