(Version 0.2.2) R package to support the specification of generalised linear mixed models using the R6 objectorientated class system.
A generalised linear mixed model (GLMM) has a mean function for observation \(i\) is
\[ \mu_i = \mathbf{x}_i\beta + \mathbf{z}_i \mathbf{u} \]
where \(\mathbf{x}_i\) is the \(i\)th row of matrix \(X\), which is a \(n \times P\) matrix of covariates, \(\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, \(\mathbf{u} \sim N(0,D)\), where \(D\) is the \(Q \times Q\) covariance matrix of the random effects terms that depends on parameters \(\theta\), and \(\mathbf{\mu}\) is the \(n\)length vector of mean values. The assumed data generating process for the study is
\[ y_i \sim G(h(\mu_i);\phi) \]
where \(\mathbf{y}\) is a \(n\)length vector of outcomes \(y_i\), \(G\) is a distribution, \(h(.)\) is the link function, and \(\phi\) additional scale parameters to complete the specification.
The package includes the function nelder()
, which we use
to generate data for the examples below. Nelder (1965) suggested a
simple notation that could express a large variety of different blocked
designs. The notation was proposed in the context of splitplot
experiments for agricultural research, where researchers often split
areas of land into blocks, subblocks, 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
that generates a data frame of a design using the notation.
There are two operations: * >
(or \(\to\) in Nelder’s notation) indicates
“clustered in”. * *
(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
abc(12)
. So for example ~cl(4) > ind(5)
means in each of five levels of cl
there are five levels of
ind
, and the individuals are different between clusters.
The formula ~cl(4) * t(3)
indicates that each of the four
levels of cl
are observed for each of the three levels of
t
. Brackets are used to indicate the order of
evaluation.
The specification of a covariance object requires three inputs: a
formula, data, and parameters. A new instance of each class can be
generated with the $new()
function, for example
Covariance$new(...)
.
A covariance function is specified as an additive formula made up of
components with structure (1f(j))
. The left side of the
vertical bar specifies the covariates in the model that have a random
effects structure. The right side of the vertical bar specify the
covariance function f
for that term using variable named in
the data j
. Covariance functions on the right side of the
vertical bar are multiplied together, i.e., (1f(j)*g(t))
.
The table below shows the currently implemented covariance functions
Function  \(Cov(x_i,x_{i'})\)  \(\theta\)  Code 

Identity/ Group membership  \(\theta_1^2 \mathbf{1}(x_i = x_{i'})\)  \(\theta_1 > 0\)  gr(x) 
Exponential  \(\theta_1 \text{exp}( \vert x_i  x_{i'}\vert / \theta_2 )\)  \(\theta_1,\theta_2 > 0\)  fexp(x) 
\(\text{exp}( \vert x_i  x_{i'}\vert /\theta_1)\)  \(\theta_1 > 0\)  fexp0(x) 

Squared Exponential  \(\theta_1 \text{exp}( (\vert x_i  x_{i'}\vert / \theta_2)^2)\)  \(\theta_1,\theta_2 > 0\)  sqexp(x) 
\(\text{exp}(( \vert x_i  x_{i'}\vert/\theta_1)^2 )\)  \(\theta_1 > 0\)  sqexp0(x) 

Autoregressive order 1  \(\theta_1^{\vert x_i  x_{i'} \vert}\)  \(0 < \theta_1 < 1\)  ar1(x) 
Bessel  \(K_{\theta_1}(x)\)  \(\theta_1\) > 0  bessel(x) 
Matern  \(\frac{2^{1\theta_1}}{\Gamma(\theta_1)}\left( \sqrt{2\theta_1}\frac{x}{\theta_2} \right)^{\theta_1} K_{\theta_1}\left(\sqrt{2\theta_1}\frac{x}{\theta_2})\right)\)  \(\theta_1,\theta_2 > 0\)  matern(x) 
Compactly supported*  
Wendland 0  \(\theta_1(1y)^{\theta_2}, 0 \leq y \leq 1; 0, y \geq 1\)  \(\theta_1>0, \theta_2 \geq (d+1)/2\)  wend0(x) 
Wendland 1  \(\theta_1(1+(\theta_2+1)y)(1y)^{\theta_2+1}, 0 \leq y \leq 1; 0, y \geq 1\)  \(\theta_1>0, \theta_2 \geq (d+3)/2\)  wend1(x) 
Wendland 2  $_1(1+(_2+2)y + ((_2+2)^2  1)y^{2)(1y)}{_2+2}, 0 y $  \(\theta_1>0,\theta_2 \geq (d+5)/2\)  wend1(x) 
\(0, y \geq 1\)  
WhittleMatern \(\times\) Wendland**  \(\theta_1\frac{2^{1\theta_2}}{\Gamma(\theta_2)}y^{\theta_2}K_{\theta_2}(y)(1+\frac{11}{2}y + \frac{117}{12}y^2)(1y), 0 \leq y \leq 1; 0, y \geq 1\)  \(\theta_1,\theta_2 > 0\)  prodwm(x) 
Cauchy \(\times\) Bohman***  \(\theta_1(1y^{\theta_2})^{3}\left( (1y)\text{cos}(\pi y)+\frac{1}{\pi}\text{sin}(\pi y) \right), 0 \leq y \leq 1; 0, y \geq 1\)  \(\theta_1>0, 0 \leq \theta_2 \leq 2\)  prodcb(x) 
Exponential \(\times\) Kantar****  \(\theta_1\exp{(y^{\theta_2})}\left( (1y)\frac{\sin{(2\pi y)}}{2\pi y} + \frac{1}{\pi}\frac{1\cos{(2\pi y)}}{2\pi y} \right), 0 \leq y \leq 1\)  \(\theta_1,\theta_2 > 0\)  prodek(x) 
\(0, y \geq 1\) 
\(\vert . \vert\) is the Euclidean
distance. \(K_a\) is the modified
Bessel function of the second kind. Variable \(y\) is defined as \(x/r\) where \(r\) is the effective range. For the
compactly supported functions \(d\) is
the number of dimensions in x
. Permissible in one
or two dimensions. Only permissible in one dimension.
****Permissible in up to three dimensions.
One combines functions to provide the desired covariance function.
For example, for a steppedwedge cluster trial we could consider the
standard specification with an exchangeable random effect for the
cluster level, and a separate exchangeable random effect for the
clusterperiod, which would be ~(1gr(j))+(1gr(j,t))
or
~(1gr(j))+(1gr(j)*gr(t))
. Alternatively, we could
consider an autoregressive clusterlevel random effect that decays
exponentially over time so that, for persons \(i\) in cluster \(j\) at time \(t\), \(Cov(y_{ijt},y_{i'jt}) = \theta_1\), for
\(i\neq i'\), \(Cov(y_{ijt},y_{i'jt'}) = \theta_1
\theta_2^{\vert tt' \vert}\) for \(t \neq t'\), and \(Cov(y_{ijt},y_{i'j't}) = 0\) for
\(j \neq j'\). This function would
be specified as ~(1gr(j)*ar1(t))
.
Parameters are provided to the covariance function as a vector. The covariance functions described in the Table have different parameters \(\theta\), and a value is required to be provided to generate the matrix \(D\) and related objects for analyses and which serve as starting values for model fitting. The elements of the vector correspond to each of the functions in the covariance formula in the order they are written.
A full call to create a new covariance object is:
R> df < nelder(~ (j(10)* t(5)) > ind(10))
R> cov < Covariance$new(formula = ~(1gr(j)*ar1(t)),
R> parameters = c(0.05,0.8),
R> data= df)
in this call, the parameters
are optional, and if
provided as a list of arguments to a Model
object (see
below), then the data
argument is also optional. A
compactly supported function is used, then the effective range
parameters should be provided in the order the function appears in the
formula.
R> cov < Covariance$new(formula = ~(1prodwm(x,y)),
R> parameters = c(0.25,0.5),
R> eff_range = 0.5,
R> data= df)
Specification of the mean function follows standard model formulae in R. A vector of values of the mean function parameters is required to complete the model specification along with the distribution as a standard R family object. A complete specification is thus:
R> mf < MeanFunction$new(formula = ~ factor(t)+ int  1,
R> data = df,
R> parameters = rep(0,6),
R> family = gaussian())
As before, the parameters
, data
, and
family
are optional and can instead be provided directly to
the Model
call below. Note that factor
in this
function does not drop one level, unlike standard R formulae, so
removing the intercept is required to prevent a collinearity
problem.
A model can be created by specifying a Covariance
and
MeanFunction
object:
R> model < Model$new(covariance = cov,
R> mean = mf,
R> var_par = 1)
Alternatively, we can provide a list of arguments to the
covariance
and mean
arguments:
R> model < Model$new(covariance = list(formula = ~(1gr(j)*ar1(t))),
R> mean = list(formula = ~ factor(t)+ int  1),
R> data = df,
R> family = gaussian(),
R> var_par = 1)
where, as required, parameters can be supplied to covariance and mean function argument lists.
For Gaussian models, and other distributions requiring an additional
scale parameter \(\phi\), one must also
specify the option var_par
which is the conditional
variance \(\phi = \sigma\) at the
individual level. The default value is 1. Alternatively, one can specify
a design by providing the list of arguments directly to
covariance
and mean.function
instead of model
objects.
The package and associated packages (glmmrMCML
and
glmmrOptim
) currently support the following families and
link functions  Family  Link functions  ——–————————  Gaussian 
Identity, log   Binomial  Logit, log, identity   Poisson  Log,
Identity   Gamma  Log, Inverse, Identity  Beta  Logit 
The Beta family is provided by the package function
Beta()
, which generates a barebones list specifying the
family and link. We use a meanvariance parameterisation of the Beta
family. The likelihood is:
\[ f(y_i  \mu_i, \phi) = \frac{y_i^{\mu_i\phi  1}(1y_i)^{(1\mu_i)\phi  1}}{B(\mu_i\phi, (1\mu_i)\phi)} \]
where \(B()\) is the Beta function, and we use logit link
\[ \log\left( \frac{\mu_i}{1\mu_i} \right) = \mathbf{x}_i\beta + \mathbf{z}_i \mathbf{u} \]
We similarly use a meanvariance parameterisation for the Gamma regression function:
\[ f(y_i  \mu_i, \nu) = \frac{1}{\Gamma(\nu)}\left( \frac{\nu y_i}{\mu_i} \right)^\nu \exp \left( \frac{\nu y_i}{\mu_i} \right) \frac{1}{y} \]
where we also provide logit, inverse, and identity link functions for the specification of \(\mu_i\).
Each class holds associated matrices and has member functions to
compute basic summaries and analyses. The Matrix
package is
used for matrix operations in R. For example, a Covariance
object holds the matrix \(D\)
R> cov$D[1:10,1:10]
10 x 10 sparse Matrix of class "dsCMatrix"
[1,] 0.002500 0.00200 0.0016 0.00128 0.001024 . . . . .
[2,] 0.002000 0.00250 0.0020 0.00160 0.001280 . . . . .
[3,] 0.001600 0.00200 0.0025 0.00200 0.001600 . . . . .
[4,] 0.001280 0.00160 0.0020 0.00250 0.002000 . . . . .
[5,] 0.001024 0.00128 0.0016 0.00200 0.002500 . . . . .
[6,] . . . . . 0.002500 0.00200 0.0016 0.00128 0.001024
[7,] . . . . . 0.002000 0.00250 0.0020 0.00160 0.001280
[8,] . . . . . 0.001600 0.00200 0.0025 0.00200 0.001600
[9,] . . . . . 0.001280 0.00160 0.0020 0.00250 0.002000
[10,] . . . . . 0.001024 0.00128 0.0016 0.00200 0.002500
glmmrBase
in other packagesThis package provides a foundation for other packages providing
analysis or estimation of generalised linear mixed models. For example,
we have the glmmrMCML
package, which provides Markov Chain
Monte Carlo Maximum Likelihood estimations for these models, and
glmmrOptim
, which finds approximate optimal designs based
on these models. New classes can be defined that inherit from the
classes included in this package. glmmrMCML
defines the
modelMCML
class that adds the member function
MCML
. Then the new functions can access the model elements,
such as covariance matrices, and benefit from automatic updating of
these elements when specifications or parameters change. As an example
we can define a new class that has a member function that returns the
determinant of the matrix \(D\):
R> CovDet < R6::R6Class("CovDet",
R> inherit = Covariance,
R> public = list(
R> det = function(){
R> return(Matrix::determinant(self$D))
R> }))
R> cov < CovDet$new(formula = ~(1gr(j)*ar1(t)),
R> parameters = c(0.05,0.8),
R> data= df)
R> cov$det()
$modulus
[1] 340.4393
attr(,"logarithm")
[1] TRUE
$sign
[1] 1
attr(,"class")
[1] "det"