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"),
       ncpus = getOption("allFit.ncpus", 1L), cl = NULL,
       catch.errs = 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 <- detectCores()-1
  optCls <- makeCluster(nc, type = "SOCK")
  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

From Josh Banta

References

Joshua A. Banta, Martin H. H Stevens, and Massimo Pigliucci (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"),
	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, A.C. and Hinkley, D.V. (1997) Bootstrap Methods and Their Application. Cambridge University Press.

Morris, J. S. (2002). The BLUPs Are Not ‘best’ When It Comes to Bootstrapping. Statistics & Probability Letters 56(4): 425–430. 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, F. E. (1938) Chocolate cake, I. Optimum baking temperature. Master's Thesis, Iowa State College.

Cochran, W. G., and Cox, G. M. (1957) Experimental designs, 2nd Ed. New York, John Wiley & Sons.

Lee, Y., Nelder, J. A., and Pawitan, Y. (2006) Generalized linear models with random effects. Unified analysis via H-likelihood. Boca Raton, Chapman and Hall/CRC.

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).

Format

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.

Details

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

Source

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, 27–40.

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))

Description

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

Usage

checkConv(derivs, coefs, ctrl, lbound, debug = FALSE)

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 = TRUE, ...)
## 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, signames=oldNames, ...),
            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
                              oldNames = FALSE))
  ## --> ~ 2.2 seconds (2022)
  set.seed(123) # (reproducibility when using bootstrap)
  system.time(fm1B <- confint(fm1, method="boot", oldNames=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

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,
        transfuns = list(from.chol = Cv_to_Sv,
                           to.chol = Sv_to_Cv,
                             to.sd = identity), ...)

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

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.”

Source

O.L. Davies and P.L. Goldsmith (eds), Statistical Methods in Research and Production, 4th ed., Oliver and Boyd, (1972), section 6.4

G.E.P. Box and G.C. Tiao, Bayesian Inference in Statistical Analysis, Addison-Wesley, (1973), 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

From the right hand side of a formula for a mixed-effects model, expand terms with the double vertical bar operator into separate, independent random effect terms.

Usage

expandDoubleVerts(term)

Arguments

Value

the modified term

Note

Because || works at the level of formula parsing, it has no way of knowing whether a variable is a factor. It just takes the terms within a random-effects term and literally splits them into the intercept and separate no-intercept terms, e.g. (1+x+y|f) would be split into (1|f) + (0+x|f) + (0+y|f). However, || will fail to break up factors into separate terms; the dummy function can be useful in this case, although it is not as convenient as ||.

See Also

formula, model.frame, model.matrix, dummy.

Other utilities: mkRespMod, mkReTrms, nlformula, nobars, subbars

Examples

m <- ~ x + (x || g)
  expandDoubleVerts(m)
  set.seed(101)
  dd <- expand.grid(f=factor(letters[1:3]),g=factor(1:200),rep=1:3)
  dd$y <- simulate(~f + (1|g) + (0+dummy(f,"b")|g) + (0+dummy(f,"c")|g),
          newdata=dd,
          newparams=list(beta=rep(0,3),
                         theta=c(1,2,1),
                         sigma=1),
          family=gaussian)[[1]]
  m1 <- lmer(y~f+(f|g),data=dd)
  VarCorr(m1)
  m2 <- lmer(y~f+(1|g) + (0+dummy(f,"b")|g) + (0+dummy(f,"c")|g),
               data=dd)
  VarCorr(m2)

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

From the right hand side of a formula for a mixed-effects model, determine the pairs of expressions that are separated by the vertical bar operator. Also expand the slash operator in grouping factor expressions and expand terms with the double vertical bar operator into separate, independent random effect terms.

Usage

findbars(term)

Arguments

Value

pairs of expressions that were separated by vertical bars

Note

This function is called recursively on individual terms in the model, which is why the argument is called term and not a name like form, indicating a formula.

See Also

formula, model.frame, model.matrix.

Other utilities: mkRespMod, mkReTrms, nlformula, nobars, subbars

Examples

findbars(f1 <- Reaction ~ Days + (Days | Subject))
## => list( Days | Subject )
## These two are equivalent:% tests in ../inst/tests/test-doubleVertNotation.R
findbars(y ~ Days + (1 | Subject) + (0 + Days | Subject))
findbars(y ~ Days + (Days || Subject))
## => list of length 2:  list ( 1 | Subject ,  0 + Days | Subject)
findbars(~ 1 + (1 | batch / cask))
## => list of length 2:  list ( 1 | cask:batch ,  1 | batch)

Description

Extract the fixed-effects estimates

Usage

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

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

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", "REML", "is_REML",
               "n_rtrms", "n_rfacs", "N", "n", "p", "q",
               "p_i", "l_i", "q_i", "k", "m_i", "m",
               "cnms", "devcomp", "offset", "lower", "devfun", "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

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

Qing Liu and Donald A. Pierce (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.

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

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, D. A., R. Moss, T. Boulinier, C. Arrowsmith, and X. Lambin. 2001. "Analysis of Aggregation, a Worked Example: Numbers of Ticks on Red Grouse Chicks." Parasitology 122 (05): 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, ncores = getOption("mc.cores",1), start, ...)
## 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 [email protected]

References

Fox, J. and Weisberg, S. (2019) An R Companion to Applied Regression, Third 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

Does every level of f1 occur in conjunction with exactly one level of f2? The function is based on converting a triplet sparse matrix to a compressed column-oriented form in which the nesting can be quickly evaluated.

Usage

isNested(f1, f2)

Arguments

Value

TRUE if factor 1 is nested within factor 2

Examples

with(Pastes, isNested(cask, batch))   ## => FALSE
with(Pastes, isNested(sample, batch))  ## => TRUE

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(x, tol = 1e-4)

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)

