Ting WangBenjamin GravesYves RosseelEdgar C. Merkle American Board of Family
MedicineUniversity of MissouriGhent UniversityUniversity of Missouri

# Computation and application of generalized linear mixed model derivatives
using lme4

###### Abstract

Maximum likelihood estimation of generalized linear mixed models (GLMMs) is
difficult due to marginalization of the random effects. Derivative
computations of a fitted GLMM’s likelihood is also difficult, especially
because the derivatives are not by-products of popular estimation algorithms.
In this paper, we first describe theoretical results related to GLMM
derivatives along with a quadrature method to efficiently compute the
derivatives, focusing on fitted lme4 models with a single clustering variable.
We describe how psychometric results related to item response models are
helpful for obtaining the derivatives, as well as for verifying the
derivatives’ accuracies. We then provide a tutorial on the many possible uses
of these derivatives, including robust standard errors, score tests of fixed
effect parameters, and likelihood ratio tests of non-nested models. The
derivative computation methods and applications described in the paper are all
available in easily-obtained R packages.

This work was supported by NSF grant 1460719. Correspondence to Edgar Merkle,
University of Missouri. Email<EMAIL_ADDRESS>
##  1 Introduction

Maximum likelihood estimation of generalized linear mixed models [<]GLMMs;
e.g.,¿stroup12 is notoriously complicated due to the fact that random effects
are integrated out of the model likelihood. In general, the integrals cannot
be solved analytically, which means that we must use numerical methods to
approximate the integrals. Along with model estimation, these issues make it
difficult to apply other statistical methods to estimated GLMMs, because the
required pieces of the estimated model are not generally available. For
example, consider the computation of “robust” (Huber-White) standard errors
[<]e.g.,¿white80, huber67, as applied to GLMM. In addition to the model’s
maximum likelihood estimates, we require first and second partial derivatives
of the model’s likelihood function. These derivatives also require integral
approximations, which do not necessarily arise as by-products of the model
estimation algorithm.

Of primary importance for this paper, the partial derivatives do not arise as
by-products of model estimation via the _lme4_ package [Bates . (2015)]. This
package uses a penalized, iteratively re-weighted least squares (PIRLS)
algorithm that indirectly maximizes the marginal likelihood by optimizing a
second function that involves conditional random effects [<]conditional on
random effect (co-)variances;¿bates2021. Although this conditional approach
bypasses the difficult integration, it also loses the ability to produce the
likelihood derivatives of interest. This makes it difficult to apply many
relevant methods that are already implemented within the R ecosystem,
including sandwich estimators from package _sandwich_ [Zeileis (2004), Zeileis
(2006), Zeileis . (2020)], score-based tests from _strucchange_ [Zeileis .
(2002)], model-based recursive partitioning from _partykit_ [Hothorn Zeileis
(2015)], and Vuong tests from _nonnest2_ [Merkle You (2018)]. These packages
all rely on partial derivatives of the model likelihood function (evaluated at
the maximum likelihood estimates, after model estimation), which to date have
not been available for GLMMs estimated by _lme4_. So the overall goal of this
work is to connect existing statistical methods with GLMMs estimated by
_lme4_. The paper’s contributions towards this goal include (i) theoretical
background on GLMM derivatives, as well as a quadrature method that
capitalizes on the fact that we are dealing with estimated models; (ii) a
general-purpose implementation of the methods via the _merDeriv_ package; and
(iii) a tutorial on how these derivatives can be used in applied research
settings, including a variety of R examples.

Our derivations are informed by previous results from both statistics and
psychometrics, which include diverse motivations for the GLMM. In particular,
the statistics community often views the GLMM as an extension of the linear
mixed model, whereas the psychometrics community additionally considers
connections between the GLMM and item response theory (IRT) models
[<]e.g.,¿de2011,dorbat07. The latter connections are seldom noticed in the
statistics literature, though skrrab04 is noteworthy in that LMMs, GLMMs, and
IRT models are included within a larger latent variable framework. We describe
below how IRT results can help us obtain derivatives of the GLMM likelihood
function with respect to both fixed parameters and random effect
hyperparameters (e.g., random effect variances) after model estimation.

In the following sections, we first fix notation and define the GLMM. We then
present theoretical results related to derivatives of the GLMM likelihood
function, including a quadrature method that can be applied to estimated
models. Next, we provide a tutorial on the application of these results to
GLMMs estimated via _lme4_. This is accomplished with the help of R package
_merDeriv_ [Wang Merkle (2018)], which implements the methods described here,
combined with other packages like _mirt_ [Chalmers (2012)], _sandwich_ ,
_nonnest2_ , and _strucchange_. Finally, we discuss potential future
extensions of our work.

##  2 Theoretical Background

Our presentation of the GLMM follows the _lme4_ framework of lme4, which
facilitates the R applications presented later. This framework encompasses a
variety of GLMMs from the exponential family, with binomial models being
especially popular. The framework does not allow for products between free
parameters and random effects, which becomes important when we discuss
relationships between GLMMs and IRT models below [<]also see¿de2011,dorbat07.

###  2.1 Model and Notation

Let $\bm{y}_{i}$ be a vector containing the response variable for the $i$th
cluster, each entry of which is assumed to follow a specific probability
distribution (e.g., binomial or Poisson). The sample size of cluster $i$ is
denoted as $n_{i}$, so the total sample size across all $I$ clusters is given
as $N=\sum_{i=1}^{I}n_{i}$. Let $\bm{X}_{i}$ be the $n_{i}\times p$ design
matrix corresponding to fixed effects for cluster $i$; $\bm{\beta}$ is the
fixed effect vector of length $p$; $\bm{Z}_{i}$ is the $n_{i}\times q$ design
matrix corresponding to random effects for cluster $i$; and $\bm{u}_{i}$ is
the random effect vector of length $q$. Then the model can be written as

$\displaystyle E(\bm{y}_{i}|\bm{u}_{i},\bm{\Lambda}_{\bm{\theta}})$
$\displaystyle=\bm{\mu}_{i}|\bm{\Lambda}_{\bm{\theta}},\bm{u}_{i}$ (1)
$\displaystyle\bm{\mu}_{i}$
$\displaystyle=g^{-1}(\bm{\eta}_{i}|\bm{\Lambda}_{\bm{\theta}},\bm{u}_{i})$
(2) $\displaystyle\bm{\eta}_{i}$
$\displaystyle=\bm{X}_{i}\bm{\beta}+\bm{Z}_{i}\bm{b}_{i}$ (3)
$\displaystyle\bm{b}_{i}$ $\displaystyle=\bm{\Lambda}_{\bm{\theta}}\bm{u}_{i}$
(4) $\displaystyle\bm{u}_{i}$ $\displaystyle\sim N(\bm{0},\bm{I}_{q}).$ (5)

The above equations express the idea that the bounded support of the expected
value of $\bm{y}_{i}$ can be transformed to an unbounded support of the linear
combination $\bm{X}_{i}\bm{\beta}+\bm{Z}_{i}\bm{b}_{i}$ through the link
function $g()$. The random effects are in $\bm{b}_{i}$, which equals
$\bm{\Lambda}_{\bm{\theta}}\bm{u}_{i}$. The vector $\bm{u}_{i}$ follows the
standard normal distribution $N(\bm{0},\bm{I}_{q})$, with
$\bm{\Lambda}_{\bm{\theta}}$ being the relative covariance factor, which can
be seen as the Cholesky decomposition of the usual random effect covariance
matrix $\bm{G}$. Reparameterizing $\bm{b}_{i}$ as the product of the relative
covariance factor and standard normal distribution makes it easier to compare
GLMM to IRT. We provide further discussion of this comparison in the next
section.

Following the above notation, the model’s log-likelihood (marginal over random
effects) can be expressed as

$\ell=\displaystyle\sum_{i=1}^{I}\ell_{i}=\displaystyle\sum_{i=1}^{I}\log\int
f_{\bm{y}_{i}|\bm{u}_{i}}(\bm{y}_{i}|\bm{u}_{i})f_{\bm{u}_{i}}(\bm{u}_{i})d\bm{u}_{i},$
(6)

where $I$ represents the number of total clusters. We further define the
following “across-cluster” matrices:

$\displaystyle\bm{y}$ $\displaystyle=$
$\displaystyle\\{\bm{y}_{1},\bm{y}_{2},\ldots,\bm{y}_{i},\ldots,\bm{y}_{I}\\}$
(7) $\displaystyle\bm{X}$ $\displaystyle=$
$\displaystyle\\{\bm{X}_{1},\bm{X}_{2},\ldots,\bm{X}_{i},\ldots,\bm{X}_{I}\\}$
(8) $\displaystyle\bm{Z}$ $\displaystyle=$
$\displaystyle\\{\bm{Z}_{1},\bm{Z}_{2},\ldots,\bm{Z}_{i},\ldots,\bm{Z}_{I}\\}$
(9) $\displaystyle\bm{b}$ $\displaystyle=$
$\displaystyle\\{\bm{b}_{1},\bm{b}_{2},\ldots,\bm{b}_{i},\ldots,\bm{b}_{I}\\}.$
(10)

