Package 'lme4'

Description

lme4 provides functions for fitting and analyzing mixed models: linear (lmer), generalized linear (glmer) and nonlinear (nlmer.)

Differences between nlme and lme4

lme4 covers approximately the same ground as the earlier nlme package. The most important differences are:

• lme4 uses modern, efficient linear algebra methods as implemented in the Eigen package, and uses reference classes to avoid undue copying of large objects; it is therefore likely to be faster and more memory-efficient than nlme.

• lme4 includes generalized linear mixed model (GLMM) capabilities, via the glmer function.

• lme4 does not currently implement nlme's features for modeling heteroscedasticity and correlation of residuals.

• lme4 does not currently offer the same flexibility as nlme for composing complex variance-covariance structures, but it does implement crossed random effects in a way that is both easier for the user and much faster.

• lme4 offers built-in facilities for likelihood profiling and parametric bootstrapping.

• lme4 is designed to be more modular than nlme, making it easier for downstream package developers and end-users to re-use its components for extensions of the basic mixed model framework. It also allows more flexibility for specifying different functions for optimizing over the random-effects variance-covariance parameters.

• lme4 is not (yet) as well-documented as nlme.

Differences between current (1.0.+) and previous versions of lme4

• [gn]lmer now produces objects of class merMod rather than class mer as before

• the new version uses a combination of S3 and reference classes (see ReferenceClasses, merPredD-class, and lmResp-class) as well as S4 classes; partly for this reason it is more interoperable with nlme

• The internal structure of [gn]lmer is now more modular, allowing finer control of the different steps of argument checking; construction of design matrices and data structures; parameter estimation; and construction of the final merMod object (see modular)

• profiling and parametric bootstrapping are new in the current version

• the new version of lme4 does not provide an mcmcsamp (post-hoc MCMC sampling) method, because this was deemed to be unreliable. Alternatives for computing p-values include parametric bootstrapping (bootMer) or methods implemented in the pbkrtest package and leveraged by the lmerTest package and the Anova function in the car package (see pvalues for more details).

Caveats and trouble-shooting

• Some users who have previously installed versions of the RcppEigen and minqa packages may encounter segmentation faults (!!); the solution is to make sure to re-install these packages before installing lme4. (Because the problem is not with the explicit version of the packages, but with running packages that were built with different versions of Rcpp in conjunction with each other, simply making sure you have the latest version, or using update.packages, will not necessarily solve the problem; you must actually re-install the packages. The problem is most likely with minqa.)

Description

Attempt to re-fit a [g]lmer model with a range of optimizers. The default is to use all known optimizers for R that satisfy the requirements (i.e. they do not require functions and allow box constraints: see ‘optimizer’ in lmerControl). These optimizers fall in four categories; (i) built-in (minqa::bobyqa, lme4::Nelder_Mead, nlminbwrap), (ii) wrapped via optimx (most of optimx's optimizers that allow box constraints require an explicit gradient function to be specified; the two provided here are the base R functions that can be accessed via optimx), (iii) wrapped via nloptr (see examples for the list of options), (iv) ‘dfoptim::nmkb’ (via the (unexported) nmkbw wrapper: this appears as ‘nmkbw’ in meth.tab)

Usage

allFit(object, meth.tab = NULL, data = NULL,
       verbose = TRUE,
       show.meth.tab = FALSE,
       maxfun = 1e5,
       parallel = c("no", "multicore", "snow", "future"),
       ncpus = getOption("allFit.ncpus", 1L), cl = NULL,
       catch.errs = TRUE,
       start_from_mle = TRUE)

Arguments

Details

• Needs packages optimx, and dfoptim to use all optimizers.

• If you are using parallel="snow" (e.g. when running in parallel on Windows), you will need to set up a cluster yourself and run clusterEvalQ(cl,library("lme4")) before calling allFit to make sure that the lme4 package is loaded on all of the workers.

• Control arguments in control$optCtrl that are unused by a particular optimizer will be silently ignored (in particular, the maxfun specification is only respected by bobyqa, Nelder_Mead, and nmkbw).

• Because allFit works by calling update, it may be fragile if the original model call contains references to variables, especially if they were originally defined in other environments or no longer exist when allFit is called.

Value

an object of type allFit, which is a list of fitted merMod objects (unless show.meth.tab is specified, in which case a data frame of methods is returned). The summary method for this class extracts tables with a variety of useful information about the different fits (see examples).

See Also

slice,slice2D from the bbmle package

Examples

if (interactive()) {
library(lme4)
  gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
             data = cbpp, family = binomial)
  ## show available methods
  allFit(show.meth.tab=TRUE) 
  gm_all <- allFit(gm1)
  ss <- summary(gm_all)
  ss$which.OK            ## logical vector: which optimizers worked?
  ## the other components only contain values for the optimizers that worked
  ss$llik                ## vector of log-likelihoods
  ss$fixef               ## table of fixed effects
  ss$sdcor               ## table of random effect SDs and correlations
  ss$theta               ## table of random effects parameters, Cholesky scale
} 
## Not run: 
  ## Parallel examples for Windows
  nc <- max(detectCores() - 1L, 1L, na.rm = TRUE)
  optCls <- makeCluster(nc)
  clusterEvalQ(optCls,library("lme4"))
  ### not necessary here because using a built-in
  ## data set, but in general you should clusterExport() your data
  clusterExport(optCls, "cbpp")
  system.time(af1 <- allFit(m0, parallel = 'snow', 
                          ncpus = nc, cl=optCls))
  stopCluster(optCls)

## End(Not run)

Description

Data on genetic variation in responses to fertilization and simulated herbivory in Arabidopsis

Usage

data("Arabidopsis")

Format

A data frame with 625 observations on the following 8 variables.

reg

region: a factor with 3 levels NL (Netherlands), SP (Spain), SW (Sweden)

popu

population: a factor with the form n.R representing a population in region R

gen

genotype: a factor with 24 (numeric-valued) levels

rack

a nuisance factor with 2 levels, one for each of two greenhouse racks

nutrient

fertilization treatment/nutrient level (1, minimal nutrients or 8, added nutrients)

amd

simulated herbivory or "clipping" (apical meristem damage): unclipped (baseline) or clipped

status

a nuisance factor for germination method (Normal, Petri.Plate, or Transplant)

total.fruits

total fruit set per plant (integer)

Source

(Banta et al. 2010)

References

Banta JA, Stevens MHH, Pigliucci M (2010). “A comprehensive test of the 'limiting resources' framework applied to plant tolerance to apical meristem damage.” Oikos, 119(2), 359-369. doi:10.1111/j.1600-0706.2009.17726.x.

Examples

data(Arabidopsis)
summary(Arabidopsis[,"total.fruits"])
table(gsub("[0-9].","",levels(Arabidopsis[,"popu"])))
library(lattice)
stripplot(log(total.fruits+1) ~ amd|nutrient, data = Arabidopsis,
          groups = gen,
          strip=strip.custom(strip.names=c(TRUE,TRUE)),
          type=c('p','a'), ## points and panel-average value --
          ## see ?panel.xyplot
          scales=list(x=list(rot=90)),
          main="Panel: nutrient, Color: genotype")

Description

Perform model-based (Semi-)parametric bootstrap for mixed models.

Usage

bootMer(x, FUN, nsim = 1, seed = NULL, use.u = FALSE, re.form=NA,
	type = c("parametric", "semiparametric"),
	verbose = FALSE, .progress = "none", PBargs = list(),
	parallel = c("no", "multicore", "snow", "future"),
	ncpus = getOption("boot.ncpus", 1L), cl = NULL)

Arguments

Details

The semi-parametric variant is only partially implemented, and we only provide a method for lmer and glmer results.

Information about warning and error messages incurred during the bootstrap returns can be retrieved via the attributes

bootFail

number of failures (errors)

boot.fail.msgs

error messages

boot.all.msgs

messages, warnings, and error messages

e.g. attr("boot.fail.msgs") to retrieve error messages

The working name for bootMer() was “simulestimate()”, as it is an extension of simulate (see simulate.merMod), but we want to emphasize its potential for valid inference.

• If use.u is FALSE and type is "parametric", each simulation generates new values of both the “spherical” random effects $u$ and the i.i.d. errors $\epsilon$, using rnorm() with parameters corresponding to the fitted model x.

• If use.u is TRUE and type=="parametric", only the i.i.d. errors (or, for GLMMs, response values drawn from the appropriate distributions) are resampled, with the values of $u$ staying fixed at their estimated values.

• If use.u is TRUE and type=="semiparametric", the i.i.d. errors are sampled from the distribution of (response) residuals. (For GLMMs, the resulting sample will no longer have the same properties as the original sample, and the method may not make sense; a warning is generated.) The semiparametric bootstrap is currently an experimental feature, and therefore may not be stable.

• The case where use.u is FALSE and type=="semiparametric" is not implemented; Morris (2002) suggests that resampling from the estimated values of $u$ is not good practice.

Value

an object of S3 class "boot", compatible with boot package's boot() result. (See Details for information on how to retrieve information about errors during bootstrapping.)

Note

If you are using parallel="snow", you will need to run clusterEvalQ(cl,library("lme4")) before calling bootMer to make sure that the lme4 package is loaded on all of the workers; you may additionally need to use clusterExport if you are using a summary function that calls any objects from the environment.

References

• Davison AC, Hinkley DV (1997). Bootstrap methods and their application. Cambridge University Press.

• Morris JS (2002). “The BLUPs are not "best" when it comes to bootstrapping.” Statistics & Probability Letters, 56(4), 425-430. ISSN 0167-7152. doi:10.1016/S0167-7152(02)00041-X.

See Also

• confint.merMod, for a more specific approach to bootstrap confidence intervals on parameters.

• refit(), or PBmodcomp() from the pbkrtest package, for parametric bootstrap comparison of models.

