∎

11institutetext: Thomas Maullin-Sapey 22institutetext: Big Data Institute, Li
Ka Shing Centre for Health Information and Discovery, Old Road Campus, Oxford,
OX3 7LF, UK
22email<EMAIL_ADDRESS>
ORCID ID: 0000-0002-1890-330X33institutetext: Thomas E. Nichols
44institutetext: Big Data Institute, Li Ka Shing Centre for Health Information
and Discovery, Old Road Campus, Oxford, OX3 7LF, UK
44email<EMAIL_ADDRESS>
ORCID ID: 0000-0002-4516-5103

# Fisher Scoring for crossed factor Linear Mixed Models

Thomas Maullin-Sapey  Thomas E. Nichols

(Received: date / Accepted: date)

###### Abstract

The analysis of longitudinal, heterogeneous or unbalanced clustered data is of
primary importance to a wide range of applications. The Linear Mixed Model
(LMM) is a popular and flexible extension of the linear model specifically
designed for such purposes. Historically, a large proportion of material
published on the LMM concerns the application of popular numerical
optimization algorithms, such as Newton-Raphson, Fisher Scoring and
Expectation Maximization to single-factor LMMs (i.e. LMMs that only contain
one “factor” by which observations are grouped). However, in recent years, the
focus of the LMM literature has moved towards the development of estimation
and inference methods for more complex, multi-factored designs. In this paper,
we present and derive new expressions for the extension of an algorithm
classically used for single-factor LMM parameter estimation, Fisher Scoring,
to multiple, crossed-factor designs. Through simulation and real data
examples, we compare five variants of the Fisher Scoring algorithm with one
another, as well as against a baseline established by the R package lmer, and
find evidence of correctness and strong computational efficiency for four of
the five proposed approaches. Additionally, we provide a new method for LMM
Satterthwaite degrees of freedom estimation based on analytical results, which
does not require iterative gradient estimation. Via simulation, we find that
this approach produces estimates with both lower bias and lower variance than
the existing methods.

###### Keywords:

Fisher Scoring Linear Mixed Model Crossed Factors

##  1 Introduction

###  1.1 Background

Since its conception in the seminal work of Laird and Ware (1982), the
literature on Linear Mixed Model (LMM) estimation and inference has evolved
rapidly. At present, many software packages exist which are capable of
performing LMM estimation and inference for large and complex LMMs in an
incredibly quick and memory-efficient manner. For some packages, this
exceptional speed and efficiency arises from simplifying model assumptions,
whilst for others, complex mathematical operations such as sparse matrix
methodology and sweep operators are utilized to improve performance (Wolfinger
et al. (1994), Bates et al. (2015)).

However, situations exist in which the current methodology cannot be applied
due to practical implementation considerations. One such example is given in
the field of medical imaging, in which mass-univariate analyses, where
hundreds of thousands of models must be estimated concurrently, are standard
practice. To efficiently perform a mass-univariate analysis within a practical
time-frame, the use of vectorized computation, which exploits the repetitive
nature of simplistic operations to streamline calculation, must be employed
(Smith and Nichols (2018), Li et al. (2019)). Unfortunately, many existing LMM
tools utilize complex operations, for which vectorized support does not
currently exist. As a result, alternative methodology, using more conceptually
simplistic mathematical operations for which vectorized support exists, is
required.

The complex nature of LMM computation has partly arisen from the gradual
expansion of the definition of “Linear Mixed Model”. Previously, the term
“Linear Mixed Model” was primarily used to refer to an extension of the linear
regression model which contains random effects grouped by a single random
factor. Examples of this definition can be seen in Laird and Ware (1982), in
which “Linear Mixed Model” refers only to “single-factor” longitudinal models,
and in Lindstrom and Bates (1988), where more complex, multi-factor models are
described as an “extension” of the “Linear Mixed Model”. Consequently,
throughout the late 1970s and 1980s, one of the main focuses of the LMM
literature was to provide parameter estimation methods, such as Fisher
Scoring, Newton-Raphson and Expectation Maximization, for the single-factor
LMM (e.g. Dempster et al. (1977), Jennrich and Schluchter (1986) and Laird et
al. (1987)). By exploiting structural features of the single-factor model,
implementation of these methods required only conceptually simplistic
mathematical operations. For instance, the Fisher Scoring algorithm proposed
by Jennrich and Schluchter (1986) relies only upon vector addition and matrix
multiplication, inversion, and reshaping operations.

In recent times, usage of the term “Linear Mixed Models” has grown
substantially to include models which contain random effects grouped by
multiple random factors. Examples of this more general definition are found in
Pinheiro and Bates (2009) and Tibaldi et al. (2007). For this more general
“multi-factor” LMM definition, models can be described as exhibiting either a
hierarchical factor structure (i.e. factor groupings are nested inside one
another) or crossed factor structure (i.e. factor groupings are not nested
inside one another). For instance, a study involving students grouped by the
factors “school” and “class” contains hierarchical factors (as every class
belongs to a specific school). In contrast, a study that involves observations
grouped by ”subject” and ”location,” where subjects change location between
recorded measurements, contains crossed factors (as each subject is not
generally associated with a specific location or vice versa). In either case,
due to the more complex nature of the multi-factor LMM, the computational
approaches used in the single-factor setting can not be directly applied to
the multi-factor model. For this reason, parameter estimation of the multi-
factor LMM has often been viewed as a much more “difficult” problem than its
single-factor counterpart (see, for example, the discussion in chapters $2$
and $8$ of West et al. (2014)).

Several authors have proposed and implemented methods for multi-factor LMM
estimation. However, such methods typically require conceptually complex
mathematical operations, which are not naturally amenable to vectorized
computation, or restrictive simplifying model assumptions, which prevent some
designs from being estimated. For example, the popular R package lmer performs
LMM estimation via minimization of a penalized least squares cost function
based on a variant of Henderson’s mixed model equations (Henderson et al.
(1959), Bates et al. (2015)). However, sparse matrix methods are required to
achieve fast evaluation of the cost function, and advanced numerical
approximation methods are needed for optimization (e.g. the Bound Optimization
BY Quadratic Approximation, BOBYQA, algorithm, Powell (2009)). The commonly
used SAS and SPSS packages, PROC-MIXED and MIXED, employ a Newton-Raphson
algorithm proposed by Wolfinger, Tobias, and Sall (1994). In this approach,
though, derivatives and Hessians must be computed using the sweep operator,
W-transformation, and Cholesky decomposition updating methodology (SAS
Institute Inc. (2015), IBM Corp (2015)). The Hierarchical Linear Models (HLM)
package takes an alternative option and restricts its inputs to only LMMs
which contain hierarchical factors (Raudenbush and Bryk (2002)); although, the
Hierarchical Cross-classified Models (HCM) sub-module does make allowances for
specific use-case crossed LMMs. To perform parameter estimation, HLM employs a
range of different methods, each tailored to a particular model design. As a
result, there are many models for which HLM does not provide support.

Methods have also been proposed for evaluating the score vector (the
derivative of the log-likelihood function) and the Fisher Information matrix
(the expected Hessian of the negative log-likelihood function) required to
perform Fisher Scoring for the multi-factor LMM. For example, alongside the
Newton-Raphson approach adopted by PROC-MIXED and MIXED, Wolfinger et al.
(1994) also describe a Fisher Scoring algorithm. However, evaluation of the
expressions they provide requires use of the sweep operator, W-transformation,
and Cholesky decomposition updating methodology; operations for which
widespread vectorized support do not yet exist. More recently, expressions for
score vectors and Fisher Information matrices were provided by Zhu and Wathen
(2018). However, the approach Zhu and Wathen (2018) adopt to derive these
expressions produces an algorithm that requires independent computation for
each variance parameter in the LMM. This algorithm’s serialized nature results
in substantial overheads in terms of computation time, thus limiting the
method’s utility in practical situations where time-efficiency is a crucial
consideration.

To our knowledge, no approach has yet been provided for multi-factor LMM
parameter estimation, which utilizes only simplistic, universally available
operations which are naturally amenable to vectorized computation. In this
work we revisit a Fisher Scoring approach suggested for the single-factor LMM,
described in Demidenko (2013), extending it to the multi-factor setting, with
the intention of revisiting the motivating example of the mass-univariate
model in later work.

The novel contribution of this work is to provide new derivations and closed-
form expressions for the score vector and Fisher Information matrix of the
multi-factor LMM. We show how these expressions can be employed for LMM
parameter estimation via the Fisher Scoring algorithm and can be further
adapted for constrained optimization of the random effects covariance matrix.
Additionally, we demonstrate that such closed-form expressions have use beyond
the scope of the Fisher Scoring algorithm. An example of this is provided in
the setting of mixed model inference, where the degrees of freedom of the
approximate T-statistic are not known and must be estimated (Verbeke and
Molenberghs (2001)). We show how our methods may be combined with the
Satterthwaite-based method for approximating degrees of freedom for the LMM by
using an approach based on the work of Kuznetsova et al. (2017).

In this paper, we first propose five variants of the Fisher Scoring algorithm.
Following this, we provide a discussion of initial starting values for the
algorithm, approaches to ensuring positive definiteness of the random effects
covariance matrix, and methods for improving the algorithm’s computational
efficiency during implementation. Detail on constrained optimization, allowing
for structural assumptions to be placed on the random effects covariance
matrix, is then provided. Proceeding this, new expressions for Satterthwaite
estimation of the degrees of freedom of the approximate T-statistic for the
multi-factor LMM are given. Finally, we verify the correctness of the proposed
algorithms and degrees of freedom estimates via simulation and real data
examples, benchmarking the performance against that of the R package lmer.

###  1.2 Preliminaries

In this section, we introduce preliminary statistical concepts and notation
which will be employed throughout the remainder of this work. In Section
1.2.1, a formal definition of the LMM is provided and in Section 1.2.2,
additional notation is detailed.

####  1.2.1 The model

Both the single-factor and multi-factor LMM, for $n$ observations, take the
following form:

$Y=X\beta+Zb+\epsilon\newline $ (1) $\epsilon\sim N(0,\sigma^{2}I_{n}),\hskip
8.5359ptb\sim N(0,\sigma^{2}D),$ (2)

where the known quantities are:

  * $\bullet$

$Y$: the $(n\times 1)$ vector of responses.

  * $\bullet$

$X$: the $(n\times p)$ fixed effects design matrix.

  * $\bullet$

$Z$: the $(n\times q)$ random effects design matrix.

And the unknown parameters are:

  * $\bullet$

$\beta$: the $(p\times 1)$ fixed effects parameter vector.

  * $\bullet$

$\sigma^{2}$: the scalar fixed effects variance.

  * $\bullet$

$D$: the $(q\times q)$ random effects covariance matrix.

From (1) and (2) the marginal distribution of the response vector, $Y$ is:

$Y\sim N(X\beta,\sigma^{2}(I_{n}+ZDZ^{\prime})).$

The log-likelihood for the LMM specified by (1) and (2) is derived from the
marginal distribution of $Y$. Dropping constant terms, this log-likelihood is
given by:

$l(\theta)=-\frac{1}{2}\bigg{\\{}n\log(\sigma^{2})+\sigma^{-2}e^{\prime}V^{-1}e+\log|V|\bigg{\\}},$
(3)

where $\theta$ is shorthand for all the parameters $(\beta,\sigma^{2},D)$,
$V=I_{n}+ZDZ^{\prime}$ and $e=Y-X\beta$. Throughout the main body of this
work, we shall consider parameter estimation performed via Maximum Likelihood
(ML) estimation of (3). However, we note that the approaches we describe for
ML estimation can also be easily adapted to use a Restricted Maximum
Likelihood (ReML) criterion. Further detail on ReML estimation is provided in
Appendix 6.4.

The distinction between the multi-factor and single-factor LMM lies in the
specification of the random effects in the model. Random effects are often
described in terms of factors, categorical variables that group the random
effects, and levels, individual instances of such a categorical variable. At
this point, we highlight that the term “factor”, in this work, refers only to
categorical variables which group random effects and does not refer to
groupings of fixed effects. We denote the total number of factors in the model
as $r$ and denote the $k^{th}$ factor in the model as $f_{k}$ for
$k\in\\{1,...,r\\}$. For a given factor $f_{k}$, $l_{k}$ will be used to
denote the number of levels possessed by $f_{k}$, and $q_{k}$ the number of
random effects which $f_{k}$ groups. The single-factor LMM corresponds to the
case $r=1$, while the multi-factor setting corresponds to the case $r>1$.

An example of how this notation may be used in practice is given as follows.
Suppose an LMM contains observations that are grouped according to ‘subject’
(i.e. the participant whose observation was recorded) and ‘location’ (i.e. the
place the observation was recorded). Further, suppose that the LMM includes a
random intercept and random slope which each model subject-specific behaviour,
and a random intercept which models location-specific behaviour. Two factors
are present in this design: the factor $f_{1}$ is ‘subject’ and the factor
$f_{2}$ is ‘location’. Therefore, $r=2$. The number of subjects is $l_{1}$ and
the number of locations is $l_{2}$. The number of covariates grouped by the
first factor, $q_{1}$, is $2$ (i.e. the random intercept and the random slope)
and the number grouped by the second factor, $q_{2}$, is $1$ (i.e. the random
intercept).

The values of $r$, $\\{q_{k}\\}_{k\in\\{1,\ldots,r\\}}$ and
$\\{l_{k}\\}_{k\in\\{1,\ldots,r\\}}$ determine the structure of the random
effects design matrix, $Z$, and random effects covariance matrix, $D$. To
specify $Z$ formally is notationally cumbersome and of little relevance to the
aims of this work. For this reason, the reader is referred to the work of
Bates et al. (2015) for further detail on the construction of $Z$. Here it
will be sufficient to note that, under the assumption that it’s columns are
appropriately ordered, $Z$ will be comprised of $r$ horizontally concatenated
blocks. The $k^{th}$ block of $Z$, denoted $Z_{(k)}$, has dimension $(n\times
l_{k}q_{k})$ and describes the random effects which are grouped by the
$k^{th}$ factor. Additionally, each block, $Z_{(k)}$, can be further
partitioned column-wise into $l_{k}$ blocks of dimension $(n\times q_{k})$.
The $j^{th}$ block of $Z_{(k)}$, denoted $Z_{(k,j)}$, corresponds to the
random effects which belong to the $j^{th}$ level of the $k^{th}$ factor. In
summary,

$\displaystyle Z=[Z_{(1)},Z_{(2)},\ldots Z_{(r)}],$ $\displaystyle
Z_{(k)}=[Z_{(k,1)},Z_{(k,2)},\ldots Z_{(k,l_{k})}]\hskip 14.22636pt\text{(for
}k\in\\{1,\ldots,r\\}\text{)}$

An important property of the matrix $Z$ is that for arbitrary $k\leq r$,
$Z_{(k,i)}$ and $Z_{(k,j)}$ satisfy the below condition for all $i\neq j$:

$\displaystyle\text{ If row }w\text{ of }Z_{(k,i)}\text{ contains any non-zero
values }$ (4) $\displaystyle\implies\text{then row }w\text{ of
}Z_{(k,j)}\text{ contains only zero values. }$

A consequence of the above property is that for any arbitrary factor $f_{k}$,
the rows of $Z_{(k)}$ can be permuted in order to obtain a block diagonal
matrix. As, for the single-factor LMM, $Z\equiv Z_{(1)}$, it follows that the
observations of the single-factor LMM can be arranged such that $Z$ is block
diagonal. This feature of the single-factor LMM simplifies the derivation of
the Fisher Information matrix and score vector required for Fisher Scoring.
However, this simplification cannot be generalized to the multi-factor LMM. In
general, it is not true that the rows of $Z$ can be permuted in such a way
that the resultant matrix is block-diagonal. As emphasized in Section 1.1, due
to this, many of the results derived in the single-factor LMM have not been
generalized to the multi-factor setting.

To describe the random effects covariance matrix, $D$, it is assumed that
factors are independent from one another and that for each factor, factor
$f_{k}$, there is a $(q_{k}\times q_{k})$ unknown covariance matrix, $D_{k}$,
representing the “within-factor” covariance for the random effects grouped by
$f_{k}$. The random effects covariance matrix, $D$, appearing in (2), can now
be given as:

$D=\bigoplus_{k=1}^{r}(I_{l_{k}}\otimes D_{k}),$

where $\oplus$ represents the direct sum, and $\otimes$ the Kronecker product.
Note that, whilst $D$ is large in dimension (having dimension $(q\times q)$
where $q=\sum_{i}q_{i}l_{i}$), $D$ contains
$q_{u}=\sum_{i}\frac{1}{2}q_{i}(q_{i}+1)$ unique elements. Typically, it is
true that $q_{u}<<q^{2}$. As a result, the task of describing the random
effects covariance matrix $D$ reduces in practice to specifying only a small
number of parameters.

####  1.2.2 Notation

In this section, we introduce notation which will be used throughout the
remainder of this work. The conventional hat operator notation is adopted to
denote estimators resulting from likelihood maximisation procedures. For
example, $\beta$ represents the true (unknown) fixed effects parameter vector
for the LMM, whilst the maximum likelihood estimate of $\beta$ is denoted
$\hat{\beta}$. The conventional subscript notation, $A_{[X,Y]}$, is used to
denote the sub-matrix of matrix $A$, composed of all elements of $A$ with row
indices $x\in X$ and column indices $y\in Y$. The replacement of $X$ or $Y$
with a colon, $:$, represents all rows or all columns of $A$, respectively. If
a scalar, $x$ or $y$, replaces $X$ or $Y$, this represents the elements with
row indices $x=X$ or column indices $y=Y$, respectively. The notation $(k)$
may also replace $X$ and $Y$ where $(k)$ represents the indices of the columns
of $Z$ which correspond to factor $f_{k}$. Similarly, $X$ and $Y$ may be
substituted for the ordered pair $(k,j)$ where $(k,j)$ represents the indices
of the columns of $Z$ which correspond to level $j$ of factor $f_{k}$. We
highlight again our earlier notation, $Z_{(k,j)}$, which, due to its frequent
occurrence acts as a shorthand for $Z_{[:,(k,j)]}$, i.e. the columns of $Z$
corresponding to level $j$ of factor $f_{k}$.