###  2.2 GLMM scores

One of the most popular IRT models is the two-parameter logistic model
[<]e.g.,¿embrei00,lornov68, which can be viewed as a binomial GLMM with logit
link function. Consider an IRT model parameterized as
$\text{logit}^{-1}(p_{ij})=\alpha_{j}\theta_{i}-\beta_{j}$, with each item
$j$’s difficulty described by $\beta_{j}$ and discrimination described by
$\alpha_{j}$. The alternative parameterization as
$\alpha_{j}(\theta_{i}-\beta_{j})$ is also applicable, but less convenient for
comparison. In the former parameterization, the IRT $\beta_{j}$ parameters are
similar to the negative of the GLMM fixed parameter $\bm{\beta}$. The IRT
$\alpha_{j}$ parameters are then similar to the relative covariance factor in
the GLMM, with the _lme4_ package requiring the covariance factor to be equal
for all items. This means that we cannot fit a 2PL model in _lme4_ , though
other GLMM software such as SAS PROC NLMIXED may allow for 2PL estimation.

In the context of IRT, glas92,glas98,glas99 utilized an identity from lou82 to
obtain first derivatives of the marginal log-likelihood (marginal over person
parameters $\theta_{i}$). This identity can be used to show that the first
derivative of the marginal log-likelihood with respect to difficulty and
discrimination parameters equals an expected value involving first derivatives
of the conditional likelihood (conditioned on person proficiency). That is, we
can obtain derivatives of the marginal likelihood by taking an expected value
that involves the conditional likelihood.

The same idea can be applied to GLMM [McCulloch Neuhaus (2001)], where
conditioning on person proficiency is replaced with conditioning on random
effects. In the next sections, we will formalize these GLMM score derivations.
Please note that, throughout this paper, _scores_ refer to first derivatives
of the clusterwise log-likelihood function with respect to some model
parameters. They are different from _factor scores_ and from _scoring_ in
psychometrics, which involve prediction of a model’s random parameters.

####  2.2.1 Fixed effect scores

Drawing on derivations by Glas as well as by mcc01, the GLMM score with
respect to the fixed effect parameter $\bm{\beta}$ can be expressed in the
following form:

$\frac{\partial\ell_{i}}{\partial\bm{\beta}}=\frac{\displaystyle\int\frac{\partial\log
f_{\bm{y}_{i}|\bm{u}_{i}}(\bm{y}_{i}|\bm{u}_{i})}{\partial\bm{\beta}}f_{\bm{y}_{i}|\bm{u}_{i}}(\bm{y}_{i}|\bm{u}_{i})f_{\bm{u}_{i}}(\bm{u}_{i})d\bm{u}_{i}}{f_{\bm{y}_{i}}(\bm{y}_{i})},$
(11)

where $f_{\bm{y}_{i}}(\bm{y}_{i})=\displaystyle\int
f_{\bm{y}_{i}|\bm{u}_{i}}(\bm{y}_{i}|\bm{u}_{i})f_{\bm{u}_{i}}(\bm{u}_{i})d\bm{u}_{i}$.

The first term in the numerator of Equation (11) can be seen as the score of a
Generalized Linear Model (GLM), which can be expressed in matrix form as

$\frac{\partial\log
f_{\bm{y}_{i}|\bm{u}_{i}}(\bm{y}_{i}|\bm{u}_{i})}{\partial\bm{\beta}}=\bm{X}_{i}^{T}\bm{D}_{i}^{-1}\bm{V}_{i}^{-1}(\bm{y}_{i}-\bm{\mu}_{i}),$
(12)

where $\bm{D}_{i}$ and $\bm{V}_{i}$ are $n_{i}\times n_{i}$ diagonal matrices
with diagonal entries as
$\frac{\partial(\eta_{t}|u_{t})}{\partial(\mu_{t}|u_{t})}$ and
$a(\phi_{t})\text{Var}(\mu_{t}|u_{t})$, respectively. The $t$ subscript
indexes an observation within cluster $i$, $1,2,...,n_{i}$. Further, the
$a(\phi_{t})$ function is unique to each distribution from the exponential
family. For example, $a(\phi_{t})=1$ for the binomial distribution and for the
Poisson distribution. The value of $a(\phi_{t})$ for other exponential family
distributions can be found in, e.g., mccnel89. Many of the relevant
derivations are also supplied by the R family() function. Note that, if we use
the canonical link function,
$\frac{\partial(\eta_{t}|u_{t})}{\partial(\mu_{t}|u_{t})}$ and
$\text{Var}(\mu_{t}|u_{t})$ will cancel out. This feature creates a shortcut
for distributions using the canonical link.

The second term in the numerator of Equation (11) is the distribution of the
GLM given $\bm{u}_{i}$. We use the following matrix form to express all
distributions belonging to the exponential family:

$f_{\bm{y}_{i}|\bm{u}_{i}}(\bm{y}_{i}|\bm{u}_{i})=\text{exp}\left(\bm{y}_{i}^{T}\bm{A}_{i}\bm{\kappa}_{i}-\bm{1}^{T}\bm{A}_{i}h(\bm{\kappa}_{i})+c(\bm{y}_{i},\bm{\psi}_{i})\right),$
(13)

where $\bm{A}_{i}$ is a $n_{i}\times n_{i}$ diagonal matrix with diagonal
element as $\frac{1}{a(\phi_{t})}$; $\bm{\kappa}_{i}$ is the vector of
canonical parameters; $\bm{1}$ is a $n_{i}\times 1$ vector with each entry as
$1$; $h(\bm{\kappa}_{i})$ is an $n_{i}\times 1$ vector defined by applying the
distribution-specific function $h()$ to each element of $\bm{\kappa}_{i}$; and
$c(\bm{y}_{i},\bm{\psi}_{i})$ is an $n_{i}\times 1$ vector of remaining terms
not depending on $\bm{\kappa}_{i}$, with $\bm{\psi}_{i}$ containing scale
parameters. For exponential distributions, these terms can also be found in
mccnel89 or in the R family() functions.

The above results based on generalized linear models are straightforward,
while the difficulty involves the integration over $\bm{u}$. In the same
spirit, the denominator can be viewed as the integration of the GLM
distribution over the random variable $\bm{u}$. Both integrals have no closed
form for GLMMs. We discuss use of quadrature to approximate the integrals
below, after describing derivatives of random effect hyperparameters.

####  2.2.2 Random effect hyperparameter scores

Following the same type of derivation, the scores w.r.t. the random effect
hyperparameters can be seen as the scores w.r.t. parameters in the
$\bm{\Lambda}_{\bm{\theta}}$ matrix. The derivation can thus be expressed as:

$\frac{\partial\ell_{i}}{\partial\bm{\Lambda}_{\bm{\theta}}}=\frac{\displaystyle\int\frac{\partial\log
f_{\bm{y}_{i}|\bm{u}_{i}}(\bm{y}_{i}|\bm{u}_{i})}{\partial\bm{\Lambda}_{\bm{\theta}}}f_{\bm{y}_{i}|\bm{u}_{i}}(\bm{y}_{i}|\bm{u}_{i})f_{\bm{u}_{i}}(\bm{u}_{i})d\bm{u}_{i}}{f_{\bm{y}_{i}}(\bm{y}_{i})},$
(14)

where $\frac{\partial\log
f_{\bm{y}_{i}|\bm{u}_{i}}(\bm{y}_{i}|\bm{u}_{i})}{\partial\bm{\Lambda}_{\bm{\theta}}}$
equals
$\bm{u}_{i}^{T}\frac{\partial\bm{\Lambda}_{\bm{\theta}}}{\partial\theta}\bm{Z}_{i}^{T}(\bm{y}_{i}-\bm{\mu}_{i})$,
with $\frac{\partial\bm{\Lambda}_{\bm{\theta}}}{\partial\theta}$ as a matrix
composed of 1s (corresponding to a particular random effect hyperparameter
$\theta$) and 0s (not corresponding to a particular random effect
hyperparameter $\theta$). This derivation is similar to the score derivation
for the IRT discrimination parameter. An equivalent approach is to rearrange
terms using the trace operator [<]e.g.,¿peter08, which results in the
expression
$\text{Tr}\left((\bm{Z}_{i}^{T}(\bm{y}_{i}-\bm{\mu}_{i})\bm{u}_{i}^{T})^{T}\frac{\partial\bm{\Lambda}_{\bm{\theta}}}{\partial\theta}\right)$.

####  2.2.3 Reparameterization

As mentioned above, $\bm{\Lambda}_{\bm{\theta}}$ is a Cholesky decomposition
of the usual variance covariance matrix $\bm{G}$, so our derivations are taken
with respect to the Cholesky decomposition. In order to obtain the scores with
respect to the variance-covariance parameters contained in $\bm{G}$, we
utilize the chain rule:

$\displaystyle\frac{\partial\ell}{\partial\bm{G}}$ $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\bm{\Lambda}_{\bm{\theta}}}\frac{\partial\bm{\Lambda}_{\bm{\theta}}}{\partial\bm{G}}$
(15) $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\bm{\Lambda}_{\bm{\theta}}}\left\\{\frac{\partial\bm{\Lambda}_{\bm{\theta}}}{\partial(\bm{\Lambda}_{\bm{\theta}}\bm{\Lambda}_{\bm{\theta}}^{T})}\right\\}$
(16) $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\bm{\Lambda}_{\bm{\theta}}}\left\\{\frac{\partial(\bm{\Lambda}_{\bm{\theta}}\bm{\Lambda}_{\bm{\theta}}^{T})}{\partial\bm{\Lambda}_{\bm{\theta}}}\right\\}^{-1}.$
(17)

For the entry in row $i$ and column $j$ of $\bm{\Lambda}_{\bm{\theta}}$, we
have that

$\frac{\partial(\bm{\Lambda}_{\bm{\theta}}\bm{\Lambda}_{\bm{\theta}}^{T})}{\partial\bm{\Lambda}_{\bm{\theta}ij}}=\bm{\Lambda}_{\bm{\theta}}\bm{J}_{ji}+\bm{J}_{ij}\bm{\Lambda}_{\bm{\theta}}^{T},$
(18)

where $\bm{J}_{ij}$ is a matrix with entry $(i,j)$ equal to 1 and 0 elsewhere.
The derivatives with respect to all unique, nonzero entries of
$\bm{\Lambda}_{\bm{\theta}}$ can be computed in this manner to obtain the
desired scores.

As an alternative to variances and covariances, users may wish to parameterize
the model via standard deviations and correlations. The scores with respect to
standard deviations and correlations can be obtained by applying another chain
rule to the above scores that are taken with respect to $\bm{G}$. For example,
assume a GLMM with two correlated random effects. In the variance-covariance
parameterization, we would have parameters $\sigma_{0}^{2}$, $\sigma_{1}^{2}$,
and $\sigma_{01}$, while, in the standard deviation-correlation
parameterization, we would have parameters $\sigma_{0}$, $\sigma_{1}$, and
$\rho$. Derivatives for the latter parameterization are:

$\displaystyle\frac{\partial\ell}{\sigma_{0}}$ $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\sigma_{0}^{2}}\frac{\partial\sigma_{0}^{2}}{\partial\sigma_{0}}$
(19) $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\sigma_{0}^{2}}(2\sigma_{0})$ (20)
$\displaystyle\frac{\partial\ell}{\partial\sigma_{1}}$ $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\sigma_{1}^{2}}\frac{\partial\sigma_{1}^{2}}{\partial\sigma_{1}}$
(21) $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\sigma_{1}^{2}}(2\sigma_{1})$ (22)
$\displaystyle\frac{\partial\ell}{\partial\rho}$ $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\rho\sigma_{0}\sigma_{1}}\frac{\partial\rho\sigma_{0}\sigma_{1}}{\partial\rho}$
(23) $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\sigma_{01}}(\sigma_{0}\sigma_{1}).$
(24)

####  2.2.4 Quadrature

All the derivatives above involve integrals that marginalize over the model
random effects $\bm{u}$. These integrals do not have closed forms, requiring
numerical methods for approximation. The method implemented in R package
_merDeriv_ is a simplified version of multivariate adaptive Gauss-Hermite
quadrature [Liu Pierce (1994), Naylor Smith (1982)], with the simplifications
being based on the fact that we are computing derivatives after model
estimation. This means that we already have information about posterior modes
and variances of random effects from _lme4_ , and we can make use of this
information in place of the “adaptive” part of the algorithm. merfur19
recently used a similar method to compute marginal versions of Bayesian
information criteria (see especially their Appendix C), with that method being
based on earlier methods described by pinbat95 and rab05. While it would be
possible to simply use a traditional adaptive quadrature method here, we would
have to use it separately for each case in the data (because we seek to
compute casewise derivatives). This would be much slower and infeasible for
many datasets, as compared to our quadrature method described here.

Focusing on the GLMM framework, the integrals from Equations (11) and (14) are
both of the form

$\displaystyle\int
g(\bm{y}|\bm{u},\bm{\omega})f_{\bm{y}|\bm{u},\bm{\omega}}(\bm{y}|\bm{u},\bm{\omega})f_{\bm{u}|\bm{\omega}}(\bm{u}|\bm{\omega})d\bm{u},$
(25)

where $g()$ differs depending on the integral, and $\bm{\omega}$ is a vector
of model parameters excluding the random effects $\bm{u}$. This conditioning
on $\bm{\omega}$ is implicit in earlier sections but was excluded to simplify
notation.

For a single clustering variable with $I$ levels, the clusters $i$ are
independent. Therefore, the above equation can be written as

$\displaystyle\prod_{i=1}^{I}\displaystyle\int
g(\bm{y}_{i}|\bm{u}_{i},\bm{\omega})f_{\bm{y}_{i}|\bm{u}_{i},\bm{\omega}}(\bm{y}_{i}|\bm{u}_{i},\bm{\omega})f_{\bm{u}_{i}|\bm{\omega}}(\bm{u}_{i}|\bm{\omega})d\bm{u}_{i}.$
(26)

To compute scores, we are interested in the elements of the above product: the
integral for each cluster $i$. For $M$ quadrature points, we use Gauss-Hermite
quadrature to approximate the integral for cluster $i$ by:

$\displaystyle\sum_{m=1}^{M}w^{\ast}_{im}g(\bm{y}_{i}|\bm{a}^{\ast}_{im},\bm{\omega})f_{\bm{y}_{i}|\bm{u}_{i},\bm{\omega}}(\bm{y}_{i}|\bm{a}^{\ast}_{im},\bm{\omega}).$
(27)

That is, the integral is approximated by a weighted sum of function
evaluations, where the functions are evaluated at different random effect
values represented by $\bm{a}^{\ast}_{im}$, $m=1,\ldots,M$. For a random
effect of dimension $d$, the quadrature locations and weights are computed by

$\displaystyle\bm{a}^{\ast}_{im}$
$\displaystyle=\tilde{\bm{b}}_{i}+\tilde{\bm{C}}_{i}\times\bm{a}_{m}$ (28)
$\displaystyle w^{\ast}_{im}$
$\displaystyle=w_{m}\times(2\pi)^{d/2}\times\det{(\tilde{\bm{C}}_{i})}\times\exp{(0.5\times\bm{a}_{m}\bm{a}_{m}^{\prime})}\times\bm{\phi}(\bm{a}^{\ast}_{im}|\bm{0},\hat{\bm{G}})$
(29)

where $\tilde{\bm{b}}_{i}$ are the posterior modes of random effects for
cluster $i$, $\tilde{\bm{C}}_{i}$ is the Cholesky factor of the conditional
covariance matrix of the random effects for cluster $i$ (obtained from the
_lme4_ function ranef()), $\bm{\phi}()$ is the normal density function, and
$\hat{\bm{G}}$ is the estimated covariance matrix of the random effects
(obtained from the _lme4_ function VarCorr()). Finally, $\bm{a}_{m}$ and
$w_{m}$ are the usual Gauss-Hermite locations and weights, respectively.

####  2.2.5 Second derivatives

While we have focused on first derivatives, the lou82 identity can also aid in
computation of second derivatives, leading to the model Hessian and
information matrix. We do not present the equations here because, for models
estimated via glmer() (but not lmer()), a Hessian is already computed and
stored in the resulting model object (specifically in the optinfo slot).
According to the _lme4_ documentation, this Hessian is computed using a finite
difference approach. The _merDeriv_ package provides a convenience function to
access this Hessian, and we use it in our applications later.

By default, the _lme4_ Hessian is parameterized via the Cholesky decomposition
of random effects. The Hessian based on the standard deviation/correlation
parameterization can alternatively be obtained via the devfun2() function in
_lme4_ , which uses the profile likelihood. The Hessian for the
variance/covariance parameterization is then related to the latter option,
through the chain rule mentioned earlier. The _merDeriv_ package incorporates
these computations and enables researchers to request the parameterization of
interest via the ranpar argument (taking possible values of "var", "sd", or
"theta"). This is illustrated in the tutorials below.

##  3 Tutorial on the Derivatives’ Uses in R

We now provide a tutorial on R package _merDeriv_ , which can carry out the
computations described above and which can be used to solve applied problems.
As we go, we provide snippets of code that illustrate how _merDeriv_ interacts
with other packages, which readers can adapt to other models and datasets. We
first provide some evidence that _merDeriv_ operates in the manner expected,
by comparing a Rasch model estimated via _lme4_ to a Rasch model estimate via
_mirt_ [Chalmers (2012)]. We then consider a variety of other applications.

###  3.1 Verifying the Computations

Before using the scores from _merDeriv_ in GLMM applications, we use the
relationship between GLMM and IRT to verify the correctness of the quadrature
implementation. We specifically compare the score computations to those of
package _mirt_ [Chalmers (2012)], which estimates many types of item response
models. We make use of the fact that the Rasch model can be estimated as a
generalized linear mixed model, which was illustrated by de2011. We also make
use of the fact that _mirt_ has its own, independent quadrature method for
score computation, which was used by schcha19 to apply Vuong tests to item
response models.

####  3.1.1 Method

For comparing the two score computation algorithms, we use the LSAT7 data
[Bock Lieberman (1970)] included with _mirt_. This dataset includes the item
responses (correct/incorrect) of 1,000 individuals across 5 items of the LSAT.

The code in Figure 1 shows how a Rasch model can be fit to the data using both
_mirt_ and _lme4_. For _mirt_ , we require the LSAT7 data to be arranged in
wide format, where each row is a person and each column is an item. If we then
rearrange the data to be in long format, as shown in Figure 1, we can fit the
Rasch model via _lme4_. We use the nAGQ argument to employ adaptive quadrature
during _lme4_ model estimation, avoiding the glmer() default, nAGQ=1, which
uses the Laplace approximation. The quadrature leads to a more accurate
approximation of the model log-likelihood, which in turn leads to maximum
likelihood estimates that tend to be closer to the true maximum of the
likelihood. The _mirt_ package employs a fixed quadrature method with 61
quadrature points.

Figure 1: Code to fit Rasch models using _mirt_ and _lme4_ , then calculate
scores.

    ## mirt:
    library("mirt")
    ls7 <- expand.table(LSAT7)
    mirtmod <- mirt(ls7[,1:5], 1, itemtype = "Rasch", SE = TRUE)

    ## reshape data and fit with glmer():
    library("reshape2")
    ls7$person <- 1:nrow(ls7)
    ls7long <- melt(ls7, id = "person")
    lme4mod <- glmer(value ~ -1 + variable + (1 | person), family = binomial,
                     data = ls7long, nAGQ = 5L)

    ## score calculation:
    mirtsc <- estfun.AllModelClass(mirtmod)
    lme4sc <- estfun.glmerMod(lme4mod, ranpar = "var")

####  3.1.2 Results

As shown at the bottom of Figure 1, scores for the two models are obtained via
their respective estfun() functions. The function for _mirt_ models is
included directly within the _mirt_ package, whereas the function for _lme4_
models is included in _merDeriv_. Both functions output a score matrix, where
rows index people and columns index model parameters. For the glmer model, we
use the ranpar argument so that the _merDeriv_ scores involve the variance-
covariance parameterization, which matches the _mirt_ output.

In comparing the two sets of scores, we arrive at Figure 2. The x-axis depicts
scores from _merDeriv_ , the y-axis depicts scores from _mirt_ , and each
point is a particular score. We see that the values are nearly exactly equal
for _mirt_ and for _merDeriv_ , falling directly on the identity line. One can
also compare the parameter variance-covariance matrix of _merDeriv_ and of
_mirt_ , using the vcov() method of each package. That comparison, not shown,
exhibits agreement similar to the score comparison. These provide evidence
that the _merDeriv_ code is performing as expected. Now that we have obtained
this evidence, we move on to illustrate practical uses of the scores in GLMM
applications.

Figure 2: Comparison of Rasch model scores from _mirt_ and from _merDeriv_. | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ | $\beta_{5}$ | $\sigma^{2}$
---|---|---|---|---|---|---
lme4 | 0.1004 | 0.0811 | 0.0913 | 0.0787 | 0.1037 | 0.1300
sandwich | 0.0996 | 0.0814 | 0.0898 | 0.0785 | 0.1058 | 0.1311
Table 1: Comparison between Rasch model standard errors reported by _lme4_ and
robust standard errors reported by _sandwich_. The $\beta$ columns correspond
to item difficulties, while the $\sigma^{2}$ column corresponds to person
(intercept) variance.

###  3.2 Huber-White estimator

Let $\bm{\omega}$ be the model parameter vector, which in a GLMM would contain
fixed effect parameters and random effect (co-)variances. Then the Huber-White
[<]e.g.,¿white80, huber67 sandwich estimator of the covariance matrix of
$\bm{\omega}$ is

$\bm{V}(\hat{\bm{\omega}})=(\bm{A})^{-1}\bm{B}(\bm{A})^{-1},$ (30)

where $\bm{A}$ is the negative expectation of the model Hessian and $\bm{B}$
is the covariance matrix of scores [<]see¿[for further discussion in the
context of linear mixed models]merDeriv. The score computations described in
the previous sections facilitate computation of this $\bm{B}$ matrix. The
square root of the diagonal elements of $\bm{V}$ are then typically called
“robust standard errors.”

Robust standard errors are used to address model misspecifications such as
unmodeled dependence between observations or deviations from normality. While
random effects are typically used in GLMMs to account for dependence between
observations, the Huber-White estimator can be used on top of a GLMM to
account for further model misspecifications. Further, str20 recently provided
evidence that quadrature can lead to downward-biased variance estimates in
GLMMs, resulting in inflated Type I error rates. The Huber-White estimator may
be considered in light of this result.

We can easily compute Huber-White standard errors using the scores from the
previous section, paired with the _sandwich_ package, as shown in Figure 3. In
that figure, the bread.glmerMod() and meat() functions come from _merDeriv_ ,
while sandwich() comes from the _sandwich_ package. Applying this result to
the Rasch model estimated in the previous section, we obtain the results in
Table 1. For this particular application, the _lme4_ standard errors and
_sandwich_ standard errors are virtually equal, likely due to the large sample
size (large by GLMM standards, at least).

Figure 3: Example code for calculating Huber-White standard errors.

    library("sandwich")
    sandwich(lme4mod, bread. = bread.glmerMod, meat. = meat(lme4mod, level = 2))

###  3.3 Score tests

Researchers have long been familiar with score tests, also known as Lagrange
multiplier tests, that can be used as an alternatives to the likelihood ratio
test or to the Wald test [<]e.g.,¿eng84,glas92,glas98,glas99. In typical score
test applications, a constrained model is fit to data, then first derivatives
of the likelihood function are used to test whether or not some constraint
should be relaxed. In contrast, the likelihood ratio test requires us to
estimate two models (a constrained model and an unconstrained model), and the
Wald test requires us to estimate only the unconstrained model.

This score test framework has expanded to a class of “parameter instability”
tests, where we test whether an estimated model’s parameters differ with
respect to unmodeled auxiliary variables (with different test statistics being
used for continuous, ordinal, or discrete auxiliary variables). zeihor07
summarized much previous work on this topic, developing a family of score-
based tests that can be used within an M-estimation framework (of which
maximum likelihood estimation is a special case). They also developed R
package _strucchange_ [Zeileis . (2002)], which can be used to compute the
test statistics so long as a model’s scores and Hessian are available. The
family of score-based tests has subsequently been studied in the context of
many specific types of models, including linear mixed models [Wang . (2020)],
structural equation models [Merkle Zeileis (2013), Merkle . (2014)], and item
response models [Komboz . (2018), Strobl . (2015), Wang . (2018)]. The
developments in the current paper make it possible to apply score-based tests
to GLMMs, yielding test statistics for GLMMs that have been unavailable up to
now. The score computations described above can be used to construct the
cumulative scores, which are further used to compute test statistics [Merkle
Zeileis (2013), Merkle . (2014)].

In this section, we show how scores can be used to test fixed effect
parameters that are not directly included in a GLMM model. This is potentially
useful in situations where a model with the fixed effect included does not
converge, which often happens in applied mixed modeling
[<]see¿barlev13,matkli17. In these situations, if we can get a model to
converge without some fixed effect of interest, it is possible to apply score-
based tests to the fitted model in order to test the omitted fixed effect.
While the more popular approach here is to drop random effects (as opposed to
fixed effects) from the model, dropping fixed effects may be useful in
instances where, e.g., the random effect variances are all large, yet the
model still exhibits convergence problems.

####  3.3.1 Method

We use data from 500 respondents on the Nerdy Personality Attributes Scale
(NPAS), a personality test designed for personal entertainment on the Open
Source Psychometrics Project website [Open Source Psychometrics Project ()].
The questionnaire consists of 26 items that attempt to define the concept of
“nerdiness”. Responses were originally measured on 5-point Likert scales, but
we converted them to binary responses for this example (where 0 corresponds to
3 or less and 1 corresponds to 4 or 5). The items ask about different aspects
of nerdiness, including hobbies and interests that are usually associated with
nerds, social interactions, personality traits, and academic or intellectual
endeavors. The data also include various demographic variables and other
personality measures assessing the “Big Five” personality factors.

Here, we assess whether item responses vary across extraversion, while also
accounting for inherent item differences (which would be called “item
difficulties” in an IRT context). The _lme4_ syntax for this model is shown at
the top of Figure 4, where the variable names are generally self-explanatory.
Note that inclusion of the interaction term (item*ext) automatically includes
main effects of both item and extraversion, in addition to the interaction.
This GLMM can be viewed as a person-by-covariate item response model, falling
into the class of explanatory item response models considered by debwil04 and
de2011.