Value

a logical value

References

Dale J. Barr, Roger Levy, Christoph Scheepers, and Harry J. Tily (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal; Journal of Memory and Language 68(3), 255–278.

Douglas Bates, Reinhold Kliegl, Shravan Vasishth, and Harald Baayen (2015). Parsimonious Mixed Models; preprint (https://arxiv.org/abs/1506.04967).

Yeojin Chung, Sophia Rabe-Hesketh, Vincent Dorie, Andrew Gelman, and Jingchen Liu (2013). A nondegenerate penalized likelihood estimator for variance parameters in multilevel models; Psychometrika 78, 685–709; doi:10.1007/S11336-013-9328-2.

Andrew Gelman and Jennifer Hill (2006). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.

Hannes Matuschek, Reinhold Kliegl, Shravan Vasishth, Harald Baayen, and Douglas Bates (2017). Balancing type I error and power in linear mixed models. Journal of Memory and Language 94, 305–315.

Richard McElreath (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 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 may fail.

• When handling perfectly collinear predictor variables (i.e. design matrices of less than full rank), [gn]lmer is not quite as sophisticated as some simpler 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, representing the theta vector. This vector defines the scaled variance-covariance matrices of the random effects, in the Cholesky parameterization. For models with only simple (intercept-only) random effects, theta is a vector of the standard deviations of the random effects. For more complex or multiple random effects, running 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.

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.

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

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

(vv <- vcov.merMod(fm2, corr=TRUE))
as(vv, "corMatrix")# extracts the ("hidden") 'correlation' entry in @factors

## 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 = TRUE,
    use.last.params = FALSE,
    sparseX = FALSE,
    standardize.X = FALSE,
    ## 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.grad     = .makeCC("warning", tol = 2e-3, relTol = NULL),
    check.conv.singular = .makeCC(action = "message", tol = formals(isSingular)$tol),
    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 = TRUE,
    use.last.params = FALSE,
    sparseX = FALSE,
    standardize.X = FALSE,
    ## 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.grad     = .makeCC("warning", tol = 2e-3, relTol = NULL),
    check.conv.singular = .makeCC(action = "message", tol = formals(isSingular)$tol),
    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

Details

Note that (only!) the pre-fitting “checking options” (i.e., all those starting with "check." but not including the convergence checks ("check.conv.*") or rank-checking ("check.rank*") options) may also be set globally via options. In that case, (g)lmerControl will use them rather than the default values, but will not override values that are passed as explicit arguments.

For example, options(lmerControl=list(check.nobs.vs.rankZ = "ignore")) will suppress warnings that the number of observations is less than the rank of the random effects model matrix Z.

Value

The *Control functions return a list (inheriting from class "merControl") containing

1. general control parameters, such as optimizer, restart_edge;

2. (currently not for nlmerControl:) "checkControl", a list of data-checking specifications, e.g., check.nobs.vs.rankZ;

3. parameters to be passed to the optimizer, i.e., the optCtrl list, which may contain maxiter.

.makeCC returns a list containing the check specification (action, tolerance, and optionally relative tolerance).

See Also

convergence and allFit() which fits for a couple of optimizers; nloptwrap for the lmerControl() default optimizer.

Examples

str(lmerControl())
str(glmerControl())
## fit with default algorithm [nloptr version of BOBYQA] ...
fm0 <- lmer(Reaction ~ Days +   ( 1 | Subject), sleepstudy)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
## or with "bobyqa" (default 2013 - 2019-02) ...
fm1_bobyqa <- update(fm1, control = lmerControl(optimizer="bobyqa"))
## or with "Nelder_Mead" (the default till 2013) ...
fm1_NMead <- update(fm1, control = lmerControl(optimizer="Nelder_Mead"))
## or with the nlminb function used in older (<1.0) versions of lme4;
## this will usually replicate older results
if (require(optimx)) {
    fm1_nlminb <- update(fm1,
                         control = lmerControl(optimizer= "optimx",
                                               optCtrl  = list(method="nlminb")))
    ## The other option here is method="L-BFGS-B".
}

## Or we can wrap base::optim():
optimwrap <- function(fn,par,lower,upper,control=list(),
                      ...) {
    if (is.null(control$method)) stop("must specify method in optCtrl")
    method <- control$method
    control$method <- NULL
    ## "Brent" requires finite upper values (lower bound will always
    ##  be zero in this case)
    if (method=="Brent") upper <- pmin(1e4,upper)
    res <- optim(par=par, fn=fn, lower=lower,upper=upper,
                 control=control,method=method,...)
    with(res, list(par  = par,
                   fval = value,
                   feval= counts[1],
                   conv = convergence,
                   message = message))
}
fm0_brent <- update(fm0,
                    control = lmerControl(optimizer = "optimwrap",
                                          optCtrl = list(method="Brent")))

## You can also use functions (in addition to the lmerControl() default "NLOPT_BOBYQA")
## from the 'nloptr' package, see also  '?nloptwrap' :
if (require(nloptr)) {
    fm1_nloptr_NM <- update(fm1, control=lmerControl(optimizer="nloptwrap",
                                      optCtrl=list(algorithm="NLOPT_LN_NELDERMEAD")))
    fm1_nloptr_COBYLA <- update(fm1, control=lmerControl(optimizer="nloptwrap",
                                      optCtrl=list(algorithm="NLOPT_LN_COBYLA",
                                                   xtol_rel=1e-6,
                                                   xtol_abs=1e-10,
                                                   ftol_abs=1e-10)))
}
## other algorithm options include NLOPT_LN_SBPLX

Description

Fit a list of lm or glm objects with a common model for different subgroups of the data.

Usage

lmList(formula, data, family, subset, weights, na.action,
       offset, pool = !isGLM || .hasScale(family2char(family)),
       warn = TRUE, ...)

Arguments

Details

• While data is optional, the package authors strongly recommend its use, especially when later applying methods such as update and drop1 to the fitted model (such methods are not guaranteed to work properly if data is omitted). If data is omitted, variables will be taken from the environment of formula (if specified as a formula) or from the parent frame (if specified as a character vector).

• Since lme4 version 1.1-16, if there are errors (see stop) in the single (lm() or glm()) fits, they are summarized to a warning message which is returned as attribute "warnMessage" and signalled as warning() when the warn argument is true.

  In previous lme4 versions, a general (different) warning had been signalled in this case.

Value

an object of class lmList4 (see there, notably for the methods defined).

See Also

lmList4

Examples

fm.plm  <- lmList(Reaction ~ Days | Subject, sleepstudy)
coef(fm.plm)
fm.2  <- update(fm.plm, pool = FALSE)
## coefficients are the same, "pooled or unpooled":
stopifnot( all.equal(coef(fm.2), coef(fm.plm)) )

(ci <- confint(fm.plm)) # print and rather *see* :
plot(ci)                # how widely they vary for the individuals

Description

Class "lmList4" is an S4 class with basically a list of objects of class lm with a common model (but different data); see lmList() which returns these.

Package nlme's lmList() returns objects of S3 class "lmList" and provides methods for them, on which our methods partly build.

Objects from the Class

Objects can be created by calls of the form new("lmList4", ...) or, more commonly, by a call to lmList().

Methods

A dozen methods are provided. Currently, S4 methods for show, coercion (as(.,.)) and others inherited via "list", and S3 methods for coef, confint, fitted, fixef, formula, logLik, pairs, plot, predict, print, qqnorm, ranef, residuals, sigma, summary, and update.

sigma(object)

returns the standard deviation $\hat{\sigma}$ (of the errors in the linear models), assuming a common variance $\sigma^2$ by pooling (even when pool = FALSE was used in the fit).

See Also

lmList

Examples

if(getRversion() >= "3.2.0") {
  (mm <- methods(class = "lmList4"))
  ## The S3 ("not S4") ones :
  mm[!attr(mm,"info")[,"isS4"]]
}
## For more examples:  example(lmList)  i.e., ?lmList

Description

The generator objects for the lmResp, lmerResp, glmResp and nlsResp reference classes. Such objects are primarily used through their new methods.

Usage

lmResp(...)

Arguments

Methods

new(y=y):

Create a new lmResp or lmerResp object.

new(family=family, y=y):

Create a new glmResp object.

new(y=y, nlmod=nlmod, nlenv=nlenv, pnames=pnames, gam=gam):

Create a new nlsResp object.

Note

Arguments to the new methods must be named arguments.

y

the numeric response vector

family

a family object

nlmod

the nonlinear model function

nlenv

an environment holding data objects for evaluation of nlmod

pnames

a character vector of parameter names

gam

a numeric vector - the initial linear predictor

See Also

lmResp, lmerResp, glmResp, nlsResp

Description

Reference classes for response modules, including linear models, "lmResp", generalized linear models, "glmResp", nonlinear models, "nlsResp" and linear mixed-effects models, "lmerResp". Each reference class is associated with a C++ class of the same name. As is customary, the generator object for each class has the same name as the class.

Extends

All reference classes extend and inherit methods from "envRefClass". Furthermore, "glmResp", "nlsResp" and "lmerResp" all extend the "lmResp" class.

Note

Objects from these reference classes correspond to objects in C++ classes. Methods are invoked on the C++ classes 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 are redundant fields containing enough information as R objects to be able to regenerate the C++ object. The convention is that a field whose name begins with an upper-case letter is an R object and the corresponding field whose name begins with the lower-case letter is a method. Access to the external pointer should be through the method, not through the field.

See Also

lmer, glmer, nlmer, merMod.

Examples

showClass("lmResp")
str(lmResp$new(y=1:4))
showClass("glmResp")
str(glmResp$new(family=poisson(), y=1:4))
showClass("nlsResp")
showClass("lmerResp")
str(lmerResp$new(y=1:4))