Generalised Linear Mixed Model Analysis, Maximum Likelihood Fitting, Prediction, Simulation, and Optimal Design in \proglangR

A generalised linear mixed model (GLMM) has the linear predictor for observation $i$

where $\mathbf{x}_{i}$ is the $i$th row of matrix $X$, which is a $n\times P$ matrix of covariates, $\boldsymbol{\beta}$ is a vector of parameters, $\mathbf{z}_{i}$ is the $i$th row of matrix $Z$, which is the $n\times Q$ “design matrix” for the random effects, and $\mathbf{u}\sim N(0,D)$, where $D$ is the $Q\times Q$ covariance matrix of the random effects terms that depends on parameters $\boldsymbol{\theta}$. In this article we use bold lower case characters, e.g. $\mathbf{x}_{i}$ to represent vectors, normal script lower case, e.g. $\beta_{1}$, to represent scalars, and upper case letters, e.g. $X$, to represent matrices.

The model is then

$\mathbf{y}$ is a $n$-length vector of outcomes with elements $y_{i}$, $G$ is a distribution, $h(.)$ is the link function such that $\mu_{i}=h(\eta_{i})$ where $\mu_{i}$ is the mean value, and $\phi$ is an additional scale parameter to complete the specification.

When $G$ is a distribution in the exponential family, we have:

The likelihood of this model is given by:

The likelihood (2) generally has no closed form solution, and so different algorithms and approximations have been proposed to estimate the model parameters. We discuss model fitting in Section 3. Where relevant we represent the set of all parameters as $\boldsymbol{\Theta}=(\boldsymbol{\beta},\phi,\boldsymbol{\theta})$.

We introduce methods to generate data for hierarchical models and blocked designs. As we show in subsequent examples, ‘fake’ data generation is typically needed to specify a particular data structure at the design stage of a study, but can be useful in other circumstances. Nelder (1965) suggested a simple notation that could express a large variety of different blocked designs. The notation was proposed in the context of split-plot experiments for agricultural research, where researchers often split areas of land into blocks, sub-blocks, and other smaller divisions, and apply different combinations of treatments. However, the notation is useful for expressing a large variety of experimental designs with correlation and clustering, including cluster trials, cohort studies, and spatial and temporal prevalence surveys. We have included the function \codenelder() in the package that generates a data frame of a design using Nelder’s notation.

There are two operations:

1. 1.
   \code

   > (or $\to$ in Nelder’s notation) indicates “clustered in”.

2. 2.
   \code

   * (or $\times$ in Nelder’s notation) indicates a crossing that generates all combinations of two factors.

The function takes a formula input indicating the name of the variable and a number for the number of levels, such as \codeabc(12). So for example \code cl(4) > ind(5) means in each of five levels of \codecl there are five levels of \codeind, and the individuals are different between clusters. Brackets are used to indicate the order of evaluation. Some specific examples are illustrated in Table 2.

Use of this function produces a data frame: {CodeChunk} {CodeInput} data <- nelder( (j(4) * t(5)) > i(5)) head(data) {CodeOutput} > j t i > 1 1 1 1 > 2 1 1 2 > 3 1 1 3 > 4 1 1 4 > 5 1 1 5 > 6 1 2 6

The data frame shown above may represent, for example, a cluster randomised study with cross-sectional sampling. Such an approach to study design assumes the same number of each factor for each other factor, which is not likely adequate for certain study designs. We may expect unequal cluster sizes, staggered inclusion/drop-out, and so forth, and so a user-generated data set would instead be required. Certain treatment conditions may be specified with this approach including parallel trial designs, stepped-wedge implementation, or factorial approaches by specifying a treatment arm as part of the block structure.

The package allows for flexible specification of $g(\mathbf{x_{i}},\boldsymbol{\beta})$. Standard model formulae in \proglangR can be used, for example \codeformula =   factor(t) - 1, the package also allows a wide range of specifications including non-linear functions and naming parameters in the function, including allowing the same parameter to appear in multiple places. For example, the function $\beta_{0}+\beta_{int}\texttt{int}_{i}+\beta_{1}\exp(\beta_{2}\texttt{x}_{i})$ can be specified as \codeformula =   int + b_1*exp(b_2*x). Names of data are assumed to be multiplied by a parameter as in the standard linear predictor \code  x would produce $\beta_{x}x$. To include data without a parameter, it can be wrapped in brackets such as \code (x), which can also provide a way of specifying offsets, although this can also be achieved by directly specifying the offset.