Figure 4: Models of the NPAS data. The first model has issues with non-
convergence, leading us to the simpler, second model. A score test is then
used to study the interaction.

    ## Model that has problems with convergence:
    m1 <- glmer(answer ~ -1 + item*ext + (1 | subject),
                data = npas.sampled, family = binomial)

    ## Model with only main effects, which converges:
    m2 <- glmer(answer ~ -1 + item + ext + (1 | subject),
                data = npas.sampled, family = binomial,
                control = glmerControl(optimizer=’bobyqa’))

    ## Score test:
    ext <- with(npas.sampled, as.numeric(tapply(ext, subject, head, 1)))

    sc1 <- sctest(m2, fit=NULL, scores=estfun.glmerMod, order.by=ext,
           parm=1:26, functional=’maxLMo’)

####  3.3.2 Results

The first model in Figure 4 did not converge, even after making changes to the
optimizer and its settings. We could have experimented further, perhaps
finding some combination of settings that led to a converged model and that
would render the score test unnecessary. But each attempted model estimation
took about ten minutes, so we could easily have spent hours tweaking the
settings. In contrast, the score test could be immediately applied to a
simpler model that converged more easily.

Our simpler model was the second model in Figure 4, which used the `bobyqa`
optimizer [Powell (2009)] instead of the default `Nelder_Mead`. In estimating
this second model, we capitalize on the fact that score tests require only a
“constrained” model, which here assumes that responses to items do not vary
across levels of extraversion. We can then obtain a score test statistic for
the interaction without directly including the interaction in the model.

Figure 5: M-fluctuation test for NPAS data. This graph presents item parameter
fluctuation across varying levels of extraversion. Peaks of the graph suggest
extraversion cutpoints that isolate individuals with similar item parameters.

To obtain a test statistic for this interaction, we carry out the score test
using the code at the bottom of Figure 4. This makes use of the sctest()
function found in the R package _strucchange_ as well as estfun.glmerMod()
function found in _merDeriv_. It simultaneously tests all 26 item parameters
for fluctuations with respect to extraversion, which is similar to including
an item$\times$extraversion interaction. Results of this test are visualized
in Figure 5, which shows how the scores fluctuate across different values of
extraversion (x-axis). We can see that there is significant parameter
fluctuation in items across levels of extraversion, because the black line
goes above the “critical value” that is depicted by the red line. The peaks in
the black line, around extraversion values of 2.5 and 4.5, suggest cutpoints
for subgroups of individuals that exhibit similar item parameters. The test
provides information about the nature of the interaction that was not easily
obtained by including extraversion in the model, due to model convergence
problems.

###  3.4 Vuong tests

Scores also play a role in Vuong tests [Vuong (1989)], which can be used to
compare nested and non-nested models to one another. In the nested case, the
tests can be viewed as extensions of the traditional likelihood ratio test,
which, unlike the traditional likelihood ratio test, make no assumptions about
the more complex model being correct. In the non-nested case, the tests
provide a formal way of comparing the fits of the two models. The scores
described in this paper can be used in tandem with package _nonnest2_ [Merkle
You (2018)] to compare GLMMs, providing new capabilities for comparing models
with different predictor variables and different random effects. Specifically,
our score computations are used to compute the null distribution of the test
statistic, which is a weighted sum of chi-square distributions. Further
descriptions of the tests and applications to psychometric models can be found
in merkle16 and in schcha19. An illustration involving GLMMs is provided here.

####  3.4.1 Method

The data used for this example comes from the SPISA data set, which can be
found in the R package _psychotree_ [Strobl . (2015)]. The data is a subsample
of 1,075 Bavarian university students who took an online, general knowledge
quiz called “Studentenpisa” administered by a German weekly news magazine
[Trepte Verbeet (2010)]. The quiz consists of 45 items on 5 topics, and we
focus here on a subset of nine questions dealing with natural science. The
data set includes several covariates such as age, gender, semester of
university enrollment, and elite university status.

Using a similar item response model as in the previous example, we construct
two non-nested models with different covariates. These models are based on a
common reduced model that only contains item and subject effects. The first
model uses age and gender as covariates, while the second model uses semester
of university enrollment and whether the student’s university has been granted
“elite” status or not. The code for these models is shown in Figure 6. Similar
to the previous application, the models here did not immediately converge, and
we switched optimizers in order to attain convergence. Following model
estimation, we obtained scores and compared the two models using a Vuong test
computed via the R package _nonnest2_ [Merkle You (2018)].

Figure 6: Code for the non-nested models to be compared using the Vuong test.
The first model uses age and gender as potential predictors, while the second
model uses number of semesters at the university and elite university status.

    mod1 <- glmer(response ~ -1 + item + agecent + gender + (1 | pnum),
                  data = spisa, family = binomial,
                  control = glmerControl(optimizer=’bobyqa’))

    mod2 <- glmer(response ~ -1 + item + semester + elite + (1 | pnum),
                  data = spisa, family = binomial,
                  control = glmerControl(optimizer=’bobyqa’))

####  3.4.2 Results

The _nonnest2_ code and output for the Vuong test is shown in Figure 7. First,
we create a convenience function, vcg(), to compute the full parameter
covariance matrix (including random effect variances/covariances) for each of
the models. This function, along with functions from _merDeriv_ for
calculating the likelihoods and scores, is then sent to vuongtest().