Finally, we shall also adopt the notations ‘vec’, ‘vech’, $N_{k}$, $K_{m,n}$,
$\mathcal{D}_{k}$ and $\mathcal{L}_{k}$ as used in Magnus and Neudecker
(1980), defined as follows:

  * $\bullet$

‘vec’ represents the mathematical vectorizarion operator which transforms an
arbitrary $(k\times k)$ matrix, $A$, to a $(k^{2}\times 1)$ column vector,
vec$(A)$, composed of the columns of $A$ stacked into a column vector. (Note:
this is not to be confused with the concept of computational vectorization
discussed in Section 1.1).

  * $\bullet$

‘vech’ represents the half-vectorization operator which transforms an
arbitrary square matrix, $A$, of dimension $(k\times k)$ to a $(k(k+1)/2\times
1)$ column vector, vech$(A)$, composed by stacking the elements of $A$ which
fall on and below the diagonal into a column vector.

  * $\bullet$

$N_{k}$ is defined as the unique matrix of dimension $(k^{2}\times k^{2})$
which implements symmetrization for any arbitrary square matrix $A$ of
dimension $(k\times k)$ in vectorized form. I.e. $N_{k}$ satisfies the
following relation:

$N_{k}\text{vec}(A)=\text{vec}(A+A^{\prime})/2.$

  * $\bullet$

$K_{m,n}$ is the unique “Commutation” matrix of dimension $(mn\times mn)$,
which permutes, for any arbitrary matrix $A$ of dimension $(m\times n)$, the
vectorization of $A$ to obtain the vectorization of the transpose of $A$. I.e.
$K_{m,n}$ satisfies the following relation:

$\text{vec}(A)=K_{m,n}\text{vec}(A^{\prime}).$

  * $\bullet$

$\mathcal{D}_{k}$ is the unique “Duplication” matrix of dimension
$(k^{2}\times k(k+1)/2)$, which maps the half-vectorization of any arbitrary
symmetric matrix $A$ of dimension $(k\times k)$ to it’s vectorization. I.e.
$\mathcal{D}_{k}$ satisfies the following relation:

$\text{vec}(A)=\mathcal{D}_{k}\text{vech}(A).$

  * $\bullet$

$\mathcal{L}_{k}$ is the unique $1-0$ “Elimination” matrix of dimension
$(k(k+1)/2\times k^{2})$, which maps the vectorization of any arbitrary lower
triangular matrix $A$ of dimension $(k\times k)$ to it’s half-vectorization.
I.e. $\mathcal{L}_{k}$ satisfies the following relation:

$\text{vech}(A)=\mathcal{L}_{k}\text{vec}(A).$

To help track the notational conventions employed in this work, an index of
notation is provided in the Supplementary Material, Section $S10$.

##  2 Methods

###  2.1 Fisher Scoring Algorithms

In this section, we employ the Fisher Scoring algorithm for ML estimation of
the parameters $(\beta,\sigma^{2},D)$. The Fisher Scoring algorithm update
rule takes the following form:

$\theta_{s+1}=\theta_{s}+\alpha_{s}\mathcal{I}(\theta_{s})^{-1}\frac{dl(\theta_{s})}{d\theta}_{,}$
(5)

where $\theta_{s}$ is the vector of parameter estimates given at iteration
$s$, $\alpha_{s}$ is a scalar step size, the score vector of $\theta_{s}$,
$\frac{dl(\theta_{s})}{d\theta}$, is the derivative of the log-likelihood with
respect to $\theta$ evaluated at $\theta=\theta_{s}$, and
$\mathcal{I}(\theta_{s})$ is the Fisher Information matrix of $\theta_{s}$;

$\mathcal{I}(\theta_{s})=E\bigg{[}\bigg{(}\frac{dl(\theta)}{d\theta}\bigg{)}\bigg{(}\frac{dl(\theta)}{d\theta}\bigg{)}^{\prime}\bigg{|}\theta=\theta_{s}\bigg{]}.$

A more general formulation of Fisher Scoring, which allows for low-rank Fisher
Information matrices, is given by Rao and Mitra (1972):

$\theta_{s+1}=\theta_{s}+\alpha_{s}\mathcal{I}(\theta_{s})^{+}\frac{dl(\theta_{s})}{d\theta}_{,}$
(6)

where superscript plus, +, is the Moore-Penrose (or “Pseudo”) inverse. For
notational brevity, when discussing algorithms of the form (5) and (6) in the
following sections, the subscript $s$, representing iteration number, will be
suppressed unless its inclusion is necessary for clarity.

For the LMM, several different representations of the parameters of interest,
$(\beta,\sigma^{2},D)$, can be used for numerical optimization and result in
different Fisher Scoring iteration schemes. In this section, we consider the
following three representations for $\theta$:

$\theta^{h}=\begin{bmatrix}\beta\\\ \sigma^{2}\\\ \text{vech}(D_{1})\\\
\vdots\\\ \text{vech}(D_{r})\\\
\end{bmatrix},\theta^{f}=\begin{bmatrix}\beta\\\ \sigma^{2}\\\
\text{vec}(D_{1})\\\ \vdots\\\ \text{vec}(D_{r})\\\
\end{bmatrix},\theta^{c}=\begin{bmatrix}\beta\\\ \sigma^{2}\\\
\text{vech}(\Lambda_{1})\\\ \vdots\\\ \text{vech}(\Lambda_{r})\\\
\end{bmatrix},$

where $\Lambda_{k}$ represents the lower triangular Cholesky factor of
$D_{k}$, such that $D_{k}=\Lambda_{k}\Lambda_{k}^{\prime}$. We will refer to
the representations ($\theta^{h}$, $\theta^{f}$ and $\theta^{c}$) as the
“half”, “full” and “Cholesky” representations of $(\beta,\sigma^{2},D)$,
respectively. In a slight abuse of notation, the function $l$ will be allowed
to take any representation of $\theta$ as input, with the interpretation
unchanged (i.e. $l(\theta^{f})=l(\theta^{h})=l(\theta^{c})$). For example, if
the full representation is being used, the log-likelihood will be denoted
$l(\theta^{f})$, but if the half representation is being used, the log
likelihood will be denoted $l(\theta^{h})$.

In the following sections, the score vectors and Fisher Information matrices
required to perform five variants of Fisher Scoring for the multi-factor LMM
will be stated with proofs provided in Appendices 6.1 and 6.2. For notational
convenience, we denote the sub-matrix of the Fisher Information matrix of
$\theta^{h}$ with row indices corresponding to parameter vector $a$ and column
indices corresponding to parameter vector $b$ as $\mathcal{I}^{h}_{a,b}$. In
other words, $\mathcal{I}^{h}_{a,b}$ is the sub-matrix of
$\mathcal{I}(\theta^{h})$, defined by:

$\mathcal{I}^{h}_{a,b}=E\bigg{[}\bigg{(}\frac{dl(\theta^{h})}{da}\bigg{)}\bigg{(}\frac{dl(\theta^{h})}{db}\bigg{)}^{\prime}\bigg{|}\theta=\theta_{s}\bigg{]}.$

For further simplification of notation, when $a$ and $b$ are equal, the second
subscript will be dropped and the matrix
$\mathcal{I}^{h}_{a,b}=\mathcal{I}^{h}_{a,a}$ will be denoted simply as
$\mathcal{I}^{h}_{a}$. Analogous notation is used for the full and Choleksy
representations.

####  2.1.1 Fisher Scoring

The first variant of Fisher Scoring we provide uses the “half”-representation
for $(\beta,\sigma^{2},D)$, $\theta^{h}$, and is based on the standard form of
Fisher Scoring given by (5). This may be considered the most natural approach
for a Fisher Scoring algorithm as $\theta^{h}$ is an unmodified vector of the
unique parameters of the LMM and (5) is the standard update rule. For this
approach the elements of the score vector are:

$\frac{dl(\theta^{h})}{d\beta}=\sigma^{-2}X^{\prime}V^{-1}e,$ (7)
$\frac{dl(\theta^{h})}{d\sigma^{2}}=-\frac{n}{2}\sigma^{-2}+\frac{1}{2}\sigma^{-4}e^{\prime}V^{-1}e.$
(8)

For $k\in\\{1,\ldots,r\\}$:

$\frac{dl(\theta^{h})}{d\text{vech}(D_{k})}=\frac{1}{2}\mathcal{D}_{q_{k}}^{\prime}\text{vec}\bigg{(}\sum_{j=1}^{l_{k}}Z^{\prime}_{(k,j)}V^{-1}\bigg{(}\frac{ee^{\prime}}{\sigma^{2}}-V\bigg{)}V^{-1}Z_{(k,j)}\bigg{)}.$
(9)

and the entries of the Fisher Information matrix are given by:

$\mathcal{I}^{h}_{\beta}=\sigma^{-2}X^{\prime}V^{-1}X,\hskip
11.38092pt\mathcal{I}^{h}_{\beta,\sigma^{2}}=\mathbf{0}_{p,1},\hskip
11.38092pt\mathcal{I}^{h}_{\sigma^{2}}=\frac{n}{2}\sigma^{-4}.$ (10)

For $k\in\\{1,\ldots,r\\}$:

$\displaystyle\mathcal{I}^{h}_{\beta,\text{vech}(D_{k})}=\mathbf{0}_{p,q_{k}(q_{k}+1)/2},$
(11)
$\displaystyle\mathcal{I}^{h}_{\sigma^{2},\text{vech}(D_{k})}=\frac{1}{2}\sigma^{-2}\text{vec}^{\prime}\bigg{(}\small\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}V^{-1}Z_{(k,j)}\bigg{)}\mathcal{D}_{q_{k}}.$

For $k_{1},k_{2}\in\\{1,\ldots,r\\}$:

$\displaystyle\mathcal{I}^{h}_{\text{vech}(D_{k_{1}}),\text{vech}(D_{k_{2}})}=$
(12)
$\displaystyle\frac{1}{2}\mathcal{D}^{\prime}_{q_{k_{1}}}\sum_{j=1}^{l_{k_{2}}}\sum_{i=1}^{l_{k_{1}}}(Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)}\otimes
Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)})\mathcal{D}_{q_{k_{2}}}.$

where $\mathbf{0}_{n,k}$ denotes the $(n\times k)$ dimensional matrix of
zeros. Due to its to natural approach to Fisher Scoring, this algorithm is
referred to as FS in the remainder of this text. Pseudocode for the FS
algorithm is given in Algorithm 1. Discussion of the initial estimates used in
the algorithm is deferred to Section 2.2.

1 Assign $\theta^{h}$ to an initial estimate using (25) and (26) while
_current $l(\theta^{h})$ and previous $l(\theta^{h})$ differ by more than a
predefined tolerance_ do

2

3 Calculate $\frac{dl(\theta^{h})}{d\theta^{h}}$ using (7)-(9). Calculate
$\mathcal{I}(\theta^{h})$ using (10)-(12) Assign
$\theta^{h}=\theta^{h}+\alpha\mathcal{I}(\theta^{h})^{-1}\frac{dl(\theta^{h})}{d\theta^{h}}$
Assign $\alpha=\frac{\alpha}{2}$ if $l(\theta^{h})$ has decreased in value.

4 end while

Algorithm 1 Fisher Scoring (FS)

####  2.1.2 Full Fisher Scoring

The second variant of Fisher Scoring considered in this work uses the “full”,
$\theta^{f}$, representation of the model parameters, and shall therefore be
referred to as “Full Fisher Scoring” (FFS). In this approach, for each factor,
$f_{k}$, the elements of vec$(D_{k})$ are to be treated as distinct from one
another with numerical optimization for $D_{k}$ performed over the space of
all $(q_{k}\times q_{k})$ matrices. This approach differs from the FS method
proposed in the previous section, in which optimization was performed on the
space of only those $(q_{k}\times q_{k})$ matrices that are symmetric. This
optimization procedure is realized by treating symmetric matrix elements of
$D_{k}$ as distinct and, for a given element, using the partial derivative
with respect to the element during optimization instead of the total
derivative with respect to both the element and its symmetric counterpart.
This change is reflected by denoting the elements of the score vector which
correspond to vec$(D_{k})$ using a partial derivative operator, $\partial$,
rather than the total derivative operator, $d$. The primary motivation for the
inclusion of the FFS approach is that it serves as a basis for which the
constrained covariance approaches of Section 2.5 can be built upon. However,
it should be noted that, as it does not require the construction or use of
duplication matrices, FFS also provides simplified expressions and potential
improvement in terms of computation speed. As a result, FFS is of some
theoretical and practical interest and is detailed in full here.

An immediate obstacle to this approach is that the Fisher Information matrix
of $\theta^{f}$ is rank-deficient, and, therefore, cannot be inverted.
Intuitively, this is to be expected, as repeated entries in $\theta^{f}$
result in repeated rows in $\mathcal{I}(\theta^{f})$. Mathematically, this can
be seen by noting that $\mathcal{I}^{f}_{\text{vec}(D_{k})}$ can be expressed
as a product containing the matrix $N_{q_{k}}$ (defined in Section 1.2.2),
which is low-rank by construction (see Appendix 6.3). Consequently, the Fisher
Scoring update rule for $\theta^{f}$ must be based on the Pseudo-inverse
formulation of Fisher Scoring given in (6).

As the derivatives of the log-likelihood with respect to $\beta$ and
$\sigma^{2}$ do not depend upon the parameterisation of $D$, both FFS and FS
employ the same expressions for the elements of the score vector which
correspond to $\beta$ and $\sigma^{2}$, given by (7) and (8) respectively. The
score vector for $\\{$vec($D_{k})\\}_{k\in\\{1,\dots r\\}}$ used by FFS is
given by:

$\frac{\partial
l(\theta^{f})}{\partial\text{vec}(D_{k})}=\frac{1}{2}\text{vec}\bigg{(}\sum_{j=1}^{l_{k}}Z^{\prime}_{(k,j)}V^{-1}\bigg{(}\frac{ee^{\prime}}{\sigma^{2}}-V\bigg{)}V^{-1}Z_{(k,j)}\bigg{)}.$
(13)

The entries of the Fisher Information matrix of $\theta^{f}$, based on the
derivatives given in (7), (8) and (13), are given by:

$\mathcal{I}^{f}_{\beta}=\mathcal{I}^{h}_{\beta},\hskip
19.91684pt\mathcal{I}^{f}_{\beta,\sigma^{2}}=\mathcal{I}^{h}_{\beta,\sigma^{2}},\hskip
19.91684pt\mathcal{I}^{f}_{\sigma^{2}}=\mathcal{I}^{h}_{\sigma^{2}}.$

For $k\in\\{1,\ldots,r\\}$:

$\displaystyle\mathcal{I}^{f}_{\beta,\text{vec}(D_{k})}=\mathbf{0}_{p,q_{k}^{2}},$
(14)
$\displaystyle\mathcal{I}^{f}_{\sigma^{2},\text{vec}(D_{k})}=\frac{1}{2}\sigma^{-2}\text{vec}^{\prime}\bigg{(}\small\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}V^{-1}Z_{(k,j)}\bigg{)}.$

For $k_{1},k_{2}\in\\{1,\ldots,r\\}$:

$\displaystyle\mathcal{I}^{f}_{\text{vec}(D_{k_{1}}),\text{vec}(D_{k_{2}})}=$
(15)
$\displaystyle\frac{1}{2}\sum_{j=1}^{l_{k_{2}}}\sum_{i=1}^{l_{k_{1}}}(Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)}\otimes
Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)})N_{q_{k}}.$

Derivations for the above can be found in Appendices 6.1 and 6.2. In addition,
in Appendix 6.3, we show that the Full Fisher Scoring algorithm can also be
expressed in the form:

$\theta^{f}_{s+1}=\theta^{f}_{s}+\alpha_{s}F(\theta_{s}^{f})^{-1}\frac{\partial
l(\theta_{s}^{f})}{\partial\theta}_{,}$ (16)

where, unlike in (5) and (6), $F(\theta^{f})$ is not the Fisher Information
matrix. Rather, $F(\theta^{f})$ is a matrix of equal dimensions to
$\mathcal{I}(\theta^{f})$ with all of its elements equal to those of
$\mathcal{I}(\theta^{f})$, apart from those which were specified by (15),
which instead are, for $k_{1},k_{2}\in\\{1,\ldots,r\\}$:

$\begin{aligned} &F_{\text{vec}(D_{k_{1}}),\text{vec}(D_{k_{2}})}=\\\
&\frac{1}{2}\sum_{j=1}^{l_{k_{2}}}\sum_{i=1}^{l_{k_{1}}}(Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)}\otimes
Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)})\end{aligned},$ (17)

where the same subscript notation has been adopted to index $F(\theta^{f})$ as
was adopted for $\mathcal{I}(\theta^{f})$.

The above approach, which utilizes the matrix $F(\theta^{f})$, rather than
$\mathcal{I}(\theta^{f})$, is adopted to match the formulation of the Fisher
Scoring algorithm given in Section 2.11 of Demidenko (2013) for the single-
factor LMM. Pseudocode for the FFS algorithm using the representation of the
update rule given by (16), is provided by Algorithm 2.