• boot(), and then boot.ci, from the boot package.

• profile-methods, for likelihood-based inference, including confidence intervals.

• pvalues, for more general approaches to inference and p-value computation in mixed models.

Examples

if (interactive()) {
fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)
## see ?"profile-methods"
mySumm <- function(.) { s <- sigma(.)
    c(beta =getME(., "beta"), sigma = s, sig01 = unname(s * getME(., "theta"))) }
(t0 <- mySumm(fm01ML)) # just three parameters
## alternatively:
mySumm2 <- function(.) {
    c(beta=fixef(.),sigma=sigma(.), sig01=sqrt(unlist(VarCorr(.))))
}

set.seed(101)
## 3.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
system.time( boo01 <- bootMer(fm01ML, mySumm, nsim = 100) )

## to "look" at it
if (requireNamespace("boot")) {
    boo01
    ## note large estimated bias for sig01
    ## (~30% low, decreases _slightly_ for nsim = 1000)

    ## extract the bootstrapped values as a data frame ...
    head(as.data.frame(boo01))

    ## ------ Bootstrap-based confidence intervals ------------

    ## warnings about "Some ... intervals may be unstable" go away
    ##   for larger bootstrap samples, e.g. nsim=500

    ## intercept
    (bCI.1 <- boot::boot.ci(boo01, index=1, type=c("norm", "basic", "perc")))# beta

    ## Residual standard deviation - original scale:
    (bCI.2  <- boot::boot.ci(boo01, index=2, type=c("norm", "basic", "perc")))
    ## Residual SD - transform to log scale:
    (bCI.2L <- boot::boot.ci(boo01, index=2, type=c("norm", "basic", "perc"),
                       h = log, hdot = function(.) 1/., hinv = exp))

    ## Among-batch variance:
    (bCI.3 <- boot::boot.ci(boo01, index=3, type=c("norm", "basic", "perc"))) # sig01

    
    confint(boo01)
    confint(boo01,type="norm")
    confint(boo01,type="basic")

    ## Graphical examination:
    plot(boo01,index=3)

    ## Check stored values from a longer (1000-replicate) run:
    (load(system.file("testdata","boo01L.RData", package="lme4")))# "boo01L"
    plot(boo01L, index=3)
    mean(boo01L$t[,"sig01"]==0) ## note point mass at zero!
} 
}

Description

Data on the breakage angle of chocolate cakes made with three different recipes and baked at six different temperatures. This is a split-plot design with the recipes being whole-units and the different temperatures being applied to sub-units (within replicates). The experimental notes suggest that the replicate numbering represents temporal ordering.

Format

A data frame with 270 observations on the following 5 variables.

replicate

a factor with levels 1 to 15

recipe

a factor with levels A, B and C

temperature

an ordered factor with levels 175 < 185 < 195 < 205 < 215 < 225

angle

a numeric vector giving the angle at which the cake broke.

temp

numeric value of the baking temperature (degrees F).

Details

The replicate factor is nested within the recipe factor, and temperature is nested within replicate.

Source

Original data were presented in Cook (1938), and reported in Cochran and Cox (1957, p. 300). Also cited in Lee, Nelder and Pawitan (2006).

References

• Cook FE (1938). Chocolate cake, I. Optimum baking temperature. Master's Thesis, Iowa State College.

• Cochran WG, Cox GM (1957). Experimental Designs, 2nd edition. John Wiley & Sons, New York.

• Lee Y, Nelder JA, Pawitan Y (2006). Generalized linear models with random effects: unified analysis via H-likelihood. Chapman and Hall/CRC, Boca Raton.

Examples

str(cake)
## 'temp' is continuous, 'temperature' an ordered factor with 6 levels

(fm1 <- lmer(angle ~ recipe * temperature + (1|recipe:replicate), cake, REML= FALSE))
(fm2 <- lmer(angle ~ recipe + temperature + (1|recipe:replicate), cake, REML= FALSE))
(fm3 <- lmer(angle ~ recipe + temp        + (1|recipe:replicate), cake, REML= FALSE))

## and now "choose" :
anova(fm3, fm2, fm1)

Description

Contagious bovine pleuropneumonia (CBPP) is a major disease of cattle in Africa, caused by a mycoplasma. This dataset describes the serological incidence of CBPP in zebu cattle during a follow-up survey implemented in 15 commercial herds located in the Boji district of Ethiopia. The goal of the survey was to study the within-herd spread of CBPP in newly infected herds. Blood samples were quarterly collected from all animals of these herds to determine their CBPP status. These data were used to compute the serological incidence of CBPP (new cases occurring during a given time period). Some data are missing (lost to follow-up).

There are two datasets, although their precise origins are unclear. The dataset used in the classical lme4 examples, cbpp, matches the file maintained by the package authors. The extended dataset, cbpp2, corresponds to Table 1 of (Lesnoff et al. 2004), with the exception of herd 6, which is a known typographical error in the original paper.

Format

cbpp is a data frame with 56 observations on the following 4 variables.

herd

A factor identifying the herd (1 to 15).

incidence

The number of new serological cases for a given herd and time period.

size

A numeric vector describing herd size at the beginning of a given time period.

period

A factor with levels 1 to 4.

The extended version, cbpp2, has the additional variables:

herd

A factor identifying the herd (1 to 15).

treatment

A factor referring to the control measure used to manage CBPP.

• Complete = complete isolation or antibiotic treatment,

• Partial/null = partial/null isolation and no antibiotic treatment,

• Unknown = strategy remained.

avg_size

The average number of animals housed in a della (a temporary paddock used for holding cattle on the farm).

Details

Serological status was determined using a competitive enzyme-linked immuno-sorbent assay (cELISA).

Source

(Lesnoff et al. 2004)

References

Lesnoff M, Laval G, Bonnet P, Abdicho S, Workalemahu A, Kifle D, Peyraud A, Lancelot R, Thiaucourt F (2004). “Within-herd spread of contagious bovine pleuropneumonia in Ethiopian highlands.” Preventive Veterinary Medicine, 64(1), 27-40. ISSN 0167-5877. doi:10.1016/j.prevetmed.2004.03.005.

Examples

## response as a matrix
(m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
             family = binomial, data = cbpp))
## response as a vector of probabilities and usage of argument "weights"
m1p <- glmer(incidence / size ~ period + (1 | herd), weights = size,
             family = binomial, data = cbpp)
## Confirm that these are equivalent:
stopifnot(all.equal(fixef(m1), fixef(m1p), tolerance = 1e-5),
          all.equal(ranef(m1), ranef(m1p), tolerance = 1e-5))

## GLMM with individual-level variability (accounting for overdispersion)
cbpp$obs <- 1:nrow(cbpp)
(m2 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd) +  (1|obs),
              family = binomial, data = cbpp))
              
## Fitting the model for cbpp2
gm1 <- glmer(incidence/size ~ period + treatment + avg_size + (1 | herd),
             family = binomial,
             data = cbpp2, weights = size,
             control = glmerControl(optimizer="bobyqa"))
## Adding an observation-level random effect
cbpp2 <- transform(cbpp2,obs=factor(seq(nrow(cbpp2))))
## Herd and observation-level REs (below causes singular fit issues)
gm2 <- update(gm1,.~.+(1|obs)) 
## observation-level REs only (no singular fit issue)
gm3 <- update(gm1,.~.-(1|herd)+(1|obs))

Description

Primarily internal code for checking optimization convergence, see convergence for a more detailed discussion.

Usage

checkConv(derivs, coefs, ctrl, lbound, ubound,
          debug = FALSE, nobs = NULL, ndim = NULL)

Arguments

Value

A result list containing

See Also

convergence

Description

Compute confidence intervals on the parameters of a *lmer() model fit (of class"merMod").

Usage

## S3 method for class 'merMod'
confint(object, parm, level = 0.95,
        method = c("profile", "Wald", "boot"), zeta,
        nsim = 500,
        boot.type = c("perc", "basic", "norm"),
        FUN = NULL, quiet = FALSE,
        oldNames, signames = TRUE,
        boot.scale = c("sdcor", "vcov"),
        ...)
## S3 method for class 'thpr'
confint(object, parm, level = 0.95,
        zeta, non.mono.tol = 1e-2, ...)

Arguments

Details

Depending on the method specified, confint() computes confidence intervals by

"profile":

computing a likelihood profile and finding the appropriate cutoffs based on the likelihood ratio test;

"Wald":

approximating the confidence intervals (of fixed-effect parameters only; all variance-covariance parameters CIs will be returned as NA) based on the estimated local curvature of the likelihood surface;

"boot":

performing parametric bootstrapping with confidence intervals computed from the bootstrap distribution according to boot.type (see bootMer, boot.ci).

Value

a numeric table (matrix with column and row names) of confidence intervals; the confidence intervals are computed on the standard deviation scale.

Note

The default method "profile" amounts to

    confint(profile(object, which=parm, ...),
            level, zeta)

where the profile method profile.merMod does almost all the computations. Therefore it is typically advisable to store the profile(.) result, say in pp, and then use confint(pp, level=*) e.g., for different levels.

Examples

if (interactive() || lme4_testlevel() >= 3) {
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
fm1W <- confint(fm1, method="Wald")# very fast, but not useful for "sigmas" = var-cov pars
fm1W
(fm2 <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy))
(CI2 <- confint(fm2, maxpts = 8)) # method = "profile"; 8: to be much faster

if (lme4_testlevel() >= 3) {
  system.time(fm1P <- confint(fm1, method="profile", ## <- default
                              signames = FALSE))
  ## --> ~ 2.2 seconds (2022)
  set.seed(123) # (reproducibility when using bootstrap)
  system.time(fm1B <- confint(fm1, method="boot", signames=FALSE,
                              .progress="txt", PBargs= list(style=3)))
  ## --> ~ 6.2 seconds (2022) and warning, messages
} else {
    load(system.file("testdata","confint_ex.rda",package="lme4"))
}
fm1P
fm1B
} ## if interactive && testlevel>=3

