# Modular Functions for Mixed Model Fits Modular functions for mixed model fits ```r glFormula(formula, data = NULL, family = gaussian, subset, weights, na.action, offset, contrasts = NULL, start, mustart, etastart, control = glmerControl(), ...) ``` ## Arguments - `formula`: a two-sided linear formula object describing both the fixed-effects and random-effects parts of the model, with the response on the left of a `~` operator and the terms, separated by `+` operators, on the right. Random-effects terms are distinguished by vertical bars (`"|"`) separating expressions for design matrices from grouping factors. - `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). - `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 `NA`s. 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`. - `control`: a list giving - **for `[g]lFormula`:**: all options for running the model, see `lmerControl`; - **for `mkLmerDevfun,mkGlmerDevfun`:**: options for the inner optimization step; - **for `optimizeLmer` and `optimizeGlmer`:**: control parameters for nonlinear optimizer (typically inherited from the `...` argument to `lmerControl`). - `start`: starting values (see `lmer`; for `glFormula`, should be just a numeric vector of fixed-effect coefficients) - `family`: a GLM family; see `glm` and `family`. - `mustart`: optional starting values on the scale of the conditional mean; see `glm` for details. - `etastart`: optional starting values on the scale of the unbounded predictor; see `glm` for details. - ``...``: other potential arguments; for `optimizeLmer` and `optimizeGlmer`, these are passed to internal function `optwrap`, which has relevant parameters `calc.derivs` and `use.last.params` (see `lmerControl`). ## Returns `lFormula` and `glFormula` return a list containing components: - **fr**: model frame - **X**: fixed-effect design matrix - **reTrms**: list containing information on random effects structure: result of `mkReTrms` `mkLmerDevfun` and `mkGlmerDevfun` return a function to calculate deviance (or restricted deviance) as a function of the theta (random-effect) parameters. `updateGlmerDevfun` returns a function to calculate the deviance as a function of a concatenation of theta and beta (fixed-effect) parameters. These deviance functions have an environment containing objects required for their evaluation. CAUTION: The `environment` of functions returned by `mk(Gl|L)merDevfun` contains reference class objects (see `ReferenceClasses`, `merPredD-class`, `lmResp-class`), which behave in ways that may surprise many users. For example, if the output of `mk(Gl|L)merDevfun` is naively copied, then modifications to the original will also appear in the copy (and vice versa). To avoid this behavior one must make a deep copy (see `ReferenceClasses` for details). `optimizeLmer` and `optimizeGlmer` return the results of an optimization. ## Details These functions make up the internal components of an [gn]lmer fit. * `[g]lFormula` takes the arguments that would normally be passed to `[g]lmer`, checking for errors and processing the formula and data input to create a list of objects required to fit a mixed model. * `mk(Gl|L)merDevfun` takes the output of the previous step (minus the `formula` component) and creates a deviance function * `optimize(Gl|L)mer` takes a deviance function and optimizes over `theta` (or over `theta` and `beta`, if `stage` is set to 2 for `optimizeGlmer` * `updateGlmerDevfun` takes the first stage of a GLMM optimization (with `nAGQ=0`, optimizing over `theta` only) and produces a second-stage deviance function * `mkMerMod` takes the **environment** of a deviance function, the results of an optimization, a list of random-effect terms, a model frame, and a model all and produces a `[g]lmerMod` object. ## Examples ```r library(lme4) ### Fitting a linear mixed model in 4 modularized steps ## 1. Parse the data and formula: lmod <- lFormula(Reaction ~ Days + (Days|Subject), sleepstudy) names(lmod) ## 2. Create the deviance function to be optimized: (devfun <- do.call(mkLmerDevfun, lmod)) ls(environment(devfun)) # the environment of 'devfun' contains objects # required for its evaluation ## 3. Optimize the deviance function: opt <- optimizeLmer(devfun) opt[1:3] ## 4. Package up the results: mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr) ### Same model in one line lmer(Reaction ~ Days + (Days|Subject), sleepstudy) ### Fitting a generalized linear mixed model in six modularized steps ## 1. Parse the data and formula: glmod <- glFormula(cbind(incidence, size - incidence) ~ period + (1 | herd), data = cbpp, family = binomial) #.... see what've got : str(glmod, max=1, give.attr=FALSE) ## 2. Create the deviance function for optimizing over theta: (devfun <- do.call(mkGlmerDevfun, glmod)) ls(environment(devfun)) # the environment of devfun contains lots of info ## 3. Optimize over theta using a rough approximation (i.e. nAGQ = 0): (opt <- optimizeGlmer(devfun)) ## 4. Update the deviance function for optimizing over theta and beta: (devfun <- updateGlmerDevfun(devfun, glmod$reTrms)) ## 5. Optimize over theta and beta: opt <- optimizeGlmer(devfun, stage=2) str(opt, max=1) # seeing what we'got ## 6. Package up the results: (fMod <- mkMerMod(environment(devfun), opt, glmod$reTrms, fr = glmod$fr)) ### Same model in one line fM <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), data = cbpp, family = binomial) all.equal(fMod, fM, check.attributes=FALSE, tolerance = 1e-12) # ---- -- even tolerance = 0 may work ```

modular function