prt-utilities {lme4}

Print and Summary Method Utilities for Mixed Effects

Description

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.

Usage

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))
(lAM <- llikAIC(fmM))
(lAg <- llikAIC(gmA))

(m1 <- methTitle(fm1 @ devcomp $ dims))
(mM <- methTitle(fmM @ devcomp $ dims))
(mG <- methTitle(gm1 @ devcomp $ dims))
(mA <- methTitle(gmA @ devcomp $ dims))

.prt.methTit(m1, class(fm1))
.prt.methTit(mA, class(gmA))

.prt.family(gaussian())
.prt.family(binomial())
.prt.family( poisson())

.prt.resids(residuals(fm1), digits = 4)
.prt.resids(residuals(fmM), digits = 2)

.prt.call(getCall(fm1))
.prt.call(getCall(gm1))

.prt.aictab ( lA1 $ AICtab ) # REML
.prt.aictab ( lAM $ AICtab ) # ML --> AIC, BIC, ...

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:
.prt.VC( V1, digits = 4)
## 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()
print(Vq, comp = c("Std.Dev.", "Variance"))
print(Vq, comp = c("Std.Dev.", "Variance"), corr=FALSE)
print(Vq, comp = "Variance")

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

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