Description

[g]lmer fits may produce convergence warnings; these do not necessarily mean the fit is incorrect (see “Theoretical details” below). The following steps are recommended assessing and resolving convergence warnings (also see examples below):

• double-check the model specification and the data

• adjust stopping (convergence) tolerances for the nonlinear optimizer, using the optCtrl argument to [g]lmerControl (see “Convergence controls” below)

• center and scale continuous predictor variables (e.g. with scale)

• double-check the Hessian calculation with the more expensive Richardson extrapolation method (see examples)

• restart the fit from the reported optimum, or from a point perturbed slightly away from the reported optimum

• use allFit to try the fit with all available optimizers (e.g. several different implementations of BOBYQA and Nelder-Mead, L-BFGS-B from optim, nlminb, ...). While this will of course be slow for large fits, we consider it the gold standard; if all optimizers converge to values that are practically equivalent, then we would consider the convergence warnings to be false positives.

Details

Convergence controls

• the controls for the nloptwrap optimizer (the default for lmer) are

  ftol_abs

  (default 1e-6) stop on small change in deviance

  ftol_rel

  (default 0) stop on small relative change in deviance

  xtol_abs

  (default 1e-6) stop on small change of parameter values

  xtol_rel

  (default 0) stop on small relative change of parameter values

  maxeval

  (default 1000) maximum number of function evaluations

  Changing ftol_abs and xtol_abs to stricter values (e.g. 1e-8) is a good first step for resolving convergence problems, at the cost of slowing down model fits.

• the controls for minqa::bobyqa (default for glmer first-stage optimization) are

  rhobeg

  (default 2e-3) initial radius of the trust region

  rhoend

  (default 2e-7) final radius of the trust region

  maxfun

  (default 10000) maximum number of function evaluations

  rhoend, which describes the scale of parameter uncertainty on convergence, is approximately analogous to xtol_abs.

• the controls for Nelder_Mead (default for glmer second-stage optimization) are

  FtolAbs

  (default 1e-5) stop on small change in deviance

  FtolRel

  (default 1e-15) stop on small relative change in deviance

  XtolRel

  (default 1e-7) stop on small change of parameter values

  maxfun

  (default 10000) maximum number of function evaluations

Theoretical issues