We use a recursive, operator precedence algorithm to parse the formula and generate a reverse Polish notation representing the formula. The algorithm first splits the formula at the first ‘+’ or ‘-’ symbol outside of a bracket, then if there are none, at a multiply (‘*’) or divide (‘/’) symbol outside of a bracket, then a power (‘’̂) symbol, and finally brackets, until a token is produced (e.g. \codeb_2 or \codeexp). The resulting token is first checked against available functions, then if it is not a function it is checked against the column names of the data, and if it is not in the data it is considered the name of a parameter. The resulting tree structure can then be parsed up each of its branches to produce an instruction set. Figure 1 shows an example of the formula parsing algorithm, and calculation of the value of the formula. This scheme also allows for an auto-differention approach to obtain the first and second derivatives where required.

We adapt and extend the random effects specification approach of the popular \proglangR package \pkglme4 (Bates et al., 2015) and other related packages, building on the flexible formula parsing and calculation described above, to allow for relatively complex structures through the specification of a covariance function (see Section 3 for a summary of related packages). A covariance function is specified as an additive formula made up of components with structure \code(z|f(j)). The left side of the vertical bar specifies the covariates in the model that have a random effects structure; using \code(1|f(j)) specifies a ‘random intercept’ model. The right side of the vertical bar specifies the covariance function \codef for that term using variable \codej named in the data. Multiple covariance functions on the right side of the vertical bar are multiplied together, i.e., \code(1|f(j)*g(t)). The currently implemented functions (as of version 0.6.1) are listed in Table 3.