1 Assign $\theta^{f}$ to an initial estimate using (25) and (26) while
_current $l(\theta^{f})$ and previous $l(\theta^{f})$ differ by more than a
predefined tolerance_ do

2

3 Calculate $\frac{\partial l(\theta^{f})}{\partial\theta^{f}}$ using (7),(8)
and (13). Calculate $F(\theta^{f})$ using (14) and (17). Assign
$\theta^{f}=\theta^{f}+\alpha F(\theta^{f})^{-1}\frac{\partial
l(\theta^{f})}{\partial\theta^{f}}$ Assign $\alpha=\frac{\alpha}{2}$ if
$l(\theta^{f})$ has decreased in value.

4 end while

Algorithm 2 Full Fisher Scoring (FFS)

####  2.1.3 Simplified Fisher Scoring

The third Fisher Scoring algorithm proposed in this work relies on the half-
representation of the parameters $(\beta,\sigma^{2},D)$ and takes an approach,
similar to that of coordinate ascent, which is commonly adopted in the single-
factor setting (c.f. Demidenko (2013)). In this approach, instead of
performing a single update step upon the entire vector $\theta^{h}$ in the
form of $(\ref{FS})$, updates for $\beta,\sigma^{2}$ and
$\\{D_{k}\\}_{k\in\\{1,\ldots,r\\}}$ are performed individually in turn. For
$\beta$ and $\sigma^{2}$ each iteration uses the standard Generalized Least
Squares (GLS) estimators given by:

$\beta_{s+1}=(X^{\prime}V_{s}^{-1}X)^{-1}X^{\prime}V_{s}^{-1}Y,$ (18)
$\sigma^{2}_{s+1}=\frac{e_{s+1}^{\prime}V^{-1}_{s}e_{s+1}}{n}.$ (19)

To update the random effects covariance matrix, $D_{k}$, for each factor,
$f_{k}$, individual Fisher Scoring updates are applied to vech$(D_{k})$. These
updates are performed using the block of the Fisher Information matrix
corresponding to vech$(D_{k})$, given by (12), and take the following form for
$k\in\\{1,\ldots,r\\}$:

$\text{vech}(D_{k,s+1})=\text{vech}(D_{k,s})+\alpha_{s}\big{(}\mathcal{I}^{h}_{\text{vech}(D_{k,s})}\big{)}^{-1}\frac{dl(\theta^{h}_{s})}{d\text{vech}(D_{k,s})}.$
(20)

In line with the naming convention used in Demidenko (2013), this method shall
be referred to as “Simplified” Fisher Scoring (SFS). This is due to the
relative simplicity, both in terms of notational and computational complexity,
of the updates (18)-(20) used in the SFS algorithm in comparison to those used
in the FS algorithm of Section 2.1.1, given by (7)-(12). Pseudocode for the
SFS algorithm is given by Algorithm 3.

1 Assign $\theta^{f}$ to an initial estimate using (25) and (26) while
_current $l(\theta^{h})$ and previous $l(\theta^{h})$ differ by more than a
predefined tolerance_ do

2

3 Update $\beta$ using (18) Update $\sigma^{2}$ using (19) for
_$k\in\\{1,...r\\}$_ do

4 Update vech($D_{k}$) using (20)

5 end for

6 Assign $\alpha=\frac{\alpha}{2}$ if $l(\theta^{h})$ has decreased in value.

7 end while

Algorithm 3 Simplified Fisher Scoring (SFS)

####  2.1.4 Full Simplified Fisher Scoring

The Full Simplified Fisher Scoring algorithm (FSFS) combines the “Full” and
“Simplified” approaches described in Sections 2.1.2 and 2.1.3. In the FSFS
algorithm individual updates are applied to $\beta$ and $\sigma^{2}$ using
(18) and (19) respectively and to
$\\{\text{vec}(D_{k})\\}_{k\in\\{1,\ldots,r\\}}$ using a Fisher Scoring update
step, based on the matrix $F_{\text{vec}(D_{k})}$ given by (17). The update
rule for $\\{\text{vec}(D_{k})\\}_{k\in\\{1,\ldots,r\\}}$ takes the following
form:

$\text{vec}(D_{k,s+1})=\text{vec}(D_{k,s})+\alpha_{s}F^{-1}_{\text{vec}(D_{k,s})}\frac{\partial
l(\theta^{f}_{s})}{\partial\text{vec}(D_{k,s})}.$ (21)

In Appendix 6.3, we justify the above update rule by demonstrating that it is
equivalent to using an update rule of the form (6), given by:

$\text{vec}(D_{k,s+1})=\text{vec}(D_{k,s})+\alpha_{s}\big{(}\mathcal{I}^{f}_{\text{vec}(D_{k,s})}\big{)}^{+}\frac{\partial
l(\theta^{f}_{s})}{\partial\text{vec}(D_{k,s})}.$

Pseudocode for the FSFS algorithm is given by Algorithm 4.

1 Assign $\theta^{f}$ to an initial estimate using (25) and (26) while
_current $l(\theta^{f})$ and previous $l(\theta^{f})$ differ by more than a
predefined tolerance_ do

2

3 Update $\beta$ using (18) Update $\sigma^{2}$ using (19) for
_$k\in\\{1,...r\\}$_ do

4 Update vec($D_{k}$) using (21)

5 end for

6 Assign $\alpha=\frac{\alpha}{2}$ if $l(\theta^{f})$ has decreased in value.

7 end while

Algorithm 4 Full Simplified Fisher Scoring (FSFS)

####  2.1.5 Cholesky Simplified Fisher Scoring

The final variant of the Fisher Scoring algorithm we consider is based on the
“simplified” approach described in Section 2.1.3 and uses the Cholesky
parameterisation of $(\beta,\sigma^{2},D)$, $\theta^{c}$. This approach
follows directly from the below application of the chain rule of
differentiation for vector-valued functions,

$\frac{dl(\theta^{c})}{d\text{vech}(\Lambda_{k})}=\frac{\partial\text{vech}(D_{k})}{\partial\text{vech}(\Lambda_{k})}\frac{\partial
l(\theta^{c})}{\partial\text{vech}(D_{k})}.$

An expression for the derivative which appears second in the above product was
given by (9). It therefore follows that in order to evaluate the score vector
of $\text{vech}(\Lambda_{k})$ (i.e. the derivative of $l$ with respect to
$\text{vech}(\Lambda_{k})$), only an expression for the first term of the
above product is required. This term can be evaluated to
$\mathcal{L}_{q_{k}}(\Lambda_{k}^{\prime}\otimes
I_{q_{k}})(I_{q_{k}^{2}}+K_{q_{k}})\mathcal{D}_{q_{k}}$, proof of which is
provided in Appendix 6.6.

Through similar arguments to those used to prove Corollaries 4-6 of Appendix
6.2, it can be shown that the Fisher Information matrix for $\theta^{c}$ is
given by:

$\mathcal{I}^{c}_{\beta}=\mathcal{I}^{h}_{\beta},\hskip
11.38092pt\mathcal{I}^{c}_{\beta,\sigma^{2}}=\mathcal{I}^{h}_{\beta,\sigma^{2}},\hskip
11.38092pt\mathcal{I}^{c}_{\sigma^{2}}=\mathcal{I}^{h}_{\sigma^{2}}.$

For $k\in\\{1,\ldots,r\\}$:

$\displaystyle\mathcal{I}^{c}_{\beta,\text{vech}(\Lambda_{k})}=\mathbf{0}_{p,q_{k}(q_{k}+1)/2},$
(22)
$\displaystyle\mathcal{I}^{c}_{\sigma^{2},\text{vech}(\Lambda_{k})}=\mathcal{I}^{h}_{\sigma^{2},\text{vech}(D_{k})}\bigg{(}\frac{\partial\text{vech}(D_{k})}{\partial\text{vech}(\Lambda_{k})}\bigg{)}^{\prime}.$

For $k_{1},k_{2}\in\\{1,\ldots,r\\}$:

$\displaystyle\mathcal{I}^{c}_{\text{vech}(\Lambda_{k_{1}}),\text{vech}(\Lambda_{k_{2}})}=$
(23)
$\displaystyle\bigg{(}\frac{\partial\text{vech}(D_{k_{1}})}{\partial\text{vech}(\Lambda_{k_{1}})}\bigg{)}\mathcal{I}^{h}_{\text{vech}(D_{k_{1}}),\text{vech}(D_{k_{2}})}\bigg{(}\frac{\partial\text{vech}(D_{k_{2}})}{\partial\text{vech}(\Lambda_{k_{2}})}\bigg{)}^{\prime}.$

From the above, it can be seen that a non-“simplified” Cholesky-based variant
of the Fisher Scoring algorithm, akin to the FS and FFS algorithms described
in Sections 2.1.1 and 2.1.2, may be constructed. This would involve
constructing the Fisher Information matrix of $\theta^{c}$,
$\mathcal{I}(\theta^{c})$, using (22) and (23), and employing a Fisher Scoring
procedure similar to that specified by Algorithms 1 and 2. Whilst this
approach is feasible, in terms of implementation, preliminary tests have
indicated that the performance of this approach, in terms of computation time,
is significantly worse than the previously proposed algorithms. For this
reason, we only consider the “simplified” version of the Cholesky Fisher
Scoring algorithm, analogous to the “simplified” approaches described in
Sections 2.1.3 and 2.1.4, is considered here. The Cholesky Simplified Fisher
Scoring (CSFS) algorithm adopts (18) and (19) as update rules for $\beta$ and
$\sigma^{2}$, respectively, and the following update rule for
$\\{\text{vech}(\Lambda_{k})\\}_{k\in\\{1,...,r\\}}$.

$\small{\text{vech}(\Lambda_{k,s+1})=\text{vech}(\Lambda_{k,s})+\alpha_{s}\big{(}\mathcal{I}^{c}_{\text{vech}(\Lambda_{k,s})}\big{)}^{-1}\frac{dl(\theta^{c}_{s})}{d\text{vech}(\Lambda_{k,s})}}.$
(24)

Pseudocode summarizing the CSFS approach is given by Algorithm 5.

1 Assign $\theta^{c}$ to an initial estimate using (25) and (26) while
_current $l(\theta^{c})$ and previous $l(\theta^{c})$ differ by more than a
predefined tolerance_ do

2

3 Update $\beta$ using (18) Update $\sigma^{2}$ using (19) for
_$k\in\\{1,...r\\}$_ do

4 Update vec($\Lambda_{k}$) using (24)

5 end for

6 Assign $\alpha=\frac{\alpha}{2}$ if $l(\theta^{c})$ has decreased in value.

7 end while

Algorithm 5 Cholesky Simplified Fisher Scoring (CSFS)

###  2.2 Initial values

Choosing which initial values of $\beta$, $\sigma^{2}$ and $D$ will be used as
starting points for optimization is an important consideration for the Fisher
Scoring algorithm. Denoting these initial values as $\beta_{0}$,
$\sigma^{2}_{0}$ and $D_{0}$, respectively, this work follows the
recommendations of Demidenko (2013) and evaluates $\beta_{0}$ and
$\sigma^{2}_{0}$ using the OLS estimates given by,

$\beta_{0}=(X^{\prime}X)^{-1}X^{\prime}Y,\hskip
17.07182pt\sigma^{2}_{0}=\frac{e_{0}^{\prime}e_{0}}{n}.$ (25)

where $e_{0}$ is defined as $e_{0}=Y-X\beta_{0}$. For the initial estimate of
$\\{D_{k}\\}_{k\in\\{1,\dots,r\\}}$, an approach similar to that suggested in
Demidenko (2013) is also adopted, which substitutes $V$ for $I_{n}$ in the
update rule for vec($D_{k}$), equation (21). The resulting initial estimate
for $\\{D_{k}\\}_{k\in\\{1,\dots,r\\}}$ is given by

$\begin{aligned}
&\text{vec}(D_{k,0})=&\bigg{(}\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}Z_{(k,j)}\otimes
Z_{(k,j)}^{\prime}Z_{(k,j)}\bigg{)}^{-1}\times&\\\
&&\text{vec}\bigg{(}\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}\bigg{(}\frac{e_{0}e_{0}^{\prime}}{\sigma_{0}^{2}}-I_{n}\bigg{)}Z_{(k,j)}\bigg{)}&\end{aligned}.$
(26)

###  2.3 Ensuring non-negative definiteness

In order to ensure numerical optimization produces a valid covariance matrix
when performing LMM parameter estimation, optimization must be constrained to
ensure that for each factor, $f_{k}$, $D_{k}$ is non-negative definite. This
is a well-documented obstacle for numerical approaches to LMM parameter
estimation for which many solutions have been suggested. A common approach,
utilized in Bates et al. (2015) for example, is to reparameterize optimization
in terms of the Cholesky factor of the covariance matrix, $\Lambda_{k}$, as
oppose to the covariance matrix itself, $D_{k}$. This approach is adopted by
the CSFS method, described in Section 2.1.5. However, the methods described in
Sections 2.1.1-2.1.4 require an alternative approach as optimization does not
account for the constraint that $D_{k}$ must be non-negative definite. For the
FS, FFS, SFS and FSFS algorithms, we use an approach described by Demidenko
(2013), for the single factor LMM. Denoting the eigendecomposition of
$\\{D_{k}\\}_{k\in\\{1,\ldots,r\\}}$ as $D_{k}=U_{k}\Sigma_{k}U_{k}^{\prime}$,
we project $D_{k}$ onto the space of all $(q_{k}\times q_{k})$ non-negative
definite matrices by, following each update of the covariance matrix, $D$,
replacing $\\{D_{k}\\}_{k\in\\{1,..r\\}}$ with
$\\{D_{+,k}\\}_{k\in\\{1,..r\\}}$,

$D_{+,k}=U_{k}\Sigma_{+,k}U_{k}^{\prime},$

where $\Sigma_{+,k}=\max(\Sigma_{k},0)$ with ‘max’ representing the element-
wise maximum. For each factor, $f_{k}$, this ensures $D_{k}$, and as a result
$D$, is non-negative definite.

###  2.4 Computational efficiency

This section provides discussion on the computational efficiency of evaluating
the Fisher Information matrices and score vectors of Sections 2.1.1-2.1.5.
This discussion is, in large part, motivated by the mass-univariate setting
(c.f. Section 1.1), in which not one, but rather hundreds of thousands of
models must be estimated concurrently. As a result, this discussion
prioritizes both time efficiency and memory consumption concerns and, further,
operates under the assumption that sparse matrix methodology, such as that
employed by the R package ‘lmer’, cannot be employed. This assumption has been
made primarily as a result of the current lack of support for computationally
vectorized sparse matrix operations.

A primary concern, for both memory consumption and time efficiency, stems from
the fact that many of the matrices used in the evaluation of the score vectors
and Fisher Information matrices possess dimensions which scale with $n$, the
number of observations. For example, $V$ has dimensions $(n\times n)$ and is
frequently inverted in the FSFS algorithm. In practice, it is not uncommon for
studies to have $n$ ranging into the thousands. As such, inverting $V$
directly may not be computationally feasible. To address this issue, we define
the “product forms” as;

$P=X^{\prime}X,\hskip 3.98337ptQ=X^{\prime}Y,\hskip
3.98337ptR=X^{\prime}Z,\hskip 3.98337ptS=Y^{\prime}Y,\hskip
3.98337ptT=Y^{\prime}Z,\hskip 3.98337ptU=Z^{\prime}Z.$

Working with the product forms is preferable to using the original matrices
$X,Y$ and $Z$, as the dimensions of the product forms do not scale with $n$.
Instead, the dimensions of the product forms depend only on the number of
fixed effects, $p$, and the second dimension of the random effects design
matrix, $q$. As an example, consider the expressions below, which appear
frequently in (13) and (17) and have been reformulated in terms of the product
forms:

$\displaystyle
Z_{(k_{1},i)}^{\prime}V^{-1}Z_{(k_{2},j)}=(U-UD(I_{q}+DU)^{-1}U)_{[(k_{1},i),(k_{2},j)]},$
$\displaystyle
Z_{(k,j)}^{\prime}V^{-1}e=((I_{q}-UD)(I_{q}+DU)^{-1}(T^{\prime}-R^{\prime}\beta))_{[(k,j),:]}.$

For computational purposes, the right-hand side of the above expressions are
much more convenient than the left-hand side. In order to evaluate the left-
hand side in both cases, an $(n\times n)$ inversion of the matrix $V$ must be
performed. In contrast, the right-hand side is expressible purely in terms of
the product forms and $(\beta,\sigma^{2},D)$, with the only inversion required
being that of the $(q\times q)$ matrix $(I_{q}+DU)$. These examples can be
generalized further. In fact, all of the previous expressions (7)-(26) can be
rewritten in terms of only the product forms and $(\beta,\sigma^{2},D)$. This
observation is important as it implies that an algorithm for mixed model
parameter estimation may begin by taking the matrices $X,Y$ and $Z$ as inputs,
but discard them entirely once the product forms have been constructed. As a
result, both computation time and memory consumption no longer scale with $n$.

Another potential source of concern regarding computation speed arises from
noting that the algorithms we have presented contain many summations of the
following two forms:

$\sum_{i=1}^{c_{0}}A_{i}B_{i}^{\prime}\hskip 5.69046pt\text{ and }\hskip
5.69046pt\sum_{i=1}^{c_{1}}\sum_{i=1}^{c_{2}}G_{i,j}\otimes H_{i,j}.$ (27)

where matrices $\\{A_{i}\\}$ and $\\{B_{i}\\}$ are of dimension $(m_{1}\times
m_{2})$, and matrices $\\{G_{i,j}\\}$ and $\\{H_{i,j}\\}$ are of dimension
$(n_{1}\times n_{2})$. We denote the matrices formed from vertical
concatenation of the $\\{A_{i}\\}$ and $\\{B_{i}\\}$ matrices as $A$ and $B$,
respectively, and $G$ and $H$ the matrices formed from block-wise
concatenation of the $\\{G_{i,j}\\}$ and $\\{H_{i,j}\\}$, respectively.
Instances of such summations can be found, for example, in equations (9) and
(12).

From a computational standpoint, summations of the forms shown in (27) are
typically realized by ‘for’ loops. This cumulative approach to computation can
cause a potential issue for LMM computation since typically the number of
summands corresponds to the number of levels, $l_{k}$, of some factor,
$f_{k}$. In particular applications of the LMM, such as repeated measures and
longitudinal studies, some factors may possess large quantities of levels. As
this means ‘for’ loops of this kind could hinder computation and result in
slow performance, we provide alternative methods for calculating summations of
the forms shown in (27).

For the former summation shown in (27), we utilize the “generalized
vectorisation”, or “vecm” operator, defined by Turkington (2013) as the
operator which performs the below mapping for a horizontally partitioned
matrix $M$:

$M=\begin{bmatrix}M_{1}&M_{2}&\ldots
M_{c_{0}}\end{bmatrix}\rightarrow\text{vec}_{m}(M)=\begin{bmatrix}M_{1}\\\
M_{2}\\\ \vdots\\\ M_{c_{0}}\end{bmatrix},$

where the partitions $\\{M_{i}\\}$ are evenly sized and contain $m$ columns.
Using the definition of the generalized vectorisation operator, the former
summation in (27) can be reformulated as:

$\sum_{i=1}^{l}A_{i}B_{i}^{\prime}=\text{vec}_{m_{2}}(A^{\prime})^{\prime}\text{vec}_{m_{2}}(B^{\prime}).$

The right-hand side of the expression above is of practical utility as the
‘vecm’ operator can be implemented efficiently. The ‘vecm’ operation can be
performed, for instance, using matrix resizing operations such as the
‘reshape’ operators commonly available in many programming languages such as
MATLAB and Python. The computational costs associated with this approach are
significantly lesser than those experienced when evaluating the summation
directly using ‘for’ loops.

For the latter summation in (27), we first define the notation $\tilde{M}$ to
denote the transformation below, for the block-wise partitioned matrix $M$:

$M=\begin{bmatrix}M_{1,1}&M_{1,2}&...&M_{1,c_{2}}\\\
M_{2,1}&M_{2,2}&...&M_{2,c_{2}}\\\ \vdots&\vdots&\ddots&\vdots\\\
M_{c_{1},1}&M_{c_{1},2}&...&M_{c_{1},c_{2}}\end{bmatrix}\rightarrow\tilde{M}=\begin{bmatrix}\text{vec}(M_{1,1})^{\prime}\\\
\text{vec}(M_{1,2})^{\prime}\\\ \vdots\\\ \text{vec}(M_{1,c_{2}})^{\prime}\\\
\text{vec}(M_{2,1})^{\prime}\\\ \vdots\\\
\text{vec}(M_{c_{1},c_{2}})^{\prime}\end{bmatrix}.$ (28)

Using a modification of Lemma 3.1 (i) of Neudecker and Wansbeek (1983), we
obtain the below identity:

$\text{vec}\bigg{(}\sum_{i,j}G_{i,j}\otimes H_{i,j}\bigg{)}=(I_{n_{2}}\otimes
K_{n_{1},n_{2}}\otimes I_{n_{1}})\text{vec}(\tilde{H}^{\prime}\tilde{G}),$
(29)

where $K_{n_{1},n_{2}}$ is the $(n_{1}n_{2}\times n_{1}n_{2})$ Commutation
matrix (c.f. Section 1.2.2). The matrices $\tilde{H}$ and $\tilde{G}$ can be
obtained from $H$ and $G$ using resizing operators with little computation
time. The matrix $(I_{n_{2}}\otimes K_{n_{1},n_{2}}\otimes I_{n_{1}})$ can be
calculated via simple means and is a permutation matrix depending only on the
dimensions of $H$ and $G$. As a result, this matrix can be calculated once and
stored as a vector. Therefore, the matrix multiplication in the above
expression does not need to be performed directly, but instead can be
evaluated by permuting the elements of vec$(\tilde{H}^{\prime}\tilde{G})$
according to the permutation vector representing $(I_{n_{2}}\otimes
K_{n_{1},n_{2}}\otimes I_{n_{1}})$. To summarize, evaluation of expressions of
the latter form shown in (27) can be performed by using only reshape
operations, a matrix multiplication, and a permutation. This method of
evaluation provides notably faster computation time than the corresponding
‘for’ loop evaluated over all values of $i$ and $j$.

The above method, combined with the product form approach, was used to obtain
the results of Sections 2.7-3.2.

###  2.5 Constrained covariance structure

In many applications involving the LMM, it is often desirable to use a
constrained parameterisation for $D_{k}$. Examples include compound symmetry,
first-order auto-regression and a Toeplitz structure. A more comprehensive
list of commonly employed covariance structures in LMM analyses can be found,
for example, in Wolfinger (1996). In this section, we describe how the Fisher
Scoring algorithms of the previous sections can be adjusted to model
dependence between covariance elements.

When a constraint is placed on the covariance matrix $D_{k}$, it is assumed
that the elements of $D_{k}$ can be defined as continuous, differentiable
functions of some smaller parameter vector, vecu$(D_{k})$. Colloquially,
vecu$(D_{k})$ may be thought of as the vector of “unique” parameters required
to specify the constrained parameterization of $D_{k}$. To perform constrained
optimization, we adopt a similar approach to the previous sections and define
the constrained representation of $\theta$ as
$\theta^{con}=[\beta^{\prime},\sigma^{2},\text{vecu}(D_{1})^{\prime},...\text{vecu}(D_{k})^{\prime}]^{\prime}$.
Denoting the Jacobian matrix $\partial$vec$(D_{k})/\partial$vecu$(D_{k})$ as
$\mathcal{C}_{k}$, the score vector and Fisher Information matrix of
$\theta^{con}$ can be constructed as follows; for $k\in\\{1,\ldots,r\\}$,

$\displaystyle\frac{dl(\theta^{con})}{d\text{vecu}(D_{k})}=\mathcal{C}_{k}\frac{\partial
l(\theta^{con})}{\partial\text{vec}(D_{k})},$ (30)
$\displaystyle\mathcal{I}^{con}_{\beta,\text{vecu}(\tilde{D}_{k})}=\mathbf{0}_{p,\tilde{q}_{k}},\hskip
14.22636pt\mathcal{I}^{con}_{\sigma^{2},\text{vecu}(\tilde{D}_{k})}=\mathcal{I}^{f}_{\sigma^{2},\text{vec}(D_{k})}\mathcal{C}_{k}^{\prime},$

and, for $k_{1},k_{2}\in\\{1,\ldots,r\\}$,

$\mathcal{I}^{con}_{\text{vecu}(\tilde{D}_{k_{1}}),\text{vecu}(\tilde{D}_{k_{2}})}=\mathcal{C}_{k_{1}}\mathcal{I}^{f}_{\text{vec}(D_{k_{1}}),\text{vec}(D_{k_{2}})}\mathcal{C}_{k_{2}}^{\prime},$

where $\tilde{q}_{k}$ is the length of vecu$(D_{k})$ and $\mathcal{I}^{f}$ is
the “full” Fisher Information matrix defined in Section 2.1.2. The above
expressions can be derived trivially using the definition of the Fisher
Information matrix and the chain rule. In the remainder of this work, the
matrix $\mathcal{C}_{k}$ is referred to as a “constraint matrix” due to the
fact it “imposes” constraints during optimization. A fuller discussion of
constraint matrices, alongside examples, is provided in Appendix 6.5.

We note here that this approach can be extended to situations in which all
covariance parameters can be expressed in terms of one set of parameters,
$\rho_{D}$, common to all $\\{\text{vec}(D_{k})\\}_{k\in\\{1,\ldots,r\\}}$. In
such situations, a constraint matrix, $\mathcal{C}$, may be defined as the
Jacobian matrix
$\partial[\text{vec}(D_{1})^{\prime},...\text{vec}(D_{1})^{\prime}]^{\prime}/\partial\rho_{D}$
and Fisher Information matrices and score vectors may be derived in a similar
manner to the above. Models requiring this type of commonality between factors
are rare in the LMM literature, since the covariance matrices
$\\{D_{k}\\}_{k\in\\{1,...r\\}}$ are typically defined independently of one
another. However, one such example is provided by the ACE model employed for
twin studies, in which the elements of $v(D)$ can all be expressed in terms of
three variance components;
$\rho_{D}=[\sigma^{2}_{a},\sigma^{2}_{c},\sigma^{2}_{e}]^{\prime}$. Further
information on the ACE model is presented in Section 2.8.2.

###  2.6 Degrees of freedom estimation

Several methods exist for drawing inference on the fixed effects parameter
vector, $\beta$. Often, research questions can be expressed as hypotheses of
the below form:

$H_{0}:L\beta=0,\hskip 14.22636ptH_{1}:L\beta\neq 0,$

where $L$ is a fixed and known $(1\times p)$-sized contrast matrix specifying
a hypothesis, or prior belief, about linear relationships between the elements
of $\beta$, upon which an inference is to be made. In the setting of the LMM,
a commonly employed statistic for testing hypotheses of this form is the
approximate T-statistic, given by:

$T=\frac{L\hat{\beta}}{\sqrt{\hat{\sigma}^{2}L(X^{\prime}\hat{V}^{-1}X)^{-1}L^{\prime}}},$

where $\hat{V}=I_{n}+Z\hat{D}Z^{\prime}$. As noted in Dempster et al. (1981),
using an estimate of $D$ in this manner results in an underestimation of the
true variability in $\hat{\beta}$. For this reason, the relation $T\sim t_{v}$
is not exact and is treated only as an approximating distribution (i.e. an
“approximate T statistic”) where the degrees of freedom, $v$, must be
estimated empirically. A standard method for estimating the degrees of
freedom, $v$, utilizes the Welch-Satterthwaite equation, originally described
in Satterthwaite (1946) and Welch (1947), given by:

$v(\hat{\eta})=\frac{2(S^{2}(\hat{\eta}))^{2}}{\text{Var}(S^{2}(\hat{\eta}))},$
(31)

where $\hat{\eta}$ represents an estimate of the variance parameters
$\eta=(\sigma^{2},D_{1},\ldots D_{r})$ and $S^{2}(\hat{\eta})$ is given by:

$S^{2}(\hat{\eta})=\hat{\sigma}^{2}L(X^{\prime}\hat{V}^{-1}X)^{-1}L^{\prime}.$

Typically, as the variance estimates obtained under ML estimation are biased
downwards, ReML estimation is employed to obtain the estimate of $\eta$
employed in the above expression. If ML were used to estimate $\eta$ instead,
the degrees of freedom, $v(\hat{\eta})$, would be underestimated and,
consequently, conservative p-values and a reduction in statistical power for
resulting hypothesis tests would be observed.

To obtain an approximation for $v(\hat{\eta})$, a second order Taylor
expansion is applied to the unknown variance on the denominator of (31):

$\text{Var}(S^{2}(\hat{\eta}))\approx\bigg{(}\frac{dS^{2}(\hat{\eta})}{d\hat{\eta}}\bigg{)}^{\prime}\text{Var}(\hat{\eta})\bigg{(}\frac{dS^{2}(\hat{\eta})}{d\hat{\eta}}\bigg{)}.$
(32)

This approach is well-documented, and has been adopted, most notably, by the R
package lmerTest, first presented in Kuznetsova et al. (2017). To obtain the
derivative of $S^{2}$ and asymptotic variance covariance matrix, lmerTest
utilizes numerical estimation methods. As an alternative, we define the “half”
representation of $\hat{\eta}$, $\hat{\eta}^{h}$, by
$\hat{\eta}^{h}=[\hat{\sigma}^{2},\text{vech}(\hat{D}_{1})^{\prime},\ldots,\text{vech}(\hat{D}_{r})^{\prime}]^{\prime}$
and present the below exact closed-form expression for the derivative of
$S^{2}$ in terms of $\hat{\eta}^{h}$;

$\displaystyle\frac{dS^{2}(\hat{\eta}^{h})}{d\hat{\sigma}^{2}}=L(X^{\prime}\hat{V}^{-1}X)^{-1}L^{\prime},$
$\displaystyle\frac{dS^{2}(\hat{\eta}^{h})}{d\text{vech}(\hat{D}_{k})}=\hat{\sigma}^{2}\mathcal{D}_{q_{k}}^{\prime}\bigg{(}\sum_{j=1}^{l_{k}}\hat{B}_{(k,j)}\otimes\hat{B}_{(k,j)}\bigg{)},$

where $\hat{B}_{(k,j)}$ is given by:

$\hat{B}_{(k,j)}=Z_{(k,j)}^{\prime}\hat{V}^{-1}X(X^{\prime}\hat{V}^{-1}X)^{-1}L^{\prime}.$

We obtain an expression for var$(\hat{\eta}^{h})$ by noting that the
asymptotic variance of $\hat{\eta}^{h}$ is given by
$\mathcal{I}(\hat{\eta}^{h})^{-1}$ where $\mathcal{I}(\hat{\eta}^{h})$ is a
sub-matrix of $\mathcal{I}(\hat{\theta}^{h})$, given by equations (10)-(12).

In summary, we have provided all the closed-form expressions necessary to
perform the Satterthwaite degrees of freedom estimation method for any LMM
described by (1) and (2). For the remainder of this work, estimation of $v$
using the above expressions is referred to as the “direct-SW” method. This
name reflects the direct approach taken for the evaluation of the right-hand
side of (32). We note that this approach may also be extended to models using
the constrained covariance structures of Section 2.5 by employing the Fisher
Information matrix given by equation 30 and transforming the derivative of
$S^{2}(\eta)$ appropriately (details of which are provided by Theorem 6.9 of
Appendix 6.7.3). We conclude this section by noting that this method can also
be used in a similar manner for estimating the degrees of freedom of an
approximate F-statistic based on the multi-factor LMM, as described in
Appendix 6.7.2.

###  2.7 Simulation methods

To assess the accuracy and efficiency of each of the proposed LMM parameter
estimation methods described in Sections 2.1.1-2.1.5 and the direct-SW degrees
of freedom estimation method described in Section 2.6, extensive simulations
were conducted. The methods used to perform these simulations are described in
this section, with the results of the simulations presented in Section 3.1.
All reported results were obtained using an Intel(R) Xeon(R) Gold 6126 2.60GHz
processor with 16GB RAM.

####  2.7.1 Parameter estimation

The algorithms of Sections 2.1.1-2.1.5 have been implemented in the
programming language Python. Results taken across three simulation settings,
each with a different design structure, are provided for parameter estimation
performed using each of the algorithms. The Fisher Scoring algorithms
presented in this work are compared against one another, the R package lmer,
and the baseline truth used to generate the simulations. All methods are
contrasted in terms of output, computation time and, for the methods presented
in this paper, the number of iterations until convergence. Convergence of each
method was assessed by verifying whether successive log-likelihood estimates
differed by more than a predefined tolerance of $10^{-6}$. $1000$ individual
simulation instances were run for each simulation setting.

All model parameters were held fixed across all runs. In every simulation
setting, test data was generated according to model (1). Each of the three
simulation settings imposed a different structure on the random effects design
and covariance matrices, $Z$ and $D$. The first simulation setting employed a
single factor design ($r=1$) with two random effects (i.e. $q_{1}=2$) and $50$
levels (i.e. $l_{1}=50$). The second simulation setting employed two crossed
factors ($r=2$) where the number of random effects and numbers of levels for
each factor were given by $q_{1}=3$, $q_{2}=2$, $l_{1}=100$ and $l_{2}=50$,
respectively. The third simulation setting used three crossed factors (i.e.
$r=3$) with the number of random effects and levels for each of the factors
given by $q_{1}=4$, $q_{2}=3$, $q_{3}=2$, $l_{1}=100$, $l_{2}=50$ and
$l_{3}=10$, respectively. In all simulations, the number of observations, $n$,
was held fixed at $1000$. In each simulated design, the first column of the
fixed effects design matrix, $X$, and the first random effect in the random
effects design matrix, $Z$, were treated as intercepts. Within each simulation
instance, the remaining (non-zero) elements of the variables $X$ and $Z$, as
well as those of the error vector $\epsilon$, were generated at random
according to the standard univariate normal distribution. The random effects
vector $b$ was simulated according to a normal distribution with covariance
$D$, where $D$ was predefined, exhibited no particular constrained structure
and contained a mixture of both zero and non-zero off-diagonal elements. The
assignment of observations to levels for each factor was performed at random
with the probability of an observation belonging to any specific level held
constant and uniform across all levels.

Fisher Scoring methods for the single-factor design have been well-studied
(c.f. Demidenko (2013)) and the inclusion of the first simulation setting is
only for comparison purposes. The second and third simulation settings
represent more complex crossed factor designs for which the proposed Fisher
Scoring based parameter estimation methods did not previously exist. To assess
the performance of the proposed algorithms, the Mean Absolute Error ($MAE$)
and Mean Relative Difference ($MRD$) were used as performance metrics. The
methods considered for parameter estimation were those described in Sections
2.1.1-2.1.5 and the baseline truth used for comparison was either the baseline
truth used to generate the simulated data or the lmer computed estimates.