lme4 uses general-purpose nonlinear optimizers (e.g. Nelder-Mead or Powell's BOBYQA method) to estimate the variance-covariance matrices of the random effects. Assessing the convergence of such algorithms reliably is difficult. For example, evaluating the Karush-Kuhn-Tucker conditions (convergence criteria which reduce in simple cases to showing that the gradient is zero and the Hessian is positive definite) is challenging because of the difficulty of evaluating the gradient and Hessian.

We (the lme4 authors and maintainers) are still in the process of finding the best strategies for testing convergence. Some of the relevant issues are

• the gradient and Hessian are the basic ingredients of KKT-style testing, but (at least for now) lme4 estimates them by finite-difference approximations which are sometimes unreliable.

• The Hessian computation in particular represents a difficult tradeoff between computational expense and accuracy. At present the Hessian computations used for convergence checking (and for estimating standard errors of fixed-effect parameters for GLMMs) follow the ordinal package in using a naive but computationally cheap centered finite difference computation (with a fixed step size of $10^{-4}$). A more reliable but more expensive approach is to use Richardson extrapolation, as implemented in the numDeriv package.

• it is important to scale the estimated gradient at the estimate appropriately; two reasonable approaches are

  1. scale gradients by the inverse Cholesky factor of the Hessian, equivalent to scaling gradients by the estimated Wald standard error of the estimated parameters. lme4 uses this approach; it requires the Hessian to be estimated (although the Hessian is required for reliable estimation of the fixed-effect standard errors for GLMMs in any case).

  2. use unscaled gradients on the random-effects parameters, since these are essentially already unitless (for LMMs they are scaled relative to the residual variance; for GLMMs they are scaled relative to the sampling variance of the conditional distribution); for GLMMs, scale fixed-effect gradients by the standard deviations of the corresponding input variable

• Exploratory analyses suggest that (1) the naive estimation of the Hessian may fail for large data sets (number of observations greater than approximately $10^{5}$); (2) the magnitude of the scaled gradient increases with sample size, so that warnings will occur even for apparently well-behaved fits with large data sets.

See Also

lmerControl, isSingular

Examples

if (interactive()) {
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

## 1. decrease stopping tolerances
strict_tol <- lmerControl(optCtrl=list(xtol_abs=1e-8, ftol_abs=1e-8))
if (all(fm1@optinfo$optimizer=="nloptwrap")) {
    fm1.tol <- update(fm1, control=strict_tol)
}

## 2. center and scale predictors:
ss.CS <- transform(sleepstudy, Days=scale(Days))
fm1.CS <- update(fm1, data=ss.CS)

## 3. recompute gradient and Hessian with Richardson extrapolation
devfun <- update(fm1, devFunOnly=TRUE)
if (isLMM(fm1)) {
    pars <- getME(fm1,"theta")
} else {
    ## GLMM: requires both random and fixed parameters
    pars <- getME(fm1, c("theta","fixef"))
}
if (require("numDeriv")) {
    cat("hess:\n"); print(hess <- hessian(devfun, unlist(pars)))
    cat("grad:\n"); print(grad <- grad(devfun, unlist(pars)))
    cat("scaled gradient:\n")
    print(scgrad <- solve(chol(hess), grad))
}
## compare with internal calculations:
fm1@optinfo$derivs

## compute reciprocal condition number of Hessian
H <- fm1@optinfo$derivs$Hessian
Matrix::rcond(H)

## 4. restart the fit from the original value (or
## a slightly perturbed value):
fm1.restart <- update(fm1, start=pars)
set.seed(101)
pars_x <- runif(length(pars),pars/1.01,pars*1.01)
fm1.restart2 <- update(fm1, start=pars_x,
                       control=strict_tol)

## 5. try all available optimizers

  fm1.all <- allFit(fm1)
  ss <- summary(fm1.all)
  ss$ fixef               ## fixed effects
  ss$ llik                ## log-likelihoods
  ss$ sdcor               ## SDs and correlations
  ss$ theta               ## Cholesky factors
  ss$ which.OK            ## which fits worked

}

Description

This data shows the results from a block design experiment to record the interactions of ‘guard’ crustaceans, shrimp (Alpheus lottini) and crab (Trapezia serenei), that defend their coral host (Pocillopora cf. meandrina) from seastar predators (Culcita novaeguineae) (McKeon et al. 2012). The coral and their seastar predator were placed in large, octagonal, flow-through seawater tanks measuring approximately 0.5m deep and 2m in diameter. There were 10 tanks in total, each serving as an experimental block. Four treatment groups were tested: (1) no exosymbionts, (2) shrimp only, (3) crab only, and (4) both crab and shrimp present. Each treatment was replicated twice within each block.

Usage

data("culcitalogreg")

Format

A data frame with 80 observations on the following 9 variables.

block

a numeric variable indicating the experimental block. There are 10 blocks in total, each corresponding to a large, octagonal, flow-through seawater tank approximately 0.5m deep and 2m in diameter.

ttt

a number which represents the combinations of different symbionts.

• 1: no exosymbionts,

• 2: pair of Alpheus lottini only (Alpheus; shrimp),

• 3: pair of Trapezia serenei only (Trapezia; crab),

• 4: pair of Alpheus lottini and pair of Trapezia serenei (‘Alpheus and Trapezia’).

predation

(binary) indicates whether the seastar predator consumed the coral. 0 = did not eat; 1 = ate.

ttt.1

a relabelled version of ttt that explicitly indicates the corresponding treatment condition in plain language.

crab

describes whether the crab was present in the experiment. C-: not present, C+: present.

shrimp

describes whether the shrimp was present in the experiment. S-: not present, S+: present.

ttt2

a relabelled version of ttt that uses letters (a, b, c, d) instead of numeric values (1, 2, 3, 4), respectively, to represent the treatment groups.

crab2

a relabelling of crab where n corresponds to C- (no crab), and y corresponds to C+ (crab present).

shrimp2

a relabelling of shrimp where n corresponds to S- (no shrimp), and y corresponds to S+ (shrimp present).

Details

Full details of the experiment (McKeon et al. 2012), is quoted below: "Twenty replicates in total were conducted for each exosymbiont treatment in a temporally blocked design with two replicates in each of ten temporal blocks. Trials were conducted in a large octagonal flow-through seawater tank approximately 0.5 m deep and 2 m across. The tank was divided into eight equal sections using plastic screening. Each section was provisioned with a seastar refugium constructed from concrete blocks. We placed Pocillopora colonies into the tank from the field during mid- to late afternoon. To minimize variation driven by search time, a single Culcita was placed directly on top of the coral colony. The following morning, we measured the coral size (length, width, height) and the feeding scars left by the Culcita (length, width, depth)."

Source

(McKeon et al. 2012)

References

McKeon CS, Stier AC, McIlroy SE, Bolker BM (2012). “Multiple defender effects: synergistic coral defense by mutualist crustaceans.” Oecologia, 169(4), 1095–1103. doi:10.1007/s00442-012-2275-2.

See Also

The version which only contains predation counts but includes the amount of volume of the coral and the amount of volume lost, culcitalvolume.

Examples

cul_mod <- glmer(predation ~ ttt2 + (1|block), data=culcitalogreg,
                     family = binomial(link = "logit"))

Description

This is the same experiment outlined in culcitalogreg except this data set only contains values where predation occurred, and the original volume as well as volume lost was recorded.

Usage

data("culcitalvolume")

Format

A data frame with 50 observations on the following 7 variables.

ttt

represents the combinations of different symbionts, which the treatment condition is explicitly written in plain language.

volume

describes the size of the coral, measured in cm$^3$.

predvolume

the amount of volume loss from the coral, measured in cm$^3$.

crab

describes whether the crab was present in the experiment. n: not present, y: present.

shrimp

describes whether the shrimp was present in the experiment. n: not present, y: present.

block

a numeric variable indicating the experimental block. There are 10 blocks in total, each corresponding to a large, octagonal, flow-through seawater tank approximately 0.5m deep and 2m in diameter.

ttt2

a relabelled version of ttt.

• 1: no exosymbionts,

• 2: pair of Alpheus lottini only (Alpheus; shrimp),

• 3: pair of Trapezia serenei only (Trapezia; crab),

• 4: pair of Alpheus lottini and pair of Trapezia serenei (‘Alpheus and Trapezia’).

Source

(McKeon et al. 2012)

References

McKeon CS, Stier AC, McIlroy SE, Bolker BM (2012). “Multiple defender effects: synergistic coral defense by mutualist crustaceans.” Oecologia, 169(4), 1095–1103. doi:10.1007/s00442-012-2275-2.

See Also

The version which contains whether or not predation occurred, culcitalogreg.

Examples

## Modifying to create a new response variable
vdata <- transform(culcitalvolume, 
                   propeaten = predvolume/volume,
                   tvol = log(predvolume))
## One-way analysis
(cvm1 <- lmer(tvol ~ ttt2 + (1|block), data = vdata))
(cvm2 <- lmer(propeaten ~ ttt2 + (1|block), data = vdata))
## Two-way analysis 
(cvm3 <- lmer(tvol ~ crab*shrimp + (1|block), data = vdata))
(cvm4 <- lmer(propeaten ~ crab*shrimp + (1|block), data = vdata))

Description

Return the deviance component list

Usage

devcomp(x)

Arguments

Details

A fitted model of class merMod has a devcomp slot as described in the value section.

Value

a list with components

Note

This function is deprecated, use getME(., "devcomp")

Description

The deviance is profiled with respect to the fixed-effects parameters but not with respect to sigma; that is, the function takes parameters for the variance-covariance parameters and for the residual standard deviation. The random-effects variance-covariance parameters are on the standard deviation/correlation scale, not the theta (Cholesky factor) scale.

Usage

devfun2(fm, useSc = if(isLMM(fm)) TRUE else NA,
        scale = c("sdcor", "varcov"), ...)

Arguments

Value

Returns a function that takes a vector of standard deviations and correlations and returns the deviance (or REML criterion). The function has additional attributes

optimum

a named vector giving the parameter values at the optimum

basedev

the deviance at the optimum, (i.e., not the REML criterion).

thopt

the optimal variance-covariance parameters on the “theta” (Cholesky factor) scale

stderr

standard errors of fixed effect parameters

Note

Even if the original model was fitted using REML=TRUE as by default with lmer(), this returns the deviance, i.e., the objective function for maximum (log) likelihood (ML).

For the REML objective function, use getME(fm, "devfun") instead.

Examples

m1 <- lmer(Reaction~Days+(Days|Subject),sleepstudy)
dd <- devfun2(m1, useSc=TRUE)
pp <- attr(dd, "optimum")
## extract variance-covariance and residual std dev parameters
sigpars <- pp[grepl("^\\.sig", names(pp))]
all.equal(unname(dd(sigpars)),deviance(refitML(m1)))

Description

Calculates the deviance of a model of class merMod. Some important details regarding the calculation of the deviance is found below.

Usage

## S3 method for class 'merMod'
deviance(object, REML = NULL, ...)

Arguments

Deviance and log-likelihood of GLMMs

One must be careful when defining the deviance of a GLM. For example, should the deviance be defined as minus twice the log-likelihood or does it involve subtracting the deviance for a saturated model? To distinguish these two possibilities we refer to absolute deviance (minus twice the log-likelihood) and relative deviance (relative to a saturated model, e.g. Section 2.3.1 in McCullagh and Nelder 1989).

With GLMMs however, there is an additional complication involving the distinction between the likelihood and the conditional likelihood. The latter is the likelihood obtained by conditioning on the estimates of the conditional modes of the spherical random effects coefficients, whereas the likelihood itself (i.e. the unconditional likelihood) involves integrating out these coefficients. The following table summarizes how to extract the various types of deviance for a glmerMod object:

This table requires two caveats:

• If the link function involves a scale parameter (e.g. Gamma) then object@resp$aic() - 2 * getME(object, "devcomp")$dims["useSc"] is required for the absolute-conditional case.

• If adaptive Gauss-Hermite quadrature is used, then logLik(object) is currently only proportional to the absolute-unconditional log-likelihood.

Note that summary() reports the unconditional absolute deviance, -2*logLik(object).

For more information about this topic see the misc/logLikGLMM directory in the package source.

Description

Drop allowable single terms from the model: see drop1 for details of how the appropriate scope for dropping terms is determined.

Usage

## S3 method for class 'merMod'
drop1(object, scope, scale = 0,
      test = c("none", "Chisq", "user"),
      k = 2, trace = FALSE, sumFun, ...)

Arguments

Details

drop1 relies on being able to find the appropriate information within the environment of the formula of the original model. If the formula is created in an environment that does not contain the data, or other variables passed to the original model (for example, if a separate function is called to define the formula), then drop1 will fail. A workaround (see example below) is to manually specify an appropriate environment for the formula.

Value

An object of class anova summarizing the differences in fit between the models.

Examples

fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
## likelihood ratio tests
drop1(fm1,test="Chisq")
## use Kenward-Roger corrected F test, or parametric bootstrap,
## to test the significance of each dropped predictor
if (require(pbkrtest) && packageVersion("pbkrtest")>="0.3.8") {
   KRSumFun <- function(object, objectDrop, ...) {
      krnames <- c("ndf","ddf","Fstat","p.value","F.scaling")
      r <- if (missing(objectDrop)) {
          setNames(rep(NA,length(krnames)),krnames)
      } else {
         krtest <- KRmodcomp(object,objectDrop)
         unlist(krtest$stats[krnames])
      }
      attr(r,"method") <- c("Kenward-Roger via pbkrtest package")
      r
   }
   drop1(fm1, test="user", sumFun=KRSumFun)

   if(lme4:::testLevel() >= 3) { ## takes about 16 sec
     nsim <- 100
     PBSumFun <- function(object, objectDrop, ...) {
	pbnames <- c("stat","p.value")
	r <- if (missing(objectDrop)) {
	    setNames(rep(NA,length(pbnames)),pbnames)
	} else {
	   pbtest <- PBmodcomp(object,objectDrop,nsim=nsim)
	   unlist(pbtest$test[2,pbnames])
	}
	attr(r,"method") <- c("Parametric bootstrap via pbkrtest package")
	r
     }
     system.time(drop1(fm1, test="user", sumFun=PBSumFun))
   }
}
## workaround for creating a formula in a separate environment
createFormula <- function(resp, fixed, rand) {  
    f <- reformulate(c(fixed,rand),response=resp)
    ## use the parent (createModel) environment, not the
    ## environment of this function (which does not contain 'data')
    environment(f) <- parent.frame()
    f
}
createModel <- function(data) {
    mf.final <- createFormula("Reaction", "Days", "(Days|Subject)")
    lmer(mf.final, data=data)
}
drop1(createModel(data=sleepstudy))

Description

Largely a wrapper for model.matrix that accepts a factor, f, and returns a dummy matrix with nlevels(f)-1 columns (the first column is dropped by default). Useful whenever one wishes to avoid the behaviour of model.matrix of always returning an nlevels(f)-column matrix, either by the addition of an intercept column, or by keeping one column for all levels.

Usage

dummy(f, levelsToKeep)

Arguments

Value

A model.matrix with dummy variables as columns.

Examples

data(Orthodont,package="nlme")
lmer(distance ~ age + (age|Subject) +
     (0+dummy(Sex, "Female")|Subject), data = Orthodont)

Description

The Dyestuff data frame provides the yield of dyestuff (Naphthalene Black 12B) from 5 different preparations from each of 6 different batchs of an intermediate product (H-acid). The Dyestuff2 data were generated data in the same structure but with a large residual variance relative to the batch variance.

Format

Data frames, each with 30 observations on the following 2 variables.

Batch

a factor indicating the batch of the intermediate product from which the preparation was created.

Yield

the yield of dyestuff from the preparation (grams of standard color).

Details

The Dyestuff data are described in Davies and Goldsmith (1972) as coming from “an investigation to find out how much the variation from batch to batch in the quality of an intermediate product (H-acid) contributes to the variation in the yield of the dyestuff (Naphthalene Black 12B) made from it. In the experiment six samples of the intermediate, representing different batches of works manufacture, were obtained, and five preparations of the dyestuff were made in the laboratory from each sample. The equivalent yield of each preparation as grams of standard colour was determined by dye-trial.”

The Dyestuff2 data are described in Box and Tiao (1973) as illustrating “ the case where between-batches mean square is less than the within-batches mean square. These data had to be constructed for although examples of this sort undoubtably occur in practice, they seem to be rarely published.”

References

• Davies OL, Goldsmith PL (eds.) (1972). Statistical Methods in Research and Production, 4th edition. Oliver and Boyd. Section 6.4.

• Box GEP, Tiao GC (1973). Bayesian Inference in Statistical Analysis. Addison-Wesley. Section 5.1.2.

Examples

require(lattice)
str(Dyestuff)
dotplot(reorder(Batch, Yield) ~ Yield, Dyestuff,
        ylab = "Batch", jitter.y = TRUE, aspect = 0.3,
        type = c("p", "a"))
dotplot(reorder(Batch, Yield) ~ Yield, Dyestuff2,
        ylab = "Batch", jitter.y = TRUE, aspect = 0.3,
        type = c("p", "a"))
(fm1 <- lmer(Yield ~ 1|Batch, Dyestuff))
(fm2 <- lmer(Yield ~ 1|Batch, Dyestuff2))

Description

If variables within a data frame are not factors, try to convert them. Not intended for end-user use; this is a utility function that needs to be exported, for technical reasons.

Usage

factorize(x,frloc,char.only=FALSE)

Arguments

Value

a copy of the data frame with factors converted

Description

Extract the fixed-effects estimates

Usage

## S3 method for class 'merMod'
 fixef(object, add.dropped=FALSE, noScale = NULL, ...)

Arguments

Details

Extract the estimates of the fixed-effects parameters from a fitted model.

Value

a named, numeric vector of fixed-effects estimates.

Examples

fixef(lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy))
fm2 <- lmer(Reaction ~ Days + Days2 + (1|Subject),
            data=transform(sleepstudy,Days2=Days))