One combines smaller functions to provide the desired overall covariance function. For example, for a stepped-wedge cluster randomised trial we could consider the standard specification (see, for example, Li et al. (2021)) with an exchangeable random effect for the cluster level (\codej), and a separate exchangeable random effect for the cluster-period (with \codet representing the discrete time period), which would be \code (1|gr(j))+(1|gr(j,t)). Alternatively, we could consider an autoregressive cluster-level random effect that decays exponentially over time so that in a linear mixed model, for person $i$ in cluster $j$ at time $t$, $Cov(y_{ijt},y_{i^{\prime}jt})=\theta_{1}^{2}$, for $i\neq i^{\prime}$, $Cov(y_{ijt},y_{i^{\prime}jt^{\prime}})=\theta_{1}^{2}\theta_{2}^{|t-t^{\prime}|}$ for $t\neq t^{\prime}$, and $Cov(y_{ijt},y_{i^{\prime}j^{\prime}t})=0$ for $j\neq j^{\prime}$. This function would be specified as \code (1|gr(j)*ar1(t)). We could also supplement the autoregressive cluster-level effect with an individual-level random intercept, where \codeind specified the identifier of the individual, as \code (1|gr(j)*ar1(t)) + (1|gr(ind)), and so forth. To add a further random parameter on some covariate \codex we must add an additional term, for example, \code (1|gr(j)*ar1(t)) + (1|gr(ind)) + (x|gr(ind)). As another example, to implement a spatial Gaussian process model with exponential covariance function and two Cartesian coordinates \codex and \codey with autoregressive temporal decay we can specify \code (1|ar(t)*fexp0(x,y).

The formula specifying the random effects is translated into a form that can be used to efficiently generate the value at any position in the matrix $D$ so that it can be quickly updated when the parameter values change, which facilitates functionality such as model fitting when using an arbitrary combination of covariance functions.

The matrices $D$ and $L$ are stored using sparse compressed row storage (CRS) format since for a large number of cases $D$ and its Cholesky factorisation are sparse. For example, a large number of covariance function specifications and their related study designs lead to a block diagonal structure for $D$:

Internally, the formula into parsed into blocks, both determined by the presence of a \codegr function in any particular random effect specification, and for each additive component of the specification. These blocks may be of differing dimensions and have different formulae and parameters. The data defining each element is a Euclidean distance and within each (likely dense) block, only the data for the strictly lower triangular portion needs to be stored. These distances are calculated once and stored in a flat data structure. We use a reverse Polish notation for to specify the instructions for a given block, as with the linear predictor, along with parameter indices and other necessary information.

By default, the Cholesky decomposition $L$, such that $D=LL^{T}$, is calculated using a sparse matrix LDL algorithm (Davis, 2005), which also provides an efficient forward substitution algorithm for computing quadratic forms in the multivariate Gaussian likelihood (see Section 3). Additionally, an efficient permutation is calculated on instantiation of the class using the approximation minimum degree algorithm (Furrer et al., 2006; Amestoy et al., 2004), which may provide a performance gain when factorising the matrix $PDP^{T}$ by ensuring the resulting factorisation is sparse. We provide a basic \proglangC++ sparse matrix class, the LDL decomposition, and basic operations in the \pkgSparseChol package for \proglangR.

As an alternative we can directly exploit the block diagonal structure of $D$. We can directly generate the decomposition of each block $L_{b}$ separately. We use the Cholesky–Banachiewicz algorithm for this purpose. Many calculations can also be broken down by block, for example $X^{T}\Sigma^{-1}X=\sum_{b}X^{T}_{b}L_{b}^{-T}L_{b}^{-1}X_{b}$. This example can also be further simplified using a forward substitution algorithm (see Section 3.5). To switch between sparse and dense matrix methods, one can use \codemodel$sparse(TRUE) and \codemodel$sparse(FALSE), respectively. The default is to use sparse methods.

For many types of GLMM, the matrix $D$ is dense and not blocked or sparse. For example, geospatial statistical models often specify spatial or spatio-temporal Gaussian process models for the random effects (Diggle et al., 1998). For even moderately sized data sets, this can result in very slow model fitting for the multivariate Gaussian likelihood. There are several approaches to reducing model fitting time. One approach is to use a compactly supported covariance function that leads to a sparse matrix $D$ while still approximating a Gaussian process. This method can be viewed as ‘covariance tapering’ (Kaufman et al., 2008), although the approach implemented here is what they describe as the ‘one taper’ method, which may be biased in some circumstances. As shown in Table 3, we include several compactly supported and parameterised covariance functions, which are described in Gneiting (2002). A covariance function with compact support has value zero beyond some ‘effective range’. Using these functions can result in a sparse matrix $D$, which is then amenable to the sparse matrix methods in this package. To implement these functions with an effective range of $r$, beyond which the covariance is zero, one must divide each variable by $r$. For example, if there are two spatial dimensions, \codex and \codey, in the data frame \codedata then one would create \codedata$xr <- data$x/r, and equivelently for \codey. Then, one could use the covariance function, for example, \code  truncpow3(xr,yr).

As an alternative to compactly supported covariance functions, one can instead use a Gaussian process approximation for large, dense covariance matrices. We provide two such approximations in this package.

The matrix $Z$ is constructed in a similar way as other packages, such as described for \pkglme4 (Bates et al., 2015). $Z$ is comprised of the matrices corresponding to each block $Z=[Z_{1},Z_{2},...,Z_{B}]$. For a model formula \code(1|f(x1,x2,…)), the dimension of the corresponding matrix $Z_{b}$ is $n$ rows and number of columns equal to the number of unique rows of the combination of the variables in \codex1,x2,…. The $i$th row and $j$th column of $Z_{b}$ is then equal to 1 if the $i$th individual has the value corresponding to the $j$th unique combination of the variables and zero otherwise. For formulae specifying random effects on covariates (“variable slopes models”), e.g. \code(z|f(x1,x2,…)), then the $i$th row and $j$th column of $Z_{b}$ is further multiplied by $x_{i}$. $Z$ is also stored in CRS format and uses sparse matrix methods for multiplication and related calculations.

For Gaussian models, and other distributions requiring an additional scale parameter $\phi$, one can also specify the option \codevar_par which is the conditional variance $\phi=\sigma$ at the individual level. The default value is 1. Currently (version 0.6.1), the package supports the following families (link functions): Gaussian (identity, log), Poisson (log, identity), Binomial (logit, log, probit, identity), Gamma (log, inverse, identity), and Beta (logit). For the beta distribution one must provide \codeBeta() to the \codefamily argument.

The \codeModel class can provide an approximation to the covariance matrix $\Sigma$ with the member function \code$Sigma(). This approximation is also used to calculate the expected information matrix with \code$information_matrix() and estimate power in the member function \code$power() (see Section 2.12). We use the first-order approximation based on the marginal quasi-likelihood proposed by Breslow and Clayton (1993):

where $W^{-1}=\text{diag}\left(\left(\frac{\partial h^{-1}(\boldsymbol{\eta})}{\partial\boldsymbol{\eta}}\right)^{2}\text{Var}(\mathbf{y}|\mathbf{u})\right)$, which are recognisable as the GLM iterated weights (Breslow and Clayton, 1993; McCullagh and Nelder, 1989). For Gaussian-identity mixed models this approximation is exact. The diagonal of the matrix $W$ can be obtained using the member function \code$w_matrix(). The information matrix for $\beta$ is:

Zeger et al. (1988) suggest that when using the marginal quasilikelihood, one can improve the approximation to the marginal mean by “attenuating” the linear predictor for non-linear models. For example, with the log link the “attenuated” mean is $E(y_{i})\approx h^{-1}(x_{i}\beta+z_{i}Dz_{i}^{T}/2)$ and for the logit link $E(y_{i})\approx h^{-1}(x_{i}\beta\text{det}(aDz_{i}^{T}z_{i}+I)^{-1/2})$ with $a=16\sqrt{3}/(15\pi)$. To use “attenuation” one can set \codemod$use_attenutation(TRUE), the default is not to use attenutation.

The parameters of the covariance function and linear predictor can updated using the function \codeupdate_parameters(), which then triggers re-generation of the matrices and updates the data in the underlying \proglangC++ classes: \codemodel$update_parameters( cov.pars = c(0.1,0.9)).

The \pkgR6 class system provides many of the standard features of other object orientated systems, including class inheritance. The classes specified in the \pkgglmmrBase package can be inherited from to provide the full range of calculations to other applications. As an example we can define a new class that has a member function that returns the log determinant of the matrix $D$: {CodeChunk} {CodeInput} CovDet <- R6::R6Class("CovDet", inherit = Covariance, public = list( det = function() return(Matrix::determinant(self$D)$modulus) )) cov <- CovDet$new(formula=~(1|gr(j)*ar1(t)),parameters=c(0.05,0.8),data=data)cov$det() [1] -72.26107 More complex applications may include classes to implement simulation-based analyses that simulate new data and fit a GLMM using one of the package’s model fitting methods.

Power and sample size calculations are an important component of the design of many studies, particularly randomised trials, which we use as an example. Cluster randomised trials frequently use GLMMs for data analysis, and hence are the basis for estimating the statistical power of a trial design (Hooper et al., 2016). Different approaches are used to estimate power or calculate sample size given an assumed correlation structure, most frequently using “design effects” (Hemming et al., 2020). However, given the large range of possible models and covariance structures, many software packages implement a narrow range of specific models and provide wrappers to other functions. For example, we identified eight different \proglangR packages available on CRAN or via a R Shiny web interface for calculating cluster trial power or sample size. These packages and their features are listed in Tables 4 and 5. As is evident, the range of models is relatively limited. Beyond this specific study type the R package \pkgsimr provides functions to simulate data from GLMMs and estimate statistics like power using Monte Carlo simulation for specific designs.

As an alternative, the flexible model specification and range of functionality of \pkgglmmrBase can be used to quickly estimate the power for a particular design and model. The \codeModel class member function \code$power() estimates power by calculating the information matrix $M$ described above. The square root of the diagonal of this matrix provides the (approximate) standard errors of $\beta$. The power of a two-sided test for the $i$th parameter at a type I error level of $\alpha$ is given by $\Phi(|\beta_{i}|/\sqrt{M_{ii}}-\Phi^{-1}(1-\alpha/2))$ where $\Phi$ is the Gaussian cumulative distribution function. As an example, to estimate power for a stepped-wedge parallel cluster randomised trial (see Hooper et al. (2016)) with 10 clusters and 11 time periods and ten individuals per cluster period, where we use a Binomial-logit model with time period fixed effects set to zero, and use a auto-regressive covariance function with parameters 0.25 and 0.7: {CodeChunk} {CodeInput} data <- nelder( (cl(10) * t(11)) > i(10)) data$int<-0#intistheinterventioneffectdata[data$t > data$cl,^{\prime}int^{\prime}]<-1model<-Model$new(formula =   factor(t) + int - 1 + (1|gr(cl)*ar1(t)), data = data, covariance = c(0.05,0.7), mean = c(rep(0,12),0.5), family = binomial()) model$power()\CodeOutput ParameterValueSEPower...12int0.50.18161360.7861501Onecanusethefunctionalityofthe\code{Model}classtoinvestigatepowerforarangeofmodelparameters.Forexample,toestimatepowerforcovarianceparametersintheranges(0.05,0.5)and(0.2,0.9)wecancreateadataframeholdingtherelevantcombinationofvaluesanditerativelyupdatethemodel:\CodeChunk\CodeInput param_{d}ata<-expand.grid(par1=seq(0.05,0.5,by=0.05),par2=c(0.2,0.9,by=0.1),power=NA)for(iin1:nrow(param_{d}ata)){model$update_{p}arameters(cov.pars=c(param_{d}ata$par1[i],param_{d}ata$par2[i]))param_{d}ata$power[i]<-model$power()[12,3]}head(param_{d}ata)\CodeOutput par1par2power10.050.20.834886320.100.20.785229830.150.20.738165840.200.20.694513750.250.20.654497060.300.20.6180314\par$

The \codeglmmrBase package also simplifies data simulation for GLMMs. Data simulation is a commonly used approach to evaluate the properties of estimators and statistical models. The \codeModel class has member function \codesim_data(), which will generate either vector of simulated output $y$ (argument \codetype="y"), a data frame that combines the linked dataset with the newly simulated data (argument \codetype="data"), or a list containing the simulated data, random effects, and matrices $Z$ and $X$ (argument \codetype="all"). To quickly simulate data new random effects are drawn as $\mathbf{v}\sim N(0,1)$ and then the random effects component of the model simulated as $ZL\mathbf{v}$. Simulation only of the random effects terms can be done with the relevant function in the \codeCovariance class, \code$covariance$simulate_re().

Software for fitting generalised linear mixed models has been available in \proglangR environment (R Core Team, 2021) for many years. The widely-used package \pkglme4 provides maximum likelihood and restricted maximum likelihood (REML) estimators (Bates et al., 2015). \pkglme4 builds on similar methods proposed by e.g. Bates and DebRoy (2004) and Henderson (1982). For non-linear models it uses Laplacian approximation or adaptive Gaussian quadrature. However, \pkglme4 only allows for compound symmetric (or exchangable) covariance structures, i.e. group-membership type random effects, and linear, additive linear predictors. Function \codeglmmPQL() in package \pkgMASS (Venables and Ripley, 2002) implements Penalized quasi-likelihood (PQL) methods for GLMMs, but is again limited to the same covariance structures, and linear, additive model structure. The package \pkgnlme (Pinheiro et al., 2022) also uses REML approaches for linear mixed models and approximations for models with mean specifications with non-linear functions of covariates discussed in Lindstrom and Bates (1990). \pkgnlme provides a wider range of covariance structures, such as autoregressive functions, which are also available nested within a broader group structure for linear, additive linear predictors. The package providing Frequentist model estimation of GLMMs with the broadest support for different specifications, including many different covariance functions, is \pkgglmmTMB, which uses the \pkgTMB (Template Model Builder) package to support model specification and implementation of Laplacian Approximation to estimate models, although also only permits linear, additive linear predictors.

Outside of R, GLMM model fitting is also provided in other popular statistical software packages. Stata offers a range of functions for fitting mixed models. The \codextmixed function offers both maximum likelihood and restricted maximum likelihood methods to estimate linear mixed models with multi-level group membership random effect structures and heteroscedasticity. For GLMMs, there is a wide range of models that allow for random effects, many of which are accessed through the \codemeglm suite of commands. Estimation can use Gauss-Hermite quadrature or Laplace approximation, but random effect structures are generally limited to group membership type effects. For the software SAS, \codePROC MIXED, and \codePROC GLIMMIX provide GLMM model fitting using either quadrature or Laplace approximation methods. These commands provide a relatively broad random of random effect coviariance structures, including autoregressive, Toeplitz, and exponential structures. \codePROC NLMIXED allows for non-linear functions of fixed and random effects; it uses a Gaussian quadrature method to approximate the integral in the log-likelihood, and a Newton-Raphson step for empirical Bayes estimation of the random effects. However, \codePROC NLMIXED only allows for group membership random effects. Both Stata and SAS also provide a range of options for standard error estimation and correction. Both software packages are proprietary.

We also provide in-built robust and bias-corrected standard error estimates where required. In R, the packages \pkglmerTest and \pkgpbkrtest provide a range of corrections to standard errors for models fit with \codelme4. The package architecture means that \pkgglmmrBase can provide these standard errors for the range of models available in the package, and generally in less time as they do not require refitting of the model. Stata and SAS can calculate bias-corrected and robust standard error estimates with the commands described above.

We also note that there are also several packages for Bayesian model fitting in R, including \pkgbrms (Bürkner, 2017) and \pkgrstanarm, which interface with the Stan statistical software (described below), and \pkgMCMCglmm. Stan is a fully-fledged probabilistic programming language that implements MCMC, and allows for a completely flexible approach to model specification (Carpenter et al., 2017). Stan can also provide limited maximum likelihood fitting functionalist and related tools for calculating gradients and other model operations.

Evidently, there already exists a good range of support for GLMM model fitting available in \proglangR. \pkgglmmrBase may be seen to provide a complement to existing resources rather than another substitute. As we discuss below, in a maximum likelihood context, the existing software all provide approximations rather than full likelihood model fitting. These approximations are typically accurate and fast, but they can also frequently fail to converge and in more complex models the quality of the approximation may degrade, which has lead to methods for bias corrections and improvements to these approximations (e.g. Breslow and Lin (1995); Lin and Breslow (1996); Shun and McCullagh (1995); Capanu et al. (2013)). \pkgglmmrBase provides an alternative that can be used to support more complex analyses, flexible covariance specification, and provide a comparison to ensure approximations are reliable, while also providing an array of other functionality.

Markov Chain Monte Carlo Maximum Likelihood (MCML) are a family of algorithms for estimating GLMMs using the full likelihood that treat the random effect terms as missing data, which are simulated on each iteration of the algorithm. Estimation can then be achieved to an arbitrary level of precision depending on the number of samples used. Approximate methods such as Laplace approximation or quadrature provide significant computational advantages over such approaches, however they may trade-off the quality of parameter and variance estimates versus full likelihood approaches. Indeed, the quality of the Laplace approximation method for GLMMs can deteriorate if there are smaller numbers of clustering groups, smaller numbers of observations to estimate the variance components, or if the random effect variance is relatively large.

The package provides full likelihood model fitting for the range of GLMMs that can be specified using the package via the MCML algorithms described by McCulloch (1997) or stochastic approximation expectation maximisation (SAEM) (Jank, 2006). Each algorithm has three main steps per iteration: draw samples of the random effects using Markov Chain Monte Carlo (MCMC), the mean function parameters are then estimated conditional on the simulated random effects, and then the parameters of the multivariate Gaussian likelihood of the random effects are then estimated. The process is repeated until convergence. We describe each step in turn and optimisations.

Reliable convergence of the algorithm requires a set of $m_{i}$ independent samples of the random effects on the $i$th iteration, i.e. realisations from the distribution of the random effect conditional on the observed data, and with fixed parameters. MCMC methods can be used to sample from the ‘posterior density’ $f_{\mathbf{u}|\mathbf{y}}(\mathbf{u}|\mathbf{y},\beta^{(i)},\phi^{(i)},\theta^{(i)})\propto f_{\mathbf{y}|\mathbf{u}}(\mathbf{y}|\mathbf{u},\beta^{(i)},\phi^{(i)})f_{\mathbf{u}}(\mathbf{u}|\theta^{(i)})$. To improve sampling, instead of generating samples of $\mathbf{u}$ directly, we instead sample $\mathbf{v}$ from the model:

where $\tilde{Z}=ZL$, and $I$ is the identity matrix. Once $m_{i}$ samples of $\mathbf{v}$ are returned, we then transform the samples as $\mathbf{u}=L\mathbf{v}$.

We use Stan for MCMC sampling. Stan implements Hamiltonian Monte Carlo (HMC) with advances such as the No-U-Turn Sampler (NUTS) (Carpenter et al., 2017; Homan and Gelman, 2014), which automatically ‘tunes‘ the HMC parameters. Stan can generally achieve a good effective sample size per unit time and so lead to better convergence of the MCMCML algorithm than standard MCMC algorithms.

Each iteration of the MCMC sampling must evaluate the log-likelihood (and its gradient). Since the elements of both $\mathbf{y}$ and $\mathbf{v}$ are all independent in this version of the model, evaluating the log-likelihood is parallelisable within a single chain, which significantly improves running time.

There are several methods to generate new estimates of $\beta$ and $\phi$ conditional on $\mathbf{y}$, $\mathbf{u}^{(i)}$, and $\phi^{(i)}$. We aim to minimise the negative log-likelihood:

There are three different approaches:

• •

  Monte Carlo Newton Raphson (MCNR): If the linear predictor is linear and additive in the data and parameters, then we can use the explicit formulation:

  $$\beta^{(i+1)}=\beta^{(i)}+E_{u}\left[X^{T}W(\theta^{(i)},\mathbf{u})X\right]^{-1}X^{T}\left(E_{u}\left[W(\theta^{(i)},\mathbf{u})\frac{\partial h^{-1}(\boldsymbol{\eta})}{\partial\boldsymbol{\eta}}(\mathbf{y-\boldsymbol{\mu}}(\beta^{(i)},\mathbf{u}))|\mathbf{y}\right]\right)$$

  where the expectations are with respect to the random effects. We estimate the expectations using the realised samples from the MCMC step of the algorithm. For many GLMMs, the likelihoods can be multimodal Searle et al. (1992), which means the MCNR algorithm may fail to converge at all.

• •

  Monte Carlo Quasi-Newton (MCQN): As the gradient is relatively inexpensive to calculate, one can make use of quasi-Newton algorithms, including L-BFGS.

• •

  Monte Carlo Expectation Maximisation (MCEM): We use the BOBYQA algorithm, a numerical, derivative free optimser that allows bounds for expectation maximisation (Powell, 2009).

The final step of each iteration is to generate new estimates of $\theta$ given the samples of the random effects. We aim to minimise

The multivariate Gaussian density is:

where $L_{jj}$ is the $j$th diagonal element of $L$.

As described in Section 2.5, we can take two approaches to the Cholesky factorisation of matrix $D$, either exploiting the block diagonal structure or using sparse matrix methods. In both cases the quadratic form can be rewritten as $\mathbf{u}^{T}D^{-1}\mathbf{u}=\mathbf{u}^{T}L^{-T}L^{-1}\mathbf{u}=\mathbf{z}^{T}\mathbf{z}$, so we can obtain $\mathbf{z}$ by solving the linear system $L\mathbf{z}=\mathbf{u}$ using forward substitution, which, if $L$ is also sparse provides efficiency gains. The log determinant is $\log(|D|)=\sum\log(\text{diag}(L))$. In the case of using a block diagonal formulation, we calculate the form:

where $r_{b}$ and $s_{b}$ indicates the start and end indexes of the vector, respectively, corresponding to block $b$.

The gradient of the log-likelihood with respect to $\theta$ is:

The quadratic form can similarly be solved by applying forward and backward substitution. The L-BFGS-B algorithm can be used then to minimise the log-likelihood. However, the inversion of $D$ and calculation of its partial derivatives may make the gradient computationally demanding. A derivative-free algorithm, BOBYQA, while requiring more function evaluations may be more efficient at solving for the parameter values. Both methods are available.

The efficiency of the MCML algorithms can be improved in two ways: first, one can dynamically alter the number of MCMC samples per step, which is discussed in the next section, and second, one can make use of MCMC sample information from previous iterations, reducing the number of samples per iteration. Stochastic Approximation Expectation Maximisation (SAEM) is a Ruppert-Monroe algorithm that can be used to approximate the log likelihoods (8) and (9) (Jank, 2006). For $\beta$ the log-likelihood on the $i+1$th iteration is:

and equivalently for $\theta$, where $\gamma\propto\left(\frac{1}{i}\right)^{\alpha}$ for $\alpha\in[0.5,1)$. This approach makes use of the all the MCMC samples from every iteration, where the contribution of previous values diminish over the iterations depending on the parameter $\alpha$. Larger values of $\alpha$ result in faster ‘forgetting’. For some applications, Ruppert-Polyak averaging can further improve convergence of the algorithms (Polyak and Juditsky, 1992), where the estimate of the log-likelihood on iteration $i$ is:

We provide two types of rule for terminating the algorithm. First, one can monitor the differences in the parameter estimates between successive iterations, and terminate the algorithm when the largest difference falls below some tolerance. However, this method may lead to premature termination. A running mean could prevent such issues, however, a perhaps preferable method is to assess the changes in the estimates of the log-likelihood:

An expectation maximisation algorithm is guaranteed to improve the objective function value over time. However, the estimates of the log-likelihood are subject to Monte Carlo error and so $\Delta\hat{\mathcal{L}}^{(i)}$ may increase or decrease at convergence. Caffo et al. (2005) proposed a stopping criteria that monitors the probability of convergence. The estimated upper bound for the difference is:

where $z_{\gamma}$ is the $\gamma$ quantile of a standard Gaussian distribution and $\hat{\sigma}_{\Delta}$ is the estimated standard deviation of $\Delta\hat{\mathcal{L}}^{(i)}$. This upper bound will be negative with probability $\gamma$ at convergence.

Caffo et al. (2005) use a similar argument to dynamically alter the number of MCMC samples on each iteration to acheive convergence with probability $\gamma$. The propose:

We provide MCML algorithm with and without dynamic sample sizes.

The member function \codeMCML fits the model. There are multiple options to control the algorithm and set its parameters; most options are determined by the relevant parameters of the \codeModel object. The default fitting method is SAEM without Ruppert-Polyak averaging, and with the log-likelihood stopping rule with probability 0.95, as this has been found to provide the most efficient and reliable model fitting in many tests.

We also provide model fitting using a Laplace approximation to the likelihood with the \codeModel member function \codeLA(). We include this approximation as it often provides a faster means of fitting the models in this package. These estimates may be of interest in their own right, but can also serve as useful starting values for the MCML algorithm. We approximate the log-likelihood of the re-parameterised model given in Equation (7) as:

As before, we use the first-order approximation (with or without attenutation), $W^{-1}(\theta,\mathbf{u})=\text{diag}\left(\left(\frac{\partial h^{-1}(\boldsymbol{\eta})}{\partial\boldsymbol{\eta}}\right)^{2}\text{Var}(\mathbf{y}|\mathbf{u})\right)$. We iterate model fitting to obtain estimates of $\beta$, $\mathbf{v}$ and $\phi$ and then to then obtain estimates of $\theta$ until the values converge. We provide two iterative algorithms for model fitting, using a similar pattern to the MCNR, and MCEM algorithms. First, the values of $\beta$ and $\mathbf{v}$ are updated. We can either use a scoring algorithm approach, or numerical optimisation. For the former method the score system we use is: