(Version 0.2.2) R package to support the specification of generalised linear mixed models using the R6 object-orientated 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 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
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 (1|f(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., (1|f(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(1-y)^{\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)(1-y)^{\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)y2)(1-y){_2+2}, 0 y $ | \(\theta_1>0,\theta_2 \geq (d+5)/2\) | wend1(x) |
| \(0, y \geq 1\) | |||
| Whittle-Matern \(\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)(1-y), 0 \leq y \leq 1; 0, y \geq 1\) | \(\theta_1,\theta_2 > 0\) | prodwm(x) |
| Cauchy \(\times\) Bohman*** | \(\theta_1(1-y^{\theta_2})^{-3}\left( (1-y)\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( (1-y)\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 stepped-wedge 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
cluster-period, which would be ~(1|gr(j))+(1|gr(j,t)) or
~(1|gr(j))+(1|gr(j)*gr(t)). Alternatively, we could
consider an autoregressive cluster-level 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 t-t' \vert}\) for \(t \neq t'\), and \(Cov(y_{ijt},y_{i'j't}) = 0\) for
\(j \neq j'\). This function would
be specified as ~(1|gr(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 = ~(1|gr(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 = ~(1|prodwm(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 = ~(1|gr(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 mean-variance parameterisation of the Beta
family. The likelihood is:
\[ f(y_i | \mu_i, \phi) = \frac{y_i^{\mu_i\phi - 1}(1-y_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 mean-variance 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 = ~(1|gr(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"