fixef(fm2,add.dropped=TRUE)
## first two parameters are the same ...
stopifnot(all.equal(fixef(fm2,add.dropped=TRUE)[1:2],
                    fixef(fm2)))

Description

fortify adds information to data based on a fitted model; getData retrieves data as specified in the data argument

Usage

fortify.merMod(model, data = getData(model),
    ...)
## S3 method for class 'merMod'
getData(object)

Arguments

Details

• fortify is defined in the ggplot2 package, q.v. for more details. fortify is not defined here, and fortify.merMod is defined as a function rather than an S3 method, to avoid (1) inducing a dependency on ggplot2 or (2) masking methods from ggplot2. This feature is both experimental and semi-deprecated, as the help page for fortify itself says: “Rather than using this function, I now recommend using the broom package, which implements a much wider range of methods. fortify may be deprecated in the future.” The broom.mixed package is recommended for mixed models in general.

• getData is a bare-bones implementation; it relies on a data argument having been specified and the data being available in the environment of the formula. Unlike the functions in the nlme package, it does not do anything special with na.action or subset.

Examples

fm1 <- lmer(Reaction~Days+(1|Subject),sleepstudy)
  names(fortify.merMod(fm1))

Description

Retrieve information about whether random effect terms specified with double bars (||) should be split semantically into separate, intercept-free terms, or encoded with a diag() special

Usage

getDoublevertDefault()

Details

To set the global behaviour retrieved by getDoublevertDefault() to "diag_special", run options(lme4.doublevert.default = "diag_special")

The default value of "split" is backward compatible with the behaviour of lme4 < 2.0; double bars do not work as typically expected with factors as varying terms. That is, if g is a factor then (1 + g||f) will expand (with method "split") to (1|f) + (0 + g|f); the second term will not produce a covariance matrix where the contrasts generated by the factor are independent.

"diag_special" is probably the generally preferred option; by transforming (1+a+b||f) into diag(1+a+b|f), it ensures that variation across groups in the contrasts coded by the factor will be independent. This matches the results of specifying a double-bar term in afex::mixed(), or glmmTMB.

Note that if f is a factor, even with "diag_special" set, models containing (0+f||g) vs. (f||g) will give different results. (0+f||g) creates a set of dummy variables that are indicator variables for each level of the factor, so the random effect term will describe variation in the intercept for each level of the factor. With default treatment contrasts set, (f||g) (equivalent to (1+f||g) because R automatically adds the intercept term) will estimate variation first in the intercept in the baseline level of the factor, then in the differences between observations in the baseline and successive levels of the factor (see contr.treatment).

Value

a (1-element) character vector suitable for reformulas::findbars_x specifying whether to expand random-effects terms containing double bars (||) as a set of independent random effects terms ("split": i.e., (1+a+b||f) becomes (1|f) + (0+a|f) + (0+b|f)) or whether to specify that the covariance matrix is diagonal ("diag_special": i.e., diag(1+a+b|f)).

Description

Extract (or “get”) “components” – in a generalized sense – from a fitted mixed-effects model, i.e., (in this version of the package) from an object of class "merMod".

Usage

getME(object, name, ...)

## S3 method for class 'merMod'
getME(object,
      name = c("X", "Z", "Zt", "Ztlist", "mmList", "y", "mu", "u", "b",
               "Gp", "Tp", "L", "Lambda", "Lambdat", "Lind", "Tlist",
               "A", "RX", "RZX", "sigma", "flist",
               "fixef", "beta", "theta", "ST", "par", "REML", "is_REML",
               "n_rtrms", "n_rfacs", "N", "n", "p", "q", "npar",
               "p_i", "l_i", "q_i", "k", "m_i", "m",
               "cnms", "devcomp", "offset", "lower", "devfun", "devarg",
               "glmer.nb.theta"),
      ...)

Arguments

Details

The goal is to provide “everything a user may want” from a fitted "merMod" object as far as it is not available by methods, such as fixef, ranef, vcov, etc.

Value

Unspecified, as very much depending on the name.

See Also

getCall(). More standard methods for "merMod" objects, such as ranef, fixef, vcov, etc.: see methods(class="merMod")

Examples

## shows many methods you should consider *before* using getME():
methods(class = "merMod")

(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
Z <- getME(fm1, "Z")
stopifnot(is(Z, "CsparseMatrix"),
          c(180,36) == dim(Z),
	  all.equal(fixef(fm1), b1 <- getME(fm1, "beta"),
		    check.attributes=FALSE, tolerance = 0))

## A way to get *all* getME()s :
## internal consistency check ensuring that all work:
parts <- getME(fm1, "ALL")
str(parts, max=2)
stopifnot(identical(Z,  parts $ Z),
          identical(b1, parts $ beta))

Description

Returns a list with one covariance object for each random effect in the model. If the object is an old-style (pre-version 2.0.0) fitted model, it automatically converts all terms to Covariance.us objects (except for terms including double-bar ‘||’ notation, which are converged to Covariance.diag objects).

Usage

getReCovs(object)

Arguments

Value

a list of covariance objects

Examples

m1 <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy)
  sapply(getReCovs(m1), "class")

Description

Create a univariate Gauss-Hermite quadrature rule.

Usage

GHrule(ord, asMatrix = TRUE)

Arguments

Details

This version of Gauss-Hermite quadrature provides the node positions and weights for a scalar integral of a function multiplied by the standard normal density.

Originally based on package SparseGrid's hidden GQN(), then on fastGHQuad's gaussHermiteData(.).

Value

a matrix (or data frame, is asMatrix is false) with ord rows and three columns which are z the node positions, w the weights and ldnorm, the logarithm of the normal density evaluated at the nodes.

References

Liu Q, Pierce DA (1994). “A note on Gauss—Hermite quadrature.” Biometrika, 81(3), 624–629. doi:10.2307/2337136.

See Also

a different interface is available via GQdk().

Examples

(r5  <- GHrule( 5, asMatrix=FALSE))
(r12 <- GHrule(12, asMatrix=FALSE))

## second, fourth, sixth, eighth and tenth central moments of the
## standard Gaussian N(0,1) density:
ps <- seq(2, 10, by = 2)
cbind(p = ps, "E[X^p]" = with(r5,  sapply(ps, function(p) sum(w * z^p)))) # p=10 is wrong for 5-rule
p <- 1:15
GQ12 <- with(r12, sapply(p, function(p) sum(w * z^p)))
cbind(p = p, "E[X^p]" = zapsmall(GQ12))
## standard R numerical integration can do it too:
intL <- lapply(p, function(p) integrate(function(x) x^p * dnorm(x),
                                        -Inf, Inf, rel.tol=1e-11))
integR <- sapply(intL, `[[`, "value")
cbind(p, "E[X^p]" = integR)# no zapsmall() needed here
all.equal(GQ12, integR, tol=0)# => shows small difference
stopifnot(all.equal(GQ12, integR, tol = 1e-10))
(xactMom <- cumprod(seq(1,13, by=2)))
stopifnot(all.equal(xactMom, GQ12[2*(1:7)], tol=1e-14))
## mean relative errors :
mean(abs(GQ12  [2*(1:7)] / xactMom - 1)) # 3.17e-16
mean(abs(integR[2*(1:7)] / xactMom - 1)) # 9.52e-17 {even better}

Description

Fit a generalized linear mixed-effects model (GLMM). Both fixed effects and random effects are specified via the model formula.

Usage

glmer(formula, data = NULL, family = gaussian
    , control = glmerControl()
    , start = NULL
    , verbose = 0L
    , nAGQ = 1L
    , subset, weights, na.action, offset, contrasts = NULL
    , mustart, etastart
    , devFunOnly = FALSE)

Arguments

Details

Fit a generalized linear mixed model, which incorporates both fixed-effects parameters and random effects in a linear predictor, via maximum likelihood. The linear predictor is related to the conditional mean of the response through the inverse link function defined in the GLM family.

The expression for the likelihood of a mixed-effects model is an integral over the random effects space. For a linear mixed-effects model (LMM), as fit by lmer, this integral can be evaluated exactly. For a GLMM the integral must be approximated. The most reliable approximation for GLMMs is adaptive Gauss-Hermite quadrature, at present implemented only for models with a single scalar random effect. The nAGQ argument controls the number of nodes in the quadrature formula. A model with a single, scalar random-effects term could reasonably use up to 25 quadrature points per scalar integral (although the code allows for up to 100 quadrature points).

Value

An object of class merMod (more specifically, an object of subclass glmerMod) for which many methods are available (e.g. methods(class="merMod"))

Note

In earlier version of the lme4 package, a method argument was used. Its functionality has been replaced by the nAGQ argument.

See Also

lmer (for details on formulas and parameterization); glm for Generalized Linear Models (without random effects). nlmer for nonlinear mixed-effects models.

glmer.nb to fit negative binomial GLMMs.

Examples

## generalized linear mixed model
library(lattice)
xyplot(incidence/size ~ period|herd, cbpp, type=c('g','p','l'),
       layout=c(3,5), index.cond = function(x,y)max(y))
(gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
              data = cbpp, family = binomial))
## using nAGQ=0 only gets close to the optimum
(gm1a <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               cbpp, binomial, nAGQ = 0))
## using  nAGQ = 9  provides a better evaluation of the deviance
## Currently the internal calculations use the sum of deviance residuals,
## which is not directly comparable with the nAGQ=0 or nAGQ=1 result.
## 'verbose = 1' monitors iteratin a bit; (verbose = 2 does more):
(gm1a <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               cbpp, binomial, verbose = 1, nAGQ = 9))

## GLMM with individual-level variability (accounting for overdispersion)
## For this data set the model is the same as one allowing for a period:herd
## interaction, which the plot indicates could be needed.
cbpp$obs <- 1:nrow(cbpp)
(gm2 <- glmer(cbind(incidence, size - incidence) ~ period +
    (1 | herd) +  (1|obs),
              family = binomial, data = cbpp))
anova(gm1,gm2)

## glmer and glm log-likelihoods are consistent
gm1Devfun <- update(gm1,devFunOnly=TRUE)
gm0 <- glm(cbind(incidence, size - incidence) ~ period,
           family = binomial, data = cbpp)
## evaluate GLMM deviance at RE variance=theta=0, beta=(GLM coeffs)
gm1Dev0 <- gm1Devfun(c(0,coef(gm0)))
## compare
stopifnot(all.equal(gm1Dev0,c(-2*logLik(gm0))))
## the toenail oncholysis data from Backer et al 1998
## these data are notoriously difficult to fit
## Not run: 
if (require("HSAUR3")) {
    gm2 <- glmer(outcome~treatment*visit+(1|patientID),
                 data=toenail,
                 family=binomial,nAGQ=20)
}

## End(Not run)

Description

Fits a generalized linear mixed-effects model (GLMM) for the negative binomial family, building on glmer, and initializing via theta.ml from MASS.

Usage

glmer.nb(..., interval = log(th) + c(-3, 3),
         tol = 5e-5, verbose = FALSE, nb.control = NULL,
         initCtrl = list(limit = 20, eps = 2*tol, trace = verbose,
                         theta = NULL))

Arguments

Value

An object of class glmerMod, for which many methods are available (e.g. methods(class="glmerMod")), see glmer.

Note

For historical reasons, the shape parameter of the negative binomial and the random effects parameters in our (G)LMM models are both called theta ($\theta$), but are unrelated here.

The negative binomial $\theta$ can be extracted from a fit g <- glmer.nb() by getME(g, "glmer.nb.theta").

Parts of glmer.nb() are still experimental and methods are still missing or suboptimal. In particular, there is no inference available for the dispersion parameter $\theta$, yet.

To fit a negative binomial model with known overdispersion parameter (e.g. as part of a model comparison exercise, use glmer with the negative.binomial family from the MASS package, e.g. glmer(...,family=MASS::negative.binomial(theta=1.75)).

See Also

glmer; from package MASS, negative.binomial (which we re-export currently) and theta.ml, the latter for initialization of optimization.

The ‘Details’ of pnbinom for the definition of the negative binomial distribution.

Examples

set.seed(101)
dd <- expand.grid(f1 = factor(1:3),
                  f2 = LETTERS[1:2], g=1:9, rep=1:15,
          KEEP.OUT.ATTRS=FALSE)
summary(mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2))))
dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)
str(dd)
require("MASS")## and use its glm.nb() - as indeed we have zero random effect:
## Not run: 
m.glm <- glm.nb(y ~ f1*f2, data=dd, trace=TRUE)
summary(m.glm)
m.nb <- glmer.nb(y ~ f1*f2 + (1|g), data=dd, verbose=TRUE)
m.nb
## The neg.binomial theta parameter:
getME(m.nb, "glmer.nb.theta")
LL <- logLik(m.nb)
## mixed model has 1 additional parameter (RE variance)
stopifnot(attr(LL,"df")==attr(logLik(m.glm),"df")+1)
plot(m.nb, resid(.) ~ g)# works, as long as data 'dd' is found

## End(Not run)

Description

Handle for calling the glmerLaplace C++ function. Not intended for routine use.

Usage

glmerLaplaceHandle(pp, resp, nAGQ, tol, maxit, verbose)

Arguments

Value

Value of the objective function

Generator object for the glmFamily class

Description

The generator object for the glmFamily reference class. Such an object is primarily used through its new method.

Usage

glmFamily(...)

Arguments

Methods

new(family=family)

Create a new glmFamily object

Note

Arguments to the new method must be named arguments.

See Also

glmFamily

Class "glmFamily" - a reference class for family

Description

This class is a wrapper class for family objects specifying a distibution family and link function for a generalized linear model (glm). The reference class contains an external pointer to a C++ object representing the class. For common families and link functions the functions in the family are implemented in compiled code so they can be accessed from other compiled code and for a speed boost.

Extends

All reference classes extend and inherit methods from "envRefClass".

Note

Objects from this reference class correspond to objects in a C++ class. Methods are invoked on the C++ class using the external pointer in the Ptr field. When saving such an object the external pointer is converted to a null pointer, which is why there is a redundant field ptr that is an active-binding function returning the external pointer. If the Ptr field is a null pointer, the external pointer is regenerated for the stored family field.

See Also

family, glmFamily

Examples

str(glmFamily$new(family=poisson()))

Description

"golden" is a reference class for a golden search scalar optimizer, for a parameter within an interval.

golden() is the generator for the "golden" class. The optimizer uses reverse communications.

Usage

golden(...)

Arguments

Extends

All reference classes extend and inherit methods from "envRefClass".

Examples

showClass("golden")

golden(lower= -100, upper= 1e100)

Description

Data of fresh Gopher tortoise shell remains (Ozgul et al. 2009).

Usage

data("gopherdat2")

Format

A data frame with 30 observations on the following 7 variables.

Site

a factor representing the sampling site. There are 10 levels:

• BS (Big Shoals State Park).

• CB (Camp Blanding Wildlife Management Area).

• Cent (privately-owned property in central Florida).

• CF (Cecil Field/Branan Field Wildlife and Environmental Area).

• FC (Fort Cooper State Park).

• FE (Flying Eagle Wildlife Management Area).

• GH (Gold Head Branch State Park).

• Old (Perry Oldenburg Wildlife and Environmental Area).

• Ord (Ordway-Swisher Biological Station).

• TE (Tenoroc Fish Management Area).

year

a numeric vector of the sampling year.

shells

a numeric vector of the number of shells found.

Area

a numeric vector representing site area (units unknown).

density

a numeric vector representing population density.

prev

a numeric vector representing the seroprevalence (frequency of antibodies to disease) of M. agassizi.

Details

Details of the study is described in (Ozgul et al. 2009) as follows: "The fieldwork was conducted between 2003 and 2006, during late spring/summer (May–September). Systematic surveys were conducted to locate tortoise burrows and shell remains from deceased tortoises at study sites and consisted of a line of four to eight observers spaced 10 m apart walking parallel transects across the study area."

Source

(Ozgul et al. 2009)

References

Ozgul A, Oli MK, Bolker BM, Perez-Heydrich C (2009). “Upper respiratory tract disease, force of infection, and effects on survival of gopher tortoises.” Ecological Applications, 19(3), 786-798. doi:10.1890/08-0219.1.

Examples

## Simple model gives a singular fit:
gopher_glmer <- glmer(shells ~ factor(year) + prev + offset(log(Area))
                   + (1|Site), data = gopherdat2, family = "poisson")
## The site-level variance for this model is indeed zero:
VarCorr(gopher_glmer)
## So a Poisson GLM gives the same answer here:
gopher_glm <- glm(shells ~ factor(year) + prev + offset(log(Area)),
                  data = gopherdat2, family = "poisson")
all.equal(fixef(gopher_glmer), coef(gopher_glm))

Description

Generate the sparse multidimensional Gaussian quadrature grids.

Currently unused. See GHrule() for the version currently in use in package lme4.

Usage

GQdk(d = 1L, k = 1L)
  GQN

Arguments

Value

GQdk() returns a matrix with d + 1 columns. The first column is the weights and the remaining d columns are the node coordinates.

GQN is a list of lists, containing the non-redundant quadrature nodes and weights for integration of a scalar function of a d-dimensional argument with respect to the density function of the d-dimensional Gaussian density function.
The outer list is indexed by the dimension, d, in the range of 1 to 20. The inner list is indexed by k, the order of the quadrature.

Note

GQN contains only the non-redundant nodes. To regenerate the whole array of nodes, all possible permutations of axes and all possible combinations of $\pm 1$ must be applied to the axes. This entire array of nodes is exactly what GQdk() reproduces.

The number of nodes gets very large very quickly with increasing d and k. See the charts at http://www.sparse-grids.de.

Examples

GQdk(2,5) # 53 x 3

GQN[[3]][[5]] # a 14 x 4 matrix

Description

Number of ticks on the heads of red grouse chicks sampled in the field (grouseticks) and an aggregated version (grouseticks_agg); see original source for more details

Usage

data(grouseticks)

Format

INDEX

(factor) chick number (observation level)

TICKS

number of ticks sampled

BROOD

(factor) brood number

HEIGHT

height above sea level (meters)

YEAR

year (-1900)

LOCATION

(factor) geographic location code

cHEIGHT

centered height, derived from HEIGHT

meanTICKS

mean number of ticks by brood

varTICKS

variance of number of ticks by brood

Details

grouseticks_agg is just a brood-level aggregation of the data

Source

Robert Moss, via David Elston

References

Elston DA, Moss R, Boulinier T, Arrowsmith C, Lambin X (2001). “Analysis of aggregation, a worked example: numbers of ticks on red grouse chicks.” Parasitology, 122(5), 563–569. doi:10.1017/S0031182001007740.

Examples

if (interactive()) {
data(grouseticks)
## Figure 1a from Elston et al
par(las=1,bty="l")
tvec <- c(0,1,2,5,20,40,80)
pvec <- c(4,1,3)
with(grouseticks_agg,plot(1+meanTICKS~HEIGHT,
                  pch=pvec[factor(YEAR)],
                  log="y",axes=FALSE,
                  xlab="Altitude (m)",
                  ylab="Brood mean ticks"))
axis(side=1)
axis(side=2,at=tvec+1,label=tvec)
box()
abline(v=405,lty=2)
## Figure 1b
with(grouseticks_agg,plot(varTICKS~meanTICKS,
                  pch=4,
                  xlab="Brood mean ticks",
                  ylab="Within-brood variance"))
curve(1*x,from=0,to=70,add=TRUE)
## Model fitting
form <- TICKS~YEAR+HEIGHT+(1|BROOD)+(1|INDEX)+(1|LOCATION)
(full_mod1  <- glmer(form, family="poisson",data=grouseticks))
}

Description

Returns the values on the diagonal of the hat matrix, which is the matrix that transforms the response vector (minus any offset) into the fitted values (minus any offset). Note that this method should only be used for linear mixed models. It is not clear if the hat matrix concept even makes sense for generalized linear mixed models.

Usage

## S3 method for class 'merMod'
 hatvalues(model, fullHatMatrix = FALSE, ...)

Arguments

Value

The diagonal elements of the hat matrix.

Examples

m <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
hatvalues(m)

Description

These functions compute deletion influence diagnostics for linear (fit by lmer) and generalized linear mixed-effects models (fit by glmer). The main functions are methods for the influence generic function. Other functions are provided for computing dfbeta, dfbetas, cooks.distance, and influence on variance-covariance components based on the objects computed by influence.merMod

Usage

## S3 method for class 'merMod'
influence(model, groups, data, maxfun = 1000,
          do.coef = TRUE, start = NULL,
          parallel = c("no", "multicore", "snow", "future"),
          ncpus = getOption("influence.ncpus", 1L), cl = NULL, ncores, ...)
## S3 method for class 'influence.merMod'
cooks.distance(model, ...)
## S3 method for class 'influence.merMod'
dfbeta(model, which = c("fixed", "var.cov"), ...)
## S3 method for class 'influence.merMod'
dfbetas(model, ...)

Arguments

Details

influence.merMod start with the estimated variance-covariance components from model and then refit the model omitting each group in turn, not necessarily iterating to completion. For example, maxfun=20 takes up to 20 function evaluations step away from the ML or REML solution for the full data, which usually provides decent approximations to the fully iterated estimates.

The other functions are methods for the dfbeta, dfbetas, and cooks.distance generics, to be applied to the "influence.merMod" object produced by the influence function; the dfbeta methods can also return influence on the variance-covariance components.

Value

influence.merMod returns objects of class "influence.merMod", which contain the following elements:

"fixed.effects"

the estimated fixed effects for the model.

"fixed.effects[-groups]"

a matrix with columns corresponding to the fixed-effects coefficients and rows corresponding to groups, giving the estimated fixed effects with each group deleted in turn; groups is formed from the name(s) of the grouping factor(s).

"var.cov.comps"

the estimated variance-covariance parameters for the model.

"var.cov.comps[-groups]"

a matrix with the estimated covariance parameters (in columns) with each group deleted in turn.

"vcov"

The estimated covariance matrix of the fixed-effects coefficients.

"vcov[-groups]"

a list each of whose elements is the estimated covariance matrix of the fixed-effects coefficients with one group deleted.

"groups"

a character vector giving the names of the grouping factors.

"deleted"

the possibly composite grouping factor, each of whose elements is deleted in turn.

"converged"

for influence.merMod, a logical vector indicating whether the computation converged for each group.

"function.evals"

for influence.merMod, a vector of the number of function evaluations performed for each group.

For plotting "influence.merMod" objects, see infIndexPlot.

Author(s)

J. Fox; parallel extensions by H. Bengtsson

References

Fox J, Weisberg S (2019). An R companion to applied regression, 3rd edition. Sage.

See Also

infIndexPlot, influence.measures

Examples

if (interactive()) {
  fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
  inf_fm1 <- influence(fm1, "Subject")
  if (require("car")) {
    infIndexPlot(inf_fm1)
  }
  dfbeta(inf_fm1)
  dfbetas(inf_fm1)
  gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               data = cbpp, family = binomial)
  inf_gm1 <- influence(gm1, "herd", maxfun=100)
  gm1.11 <- update(gm1, subset = herd != 11) # check deleting herd 11
  if (require("car")) {
    infIndexPlot(inf_gm1)
    compareCoefs(gm1, gm1.11)
  }
  if(packageVersion("car") >= "3.0.10") {
    dfbeta(inf_gm1)
    dfbetas(inf_gm1)
  }
 }

Description

University lecture evaluations by students at ETH Zurich, anonymized for privacy protection. This is an interesting “medium” sized example of a partially nested mixed effect model.

Format

A data frame with 73421 observations on the following 7 variables.

s

a factor with levels 1:2972 denoting individual students.

d

a factor with 1128 levels from 1:2160, denoting individual professors or lecturers.

studage

an ordered factor with levels 2 < 4 < 6 < 8, denoting student's “age” measured in the semester number the student has been enrolled.

lectage

an ordered factor with 6 levels, 1 < 2 < ... < 6, measuring how many semesters back the lecture rated had taken place.

service

a binary factor with levels 0 and 1; a lecture is a “service”, if held for a different department than the lecturer's main one.

dept

a factor with 14 levels from 1:15, using a random code for the department of the lecture.

y

a numeric vector of ratings of lectures by the students, using the discrete scale 1:5, with meanings of ‘poor’ to ‘very good’.

Each observation is one student's rating for a specific lecture (of one lecturer, during one semester in the past).

Details

The main goal of the survey is to find “the best liked prof”, according to the lectures given. Statistical analysis of such data has been the basis for a (student) jury selecting the final winners.

The present data set has been anonymized and slightly simplified on purpose.

Examples

str(InstEval)

head(InstEval, 16)
xtabs(~ service + dept, InstEval)

Description

isNewMerMod tests if a merMod object has the latest representation (in other words, if the object could have been generated using the current version of lme4). The representation of a merMod object was last changed in lme4 version 2.0-0.

forceNewMerMod updates a merMod object so that it has the latest representation. It is useful in code that might handle objects generated by old versions of lme4.

anyStructured tests if a merMod object specifies random effects with structured covariance matrices.

Usage

isNewMerMod(object)
forceNewMerMod(object, reference = object)
anyStructured(object)

Arguments

Value

isNewMerMod and anyStructured return TRUE or FALSE.

forceNewMerMod returns the argument if it is “current” and a modified copy of the argument otherwise.

See Also

Virtual class Covariance.

Examples

fm1 <- lmer(Reaction ~ Days +   us(Days | Subject), sleepstudy)
fm2 <- lmer(Reaction ~ Days + diag(Days | Subject), sleepstudy)
stopifnot(isNewMerMod(fm1), identical(fm1, forceNewMerMod(fm1)),
          !anyStructured(fm1), anyStructured(fm2))

Description

Check characteristics of models: whether a model fit corresponds to a linear (LMM), generalized linear (GLMM), or nonlinear (NLMM) mixed model, and whether a linear mixed model has been fitted by REML or not (isREML(x) is always FALSE for GLMMs and NLMMs).

Usage

isREML(x, ...)

  isLMM(x, ...)

  isNLMM(x, ...)

  isGLMM(x, ...)

Arguments

Details

These are generic functions. At present the only methods are for mixed-effects models of class merMod.

Value

a logical value

See Also

getME

Examples

fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
              data = cbpp, family = binomial)
nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
             Orange, start = c(Asym = 200, xmid = 725, scal = 350))

isLMM(fm1)
isGLMM(gm1)
## check all :
is.MM <- function(x) c(LMM = isLMM(x), GLMM= isGLMM(x), NLMM= isNLMM(x))
stopifnot(cbind(is.MM(fm1), is.MM(gm1), is.MM(nm1))
	  == diag(rep(TRUE,3)))

Description

Evaluates whether a fitted mixed model is (almost / near) singular, i.e., the parameters are on the boundary of the feasible parameter space: variances of one or more linear combinations of effects are (close to) zero.

Usage

isSingular(object, tol = getSingTol())
getSingTol()

Arguments

Details

Complex mixed-effect models (i.e., those with a large number of variance-covariance parameters) frequently result in singular fits, i.e. estimated variance-covariance matrices with less than full rank. Less technically, this means that some "dimensions" of the variance-covariance matrix have been estimated as exactly zero. For scalar random effects such as intercept-only models, or 2-dimensional random effects such as intercept+slope models, singularity is relatively easy to detect because it leads to random-effect variance estimates of (nearly) zero, or estimates of correlations that are (almost) exactly -1 or 1. However, for more complex models (variance-covariance matrices of dimension >=3) singularity can be hard to detect; models can often be singular without any of their individual variances being close to zero or correlations being close to +/-1.

This function performs a simple test to determine whether any of the random effects covariance matrices of a fitted model are singular. The rePCA method provides more detail about the singularity pattern, showing the standard deviations of orthogonal variance components and the mapping from variance terms in the model to orthogonal components (i.e., eigenvector/rotation matrices).

While singular models are statistically well defined (it is theoretically sensible for the true maximum likelihood estimate to correspond to a singular fit), there are real concerns that (1) singular fits correspond to overfitted models that may have poor power; (2) chances of numerical problems and mis-convergence are higher for singular models (e.g. it may be computationally difficult to compute profile confidence intervals for such models); (3) standard inferential procedures such as Wald statistics and likelihood ratio tests may be inappropriate.

