Skip to contents

Fitting Nonlinear Mixed-Effects Models

nlmer.Rd

Fit a nonlinear mixed-effects model (NLMM) to data, via maximum likelihood.

Usage

nlmer(formula, data = NULL, control = nlmerControl(),
      start = NULL, verbose = 0L, nAGQ = 1L, subset, weights, na.action,
      offset, contrasts = NULL, devFunOnly = FALSE)

Arguments

formula

a three-part “nonlinear mixed model” formula, of the form resp ~ Nonlin(...) ~ fixed + random, where the third part is similar to the RHS formula of, e.g., lmer. Currently, the Nonlin(..) formula part must not only return a numeric vector, but also must have a "gradient" attribute, a matrix. The functions SSbiexp, SSlogis, etc, see selfStart, provide this (and more). Alternatively, you can use deriv() to automatically produce such functions or expressions.

data

an optional data frame containing the variables named in formula. By default the variables are taken from the environment from which lmer is called. 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).

control

a list (of correct class, resulting from lmerControl() or glmerControl() respectively) containing control parameters, including the nonlinear optimizer to be used and parameters to be passed through to the nonlinear optimizer, see the *lmerControl documentation for details.

start

starting estimates for the nonlinear model parameters, as a named numeric vector or as a list with components

nlpars

required numeric vector of starting values for the nonlinear model parameters

theta

optional numeric vector of starting values for the covariance parameters

verbose

integer scalar. If > 0 verbose output is generated during the optimization of the parameter estimates. If > 1 verbose output is generated during the individual PIRLS steps (PIRLS aka PRSS, e.g. in the C++ sources).

nAGQ

integer scalar - the number of points per axis for evaluating the adaptive Gauss-Hermite approximation to the log-likelihood. Defaults to 1, corresponding to the Laplace approximation. Values greater than 1 produce greater accuracy in the evaluation of the log-likelihood at the expense of speed. A value of zero uses a faster but less exact form of parameter estimation for GLMMs by optimizing the random effects and the fixed-effects coefficients in the penalized iteratively reweighted least squares (PIRLS) step.

subset

an optional expression indicating the subset of the rows of data that should be used in the fit. This can be a logical vector, or a numeric vector indicating which observation numbers are to be included, or a character vector of the row names to be included. All observations are included by default.

weights

an optional vector of ‘prior weights’ to be used in the fitting process. Should be NULL or a numeric vector.

na.action

a function that indicates what should happen when the data contain NAs. The default action (na.omit, inherited from the ‘factory fresh’ value of getOption("na.action")) strips any observations with any missing values in any variables.

offset

this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. One or more offset terms can be included in the formula instead or as well, and if more than one is specified their sum is used. See model.offset.

contrasts

an optional list. See the contrasts.arg of model.matrix.default.

devFunOnly

logical - return only the deviance evaluation function. Note that because the deviance function operates on variables stored in its environment, it may not return exactly the same values on subsequent calls (but the results should always be within machine tolerance).

Details

Fit nonlinear mixed-effects models, such as those used in population pharmacokinetics.

Note

Adaptive Gauss-Hermite quadrature (nAGQ > 1) is not currently implemented for nlmer. Several other methods, such as simulation or prediction with new data, are unimplemented or very lightly tested.

A method argument was used in earlier versions of the lme4 package. Its functionality has been replaced by the nAGQ argument.

Examples

## nonlinear mixed models --- 3-part formulas ---
## 1. basic nonlinear fit. Use stats::SSlogis for its
## implementation of the 3-parameter logistic curve.
## "SS" stands for "self-starting logistic", but the
## "self-starting" part is not currently used by nlmer ... 'start' is
## necessary
startvec <- c(Asym = 200, xmid = 725, scal = 350)
(nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
             Orange, start = startvec))
#> Nonlinear mixed model fit by maximum likelihood  ['nlmerMod']
#> Formula: circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym | Tree
#>    Data: Orange
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#>  273.1438  280.9205 -131.5719  263.1438        30 
#> Random effects:
#>  Groups   Name Std.Dev.
#>  Tree     Asym 31.646  
#>  Residual       7.843  
#> Number of obs: 35, groups:  Tree, 5
#> Fixed Effects:
#>  Asym   xmid   scal  
#> 192.1  727.9  348.1  
## 2. re-run with "quick and dirty" PIRLS step
(nm1a <- update(nm1, nAGQ = 0L))
#> Nonlinear mixed model fit by maximum likelihood  ['nlmerMod']
#> Formula: circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym | Tree
#>    Data: Orange
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#>  273.1689  280.9456 -131.5844  263.1689        30 
#> Random effects:
#>  Groups   Name Std.Dev.
#>  Tree     Asym 31.63   
#>  Residual       7.84   
#> Number of obs: 35, groups:  Tree, 5
#> Fixed Effects:
#>  Asym   xmid   scal  
#> 191.1  722.6  344.2  

## 3. Fit the same model with a user-built function:
## a. Define formula
nform <- ~Asym/(1+exp((xmid-input)/scal))
## b. Use deriv() to construct function:
nfun <- deriv(nform,namevec=c("Asym","xmid","scal"),
              function.arg=c("input","Asym","xmid","scal"))
nm1b <- update(nm1,circumference ~ nfun(age, Asym, xmid, scal)  ~ Asym | Tree)

## 4. User-built function without using deriv():
##    derivatives could be computed more efficiently
##    by pre-computing components, but these are essentially
##    the gradients as one would derive them by hand
nfun2 <- function(input, Asym, xmid, scal) {
    value <- Asym/(1+exp((xmid-input)/scal))
    grad <- cbind(Asym=1/(1+exp((xmid-input)/scal)),
              xmid=-Asym/(1+exp((xmid-input)/scal))^2*1/scal*
                    exp((xmid-input)/scal),
              scal=-Asym/(1+exp((xmid-input)/scal))^2*
                     -(xmid-input)/scal^2*exp((xmid-input)/scal))
    attr(value,"gradient") <- grad
    value
}
stopifnot(all.equal(attr(nfun(2,1,3,4),"gradient"),
                    attr(nfun(2,1,3,4),"gradient")))
nm1c <- update(nm1,circumference ~ nfun2(age, Asym, xmid, scal)  ~ Asym | Tree)