All methods were contrasted in terms of the $MAE$ and $MRD$ of both $\beta$
and the variance product $\sigma^{2}D$. The variance product $\sigma^{2}D$ was
chosen for comparison instead of the individual components $\sigma^{2}$ and
$D$ as the variance product $\sigma^{2}D$ is typically of intrinsic interest
during the inference stage of conventional statistical analyses and is often
employed for further computation. For each of the methods proposed in Sections
2.1.1-2.1.5, both ML and ReML estimation variants of the methods were
assessed, where the ReML variants of each method were performed according to
the adjustments detailed in Appendix 6.4. The computation time for each method
was also recorded and contrasted against lmer. To ensure a fair comparison of
computation time, the recorded computation times for lmer were based only on
the performance of the ‘optimizeLmer’ function, which is the module employed
for parameter estimation by lmer.

The specific values of $\beta,\sigma^{2},$ and $D$ employed for each
simulation setting can be found in Section $S1$ of the Supplementary Material.
Formal definitions of MAE and MRD measures used for comparison are given in
Section $S2$ of the Supplementary Material.

####  2.7.2 Degrees of freedom estimation

The accuracy and validity of the direct-SW degrees of freedom estimation
method proposed in Section 2.6 were assessed through further simulation. To
achieve this, as LMM degrees of freedom are assumed to be specific to the
experiment design, a single design was chosen at random from each of the three
simulation settings described in Section 2.7.1. With the fixed effects and
random effects matrices, $X$ and $Z$, now held constant across simulations,
$1000$ simulations were run for each simulation setting. The random effects
vector, $b$, and random error vector, $\epsilon$, were allowed to vary across
simulations according to normal distributions with the appropriate variances.
In each simulation, degrees of freedom were estimated via the direct-SW method
for a predefined contrast vector, corresponding to a fixed effect that was
truly zero.

The direct-SW estimated degrees of freedom were compared to a baseline truth
and the degrees of freedom estimates produced by the R package lmerTest.
Baseline truth was established in each simulation setting using $1,000,000$
simulations to empirically estimate $\text{Var}(S^{2}(\hat{\eta}))$, giving a
single value for the denominator of $v(\hat{\eta})$ from Equation (31).
Following this, in each of the $1000$ simulation instances described above,
the numerator of (31) was recorded, giving $1000$ estimates of
$v(\hat{\eta})$. The final estimate was obtained as an average of these $1000$
values. All lmerTest degrees of freedom estimates were obtained using the same
simulated data as was used by the direct-SW method and were computed using the
‘contest1D’ function from the lmerTest package. Whilst all baseline truth
measures were computed using parameter estimates obtained using the FSFS
algorithm of Section 2.1.4, we believe this has not induced bias into the
simulations as, as discussed in Section 3.1.1, all simulated FSFS-derived
parameter estimates considered in this work agreed with those provided by lmer
within a tolerance level on the scale of machine error.

###  2.8 Real data methods

To illustrate the usage of the methodology presented in Sections 2.1-2.6 in
practical situations, we provide two real data examples. These examples are
described fully in this section and the results are presented in Section 3.2.
Again, all reported results were obtained using an Intel(R) Xeon(R) Gold 6126
2.60GHz processor with 16GB RAM. For each reported computation time, averages
were taken across $50$ repeated runs of the given analysis.

####  2.8.1 The SAT score example

The first example presented here is based on data from the longitudinal
evaluation of school change and performance (LESCP) dataset (Turnbull et al.
(1999)). This dataset has notably been previously analyzed by Hong and
Raudenbush (2008) and was chosen for inclusion in this work because it
previously formed the basis for between-software comparisons of LMM software
packages by West et al. (2014). The LESCP study was conducted in $67$ American
schools in which SAT (Student Aptitude Test) math scores were recorded for
randomly selected samples of students. As in West et al. (2014), one of the
$67$ schools from this dataset was chosen as the focus for this analysis. For
each individual SAT score, unique identifiers (ID’s) for the student who took
the test and for the teacher who prepared the student for the test were
recorded. As many students were taught by multiple teachers and all teachers
taught multiple students, the grouping of SAT scores by student ID and the
grouping of SAT scores by teacher ID constitute two crossed factors. In total,
$n=234$ SAT scores were considered for analysis. The SAT scores were taken
from $122$ students taught by a total of $12$ teachers, with each student
sitting between $1$ and $3$ tests.

The research question in this example concerns how well a student’s grade
(i.e. year of schooling) predicted their mathematics SAT score (i.e. did
students improve in mathematics over the course of their education?). For this
question, between-student variance and between-teacher variance must be taken
into consideration as different students possess different aptitudes for
mathematics exams and different teachers possess different aptitudes for
teaching mathematics. In practice, this is achieved by employing an LMM which
includes random intercept terms for both the grouping of SAT scores according
to student ID and the grouping of SAT scores according to teacher ID. For the
$k^{th}$ mathematics SAT test taken by the $i^{th}$ student, under the
supervision of the $j^{th}$ teacher, such a model could be stated as follows:

$\text{MATH}_{i,j,k}=\beta_{0}+\beta_{1}\times\text{YEAR}_{i,j,k}+s_{i}+t_{j}+\epsilon_{i,j,k},$

where MATHi,j,k is the SAT score achieved and YEARi,j,k is the grade of the
student at the time the test was taken. In the above model, $\beta_{0}$ and
$\beta_{1}$ are unknown parameters and $s_{i}$, $t_{j}$ and $\epsilon_{i,j,k}$
are independent mean-zero random variables which differ only in terms of their
covariance. $s_{i}$ is the random intercept which models between-student
variance, $t_{j}$ is the random intercept which models between-teacher
variance, and $\epsilon_{i,j,k}$ is the random error term. The random
variables $s_{i}$, $t_{j}$ and $\epsilon_{i,j,k}$ are assumed to be mutually
independent and follow the below distributions:

$\displaystyle s_{i}\sim N(0,\sigma^{2}_{s}),\hskip 28.45274ptt_{j}\sim
N(0,\sigma^{2}_{t}),$ $\displaystyle\epsilon_{i,j,k}\sim N(0,\sigma^{2}),$

where the parameters $\sigma^{2}_{s}$, $\sigma^{2}_{t}$ and $\sigma^{2}$ are
the unknown student, teacher and residual variance parameters, respectively.

The random effects in the SAT score model can be described using the notation
presented in previous sections as follows; $r=2$ (i.e. observations are
grouped by two factors, student ID and teacher ID), $q_{1}=q_{2}=1$ (i.e. one
random effect is included for each factor, the random intercepts $s_{i}$ and
$t_{j}$, respectively), and $l_{1}=122,l_{2}=12$ (i.e. there are $122$
students and $12$ teachers). When the model is expressed in the form described
by equations (1) and (2), the random effects design matrix $Z$ is a $0-1$
matrix. In this setting, the positioning of the non-zero elements in $Z$
indicates the student and teacher associated with each test score. The random
effects covariance matrices for the two factors, student ID and teacher ID,
are given by $D_{0}=[\sigma^{2}_{s}]$ and $D_{1}=[\sigma^{2}_{t}]$,
respectively.

For the SAT score model, the estimated parameters obtained using each of the
methods detailed in Sections 2.1.1-2.1.5 are reported. For comparison, the
parameter estimates obtained by the R package lmer are also given. As the
estimates for the fixed effects parameters, $\beta_{0}$ and $\beta_{1}$, are
of primary interest, ML was employed to obtain all reported parameter
estimates. For this example, methods are contrasted in terms of output,
computation time and the number of iterations performed.

####  2.8.2 The twin study example

The second example presented in this work aims to demonstrate the flexibility
of the constrained covariance approaches described in Section 2.5. This
example is based on the ACE model; an LMM commonly employed to analyze the
results of twin studies by accounting for the complex covariance structure
exhibited between related individuals. The ACE model achieves this by
separating between-subject response variation into three categories: variance
due to additive genetic effects ($\sigma^{2}_{a}$), variance due to common
environmental factors ($\sigma^{2}_{c}$), and residual error
($\sigma^{2}_{e}$). For this example, we utilize data from the Wu-Minn Human
Connectome Project (HCP) (Van Essen et al. (2013)). The HCP dataset contains
brain imaging data collected from $1,200$ healthy young adults, aged between
$22$ and $35$, including data from $300$ twin pairs and their siblings. We do
not make use the imaging data in the HCP dataset but, instead, focus on the
baseline variables for cognition and alertness.

The primary research question considered in this example focuses on how well a
subject’s quality of sleep predicts their English reading ability. The
variables of primary interest used to address this question are subject scores
in the Pittsburgh Sleep Quality Index (PSQI) questionnaire and an English
language recognition test (ENG). Other variables included the subjects’ age in
years (AGE) and sex (SEX), as well as an age-sex interaction effect. A
secondary research question considered asks “How much of the between-subject
variance observed in English reading ability test scores can be explained by
additive genetic and common environmental factors?”. To address this question,
the covariance parameters $\sigma^{2}_{a}$, $\sigma^{2}_{c}$ and
$\sigma^{2}_{e}$ must be estimated.

To model the covariance components of the ACE model, a simplifying assumption
that all family units share common environmental factors was made. Following
the work of Winkler et al. (2015), family units were first categorized by
their internal structure into what shall be referred to as “family structure
types” (i.e. unique combinations of full-siblings, half-siblings and identical
twin pairs which form a family unit present in the HCP dataset). In the HCP
dataset, $19$ such family structure types were identified. In the following
text, each family structure type shall be treated as a factor in the model.
For the $i^{th}$ observation to belong to the $j^{th}$ level of the $k^{th}$
factor in the model may be interpreted as the $i^{th}$ subject belonging to
the $j^{th}$ family exhibiting family structure of type $k$. The model
employed for this example is given by:

$\displaystyle\text{ENG}_{k,j,i}=$
$\displaystyle\beta_{0}+\beta_{1}\times\text{AGE}_{i}+\beta_{2}\times\text{SEX}_{i}+\ldots$
$\displaystyle\beta_{3}\times\text{AGE}_{i}\times\text{SEX}_{i}+\beta_{4}\times\text{PSQI}_{i}+...$
$\displaystyle\gamma_{k,j,i}+\epsilon_{k,j,i},$

where both $\gamma_{k,j,i}$ and $\epsilon_{k,j,i}$ are mean-zero random
variables. The random error, $\epsilon_{k,j,i}$, is defined by;

$\displaystyle\epsilon_{k,j,i}\sim N(0,\sigma^{2}_{e}),$

and the random term $\gamma_{k,j,i}$ models the within-“family unit”
covariance. Further detail on the specification of $\gamma_{k,j,i}$ can be
found in Appendix (6.8.1).

In the notation of the previous sections, the number of factors in the ACE
model, $r$, is equal to the number of family structure types present in the
model (i.e. $19$). For each family structure type present in the model, family
structure type $k$, $l_{k}$ is the number of families who exhibit such
structure and $q_{k}$ is the number of subjects present in any given family
unit with such structure. As there is a unique random effect (i.e. a unique
random variable, $\gamma_{k,j,i}$) associated with each individual subject,
none of which are scaled by any coefficients, the random effects design
matrix, $Z$, is the $(n\times n)$ identity matrix. To describe
$\\{D_{k}\\}_{k\in\\{1,...,r\\}}$ requires the known matricies
$\mathbf{K}^{a}_{k}$ and $\mathbf{K}^{c}_{k}$, which specify the kinship
(expected genetic material) and environmental effects shared between
individuals, respectively (See Appendix 6.8.2 for more details). Given
$\mathbf{K}^{a}_{k}$ and $\mathbf{K}^{c}_{k}$, the covariance components
$\sigma^{2}$ and $\\{D_{k}\\}_{k\in\\{1,...,r\\}}$ are given as
$\sigma^{2}=\sigma^{2}_{e}$ and
$D_{k}=\sigma^{-2}_{e}(\sigma^{2}_{a}\mathbf{K}^{a}_{k}+\sigma^{2}_{c}\mathbf{K}^{c}_{k})$
respectively.

As the covariance components $\sigma_{a}$ and $\sigma_{c}$ are of practical
interest in this example, optimization is performed according to the ReML
criterion using the adjustments described in Appendix 6.4 and covariance
structure is constrained via the methods outlined in Section 2.5. Further
detail on the constrained approach for the ACE model can be found in Appendix
6.8.3. Discussion of computational efficiency for the ACE model is also
provided in Appendix 6.8.4.

Section 3.2.2 reports the maximized restricted log-likelihood values and
parameter estimates obtained using the Fisher Scoring method. Also given are
approximate T-statistics for each fixed effects parameter, alongside
corresponding degrees of freedom estimates and p-values obtained via the
methods outlined in Section 2.6. To verify correctness, the restricted log-
likelihood of the ACE model was also maximized numerically using the
implementation of Powell’s bi-directional search based optimization method
(Powell (1964)) provided by the SciPy Python package. The maximized restricted
log-likelihood values and parameter estimates produced were then contrasted
against those provided by the Fisher Scoring method. The OLS (Ordinary Least
Squares) estimates, which would be obtained had the additive genetic and
common environmental variance components not been accounted for in the LMM
analysis, are also provided for comparison.

##  3 Results

###  3.1 Simulation results

####  3.1.1 Parameter estimation results

The results of the parameter estimation simulations of 2.7.1 were identical.
In all three settings, all parameter estimates and maximised likelihood
criteria produced by the FS, FFS, SFS and FSFS methods were identical to those
produced by lmer up to a machine error level of tolerance. In terms of the
maximisation of likelihood criteria, there was no evidence to suggest that
lmer outperformed any of the FS, FFS, SFS or FSFS methods, or vice versa. The
CSFS method provided worse performance, however, with maximised likelihood
criteria which were lower than those reported by lmer in approximately $2.5\%$
of the simulations run for simulation setting $3$. The reported maximised
likelihood criteria for these simulations were all indicative of convergence
failure. For each parameter estimation method considered during simulation,
the observed performance was unaffected by the choice of likelihood estimation
criteria (i.e. ML or ReML) employed.

The $MRD$ and $MAE$ values provided evidence for strong agreement between the
parameter estimates produced by lmer and the FS, FFS, SFS and FSFS methods.
The largest $MRD$ values observed for these methods, taken relative to lmer,
across all simulations and likelihood criteria, were $1.03\times 10^{-3}$ and
$2.12\times 10^{-3}$ for $\beta$ and $\sigma^{2}D$, respectively. The largest
observed $MAE$ values for these methods, taken relative to lmer, across all
simulations, were given by$1.02\times 10^{-5}$ and $4.30\times 10^{-4}$ for
$\beta$ and $\sigma^{2}D$, respectively. These results provide clear evidence
that the discrepancies between the parameter estimates produced by lmer and
those produced by the FS, FFS, SFS and FSFS methods are extremely small in
magnitude. Due to the extremely small magnitudes of these differences, $MAE$
and $MRD$ values are not reported in further detail here. For more detail, see
Supplementary Material Sections $S3$-$S6$.

To compare the proposed methods in terms of computational performance, Tables
LABEL:perfTable and LABEL:perfTable2 present the average time and number of
iterations performed for each of the $5$ methods, using the ML and ReML
likelihood criteria, respectively. Corresponding computation times are also
provided for lmer. For both the ML and ReML likelihood criteria, all Fisher
Scoring methods demonstrated considerable efficiency in terms of computation
speed. While the improvement in speed is minor for simulation setting $1$, for
multi-factor simulation settings $2$ and $3$, the performance gains can be
seen to be considerable. Overall, the only method that consistently
demonstrated notably worse performance than the others was the CSFS algorithm.
In addition to requiring a longer computation time, the CSFS algorithm
employed many more iterations.

Within a setting where $q$ (i.e. the number number of columns in the random
effects design matrix) is small, the results of these simulations demonstrate
strong computational efficiency for the FS, FFS, SFS and FSFS methods.
However, we emphasise that no claim is made to suggest that the observed
results will generalise to larger values of $q$ or all possible designs. For
large sparse LMMs that include a large number of random effects, without
further adaptation of the methods presented in this work to employ sparse
matrix methodology, it is expected that lmer will provide superior
performance. We again emphasise here that our purpose in undertaking this work
is not to compete with the existing LMM parameter estimation tools. Instead,
the focus of this work is to develop methodology which provides performance
comparable to the existing tools but can be used in situations in which the
existing tools may not be applicable due to practical implementation
considerations.

Two potential causes for the poor performance of the CSFS method have been
identified below. The first cause is presented in Demidenko (2013), in which
it is argued that the Cholesky parameterised algorithm will provide slower
convergence than that of the unconstrained alternative methods due to
structural differences between the likelihood surfaces over which parameter
estimation is performed. This reasoning offers a likely explanation for the
higher number of iterations and computation time observed for the CSFS method
across all simulations. However, this does not explain the small number of
simulations in simulation setting $3$ in which evidence of convergence failure
was observed.

The second possible explanation for the poor performance of the CSFS algorithm
presents a potential reason for the observed convergence failure. Pinheiro and
Bates (1996) note that the Cholesky parameterisation of an arbitrary matrix is
not unique and the number of possible choices for the Cholesky factorisation
of a square positive definite matrix of dimension $(q\times q)$ increases with
the dimension $q$. This means that, for each covariance matrix $D_{k}$, there
are multiple Cholesky factors, $\Lambda_{k}$, which satisfy
$\Lambda_{k}\Lambda_{k}^{\prime}=D_{k}$. The larger $D_{k}$ is in dimension,
the greater the number of $\Lambda_{k}$ that correspond to $D_{k}$ there are.
Pinheiro and Bates (1996) argue that when optimal solutions are numerous and
close together in the parameter space, numerical problems can arise during
optimisation. We note that when compared with the other simulation settings
considered here, the design for simulation setting $3$ contains the largest
number of factors and random effects. Consequently, the covariance matrices,
$D_{k}$, are large for this simulation setting and numerous Cholesky factors
which correspond to the same optimal solution for the covariance matrix, $D$,
exist. For this reason, simulation setting $3$ is the most susceptible to
numerical problems of the kind described by Pinheiro and Bates (1996). We
suggest that this is a likely accounting for the $2.5\%$ of simulations in
which convergence failure was observed. In summary, the simulation results
offer strong evidence for the correctness and efficiency of the Fisher Scoring
methods proposed, with the exception of the CSFS method, which experiences
slower performance and convergence failure in rare cases.