There is not yet consensus about how to deal with singularity, or more generally to choose which random-effects specification (from a range of choices of varying complexity) to use. Some proposals include:

• avoid fitting overly complex models in the first place, i.e. design experiments/restrict models a priori such that the variance-covariance matrices can be estimated precisely enough to avoid singularity (Matuschek et al 2017)

• use some form of model selection to choose a model that balances predictive accuracy and overfitting/type I error (Bates et al 2015, Matuschek et al 2017)

• “keep it maximal”, i.e. fit the most complex model consistent with the experimental design, removing only terms required to allow a non-singular fit (Barr et al. 2013), or removing further terms based on p-values or AIC

• use a partially Bayesian method that produces maximum a posteriori (MAP) estimates using regularizing priors to force the estimated random-effects variance-covariance matrices away from singularity (Chung et al 2013, blme package)

• use a fully Bayesian method that both regularizes the model via informative priors and gives estimates and credible intervals for all parameters that average over the uncertainty in the random effects parameters (Gelman and Hill 2006, McElreath 2015; MCMCglmm, rstanarm and brms packages)

getSingTol() is a utility function that just returns the default tolerance for testing singularity (1e-4, in the current version); if the user has set options(lme4.singular.tolerance=...), it will retrieve this value instead.

Value

a logical value

References

• Barr DJ, Levy R, Scheepers C, Tily HJ (2013). “Random effects structure for confirmatory hypothesis testing: Keep it maximal.” Journal of Memory and Language, 68(3), 255–278. doi:10.1016/j.jml.2012.11.001.

• Bates D, Kliegl R, Vasishth S, Baayen H (2015). “Parsimonious Mixed Models.” arXiv preprint. https://arxiv.org/abs/1506.04967.

• Chung Y, Rabe-Hesketh S, Dorie V, Gelman A, Liu J (2013). “A nondegenerate penalized likelihood estimator for variance parameters in multilevel models.” Psychometrika, 78(4), 685–709. doi:10.1007/S11336-013-9328-2.

• Gelman A, Hill J (2006). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.

• Matuschek H, Kliegl R, Vasishth S, Baayen H, Bates D (2017). “Balancing Type I error and power in linear mixed models.” Journal of Memory and Language, 94, 305-315. ISSN 0749-596X. doi:10.1016/j.jml.2017.01.001.

• McElreath R (2015). Statistical Rethinking: A Bayesian Course with Examples in R and Stan. Chapman and Hall/CRC.

See Also

getME, rePCA

Description

Reads the environment variable LME4_TEST_LEVEL to determine which tests and examples to run

Usage

lme4_testlevel()

Value

a numeric value: 1 for standard/'light' testing, larger values for more testing. Defaults to 1 if the environment variable is not set.

Description

Fit a linear mixed-effects model (LMM) to data, via restricted maximum likelihood (REML) or maximum likelihood.

Usage

lmer(formula, data = NULL, REML = TRUE, control = lmerControl(),
     start = NULL, verbose = 0L, subset, weights, na.action,
     offset, contrasts = NULL, devFunOnly = FALSE)

Arguments

Details

• If the formula argument is specified as a character vector, the function will attempt to coerce it to a formula. However, this is not recommended (users who want to construct formulas by pasting together components are advised to use as.formula or reformulate); model fits will work but subsequent methods such as drop1, update, etc. may fail.

• When handling perfectly collinear predictor variables (i.e. fixed-effect design matrices of less than full rank), [gn]lmer is not as sophisticated as modeling frameworks such as lm and glm. While it does automatically drop collinear variables (with a message rather than a warning), it does not automatically fill in NA values for the dropped coefficients; these can be added via fixef(fitted.model, add.dropped=TRUE). This information can also be retrieved via attr(getME(fitted.model, "X"), "col.dropped").

• the deviance function returned when devFunOnly is TRUE takes a single numeric vector argument which defines the scaled variance-covariance matrices of the random effects.

  • In the case of unstructured covariances, this vector is directly mapped to the theta vector, which represents the unique non-zero values in the Cholesky factor of the (scaled) covariance matrix. For models with only simple (intercept-only) random effects, par (and thus theta) is a vector of the standard deviations of the random effects. For more complex or multiple random effects, running getME(.,"par") or (equivalently) getME(., "theta") to retrieve the theta vector for a fitted model and examining the names of the vector is probably the easiest way to determine the correspondence between the elements of the theta vector and elements of the lower triangles of the Cholesky factors of the random effects.

  • For structured covariances, the getTheta method translates the parameter vector to the theta (Cholesky-factor element) scale for internal use. The parameter vector is usually composed of a set of standard-deviation values (one if hom = TRUE or many if hom = FALSE), followed by one or more parameters that determine the correlation matrix.

Value

An object of class merMod (more specifically, an object of subclass lmerMod), for which many methods are available (e.g. methods(class="merMod"))

Note

In earlier version of the lme4 package, a method argument was used. Its functionality has been replaced by the REML argument.

Also, lmer(.) allowed a family argument (to effectively switch to glmer(.)). This has been deprecated in summer 2013, and been disabled in spring 2019.

See Also

lm for linear models; glmer for generalized linear; and nlmer for nonlinear mixed models.

plot.merMod for plot diagnostics.

Examples

## linear mixed models - reference values from older code
(fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy))
summary(fm1) # (with its own print method; see class?merMod % ./merMod-class.Rd
plot(fm1) # plotting the model diagnostics; see ?plot.merMod

str(terms(fm1))
stopifnot(identical(terms(fm1, fixed.only=FALSE),
                    terms(model.frame(fm1))))
attr(terms(fm1, FALSE), "dataClasses") # fixed.only=FALSE needed for dataCl.

## Maximum Likelihood (ML), and "monitor" iterations via 'verbose':
fm1_ML <- update(fm1, REML=FALSE, verbose = 1)
(fm2 <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy))
anova(fm1, fm2)
sm2 <- summary(fm2)
print(fm2, digits=7, ranef.comp="Var") # the print.merMod()         method
print(sm2, digits=3, corr=FALSE)       # the print.summary.merMod() method

## Fit sex-specific variances by constructing numeric dummy variables
## for sex and sex:age; in this case the estimated variance differences
## between groups in both intercept and slope are zero ...
data(Orthodont,package="nlme")
Orthodont$nsex <- as.numeric(Orthodont$Sex=="Male")
Orthodont$nsexage <- with(Orthodont, nsex*age)
lmer(distance ~ age + (age|Subject) + (0+nsex|Subject) +
     (0 + nsexage|Subject), data=Orthodont)

Description

Construct control structures for mixed model fitting. All arguments have defaults, and can be grouped into

• general control parameters, most importantly optimizer, further restart_edge, etc;

• model- or data-checking specifications, in short “checking options”, such as check.nobs.vs.rankZ, or check.rankX (currently not for nlmerControl);

• all the parameters to be passed to the optimizer, e.g., maximal number of iterations, passed via the optCtrl list argument.

Usage

lmerControl(optimizer = "nloptwrap",
				    
    restart_edge = TRUE,
    boundary.tol = 1e-5,
    calc.derivs = NULL,
    use.last.params = FALSE,
    sparseX = FALSE,
    standardize.X = FALSE,
    autoscale = NULL,
    ## input checking options
    check.nobs.vs.rankZ = "ignore",
    check.nobs.vs.nlev = "stop",
    check.nlev.gtreq.5 = "ignore",
    check.nlev.gtr.1 = "stop",
    check.nobs.vs.nRE= "stop",
    check.rankX = c("message+drop.cols", "silent.drop.cols", "warn+drop.cols",
                    "stop.deficient", "ignore"),
    check.scaleX = c("warning","stop","silent.rescale",
                     "message+rescale","warn+rescale","ignore"),
    check.formula.LHS = "stop",
    ## convergence checking options
    check.conv.nobsmax  = 1e4,
    check.conv.nparmax  = 10,
    check.conv.grad     = .makeCC("warning", tol = 2e-3, relTol = NULL),
    check.conv.singular = .makeCC(action = "message", tol = getSingTol()),
    check.conv.hess     = .makeCC(action = "warning", tol = 1e-6),
    ## optimizer args
    optCtrl = list(),
    mod.type = "lmer"
)

glmerControl(optimizer = c("bobyqa", "Nelder_Mead"),
    restart_edge = FALSE,
    boundary.tol = 1e-5,
    calc.derivs = NULL,
    use.last.params = FALSE,
    sparseX = FALSE,
    standardize.X = FALSE,
    autoscale = NULL,
    ## input checking options
    check.nobs.vs.rankZ = "ignore",
    check.nobs.vs.nlev = "stop",
    check.nlev.gtreq.5 = "ignore",
    check.nlev.gtr.1 = "stop",
    check.nobs.vs.nRE= "stop",
    check.rankX = c("message+drop.cols", "silent.drop.cols", "warn+drop.cols",
                    "stop.deficient", "ignore"),
    check.scaleX = c("warning","stop","silent.rescale",
                     "message+rescale","warn+rescale","ignore"),
    check.formula.LHS = "stop",
    ## convergence checking options
    check.conv.nobsmax  = 1e4, 
    check.conv.nparmax  = 20,
    check.conv.grad     = .makeCC("warning", tol = 2e-3, relTol = NULL),
    check.conv.singular = .makeCC(action = "message", tol = getSingTol()),
    check.conv.hess     = .makeCC(action = "warning", tol = 1e-6),
    ## optimizer args
    optCtrl = list(),
    mod.type = "glmer",
    tolPwrss = 1e-7,
    compDev = TRUE,
    nAGQ0initStep = TRUE,
    check.response.not.const = "stop"
 )

nlmerControl(optimizer = "Nelder_Mead", tolPwrss = 1e-10,
             optCtrl = list())

.makeCC(action, tol, relTol, ...)

Arguments