Figure 7: Code to run Vuong test for comparing two non-nested models. The
models are able to be distinguished from each other, but one model does not
have better fit over the other.

    vcg <- function(obj) vcov(obj, full = TRUE)

    vuongtest(mod1, mod2, ll1 = llcont.glmerMod, ll2 = llcont.glmerMod,
              score1 = estfun.glmerMod, score2 = estfun.glmerMod,
              vc1 = vcg, vc2 = vcg)

    Model 1
     Class: glmerMod
     Call: glmer(formula = response ~ -1 + item + agecent + gender + (1 | ...

    Model 2
     Class: glmerMod
     Call: glmer(formula = response ~ -1 + item + semester + elite + (1 | ...

    Variance test
      H0: Model 1 and Model 2 are indistinguishable
      H1: Model 1 and Model 2 are distinguishable
        w2 = 0.033,   p = 6.25e-07

    Non-nested likelihood ratio test
      H0: Model fits are equal for the focal population
      H1A: Model 1 fits better than Model 2
        z = -0.356,   p = 0.639
      H1B: Model 2 fits better than Model 1
        z = -0.356,   p = 0.3611

The output from the function first shows a variance test, which provides
information about whether the non-nested models are distinguishable from each
other via the observed dataset. From this, we reject the hypothesis that the
models are indistinguishable from one another. We then move on to the non-
nested likelihood ratio test to examine whether one model fits better than the
other. For our example, we conclude that neither model fits better than the
other.

Figure 8 shows that the _nonnest2_ functionality can also be used to test
nested models, by adding the nested = TRUE argument. We first fit a simple
Rasch model to the data, with this model being nested in the two considered
previously. We then compute test statistics comparing this model to the second
model from Figure 6. The two test statistics in the output can each be used to
compare the nested models, providing two alternatives to the traditional
likelihood ratio test. Here, we conclude that the full model including the
“semester” and “elite” predictors fits better than the simple Rasch model
without those predictors.

Figure 8: Code for testing fit of two nested models. The full model has better
fit than the reduced model.

    mod3 <- glmer(response ~ -1 + item + (1 | pnum), data = spisa,
                  family = binomial,
                  control = glmerControl(optimizer=’bobyqa’))

    vuongtest(mod2, mod3, nested = TRUE,
              ll1 = llcont.glmerMod, ll2 = llcont.glmerMod,
              score1 = estfun.glmerMod, score2 = estfun.glmerMod,
              vc1 = vcg, vc2 = vcg)

    Model 1
     Class: glmerMod
     Call: glmer(formula = response ~ -1 + item + semester + elite + (1 | ...

    Model 2
     Class: glmerMod
     Call: glmer(formula = response ~ -1 + item + (1 | pnum), data = spisa, ...

    Variance test
      H0: Model 1 and Model 2 are indistinguishable
      H1: Model 1 and Model 2 are distinguishable
        w2 = 0.017,   p = 0.000109

    Robust likelihood ratio test of distinguishable models
      H0: Model 2 fits as well as Model 1
      H1: Model 1 fits better than Model 2
        LR = 18.680,   p = 9.08e-05

###  3.5 Poisson GLMMs

Of course, the GLMM framework is not limited solely to binomial models, and
our derivations extend to other exponential family models. In this section, we
illustrate extensions to the Poisson GLMM using the epilepsy data set [Thall
Vail (1990)] found in the package _brms_ [Bürkner (2018)].

####  3.5.1 Method

The data consist of 236 observations of seizure counts from 59 people across 4
time periods. Covariates include study group (treatment vs control),
participant age, and a base rate seizure count across 8-weeks (standardized).
For our initial model, we predict number of seizures using the patient’s base
rate (zBase), treatment group indicator (Trt), and visit number (visit). We
allow the intercept and visit slope to vary by participant, with these two
random effects being correlated. The _lme4_ code for this model is at the top
of Figure 9.

Figure 9: Code to fit a Poisson GLMM predicting the number of seizures in
epileptic patients, then compute robust standard errors and a score test
statistic.

    ## linear effect of visit number:
    epilepsy$visit <- as.numeric(epilepsy$visit)

    ## Poisson model:
    poimod <- glmer(count ~ zBase * Trt * visit + (visit | patient),
                    data = epilepsy, family = poisson)

    ## Robust standard errors:
    rse <- sandwich(poimod, bread. = bread.glmerMod,
                    meat. = meat(poimod, level = 2))

    ## Score-based test with 5 quadrature points:
    age <- with(epilepsy, tapply(Age, patient, head, 1))
    efg5 <- function(...) estfun.glmerMod(..., nAGQ = 5)

    poisc <- sctest(poimod, fit = NULL, scores = efg5,
                    order.by = age, parm = 3, functional = ’maxLMo’)

Because this model includes multiple random effects, _lme4_ requires that we
use the Laplace approximation (nAGQ = 1) for estimation. We can still choose a
larger number of quadrature points for score computation after model
estimation, however, which provides more precise approximations of these
quantities. We can also use extra quadrature points to compute the model’s
log-likelihood (via the _merDeriv_ command llcont.glmerMod()), which
potentially yields a log-likelihood that is more precise than the log-
likelihood that is output by _lme4_.

####  3.5.2 Results

We first used _merDeriv_ to repeatedly compute the log-likelihood and the
standardized gradient of the estimated Poisson GLMM, using one to ten
quadrature points per dimension (the gradient is obtained by summing scores
across people). Some of those results are shown in Figure 10, where the left
panel displays results for the log-likelihood and the right panel displays
results for the standardized gradient of a single model parameter (the fixed
intercept). We see that, for small numbers of quadrature points, both of the
displayed quantities are unstable. The log-likelihood varies by about a tenth
of a point, while the standardized gradient varies by much more. Both
quantities stabilize around five quadrature points, however, suggesting that
we should use at least that many points in practice (while also considering
total computation time). We also remark that the log-likelihood reported by
_lme4_ is the value in the left panel at 1 quadrature point, which is somewhat
different from the “stabilized” value at larger numbers of quadrature points.
We can obtain a more accurate approximation of the fitted model’s log-
likelihood using the methods described here, and this approximation could
influence some likelihood ratio tests or other statistics that rely on the
model’s log-likelihood.

Figure 10: Log-likelihood and standardized gradient of the Poisson mixed model, by number of quadrature points used. The standardized gradient shown is that of the model’s fixed intercept parameter. | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ | $\beta_{5}$ | $\beta_{6}$ | $\beta_{7}$
---|---|---|---|---|---|---|---|---
lme4 | 0.1575 | 0.1465 | 0.2194 | 0.0468 | 0.1914 | 0.0404 | 0.0656 | 0.0518
sandwich | 0.2464 | 0.1588 | 0.2301 | 0.0771 | 0.1554 | 0.0452 | 0.0679 | 0.0412
Table 2: Comparison of model generated standard errors to robust standard
errors for Poisson model.

We now illustrate how methods from the previous sections can be applied to the
Poisson GLMM. We first calculate robust standard errors using the code in the
middle of Figure 9, with Table 2 showing the results. The table shows that,
for the model considered here, the Huber-White standard errors are generally
larger.

Similarly to the previous section on score-based tests, we next examine the
Poisson GLMM parameter fluctuation across an extraneous variable. In this
example, we assess the stability of the treatment main effect across patient
age, which provides information about whether the treatment efficacy varies
for patients of different ages. The score test is carried out via the code at
the bottom of Figure 9, which is similar to that used in the score test
section above. The test statistic here (not shown) indicates that the
parameter fluctuation is not significant, suggesting that the treatment effect
does not fluctuate across the range of age. Figure 11 contains the parameter
fluctuation across values age, with the critical value being the red
horizontal line.

Figure 11: Graph of M-Fluctuation test for Poisson model. Model parameters are
stable across the range of age.

As can be seen, our methods work for other exponential family GLMMs, with the
code remaining very similar. In the General Discussion below, we provide
further detail about models that our methods cannot handle, as well as future
extensions.

##  4 General Discussion

In this paper, we have provided technical details on computing derivatives of
the GLMM likelihood function, with a focus on models estimated via package
_lme4_. We then showed how the derivatives can be used in various manners: to
obtain robust standard errors, to test predictors that were not included in
the estimated model, and to carry out Vuong tests of non-nested GLMMs. All of
these applications used the GLMM derivatives in concert with other R packages,
illustrating how the R infrastructure can be combined to obtain new statistics
that were difficult or impossible to obtain previously.

###  4.1 Computational Issues

The quadrature implementation described in this paper can be uesd to obtain
derivatives of the marginal likelihood function for many models with random
effects. This method is especially relevant because the conditional random
effects $\bm{b}$ and corresponding components in the variance covariance
matrix $\bm{G}$ are often employed in the model estimation process, in place
of derivatives [<]for example,¿cai2010, cai2010b, bau2004. Therefore, the
derivatives based on the marginal distribution are often not available, or at
least not easy to obtain. Our quadrature method took advantage of the fact
that the model was already estimated, so that the predicted modes of the
random parameters were available.

Another integral approximation method is the Laplace approximation, which is
equivalent to Gauss-Hermite quadrature with one quadrature point [McCulloch
Neuhaus (2005)]. Thus, the Laplace approximation is less accurate than Gauss-
Hermite quadrature with multiple points, but also less computationally
intensive and more flexible [Stroup (2012)]. Additionally, it is possible to
use derivatives associated with the pseudo maximum likelihood function, which
is a transformation of the $y$ response variable into $y^{\star}$, which
conditions on the random effect [Stroup (2012)]. The scores are then related
to a simpler GLM, with such a procedure being implemented in SAS
[Schabenberger (2005)]. However, the scores based on this pseudo likelihood
are not always applicable because the estimates can be problematic, such as
when $y$ follows a two-parameter exponential family distribution or sparse
Bernoulli distribution [Nelder Lee (1992)]. Finally, numerical methods and
Monte Carlo can be flexibly applied to many types of derivative computations,
but they are often too slow to be practical. In all, these remarks indicate
that there is not a single, superior method for all scenarios. The quadrature
method described here is flexible and appears to work well enough for many
types of models.

###  4.2 Additional Applications

There exist other relevant applications that are worth exploring in more
detail, including use of the derivatives in GLMM trees. GLMM trees are part of
a model-based recursive partitioning framework that has been developed by
Zeileis and colleagues [Hothorn Zeileis (2015)]. The goal of the framework is
to split a dataset into homogeneous subsamples based on auxiliary variables,
where each subsample exhibits different values of model parameters. To
accomplish this, a tree is constructed via the following steps

  1. 1.

Fit the model of interest to the data in the current node of the tree.

  2. 2.

Conduct a score-based test for each auxiliary variable.

  3. 3.

Split the current node into two nodes, based on the auxiliary variable with
the largest test statistic.

  4. 4.

Repeat steps 1–3 for the two nodes that were just created.

This procedure is continued until the score-based tests indicate no parameter
instabilities with respect to any auxiliary variables (or until a minimal
subsample size is reached).

glmertree recently applied model-based recursive partitioning to GLMMs. But,
due to the difficulty of obtaining scores associated with GLMMs, they
developed an alternative procedure where only fixed effect parameters were
allowed to vary across subgroups. The developments in the current paper make
it possible to apply the original, model-based recursive partitioning
procedure to GLMMs, allowing us to detect new types of GLMM heterogeneity in a
tree-based framework.

In addition to trees, scores may be used to study heterogeneity through “on
the fly” tests of residual covariance structures in GLMMs. These developments
could reduce computation time by testing multiple covariance structures after
fitting a single model, as opposed to requiring estimation of one model per
covariance structure. Such tests can be facilitated by the coeftest() function
of package _lmtest_ [Zeileis Hothorn (2002)], though some _merDeriv_
extensions may be necessary before this works.

###  4.3 Limitations

While the derivations in this paper work for general, exponential family
models, two-parameter distributions such as the gamma and inverse Gaussian are
additionally complicated by estimation of the extra dispersion parameter. The
current _merDeriv_ implementation does not currently handle some of these
models, nor does it handle the quasi-Poisson or quasi-binomial families (which
are not based on formal likelihood functions). Additionally, the applications
in this paper take advantage of the fact that we focused on models with a
single clustering variable. Researchers often consider three-level models and
models with crossed or partially-crossed random effects, though, which utilize
multiple clustering variables. The derivations in this paper generally work
for those models, allowing us to obtain scores for each case in the data
(i.e., for each row of the data). But most of the applications in this paper
require a way to split observations into independent groups, which is often
impossible when we have multiple clustering variables. For example,
individuals in separate groups under one clustering variable may appear in the
same group under another clustering variable, leading to different forms of
dependence between different pairs of individuals’ scores. In contrast, when
there is only one clustering variable, we know that individuals in one group
are independent of individuals in other groups.

For models with multiple clustering variables, it may be possible to de-
correlate scores after the fact, using an appropriately-specified covariance
matrix [Zeileis (2004), Zeileis . (2020)] or a self-normalization technique
that is commonly used in time series research [Shao Zhang (2010), Zhang .
(2011)]. This would allow us to split observations into uncorrelated groups,
which may be sufficient for applications. Alternatively, rasgol94 describe
methods for re-specifying a model with crossed random effects to be a fully
hierarchical model, in which case it may be possible to directly use the
results described in this paper. None of these solutions is trivial, and we
hope to further study them in the future. We aspire to a future version of
_merDeriv_ that is able to handle all of the models that _lme4_ can estimate.

## Computational Note

All results were obtained using the R system for statistical computing [R Core
Team (2020)], version 3.6.1, employing the add-on package _merDeriv_ 0.2-3 for
derivative computations and _lme4_ 1.1-26 [Bates . (2015)] for fitting of the
mixed models. Code to reproduce the results in the paper is available at
https://osf.io/58ruw/.

## References

  * Barr . (2013) barlev13Barr, DJ., Levy, R., Scheepers, C. Tily, HJ.  2013\.  Random effects structure for confirmatory hypothesis testing: Keep it maximal Random effects structure for confirmatory hypothesis testing: Keep it maximal.  Journal of Memory and Language68255–278.
  * Bates (2021) bates2021Bates, D.  2021\.  Computational methods for mixed models Computational methods for mixed models.  lme4 Package Vignette.  https://cran.r-project.org/web/packages/lme4/vignettes/Theory.pdf
  * Bates . (2015) lme4Bates, D., Mächler, M., Bolker, B. Walker, S.  2015\.  Fitting Linear Mixed-Effects Models Using _lme4_ Fitting linear mixed-effects models using _lme4_.  Journal of Statistical Software6711–48.  10.18637/jss.v067.i01
  * Bauer Curran (2004) bau2004Bauer, DJ. Curran, PJ.  2004\.  The integration of continuous and discrete latent variable models: Potential problems and promising opportunities The integration of continuous and discrete latent variable models: Potential problems and promising opportunities.  Psychological Methods913–29.
  * Bock Lieberman (1970) boclie70Bock, RD. Lieberman, M.  1970\.  Fitting a response model for $n$ dichotomously scored items Fitting a response model for $n$ dichotomously scored items.  Psychometrika35179–198.
  * Bürkner (2018) brmsBürkner, PC.  2018\.  Advanced Bayesian Multilevel Modeling with the R Package brms Advanced Bayesian multilevel modeling with the R package brms.  The R Journal101395–411.  10.32614/RJ-2018-017
  * Cai (20101) cai2010bCai, L.  20101\.  High-dimensional exploratory item factor analysis by a Metropolis-Hastings Robbins-Monro algorithm High-dimensional exploratory item factor analysis by a Metropolis-Hastings Robbins-Monro algorithm.  Psychometrika75133–57.
  * Cai (20102) cai2010Cai, L.  20102\.  A two-tier full-information item factor analysis model with applications A two-tier full-information item factor analysis model with applications.  Psychometrika754581–612.
  * Chalmers (2012) mirtChalmers, RP.  2012\.  mirt: A Multidimensional Item Response Theory Package for the R Environment mirt: A multidimensional item response theory package for the R environment.  Journal of Statistical Software4861–29.  10.18637/jss.v048.i06
  * De Boeck . (2011) de2011De Boeck, P., Bakker, M., Zwitser, R., Nivard, M., Hofman, A., Tuerlinckx, F. Partchev, I.  2011\.  The estimation of item response models with the lmer function from the lme4 package in R The estimation of item response models with the lmer function from the lme4 package in R.  Journal of Statistical Software39121–28.
  * De Boeck Wilson (2004) debwil04De Boeck, P. Wilson, M.  2004\.  Explanatory item response models: A generalized linear and nonlinear approach Explanatory item response models: A generalized linear and nonlinear approach.  New York: Springer-Verlag.
  * Doran . (2007) dorbat07Doran, H., Bates, D., Bliese, P. Dowling, M.  2007\.  Estimating the Multilevel Rasch Model: With the lme4 Package Estimating the multilevel Rasch model: With the lme4 package.  Journal of Statistical Software2021–18.  10.18637/jss.v020.i02
  * Embretson Reise (2000) embrei00Embretson, SE. Reise, SP.  2000\.  Item response theory for psychologists Item response theory for psychologists.  Mahwah, NJ: Erlbaum Associates.
  * Engle (1984) eng84Engle, RF.  1984\.  Wald, likelihood ratio, and Lagrange multiplier tests in econometrics Wald, likelihood ratio, and Lagrange multiplier tests in econometrics.  Z. Griliches MD. Intriligator (), Handbook of Econometrics Handbook of econometrics ( II).  Elsevier.
  * Fokkema . (2018) glmertreeFokkema, M., Smits, N., Zeileis, A., Hothorn, T. Kelderman, H.  2018\.  Detecting Treatment-Subgroup Interactions in clustered Data With Generalized Linear Mixed-Effects Model Trees Detecting treatment-subgroup interactions in clustered data with generalized linear mixed-effects model trees.  Behavior Research Methods502016-2034.  http://link.springer.com/article/10.3758/s13428-017-0971-x
  * Glas (1992) glas92Glas, CAW.  1992\.  A Rasch model with a multivariate distribution of ability A Rasch model with a multivariate distribution of ability.  Objective measurement: Theory into practice1236–258.
  * Glas (1998) glas98Glas, CAW.  1998\.  Detection of differential item functioning using Lagrange multiplier tests Detection of differential item functioning using Lagrange multiplier tests.  Statistica Sinica83647–667.
  * Glas (1999) glas99Glas, CAW.  1999\.  Modification indices for the 2-PL and the nominal response model Modification indices for the 2-PL and the nominal response model.  Psychometrika64273–294.
  * Hothorn Zeileis (2015) partyHothorn, T. Zeileis, A.  2015\.  _partykit_ : A Modular Toolkit for Recursive Partytioning in R _partykit_ : A modular toolkit for recursive partytioning in R.  Journal of Machine Learning Research163905-3909.  http://jmlr.org/papers/v16/hothorn15a.html
  * Huber (1967) huber67Huber, PJ.  1967\.  The Behavior of Maximum Likelihood Estimates Under Nonstandard Conditions The behavior of maximum likelihood estimates under nonstandard conditions.  Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability Proceedings of the fifth Berkeley symposium on mathematical statistics and probability ( 1, 221–233).
  * Komboz . (2018) kom18Komboz, B., Strobl, C. Zeileis, A.  2018\.  Tree-based global model tests for polytomous Rasch models Tree-based global model tests for polytomous Rasch models.  Educational and Psychological Measurement781128–166.
  * Liu Pierce (1994) liupie94Liu, Q. Pierce, DA.  1994\.  A note on Gauss-Hermite quadrature A note on Gauss-Hermite quadrature.  Biometrika81624–629.
  * Lord Novick (1968) lornov68Lord, FM. Novick, MR.  1968\.  Statistical theories of mental test scores Statistical theories of mental test scores.  Reading, MA: Addison-Wesley.
  * Louis (1982) lou82Louis, TA.  1982\.  Finding the observed information matrix when using the EM algorithm Finding the observed information matrix when using the EM algorithm.  Journal of the Royal Statistical Society. Series B (Methodological)226–233.
  * Matuschek . (2017) matkli17Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H. Bates, D.  2017\.  Balancing Type I error and power in linear mixed models Balancing Type I error and power in linear mixed models.  Journal of Memory and Language94305–315.
  * McCullagh Nelder (1989) mccnel89McCullagh, P. Nelder, JA.  1989\.  Generalized linear models Generalized linear models.  Boca Raton, FL: Chapman & Hall.
  * McCulloch Neuhaus (2001) mcc01McCulloch, CE. Neuhaus, JM.  2001\.  Generalized linear mixed models Generalized linear mixed models.  Wiley Online Library.
  * McCulloch Neuhaus (2005) mcc2005McCulloch, CE. Neuhaus, JM.  2005\.  Generalized linear mixed models Generalized linear mixed models.  Encyclopedia of biostatistics4.
  * Merkle . (2014) merfanzeiMerkle, EC., Fan, J. Zeileis, A.  2014\.  Testing for Measurement Invariance with Respect to an Ordinal Variable Testing for measurement invariance with respect to an ordinal variable.  Psychometrika79569–584.
  * Merkle . (2019) merfur19Merkle, EC., Furr, D. Rabe-Hesketh, S.  2019\.  Bayesian comparison of latent variable models: Conditional versus marginal likelihoods Bayesian comparison of latent variable models: Conditional versus marginal likelihoods.  Psychometrika84802–829.
  * Merkle You (2018) nonnest2Merkle, EC. You, D.  2018\.  _nonnest2_ : Tests of Non-Nested Models _nonnest2_ : Tests of non-nested models [].  https://cran.r-project.org/package=nonnest2 R package version 0.5-2
  * Merkle . (2016) merkle16Merkle, EC., You, D. Preacher, KJ.  2016\.  Testing Nonnested Structural Equation Models. Testing nonnested structural equation models.  Psychological Methods212151–163.
  * Merkle Zeileis (2013) merzei13Merkle, EC. Zeileis, A.  2013\.  Tests of Measurement Invariance without Subgroups: A Generalization of Classical Methods Tests of measurement invariance without subgroups: A generalization of classical methods.  Psychometrika7859–82.
  * Naylor Smith (1982) naysmi82Naylor, JC. Smith, AFM.  1982\.  Applications of a method for the efficient computation of posterior distributions Applications of a method for the efficient computation of posterior distributions.  Journal of the Royal Statistical Society C31214–225.
  * Nelder Lee (1992) nel92Nelder, J. Lee, Y.  1992\.  Likelihood, quasi-likelihood and pseudolikelihood: some comparisons Likelihood, quasi-likelihood and pseudolikelihood: some comparisons.  Journal of the Royal Statistical Society: Series B (Methodological)541273–284.
  * Open Source Psychometrics Project () osppOpen Source Psychometrics Project.  .  Open psychology data: Raw data from online Open psychology data: Raw data from online [].  [2017-10-17]https://openpsychometrics.org/_rawdata/
  * Petersen Pedersen (2012) peter08Petersen, KB. Pedersen, MS.  2012\.  The Matrix Cookbook. The matrix cookbook.  Technical University of Denmark.  http://www2.imm.dtu.dk/pubdb/p.php?3274 Version 20121115
  * Pinheiro Bates (1995) pinbat95Pinheiro, JC. Bates, DM.  1995\.  Approximations to the log-likelihood function in the nonlinear mixed-effects model Approximations to the log-likelihood function in the nonlinear mixed-effects model.  Journal of Computational Graphics and Statistics412-35.
  * Powell (2009) powell2009Powell, MJ.  2009\.  The BOBYQA algorithm for bound constrained optimization without derivatives The BOBYQA algorithm for bound constrained optimization without derivatives.  Cambridge NA Report NA2009/06, University of Cambridge, Cambridge26–46.
  * R Core Team (2020) R20R Core Team.  2020\.  R: A Language and Environment for Statistical Computing R: A language and environment for statistical computing [].  Vienna, Austria.  https://www.R-project.org/
  * Rabe-Hesketh . (2005) rab05Rabe-Hesketh, S., Skrondal, A. Pickles, A.  2005\.  Maximum likelihood estimation of limited and discrete dependent variable models with nested random effects Maximum likelihood estimation of limited and discrete dependent variable models with nested random effects.  Journal of Econometrics1282301–323.
  * Rasbash Goldstein (1994) rasgol94Rasbash, J. Goldstein, H.  1994\.  Efficient analysis of mixed hierarchical and cross-classified random structures using a multilevel model Efficient analysis of mixed hierarchical and cross-classified random structures using a multilevel model.  Journal of Educational and Behavioral Statistics19337–350.
  * Schabenberger (2005) SASSchabenberger, O.  2005\.  Introducing the GLIMMIX procedure for generalized linear mixed models Introducing the GLIMMIX procedure for generalized linear mixed models.  SUGI 30 Proceedings196.
  * Schneider . (2020) schcha19Schneider, L., Chalmers, RP., Debelak, R. Merkle, EC.  2020\.  Model selection of nested and non-nested item response models using Vuong tests Model selection of nested and non-nested item response models using Vuong tests.  Multivariate Behavioral Research55664–684.
  * Shao Zhang (2010) shao10Shao, X. Zhang, X.  2010\.  Testing for change points in time series Testing for change points in time series.  Journal of the American Statistical Association1054911228–1240.
  * Skrondal Rabe-Hesketh (2004) skrrab04Skrondal, A. Rabe-Hesketh, S.  2004\.  Generalized latent variable modeling: Multilevel, longitudinal, and structural equation modeling Generalized latent variable modeling: Multilevel, longitudinal, and structural equation modeling.  Boca Raton, FL: Chapman & Hall.
  * Strobl . (2015) strkop15Strobl, C., Kopf, J. Zeileis, A.  2015\.  Rasch Trees: A New Method for Detecting Differential Item Functioning in the Rasch Model Rasch trees: A new method for detecting differential item functioning in the Rasch model.  Psychometrika802289–316.  10.1007/s11336-013-9388-3
  * Stroup (2012) stroup12Stroup, WW.  2012\.  Generalized Linear Mixed Models: Modern Concepts, Methods and Applications Generalized linear mixed models: Modern concepts, methods and applications.  CRC Press.
  * Stroup Claassen (2020) str20Stroup, WW. Claassen, E.  2020\.  Pseudo-Likelihood or Quadrature? What We Thought We Knew, What We Think We Know, and What We Are Still Trying to Figure Out Pseudo-likelihood or quadrature? What we thought we knew, what we think we know, and what we are still trying to figure out.  Journal of Agricultural, Biological and Environmental Statistics.
  * Thall Vail (1990) thall1990Thall, PF. Vail, SC.  1990\.  Some covariance models for longitudinal count data with overdispersion Some covariance models for longitudinal count data with overdispersion.  Biometrics657–671.
  * Trepte Verbeet (2010) trever10Trepte, S. Verbeet, M. ().  2010\.  Allgemeinbildung in Deutschland – Erkenntnisse aus dem SPIEGEL Studentenpisa-Test Allgemeinbildung in Deutschland – erkenntnisse aus dem SPIEGEL Studentenpisa-Test.  WiesbadenVS Verlag.
  * Vuong (1989) vuo89Vuong, QH.  1989\.  Likelihood Ratio Tests for Model Selection and Non-Nested Hypotheses Likelihood ratio tests for model selection and non-nested hypotheses.  Econometrica57307–333.
  * Wang Merkle (2018) merDerivWang, T. Merkle, EC.  2018\.  merDeriv: Derivative Computations for Linear Mixed Effects Models with Application to Robust Standard Errors merDeriv: Derivative computations for linear mixed effects models with application to robust standard errors.  Journal of Statistical Software8711–16.  10.18637/jss.v087.c01
  * Wang . (2020) wanmer20Wang, T., Merkle, EC., Anguera, JA. Turner, BM.  2020\.  Score-based tests for detecting heterogeneity in linear mixed models Score-based tests for detecting heterogeneity in linear mixed models.  Behavior Research Methods.
  * Wang . (2018) wanstr18Wang, T., Strobl, C., Zeileis, A. Merkle, EC.  2018\.  Score-based tests of differential item functioning via pairwise maximum likelihood estimation Score-based tests of differential item functioning via pairwise maximum likelihood estimation.  Psychometrika831132–155.
  * White (1980) white80White, H.  1980\.  A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity.  Econometrica: Journal of the Econometric Society484817–838.
  * Zeileis (2004) sand1Zeileis, A.  2004\.  Econometric Computing with HC and HAC Covariance Matrix Estimators Econometric computing with HC and HAC covariance matrix estimators.  Journal of Statistical Software11101–17.  http://www.jstatsoft.org/v11/i10/
  * Zeileis (2006) sandwichbZeileis, A.  2006\.  Object-Oriented Computation of Sandwich Estimators Object-oriented computation of sandwich estimators.  Journal of Statistical Software1691–16.  10.18637/jss.v016.i09
  * Zeileis Hornik (2007) zeihor07Zeileis, A. Hornik, K.  2007\.  Generalized M-Fluctuation Tests for Parameter Instability Generalized M-fluctuation tests for parameter instability.  Statistica Neerlandica61488–508.
  * Zeileis Hothorn (2002) zeihot02Zeileis, A. Hothorn, T.  2002\.  Diagnostic Checking in Regression Relationships Diagnostic checking in regression relationships.  R News237–10.  https://CRAN.R-project.org/doc/Rnews/
  * Zeileis . (2020) sandwichcZeileis, A., Köll, S. Graham, N.  2020\.  Various Versatile Variances: An Object-Oriented Implementation of Clustered Covariances in R Various versatile variances: An object-oriented implementation of clustered covariances in R.  Journal of Statistical Software951.  10.18637/jss.v095.i01
  * Zeileis . (2002) strucchangeZeileis, A., Leisch, F., Hornik, K. Kleiber, C.  2002\.  strucchange: An R Package for Testing for Structural Change in Linear Regression Models strucchange: An R package for testing for structural change in linear regression models.  Journal of Statistical Software721–38.  http://www.jstatsoft.org/v07/i02/
  * Zhang . (2011) zhang11Zhang, X., Shao, X., Hayhoe, K. Wuebbles, DJ.  2011\.  Testing the structural stability of temporally dependent functional observations and application to climate projections Testing the structural stability of temporally dependent functional observations and application to climate projections.  Electronic Journal of Statistics51765–1796.