####  3.1.2 Degrees of freedom estimation results

Method | Truth | Direct-SW |  lmerTest
---|---|---|---
Simulation $1$ |  |  |
Mean | $910.93$ | $910.68$ |  $906.62$
Standard Deviation | $0.0$ | $2.36$ |  $2.39$
Mean Squared Error | $0.0$ | $5.64$ |  $24.23$
Simulation $2$ |  |  |
Mean | $844.61$ | $842.17$ |  $837.79$
Standard Deviation | $0.0$ | $7.16$ |  $7.26$
Mean Squared Error | $0.0$ | $57.20$ |  $99.21$
Simulation $3$ |  |  |
Mean | $707.01$ | $700.96$ |  $695.08$
Standard Deviation | $0.0$ | $17.19$ |  $17.67$
Mean Squared Error | $0.0$ | $332.01$ |  $454.53$
Table 1: The mean, standard deviation and mean squared error for $1,000$
degrees of freedom estimates in each simulation setting. Results are displayed
for both the lmerTest and Direct-SW methods, alongside “true” mean values,
which were established using the moment-matching based approach outlined in
Section 2.7.2 and computed using $1,000,000$ simulation instances.

Across all degrees of freedom simulations, results indicated that the degrees
of freedom estimates produced by the direct-SW method possessed both lower
bias and lower variance than those produced by lmerTest. The degrees of
freedom estimates, for each of the three simulation settings, are summarized
in Table 1.

It can be seen from Table 1 that throughout all simulation settings both
direct-SW and lmerTest appear to underestimate the true value of the degrees
of freedom. However, the bias observed for the lmerTest estimates is notably
more severe that of the direct-SW method, suggesting that the estimates
produced by direct-SW have a higher accuracy than those produced by lmerTest.
The observed difference in the standard deviation of the degrees of freedom
estimates between the lmerTest and direct-SW methods is less pronounced.
However, in all simulation settings, lower standard deviations are reported
for direct-SW, suggesting that the estimates produced by direct-SW have a
higher precision than those produced by lmerTest.

For both lmerTest and direct-SW, the observed bias and variance increase with
simulation complexity. The simulations provided here indicate that the
severity of disagreement between direct-SW and lmerTest increases as the
complexity of the random effects design increases. This observation matches
the expectation that the accuracy of the numerical gradient estimation
employed by lmerTest will worsen as the complexity of the parameter space
increases. Reported mean squared errors are also provided, indicating further
that the direct-SW method outperformed lmerTest in terms of accuracy and
precision in all simulation settings.

###  3.2 Real data results

####  3.2.1 The SAT score results

For the SAT score example described in Section 2.8.1, the log-likelihood,
fixed effect parameters and variance components estimates produced by the
Fisher Scoring algorithms were identical to those produced by lmer. Further,
the reported computation times for this example suggested little difference
between all methods considered in terms of computational efficiency. Of the
Fisher Scoring methods considered, the FS and FFS methods took the fewest
iterations to converge while the SFS and FSFS methods took the most iterations
to converge. The results presented here exhibit strong agreement with those
reported in West et al. (2014), in which the same model was used as the basis
for a between-software comparison of LMM software packages. For completeness,
the full table of results can be found in Supplementary Material Section $S9$.

####  3.2.2 The twin study results

The results for the twin study example described in Section 2.8.2 are
presented in Table LABEL:AceTab. It can be seen from Table LABEL:AceTab that
the Powell optimizer and Fisher Scoring method attained extremely similar
optimized likelihood values, with Fisher Scoring converging notably faster.
This result offers further evidence for the correctness of the parameter
estimates produced by the Fisher Scoring method as it is unlikely that both
Powell estimation and Fisher Scoring would converge to the same solution if
the solution were suboptimal. The parameter estimates produced by Fisher
Scoring and Powell optimization can be seen to be smaller in magnitude than
those estimated by OLS, highlighting how the inclusion of additional variance
terms in the model can have a meaningful impact on the conclusion of the
analysis.

Also provided in Table LABEL:AceTab, are approximate t-tests based on the
Fisher Scoring and Powell optimizer parameter estimates, with corresponding
standard t-tests given based on the OLS estimates. At the $5\%$ significance
level, the approximate t-tests conclude that the fixed effects parameters
corresponding to the ‘Intercept’, ‘Age and Sex Interaction’ and the ‘PSQI
Score’ are non-zero in value. While it may be expected that age should affect
reading ability, no significant effect was observed for the ‘Age’ covariate.
This lack of observed effect may be explained by the narrow age range of
subjects present in the HCP dataset, with subjects ranging from 22 to 35 years
in age, and by the fact that individual observations were recorded in units of
years. In general, the OLS-based t-tests produced similar conclusions to those
produced by the FS-based approximate t-tests. A notable exception, however, is
given by the t-tests for the fixed effect associated with the ‘Sex’ covariate.
While the standard OLS t-test reported at the $5\%$ significance level that
the ‘Sex’ fixed effect was non-zero in value, the FS-based approximate t-test
concluded that there was not evidence to support this claim. This result
further highlights the importance of modeling all relevant variance terms.

##  4 Discussion

In this work, we have presented derivations for and demonstrated potential
applications of, score vector and Fisher Information matrix expressions for
the LMMs containing multiple random factors. While many of the examples
presented in this paper were benchmarked against existing software, it is not
the authors’ intention to suggest that the proposed methods are superior to
existing software packages. Instead, this work aims to complement existing LMM
parameter estimation research. This aim is realized through careful exposition
of the score vectors and Fisher Information matrices and detailed description
of methodology and algorithms.

Modern approaches to LMM parameter estimation typically depend upon
conceptually complex mathematical operations which require support from a
range of software packages and infrastructure. This work has been, in large
part, motivated by current deficiencies in vectorized support for such
operations. Vectorized computation is a crucial requirement for medical
imaging applications, in which hundreds of thousands of mixed models must be
estimated concurrently. It is with this application in mind that this work has
been undertaken. We intend to pursue the application itself in future work.

Although the methods presented in this work are well suited for the desired
application, we stress that there are many situations where current software
packages will likely provide superior performance. One such situation can be
found by observing that the Fisher Scoring method requires the storage and
inversion of a matrix of dimensions $(q\times q)$. This is problematic, both
in terms of memory and computation accuracy, for designs which involve very
large numbers of random effects, grouping factors or factor levels. While such
applications have not been considered extensively in this work, we note that
many of the expressions provided in this document may benefit from combination
with sparse matrix methodology to overcome this issue. We suggest that the
potential for improved computation time via the combination of the Fisher
Scoring approaches described in this paper with sparse matrix methodology may
also be the subject of future research.

##  5 Declarations

###  5.1 Funding

This work was supported by the Li Ka Shing Centre for Health Information and
Discovery and NIH grant [R01EB026859] (TMS, TN) and the Wellcome Trust award
[100309/Z/12/Z] (TN).

###  5.2 Conflicts of interest/Competing interests

Not applicable.

###  5.3 Availability of data and material

Data were provided in part by the Human Connectome Project, WU-Minn Consortium
(Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657)
funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for
Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at
Washington University.

The longitudinal evaluation of school change and performance (LESCP) dataset
employed in Sections 2.8.1 and 3.2.1 is publicly available and can be found,
for example, as one of the example datasets included in the Heirachical Linear
Models (HLM) software package (Raudenbush and Bryk (2002)).

###  5.4 Code availability

All code used to generate the results of the simulations and real data
examples presented in Sections 2.7-3.2.2 is publicly available and can be
found in the below GitHub repository:

https://github.com/TomMaullin/LMMPaper

Included in this repository are detailed notebooks demonstrating how the code
may be executed in order to reproduce the results presented in this work.

###  5.5 Authors’ contributions

Not applicable.

##  6 Appendix

###  6.1 Score vectors

In this appendix, we provide full derivations for the derivatives (9) and
(13). For derivatives (7) and (8), we note that the derivations for the multi-
factor LMM are identical to those given for the single-factor setting in
Demidenko (2013). As a result, we refer the interested reader to this source
for proofs. To obtain the derivatives (9) and (13), two lemmas, lemma 1 and
lemma 2, are required. Before stating the two lemmas, we briefly make explicit
the definition of derivative employed throughout this paper.

Throughout this work, for arbitrary matrices $A$ and $B$ of dimensions
$(p\times q)$ and $(m\times n)$ respectively, the following definition of
matrix (partial) derivative has been implicitly employed:

$\frac{\partial\text{vec}(A)}{\partial\text{vec}(B)}=\begin{bmatrix}\frac{\partial
A_{[1,1]}}{\partial B_{[1,1]}}&\ldots&\frac{\partial A_{[p,1]}}{\partial
B_{[1,1]}}&\ldots&\frac{\partial A_{[p,q]}}{\partial B_{[1,1]}}\\\
\vdots&\ddots&\vdots&\ddots&\vdots\\\ \frac{\partial A_{[1,1]}}{\partial
B_{[m,1]}}&\ldots&\frac{\partial A_{[p,1]}}{\partial
B_{[m,1]}}&\ldots&\frac{\partial A_{[p,q]}}{\partial B_{[m,1]}}\\\
\vdots&\ddots&\vdots&\ddots&\vdots\\\ \frac{\partial A_{[1,1]}}{\partial
B_{[m,n]}}&\ldots&\frac{\partial A_{[p,1]}}{\partial
B_{[m,n]}}&\ldots&\frac{\partial A_{[p,q]}}{\partial B_{[m,n]}}\\\
\end{bmatrix},$

with the equivalent definition holding for the total derivative. In addition,
for a scalar value, $a$, and matrix, $B$, defined as above, we employ the
below matrix (partial) derivative definition:

$\frac{\partial{a}}{\partial B}=\begin{bmatrix}\frac{\partial a}{\partial
B_{[1,1]}}&\ldots&\frac{\partial a}{\partial B_{[1,n]}}\\\
\vdots&\ddots&\vdots\\\ \frac{\partial a}{\partial
B_{[m,1]}}&\ldots&\frac{\partial a}{\partial B_{[m,n]}}\\\ \end{bmatrix},$

with, again, the equivalent definition holding for the total derivative. The
above notation can be useful since by definition, the two derivatives defined
above satisfy the below relation:

$\text{vec}\bigg{(}\frac{\partial{a}}{\partial B}\bigg{)}=\frac{\partial
a}{\partial\text{vec}(B)}.$

The above matrix derivative definitions are commonly employed in the
mathematical and statistical literature, but are not universal. For example,
sources such as Neudecker and Wansbeek (1983) and Magnus and Neudecker (1999)
define $\partial\text{vec}(A)/\partial\text{vec}(B)$ as the transpose of the
definition supplied above. This discrepancy between definitions is noteworthy
as some of the referenced results used in the following proofs are the
transpose of those found in the original source. For a full account and in-
depth discussion of the different definitions of matrix derivatives used in
the literature, we recommend the work of Turkington (2013).

###### Lemma 1

Let $g$ be a column vector, $A$ be a square matrix and $\\{B_{s}\\}$ be a set
of arbitrary matrices of equal size. Let $K$ be an unstructured matrix which
none of $g$, $A$ or any of the $\\{B_{s}\\}$ depend on. Further, assume
$A,\\{B_{s}\\},g$ and $K$ can be multiplied as required. The below matrix
derivative expression now holds;

$\displaystyle\frac{\partial}{\partial
K}\bigg{[}g^{\prime}(A+\sum_{t}B_{t}KB_{t}^{\prime})^{-1}g\bigg{]}=$ (33)
$\displaystyle-\sum_{s}B_{s}^{\prime}\big{(}A^{\prime}+\sum_{t}B_{t}KB_{t}^{\prime}\big{)}^{-1}gg^{\prime}\big{(}A^{\prime}+\sum_{t}B_{t}KB_{t}^{\prime}\big{)}^{-1}B_{s}.$

###### Proof

To begin, denote $M=A+\sum_{s}B_{s}KB_{s}^{\prime}$. For clarity and
convenience, in this proof, we shall briefly depart from the usual notation of
$H_{[i,j]}$ for the $(i,j)^{th}$ element of an arbitrary matrix $H$ and
instead employ Ricci calculus notation, in which the $(i,j)^{th}$ element of
$H$ is given by $H^{i}_{j}$ and all summations are made implicit. The left
hand side of (33) is now given, using Ricci calculus notation, as;

$\frac{\partial}{\partial
K^{m}_{n}}\bigg{[}(g^{\prime})_{i}(M^{-1})^{i}_{j}g^{j}\bigg{]}=(g^{\prime})_{i}g^{j}\frac{\partial}{\partial
K^{m}_{n}}\big{[}M^{-1}\big{]}^{i}_{j}.$ (34)

We employ a result from Giles (2008) which states, in Ricci calculus notation;

$\frac{\partial[F^{-1}]^{i}_{j}}{\partial x}=-(F^{-1})^{i}_{q}\frac{\partial
F^{q}_{p}}{\partial x}(F^{-1})^{p}_{j}.$

Applying this result to (34), it can be seen that the right hand side of (34)
is equal to:

$-(g^{\prime})_{i}g^{j}(M^{-1})^{i}_{q}\frac{\partial M^{q}_{p}}{\partial
K^{m}_{n}}(M^{-1})^{p}_{j}.$

Evaluating the derivative in the above expression and rearranging the terms
gives the below:

$\displaystyle-(B_{s}^{\prime})^{n}_{p}(M^{-1})^{p}_{j}g^{j}(g^{\prime})_{i}(M^{-1})^{i}_{q}B^{q}_{sm}$
$\displaystyle=-[B_{s}^{\prime}M^{-1}gg^{\prime}M^{-1}B_{s}]^{n}_{m}.$

By applying a transposition to the matrices in the above expression, in order
to interchange the indices $m$ and $n$, the result follows.∎

###### Lemma 2

Let $A,\\{B_{s}\\}$ and $K$ be defined as in Lemma 1. Then the following is
true:

$\frac{\partial}{\partial
K}\log|A+\sum_{t}B_{t}KB_{t}^{\prime}|=\sum_{s}B_{s}^{\prime}(A+\sum_{t}B_{t}K^{\prime}B_{t}^{\prime})^{-1}B_{s}.$
(35)

###### Proof

As in the proof of Lemma 1, the notation of
$M=A+\sum_{t}B_{t}K^{\prime}B_{t}^{\prime}$ will be adopted and the proof will
be given in Ricci calculus notation. Expanding the derivative with respect to
$K^{m}_{n}$ using the partial derivatives of the elements of $M$ gives the
following:

$\frac{\partial}{\partial K^{m}_{n}}\log|M|=\frac{\partial\log|M|}{\partial
M^{i}_{j}}\frac{\partial M^{i}_{j}}{\partial K^{m}_{n}}.$ (36)

We now make use of the following well-known result, which can be found in
Giles (2008).

$\frac{\partial\log|X|}{\partial X^{i}_{j}}=(X^{\prime-1})^{i}_{j}.$ (37)

Inserting (37) into (36) and noting that $A^{i}_{j}$ does not depend on
$K^{m}_{n}$, yields:

$\frac{\partial}{\partial
K^{m}_{n}}\log|M|=(M^{\prime-1})^{i}_{j}\frac{\partial(B_{s}K(B_{s})^{\prime})^{i}_{j}}{\partial
K^{m}_{n}}.$

Evaluating the derivative on the right hand side of the above now returns:

$\frac{\partial}{\partial
K^{m}_{n}}\log|M|=(M^{\prime-1})^{i}_{j}(B_{s})_{m}^{i}(B_{s}^{\prime})^{n}_{j}=(B^{\prime}_{s})_{i}^{m}(M^{\prime-1})^{i}_{j}(B_{s})^{j}_{n}.$

Combining the summation terms in the above expression yields the following.

$\frac{\partial}{\partial
K^{m}_{n}}\log|M|=(B^{\prime}_{s}M^{\prime-1}B_{s})^{m}_{n}.$

Finally, converting the above into matrix notation now gives (35) as desired.
∎

We are now in a position to derive the partial derivative matrix
$\frac{\partial l}{\partial D_{k}}$, the matrix derivative with respect to
$D_{k}$ which does not account for the equal elements of $D_{k}$ induced by
symmetry.

###### Theorem 6.1

The partial derivative matrix $\frac{\partial l}{\partial D_{k}}$ is given by
the following.

$\frac{\partial l(\theta)}{\partial
D_{k}}=\frac{1}{2}\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}V^{-1}\bigg{(}\frac{ee^{\prime}}{\sigma^{2}}-V\bigg{)}V^{-1}Z_{(k,j)}.$
(38)

###### Proof

To derive (38), we use the expression for the log-likelihood of the LMM, given
by (3). As the first term inside the brackets of (3) does not depend on
$D_{k}$, we need only consider the second and third term for differentiation.
We now note that, by (4) and the block diagonal structure of $D$, it can be
seen that:

$V=I+ZDZ^{\prime}=I+\sum_{k=1}^{r}\sum_{j=1}^{l_{r}}Z_{(k,j)}D_{k}Z_{(k,j)}^{\prime}.$
(39)

By substituting (39) into the second term of (3), it can be seen that:

$\displaystyle\sigma^{-2}e^{\prime}V^{-1}e=\sigma^{-2}e^{\prime}\bigg{(}I+\sum_{k=1}^{r}\sum_{j=1}^{l_{r}}Z_{(k,j)}D_{k}Z_{(k,j)}^{\prime}\bigg{)}^{-1}e.$

By applying the result of Lemma 1, it can now be seen that the partial
derivative matrix of the above, taken with respect to $D_{k}$, can be
evaluated to:

$\frac{\partial}{\partial
D_{k}}\bigg{[}\sigma^{-2}e^{\prime}V^{-1}e\bigg{]}=\sigma^{-2}\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}V^{-1}ee^{\prime}V^{-1}Z_{(k,j)}.$

Similarly, by substituting (39) into the third term of (3), it can be seen
that:

$\log|V|=\log\big{|}I+\sum_{k=1}^{r}\sum_{j=1}^{l_{r}}Z_{(k,j)}D_{k}Z_{(k,j)}^{\prime}\big{|}.$

By now applying the result of Lemma 2, the partial derivative matrix of the
above expression, taken with respect to $D_{k}$, can be seen to be given by:

$\frac{\partial}{\partial
D_{k}}\big{[}\log|V|\big{]}=\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}V^{-1}Z_{(k,j)}.$

Combining the previous two derivative expressions, we obtain (38) as desired.∎

Through applying the vectorisation operator to (38), the following corollary
is obtained. This is the result stated by (13) in Section 2.1.2.

###### Corollary 1

The partial derivative vector of the log-likelihood with respect to
vec$(D_{k})$ is given as:

$\frac{\partial
l(\theta)}{\partial\text{vec}(D_{k})}=\frac{1}{2}\text{vec}\bigg{(}\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}V^{-1}\bigg{(}\frac{ee^{\prime}}{\sigma^{2}}-V\bigg{)}V^{-1}Z_{(k,j)}\bigg{)}.$

Using Theorem 5.12 of Turkington (2013), which states that, in our notation:

$\frac{dl(\theta)}{d\text{vec}(D_{k})}=\mathcal{D}_{q_{k}}\mathcal{D}_{q_{k}}^{\prime}\frac{\partial
l(\theta)}{\partial\text{vec}(D_{k})},$

the following corollary is now obtained.

###### Corollary 2

The total derivative vector of the log-likelihood with respect to vec($D_{k}$)
is given as:

$\displaystyle\frac{dl(\theta)}{d\text{vec}(D_{k})}=$
$\displaystyle\frac{1}{2}\mathcal{D}_{q_{k}}\mathcal{D}^{\prime}_{q_{k}}\text{vec}\bigg{(}\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}V^{-1}\bigg{(}\frac{ee^{\prime}}{\sigma^{2}}-V\bigg{)}V^{-1}Z_{(k,j)}\bigg{)}.$

Finally, by noting that the vectorisation and half-vectorisation operators
satisfy the following relationship:

$\text{vec}(D_{k})=\mathcal{D}^{+}_{q_{k}}\text{vech}(D_{k}),$

the following corollary is obtained. This is the result stated by (9) in
Section 2.1.1.

###### Corollary 3

The total derivative vector of the log-likelihood with respect to
vech$(D_{k})$ is given as:

$\frac{dl(\theta)}{d\text{vech}(D_{k})}=\frac{1}{2}\mathcal{D}^{\prime}_{q_{k}}\text{vec}\bigg{(}\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}V^{-1}\bigg{(}\frac{ee^{\prime}}{\sigma^{2}}-V\bigg{)}V^{-1}Z_{(k,j)}\bigg{)}.$

This result diverges slightly from the proof given for the single-factor LMM
in Demidenko (2013), which uses the inverse duplication matrix
$\mathcal{D}^{+}_{q_{k}}$ in the place of the transposed duplication matrix,
$\mathcal{D}^{\prime}_{q_{k}}$, in the above. The authors believe that the
derivative given in Demidenko (2013) is, in fact, vech$\big{(}\frac{\partial
l(\theta)}{\partial\text{vec}(D_{k})}\big{)}$, and not
$\frac{dl(\theta)}{d\text{vech}(D_{k})}$ as stated. It should be noted,
however, that this is a minor amendment and has little impact on the
performance of the algorithm in practice.

###  6.2 Fisher Information matrix

In this appendix, we derive the Fisher Information matrix for the full
representation of the parameter vector, $\theta^{f}$. The derivation of the
elements $\mathcal{I}^{f}_{\beta}$, $\mathcal{I}^{f}_{\beta,\sigma^{2}}$ and
$\mathcal{I}^{f}_{\sigma^{2}}$ for the multi-factor LMM is identical to that
used for the single factor LMM given in Demidenko (2013) and, therefore, will
not be repeated here. Provided here are only derivations of components of the
Fisher Information matrix which relate to $D_{k}$, for some factor $f_{k}$. To
derive these results, we will follow a similar argument to that used by
Demidenko (2013) for the single-factor setting and will begin by noting that
the Fisher Information matrix elements,
$\\{\mathcal{I}_{a,b}^{f}\\}_{a,b\in\\{\beta,\sigma^{2},\text{vec}(D_{1}),\ldots,\text{vec}(D_{r})\\}}$,
can be expressed by the following formula:

$\mathcal{I}_{a,b}^{f}=\text{cov}\bigg{(}\frac{\partial
l(\theta^{f})}{\partial a},\frac{\partial l(\theta^{f})}{\partial b}\bigg{)}.$

Theorems 6.2 and 6.3 provide derivation of equation (14) and Theorem 6.4
provides derivation of equation (15). Following this, Corollaries 4-6 detail
the derivation of equations (11) and (12).

###### Theorem 6.2

For any arbitrary integer $k$ between $1$ and $r$, the covariance of the
partial derivatives of $l(\theta^{f})$ with respect to $\beta$ and
vec$(D_{k})$ is given by:

$\mathcal{I}^{f}_{\beta,\text{vec}(D_{k})}=\text{cov}\bigg{(}\frac{\partial
l(\theta^{f})}{\partial\beta},\frac{\partial
l(\theta^{f})}{\partial\text{vec}(D_{k})}\bigg{)}=\mathbf{0}_{p,q_{k}^{2}}.$

###### Proof

First, let $u$ and $T_{(k,j)}$ denote the following quantities:

$u=\sigma^{-1}V^{-\frac{1}{2}}e,\hskip
14.22636ptT_{(k,j)}=Z_{(k,j)}^{\prime}V^{-\frac{1}{2}}.$ (40)

As $e\sim N(0,\sigma^{2}V)$, it follows trivially that $u\sim N(0,I_{n})$. Let
$c$ be an arbitrary column vector of length $q_{k}$. Noting that the
derivative of the log-likelihood function has mean zero, the below can be seen
to be true:

$\text{cov}\bigg{(}\frac{\partial
l(\theta|y)}{\partial\beta},c^{\prime}\frac{\partial l(\theta|y)}{\partial
D_{k}}c\bigg{)}=\mathbb{E}\bigg{[}\frac{\partial
l(\theta|y)}{\partial\beta}c^{\prime}\frac{\partial l(\theta|y)}{\partial
D_{k}}c\bigg{]}.$

By rewriting in terms of $u$ and $T_{(k,j)}$, and noting that
$\mathbb{E}[u]=0$, the right hand side of the above can be seen to simplify
to:

$\mathbb{E}\bigg{[}\sigma^{-1}XV^{-\frac{1}{2}}uc^{\prime}\bigg{(}\frac{1}{2}\sum_{j=1}^{l_{k}}(T_{(k,j)}u)(T_{(k,j)}u)^{\prime}\bigg{)}c\bigg{]}$
$=\frac{1}{2}\sigma^{-1}XV^{-\frac{1}{2}}\sum_{j=1}^{l_{k}}\mathbb{E}\bigg{[}uc^{\prime}T_{(k,j)}uu^{\prime}\bigg{]}T_{(k,j)}^{\prime}c=\mathbf{0}_{p,q_{k}^{2}}.$

That the above is equal to a matrix of zeros, follows directly from noting
that the third moment of the Normal distribution is $0$. The result of the
theorem now follows.∎

###### Theorem 6.3

For any arbitrary integer $k$ between $1$ and $r$, the covariance of the
partial derivatives of $l(\theta^{f})$ with respect to $\sigma^{2}$ and
vec$(D_{k})$ is given by:

$\displaystyle\mathcal{I}^{f}_{\sigma^{2},\text{vec}(D_{k})}=\text{cov}\bigg{(}\frac{\partial
l(\theta^{f})}{\partial\sigma^{2}},\frac{\partial
l(\theta^{f})}{\partial\text{vec}(D_{k})}\bigg{)}=$
$\displaystyle\frac{1}{2\sigma^{2}}\text{vec}^{\prime}\bigg{(}\sum_{j=1}^{l_{k}}Z^{\prime}_{(k,j)}V^{-1}Z_{(k,j)}\bigg{)}.$

###### Proof

To begin, $u$ and $T_{(k,j)}$ are defined as in (40). Rewriting the covariance
in Theorem 6.3 in terms of (40) and removing constant terms gives:

$\frac{1}{4\sigma^{2}}\text{cov}\bigg{(}u^{\prime}u,\text{vec}\bigg{(}\sum_{j=1}^{l_{k}}(T_{(k,j)}u)(T_{(k,j)}u)^{\prime}\bigg{)}\bigg{)}.$

Noting that the Kronecker product satisfies the property
vec$(aa^{\prime})=a\otimes a$, for all column vectors $a$, we obtain that the
above is equal to:

$\frac{1}{4\sigma^{2}}\text{cov}\bigg{(}u^{\prime}u,\sum_{j=1}^{l_{k}}\bigg{[}(T_{(k,j)}u)\otimes(T_{(k,j)}u)\bigg{]}\bigg{)}.$

We now note that the Kronecker product satisfies vec$(ABC)=(C^{\prime}\otimes
A)$vec$(B)$ for arbitrary matrices $A$, $B$ and $C$ of appropriate dimensions.
Utilizing this and applying the mixed product property of the Kronecker
product to the above expression yeilds:

$\frac{1}{4\sigma^{2}}\text{cov}\bigg{(}\text{vec}^{\prime}(I_{n})(u\otimes
u),\sum_{j=1}^{l_{k}}\bigg{[}(T_{(k,j)}\otimes T_{(k,j)})(u\otimes
u)\bigg{]}\bigg{)}.$

By moving constant values outside of the covariance function using standard
results the above now becomes:

$\frac{1}{4\sigma^{2}}\text{vec}^{\prime}(I_{n})\text{cov}(u\otimes
u)\sum_{j=1}^{l_{k}}\bigg{[}(T_{(k,j)}\otimes
T_{(k,j)})\bigg{]}^{{}^{\prime}}.$

Noting that $u\sim N(0,I_{n})$ we now employ a result from Magnus and
Neudecker (1986) which states that $\text{cov}(u\otimes u)=2N_{n}$.
Substituting this result into the previous expression gives the below:

$\frac{1}{2\sigma^{2}}\text{vec}^{\prime}(I_{n})N_{n}\sum_{j=1}^{l_{k}}\bigg{[}(T_{(k,j)}\otimes
T_{(k,j)})\bigg{]}^{\prime}.$

We now note a further result given in Magnus and Neudecker (1986); the matrix
$N_{n}$ satisfies $(A\otimes A)N_{k}=N_{n}(A\otimes A)$ for all $A$, $n$ and
$k$ such that the resulting matrix multiplications are well defined. Applying
this result to the above expression, the below can be obtained:

$\frac{1}{2\sigma^{2}}\text{vec}^{\prime}(I_{n})\sum_{j=1}^{l_{k}}\bigg{[}(T_{(k,j)}\otimes
T_{(k,j)})\bigg{]}^{\prime}N_{q_{k}}.$

Finally, again using the relationship vec$(ABC)=(C^{\prime}\otimes A)$vec$(B)$
which is satisfied by the Kronecker product for arbitrary matrices $A$, $B$
and $C$ of appropriate dimensions, the above can be seen to reduce to:

$\frac{1}{2\sigma^{2}}\text{vec}^{\prime}\bigg{(}\sum_{j=1}^{l_{k}}T_{(k,j)}T_{(k,j)}^{\prime}\bigg{)}N_{q_{k}}.$

Using the definition of $T_{(k,j)}$, given in (40), and noting that the matrix
$N_{q_{k}}$ satisfies the relationship
$\text{vec}^{\prime}(A)N_{q_{k}}=\text{vec}^{\prime}(A)$ for all appropriately
sized symmetric matrices $A$, the result now follows. ∎

###### Theorem 6.4

For any arbitrary integers $k_{1}$ and $k_{2}$ between $1$ and $r$, the
covariance of the partial derivatives of $l(\theta^{f})$ with respect to
vec$(D_{k_{1}})$ and vec$(D_{k_{2}})$ is given by:

$\displaystyle\mathcal{I}^{f}_{\text{vec}(D_{k_{1}}),\text{vec}(D_{k_{2}})}=\text{cov}\bigg{(}\frac{\partial
l(\theta^{f})}{\partial\text{vec}(D_{k_{1}})},\frac{\partial
l(\theta^{f})}{\partial\text{vec}(D_{k_{2}})}\bigg{)}=$
$\displaystyle\frac{1}{2}N_{q_{k_{1}}}\sum_{j=1}^{l_{k_{2}}}\sum_{i=1}^{l_{k_{1}}}(Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)}\otimes
Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)}).$

###### Proof

By Corollary 1, and using the notation introduced in (40), it can be seen that
the partial derivative vector with respect to vec$(D_{k})$ is given by:

$\frac{\partial
l(\theta|y)}{\partial\text{vec}(D_{k})}=\text{vec}\bigg{(}\frac{1}{2}\sum_{j=1}^{l_{k}}(T_{(k,j)}u)(T_{(k,j)}u)^{\prime}-\frac{1}{2}\sum_{j=1}^{l_{k}}T_{(k,j)}T_{(k,j)}^{\prime}\bigg{)}.$

By substituting the above into the covariance expression appearing in Theorem
4 and removing constant terms, the below is obtained:

$\displaystyle\frac{1}{4}\text{cov}\bigg{(}\sum_{j=1}^{l_{k_{1}}}\text{vec}\bigg{[}(T_{(k_{1},j)}u)(T_{(k_{1},j)}u)^{\prime}\bigg{]},$
$\displaystyle\hskip
51.21504pt\sum_{j=1}^{l_{k_{2}}}\text{vec}\bigg{[}(T_{(k_{2},j)}u)(T_{(k_{2},j)}u)^{\prime}\bigg{]}\bigg{)}.$

Noting that for any column vector $a$, the Kronecker product satisfies the
property vec$(aa^{\prime})=a\otimes a$, the above can be seen to be equal to:

$\displaystyle\frac{1}{4}\text{cov}\bigg{(}\sum_{j=1}^{l_{k_{1}}}\bigg{[}(T_{(k_{1},j)}u)\otimes(T_{(k_{1},j)}u)\bigg{]},$
$\displaystyle\hskip
51.21504pt\sum_{j=1}^{l_{k_{2}}}\bigg{[}(T_{(k_{2},j)}u)\otimes(T_{(k_{2},j)}u)\bigg{]}\bigg{)}.$

Using the mixed product property of the Kronecker product and moving constant
matrices outside of the covariance function, using standard results, gives:

$\frac{1}{4}\sum_{j=1}^{l_{k_{1}}}\bigg{[}(T_{(k_{1},j)}\otimes
T_{(k_{1},j)})\bigg{]}\text{cov}(u\otimes
u)\sum_{j=1}^{l_{k_{2}}}\bigg{[}(T_{(k_{2},j)}\otimes
T_{(k_{2},j)})\bigg{]}^{\prime}.$

Noting that $u\sim N(0,I_{n})$ and using a result from Magnus and Neudecker
(1986) which states that cov$(u\otimes u)=2N_{n}$, yields:

$\frac{1}{4}\sum_{j=1}^{l_{k_{1}}}\bigg{[}(T_{(k_{1},j)}\otimes
T_{(k_{1},j)})\bigg{]}2N_{n}\sum_{j=1}^{l_{k_{2}}}\bigg{[}(T_{(k_{2},j)}^{\prime}\otimes
T_{(k_{2},j)}^{\prime})\bigg{]}.$

We now utilize the fact that, as noted by Magnus and Neudecker (1986), the
matrix $N_{n}$ satisfies $(A\otimes A)N_{k}=N_{n}(A\otimes A)$ for all $A$,
$n$ and $k$ such that the resulting matrix multiplications are well defined.
Using this and the mixed product property of the Kronecker product the above
can now be seen to be equal to:

$\frac{1}{2}N_{q_{k_{1}}}\sum_{j=1}^{l_{k_{2}}}\sum_{i=1}^{l_{k_{1}}}\bigg{[}(T_{(k_{1},i)}T_{(k_{2},j)}^{\prime})\otimes(T_{(k_{1},i)}T_{(k_{2},j)}^{\prime})\bigg{]}.$

From the definition of $T_{(k,j)}$, it can be seen that the above is equal to
the result of Theorem 4. ∎

We now turn attention to the derivation of (11) and (12). As in the proofs of
Corollaries 2 and 3 of Appendix 6.1, we begin by noting that:

$\displaystyle\frac{dl(\theta)}{d\text{vech}(D_{k})}$
$\displaystyle=\mathcal{D}^{+}_{q_{k}}\frac{dl(\theta)}{d\text{vec}(D_{k})}$
$\displaystyle=\mathcal{D}^{+}_{q_{k}}\mathcal{D}_{q_{k}}\mathcal{D}_{q_{k}}^{\prime}\frac{\partial
l(\theta)}{\partial\text{vec}(D_{k})}$
$\displaystyle=\mathcal{D}_{q_{k}}^{\prime}\frac{\partial
l(\theta)}{\partial\text{vec}(D_{k})},$

where the first equality follows from the definition of the duplication matrix
and the second equality follows from Theorem 5.12 of Turkington (2013).
Applying the above identity to Theorems 6.2, 6.3 and 6.4 and moving the matrix
$\mathcal{D}_{q_{k}}^{\prime}$ outside the covariance function in each, leads
to the following three corollaries which, when taken in combination, provide
equations (11) and (12).

###### Corollary 4

For any arbitrary integer $k$ between $1$ and $r$, the covariance of the total
derivatives of $l(\theta^{h})$ with respect to $\beta$ and vech$(D_{k})$ is
given by:

$\mathcal{I}^{h}_{\beta,\text{vech}(D_{k})}=\text{cov}\bigg{(}\frac{dl(\theta^{h})}{d\beta},\frac{dl(\theta^{h})}{d\text{vech}(D_{k})}\bigg{)}=\mathbf{0}_{p,q_{k}(q_{k}+1)/2}.$

###### Corollary 5

For any arbitrary integer $k$ between $1$ and $r$, the covariance of the total
derivatives of $l(\theta^{h})$ with respect to $\sigma^{2}$ and vech$(D_{k})$
is given by:

$\displaystyle\mathcal{I}^{h}_{\sigma^{2},\text{vech}(D_{k})}=\text{cov}\bigg{(}\frac{dl(\theta^{h})}{d\sigma^{2}},\frac{dl(\theta^{h})}{d\text{vech}(D_{k})}\bigg{)}=$
$\displaystyle\frac{1}{2\sigma^{2}}\text{vec}^{\prime}\bigg{(}\sum_{j=1}^{l_{k}}Z^{\prime}_{(k,j)}V^{-1}Z_{(k,j)}\bigg{)}\mathcal{D}_{q_{k}}.$

###### Corollary 6

For any arbitrary integers $k_{1}$ and $k_{2}$ between $1$ and $r$, the
covariance of the total derivatives of $l(\theta^{h})$ with respect to
vech$(D_{k_{1}})$ and vech$(D_{k_{2}})$ is given by:

$\displaystyle\mathcal{I}^{h}_{\text{vech}(D_{k_{1}}),\text{vech}(D_{k_{2}})}=\text{cov}\bigg{(}\frac{dl(\theta^{h})}{d\text{vech}(D_{k_{1}})},\frac{dl(\theta^{h})}{d\text{vech}(D_{k_{2}})}\bigg{)}=$
$\displaystyle\frac{1}{2}\mathcal{D}^{\prime}_{q_{k_{1}}}\sum_{j=1}^{l_{k_{2}}}\sum_{i=1}^{l_{k_{1}}}(Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)}\otimes
Z^{\prime}_{(k_{1},i)}V^{-1}Z_{(k_{2},j)})\mathcal{D}_{q_{k_{2}}}.$

We note that the results of Corollaries 4, 5 and 6 do not contain the matrix
$N_{q_{k}}$, which appears in the corresponding theorems (Theorems 6.2, 6.3
and 6.4). This is due to another result of Magnus and Neudecker (1986), which
states that $\mathcal{D}^{\prime}_{k}N_{k}=\mathcal{D}^{\prime}_{k}$ for any
integer $k$. This concludes the derivations of Fisher Information matrix
expressions given in Sections 2.1.1-2.1.4.

###  6.3 Full Fisher Scoring

In this appendix, we prove that equations (16) and (21) are valid Fisher
Scoring rules of the form given by equation (6). To achieve this, we first
require the following lemma.

###### Lemma 3

Let $n,p$ and $k$ be arbitrary positive integers with $p<n$. Let $\tilde{N}$
be of dimension $(n\times n)$ with rank$(\tilde{N})<n$, $\tilde{K}$ be of
dimension $(n\times p)$ with rank$(\tilde{K})=p$, $C$ be of dimension
$(n\times n)$ with rank$(\tilde{N})=n$ and $\delta$ be of size $(n\times k)$.
If the following properties are satisfied:

  1. 1.

$\tilde{K}\tilde{K}^{+}=\tilde{N}$, where
$\tilde{K}^{+}=(\tilde{K}^{\prime}\tilde{K})^{-1}\tilde{K}^{\prime}$ is the
generalized inverse of $\tilde{K}$.

  2. 2.

$\tilde{N}$ is symmetric and idempotent.

  3. 3.

$\tilde{N}C=C\tilde{N}=\tilde{N}C\tilde{N}$.

  4. 4.

$\tilde{N}\delta=\delta.$

then it follows that:

$(\tilde{N}C)^{+}\delta=C^{-1}\delta.$

###### Proof

In proving the above lemma, we make use of the following result given by
Barnett (1990);

  * $\bullet$

If $A$ is a matrix of dimension $(m\times l)$ with rank$(A)=l<m$, and $B$ is a
matrix of dimension $(l\times m)$ with rank$(B)=l$, then;

$(AB)^{+}=B^{\prime}(BB^{\prime})^{-1}(A^{\prime}A)^{-1}A^{\prime}.$ (41)

By property 1, we have that:

$(\tilde{N}C)^{+}=(\tilde{K}\tilde{K}^{+}C)^{+}.$

By applying (41) to the above, with $A=\tilde{K}$ and $B=\tilde{K}^{+}C$, and
simplifying the result, the following is obtained:

$(\tilde{K}\tilde{K}^{+}C)^{+}=C^{\prime}\tilde{K}^{+^{\prime}}(\tilde{K}^{+}CC^{\prime}\tilde{K}^{+^{\prime}})^{-1}\tilde{K}^{+}.$
(42)

We now consider the inversion in the center of the above expression, i.e. the
inversion of $(\tilde{K}^{+}CC^{\prime}\tilde{K}^{+^{\prime}})$. This
inversion is given by;

$(\tilde{K}^{+}CC^{\prime}\tilde{K}^{+^{\prime}})^{-1}=\tilde{K}^{\prime}C^{-1^{\prime}}C^{-1}\tilde{K}.$

We will demonstrate this by showing that the product of the matrix inside the
inversion on the left hand side and the matrix on the right hand side is the
identity matrix:

$(\tilde{K}^{+}CC^{\prime}\tilde{K}^{+^{\prime}})(\tilde{K}^{\prime}C^{-1^{\prime}}C^{-1}\tilde{K})=\tilde{K}^{+}CC^{\prime}\tilde{K}^{+^{\prime}}\tilde{K}^{\prime}C^{-1^{\prime}}C^{-1}\tilde{K}.$

Using properties $1$, $2$ and $3$ we see that the above is equal to:

$\tilde{K}^{+}\tilde{N}CC^{\prime}C^{-1^{\prime}}C^{-1}\tilde{K}.$

Simplifying and applying property $1$ now yields that the above is equal the
below:

$\tilde{K}^{+}\tilde{N}\tilde{K}=\tilde{K}^{+}\tilde{K}=I,$

where $I$ is the identity matrix as required. Returning to (42) and applying
properties 1 and 2, we now have that:

$(\tilde{N}C)^{+}=C^{\prime}\tilde{N}C^{-1^{\prime}}C^{-1}\tilde{N}.$

By applying properties 2 and 3, the above can now be reduced to:

$(\tilde{N}C)^{+}=\tilde{N}C^{-1}\tilde{N}.$ (43)

We now note that, by property $3$, $\tilde{N}C=C\tilde{N}$. By multiplying the
left and right side of this identity by $C^{-1}$ it can be seen that:

$C^{-1}\tilde{N}=\tilde{N}C^{-1}.$

Multiplying the right hand side of the above by $\tilde{N}$ and applying
property 2 we obtain:

$C^{-1}\tilde{N}=\tilde{N}C^{-1}\tilde{N}.$

Substituting the above into (43) gives:

$(\tilde{N}C)^{+}=C^{-1}\tilde{N}.$

The result of Lemma 3 now follows by multiplying the right-hand side by
$\delta$ and applying property 4. ∎

###### Theorem 6.5

The below two update rules are equivalent, where $F(\theta^{f}_{s})$ is the
matrix defined in Section 2.1.2 and $\mathcal{I}(\theta^{f}_{s})$ is the
Fisher Information matrix of $\theta^{f}_{s}$.

$\theta^{f}_{s+1}=\theta^{f}_{s}+\alpha_{s}F(\theta_{s}^{f})^{-1}\frac{\partial
l(\theta_{s}^{f})}{\partial\theta},$
$\theta^{f}_{s+1}=\theta^{f}_{s}+\alpha_{s}\mathcal{I}(\theta_{s}^{f})^{+}\frac{\partial
l(\theta_{s}^{f})}{\partial\theta}.$

###### Proof

To prove Theorem 6.5, it suffices to prove the below equality:

$\mathcal{I}(\theta_{s}^{f})^{+}\frac{\partial
l(\theta_{s}^{f})}{\partial\theta}=F(\theta_{s}^{f})^{-1}\frac{\partial
l(\theta_{s}^{f})}{\partial\theta}.$

To prove the above, we employ Lemma 3. To achieve this, we proceed by
defining:

  * $\bullet$

$\tilde{N}$ to be the matrix:

$\begin{bmatrix}I_{p+1}&0&0&\ldots&0\\\ 0&N_{q_{1}}&0&\ldots&0\\\
0&0&N_{q_{2}}&\ldots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\ldots&N_{q_{r}}\\\ \end{bmatrix},$

  * $\bullet$

$\tilde{K}$ to be the matrix:

$\begin{bmatrix}I_{p+1}&0&0&\ldots&0\\\ 0&K_{q_{1},q_{1}}&0&\ldots&0\\\
0&0&K_{q_{2},q_{2}}&\ldots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\ldots&K_{q_{r},q_{r}}\\\ \end{bmatrix},$

  * $\bullet$

$\delta$ to be the column vector:

$\frac{\partial l(\theta^{f})}{\partial\theta^{f}},$

  * $\bullet$

$C$ to be the matrix $F(\theta^{f})$,

and claim that $F(\theta^{f})\tilde{N}=\mathcal{I}(\theta^{f})$. To prove this
equality, we first note the following two properties of the matrix $N_{k}$,
given by Magnus and Neudecker (1986).

  1. 1.

$N_{k_{1}}(A\otimes A)=N_{k_{1}}(A\otimes A)N_{k_{2}}=(A\otimes A)N_{k_{2}}$
for any arbitrary matrix $A$ and scalars $k_{1}$ and $k_{2}$, such that the
matrix multiplications are well defined.

  2. 2.

$N_{k}\text{vec}(A)=\text{vec}(A)$ for any arbitrary symmetric matrix $A$ of
appropriate dimension.

The first of the above properties, in combination with equations (15) and
(17), can be seen to imply that, for any $k_{1}$ and $k_{2}$ between $1$ and
$r$:

$\mathcal{I}^{f}_{\text{vec}(D_{k_{1}}),\text{vec}(D_{k_{2}})}=F_{\text{vec}(D_{k_{1}}),\text{vec}(D_{k_{2}})}N_{q_{k_{2}}}.$

In a similar manner, the second of the above properties, in combination with
equation (14), can be seen to imply that, for any $k$ between $1$ and $r$:

$\mathcal{I}^{f}_{\sigma^{2},\text{vec}(D_{k})}=F_{\sigma^{2},\text{vec}(D_{k})}N_{q_{k}}.$

Expanding the multiplication $F(\theta_{s}^{f})\tilde{N}$ block-wise
demonstrates that the two equations above are sufficient to prove the equality
$F(\theta_{s}^{f})\tilde{N}=\mathcal{I}(\theta^{f}_{s})$. The matrices
$\tilde{N},\tilde{K}$ and $C$ can be seen to satisfy properties 1-3 of Lemma 3
by noting standard results concerning the matrices $\mathcal{D}_{k}$,
$K_{m,n}$ and $N_{n}$, provided in Magnus and Neudecker (1986). Property 4 of
Lemma 3 can be seen to hold for the matrix $\tilde{N}$ and column vector
$\delta$ by noting the symmetry of $\frac{\partial l(\theta^{f}_{s})}{\partial
D_{k}}$ (see Theorem 6.1). By combining the previous arguments and applying
Lemma 3, the below is obtained:

$\mathcal{I}(\theta^{f})^{+}\delta=(F(\theta^{f})\tilde{N})^{+}\delta=F(\theta^{f})^{-1}\delta.$

The result of Theorem 6.5 now follows.∎

###### Theorem 6.6

For any fixed integer $k$ between $1$ and $r$, the below two update rules are
equivalent, where $F_{\text{vec}(D_{k,s})}$ is as defined in Section 2.1.2 and
$\mathcal{I}^{f}_{\text{vec}(D_{k,s})}$ is the Fisher Information matrix of
$\text{vec}(D_{k})$.

$\text{vec}(D_{k,s+1})=\text{vec}(D_{k,s})+\alpha_{s}F^{-1}_{\text{vec}(D_{k,s})}\frac{\partial
l(\theta^{f}_{s})}{\partial\text{vec}(D_{k,s})},$
$\text{vec}(D_{k,s+1})=\text{vec}(D_{k,s})+\alpha_{s}\big{(}\mathcal{I}^{f}_{\text{vec}(D_{k,s})}\big{)}^{+}\frac{\partial
l(\theta^{f}_{s})}{\partial\text{vec}(D_{k,s})}.$

###### Proof

The proof of the above proceeds by defining $\tilde{N}=N_{q_{k}}$,
$\tilde{K}=K_{q_{k},q_{k}}$, $C=F_{\text{vec}(D_{k})}$ and $\delta$ to be the
partial derivative appearing in the statement of the theorem. Following this,
the proof follows the same structure as that of Theorem 6.5. For brevity, the
argument will not be repeated in this section. ∎

###  6.4 Restricted Maximum Likelihood Estimation

In this appendix, we describe how the methods from Section 2.1 may be adapted
to use an alternative likelihood criteria: the criteria employed by Restricted
Maximum Likelihood (ReML) estimation. A well-documented issue for ML
estimation is that the variance estimates produced using ML are biased. First
described by Patterson and Thompson (1971), and proposed for the LMM in Laird
and Ware (1982), ReML offers an alternative to ML to address this issue by
maximizing the log-likelihood function of the residual vector, $e$, instead of
the response vector, $Y$. Neglecting constant terms, the Restricted Maximum
log-likelihood function, $l_{R}$, is given by:

$l_{R}(\theta^{h})=l(\theta^{h})-\frac{1}{2}\bigg{(}-p\log(\sigma^{2})+\log|X^{\prime}V^{-1}X|\bigg{)},$
(44)

where $l(\theta^{h})$ is given in (3). To derive the ReML-based FS algorithm,
akin to that described in Section 2.1.1, the following adjustments to the
score vectors given by (8) and (9) must be used.

$\frac{dl_{R}(\theta^{h})}{d\sigma^{2}}=\frac{dl(\theta^{h})}{d\sigma^{2}}+\frac{1}{2}p\sigma^{-2},$
$\displaystyle\frac{dl_{R}(\theta^{h})}{d\text{vech}(D_{k})}=\frac{dl(\theta^{h})}{d\text{vech}(D_{k})}+$
$\displaystyle\frac{1}{2}\mathcal{D}^{\prime}_{q_{k}}\text{vec}\bigg{(}\sum_{j=1}^{l_{k}}Z_{(k,j)}^{\prime}V^{-1}X(X^{\prime}V^{-1}X)^{-1}X^{\prime}V^{-1}Z_{(k,j)}\bigg{)}.$

The latter of the above results can be derived through trivial adjustment of
the proofs of Theorem 6.1 and Corollaries 1-3 given in Appendix 6.1. As the
derivative of $l_{R}(\theta^{h})$ with respect to $\beta$ is identical to that
of $l(\theta^{h})$, given in Section 2.1.1, we do not list it again above.

For a derivation of the ReML score vectors of $\beta$ and $\sigma^{2}$, we
direct the reader to the works of Demidenko (2013) where proofs may be found
for the single-factor LMM. As these proofs do not depend on the number of
factors in the model, they can be easily seen to also apply in the multi-
factor LMM setting without further adjustment.

Due to the asymptotic equivalence of the ML and ReML estimates, the parameter
estimates produced by ML and ReML have the same asymptotic covariance matrix.
Consequently, the Fisher Information matrix for the parameter estimates
produced by ReML, which can be shown to be equal to the inverse of the
asymptotic covariance matrix, is identical to that of the parameter estimates
produced by ML. As a result, the ReML-based FS algorithm utilizes the Fisher
Information matrix specified by (10)-(12) and requires no further derivation.

To summarize, the ReML-based FS algorithm for the multi-factor LMM is almost
identical to that given in Algorithm 1. The only adaptations required occur on
line 3 of Algorithm 1 where the score vectors must be substituted for their
ReML counterparts, provided above. Analogous adjustments can be made for the
algorithms presented in Sections 2.1.2-2.1.5.

###  6.5 Constraint Matrices

In this appendix, we provide further detail on the approach to constrained
optimization which was outlined in Section 2.5. We begin by providing examples
of how the notation described in Section 2.5 may be used in practice. For a
given $k$, we introduce the notation $\\{d_{i}\\}$ and $\\{\tilde{d}_{i}\\}$