# MM Algorithms for Variance Components Models

Hua Zhou
Department of Biostatistics
University of California, Los Angeles
<EMAIL_ADDRESS>
Liuyi Hu
Department of Statistics
North Carolina State University
<EMAIL_ADDRESS>
Jin Zhou
Division of Epidemiology
and Biostatistics
University of Arizona
<EMAIL_ADDRESS>
Kenneth Lange
Departments of Biomathematics,
Human Genetics, and Statistics
University of California, Los Angeles
<EMAIL_ADDRESS>
###### Abstract

Variance components estimation and mixed model analysis are central themes in
statistics with applications in numerous scientific disciplines. Despite the
best efforts of generations of statisticians and numerical analysts, maximum
likelihood estimation and restricted maximum likelihood estimation of variance
component models remain numerically challenging. Building on the minorization-
maximization (MM) principle, this paper presents a novel iterative algorithm
for variance components estimation. MM algorithm is trivial to implement and
competitive on large data problems. The algorithm readily extends to more
complicated problems such as linear mixed models, multivariate response models
possibly with missing data, maximum a posteriori estimation, penalized
estimation, and generalized estimating equations (GEE). We establish the
global convergence of the MM algorithm to a KKT point and demonstrate, both
numerically and theoretically, that it converges faster than the classical EM
algorithm when the number of variance components is greater than two and all
covariance matrices are positive definite.

Keywords: generalized estimating equations (GEE), global convergence, matrix
convexity, linear mixed model (LMM), maximum a posteriori (MAP) estimation,
maximum likelihood estimation (MLE), minorization-maximization (MM), missing
data, multivariate response, penalized estimation, restricted maximum
likelihood (REML), variance components model

##  1 Introduction

Variance components and linear mixed models are among the most potent tools in
a statistician’s toolbox. They are essential topics in graduate-level linear
model courses and the subject of many current papers and research monographs
(Rao and Kleffe,, 1988; Searle et al.,, 1992; Rao,, 1997; Khuri et al.,, 1998;
Demidenko,, 2013). Their applications in agriculture, biology, economics,
genetics, epidemiology, and medicine are too numerous to cover here in detail.
The recommended books (Verbeke and Molenberghs,, 2000; Weiss,, 2005;
Fitzmaurice et al.,, 2011) stress longitudinal data analysis.

Given an observed $n\times 1$ response vector ${\bm{y}}$ and $n\times p$
predictor matrix ${\bm{X}}$, the simplest variance components model postulates
that ${\bm{Y}}\sim
N({\bm{X}}\mbox{\boldmath$\beta$},\mbox{\boldmath$\Omega$})$, where

$\Omega$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\sigma_{i}^{2}{\bm{V}}_{i},$

and the ${\bm{V}}_{1},\ldots,{\bm{V}}_{m}$ are $m$ fixed positive semidefinite
matrices. The parameters of the model can be divided into mean effects
$(\beta_{1},\ldots,\beta_{p})$ and variance components
$(\sigma_{1}^{2},\ldots,\sigma_{m}^{2})$, summarized by vectors $\beta$ and
$\mbox{\boldmath$\sigma$}^{2}$. Throughout we assume $\Omega$ is positive
definite. The extension to singular $\Omega$ will not be pursued here.
Estimation revolves around the log-likelihood function

$\displaystyle L(\mbox{\boldmath$\beta$},\mbox{\boldmath$\sigma$}^{2})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\ln\det\mbox{\boldmath$\Omega$}-\frac{1}{2}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})^{T}\mbox{\boldmath$\Omega$}^{-1}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}).$
(1)

Among the commonly used methods for estimating variance components, maximum
likelihood estimation (MLE) (Hartley and Rao,, 1967) and restricted (or
residual) MLE (REML) (Harville,, 1977) are the most popular. REML first
projects ${\bm{y}}$ to the null space of ${\bm{X}}$ and then estimates
variance components based on the projected responses. If the columns of the
matrix ${\bm{B}}$ span the null space of ${\bm{X}}^{T}$, then REML estimates
the $\sigma_{i}^{2}$ by maximizing the log-likelihood of the redefined
response vector ${\bm{B}}^{T}{\bm{Y}}$, which is normally distributed with
mean ${\mathbf{0}}$ and covariance
${\bm{B}}^{T}\mbox{\boldmath$\Omega$}{\bm{B}}=\sum_{i=1}^{m}\sigma_{i}^{2}{\bm{B}}^{T}{\bm{V}}_{i}{\bm{B}}$.

There exists a large literature on iterative algorithms for finding MLE and
REML (Laird and Ware,, 1982; Lindstrom and Bates,, 1988, 1990; Harville and
Callanan,, 1990; Callanan and Harville,, 1991; Bates and Pinheiro,, 1998;
Schafer and Yucel,, 2002). Fitting variance component models remains a
challenge in models with a large sample size $n$ or a large number of variance
components $m$. Newton’s method (Lindstrom and Bates,, 1988) converges quickly
but is numerically unstable owing to the non-concavity of the log-likelihood.
Fisher’s scoring algorithm replaces the observed information matrix in
Newton’s method by the expected information matrix and yields an ascent
algorithm when safeguarded by step halving. However the calculation and
inversion of expected information matrices cost $O(mn^{3})+O(m^{3})$ flops for
unstructured ${\bm{V}}_{i}$ and quickly become impractical when either $n$ or
$m$ is large. The expectation-maximization (EM) algorithm initiated by
Dempster et al. is a third alternative (Dempster et al.,, 1977; Laird and
Ware,, 1982; Laird et al.,, 1987; Lindstrom and Bates,, 1988; Bates and
Pinheiro,, 1998). Compared to Newton’s method, the EM algorithm is easy to
implement and numerically stable, but painfully slow to converge. In practice,
a strategy of priming Newton’s method by a few EM steps leverages the
stability of EM and the faster convergence of second-order methods. Quasi-
Newton methods dispense with explicit calculation of the observed information
while achieving a superlinear rate of convergence.

In this paper we derive a minorization-maximization (MM) algorithm for finding
the MLE and REML estimates of variance components. We prove global convergence
of the MM algorithm to a Karush-Kuhn-Tucker (KKT) point and explain why MM
generally converges faster than EM for models with more than two variance
components. We also sketch extensions of the MM algorithm to the multivariate
response model with possibly missing responses, the linear mixed model (LMM),
maximum a posteriori (MAP) estimation, penalized estimation, and generalized
estimating equations (GEE). The numerical efficiency of the MM algorithm is
illustrated through simulated data sets and a genomic example with more than
200 variance components.

##  2 Preliminaries

### Background on MM algorithms

Throughout we reserve Greek letters for parameters and indicate the current
iteration number by a superscript $t$. The MM principle for maximizing an
objective function $f(\mbox{\boldmath$\theta$})$ involves minorizing the
objective function $f(\mbox{\boldmath$\theta$})$ by a surrogate function
$g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$ around the
current iterate $\mbox{\boldmath$\theta$}^{(t)}$ of a search (Lange et al.,,
2000). Minorization is defined by the two conditions

$\displaystyle f(\mbox{\boldmath$\theta$}^{(t)})$ $\displaystyle=$
$\displaystyle
g(\mbox{\boldmath$\theta$}^{(t)}\mid\mbox{\boldmath$\theta$}^{(t)})$ (2)
$\displaystyle f(\mbox{\boldmath$\theta$})$ $\displaystyle\geq$ $\displaystyle
g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})\>,\quad\quad\mbox{\boldmath$\theta$}\neq\mbox{\boldmath$\theta$}^{(t)}.$

In other words, the surface $\mbox{\boldmath$\theta$}\mapsto
g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$ lies below the
surface $\mbox{\boldmath$\theta$}\mapsto f(\mbox{\boldmath$\theta$})$ and is
tangent to it at the point
$\mbox{\boldmath$\theta$}=\mbox{\boldmath$\theta$}^{(t)}$. Construction of the
minorizing function
$g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$ constitutes
the first M of the MM algorithm. The second M of the algorithm maximizes the
surrogate $g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$
rather than $f(\mbox{\boldmath$\theta$})$. The point
$\mbox{\boldmath$\theta$}^{(t+1)}$ maximizing
$g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$ satisfies the
ascent property $f(\mbox{\boldmath$\theta$}^{(t+1)})\geq
f(\mbox{\boldmath$\theta$}^{(t)})$. This fact follows from the inequalities

$\displaystyle f(\mbox{\boldmath$\theta$}^{(t+1)})$ $\displaystyle\geq$
$\displaystyle
g(\mbox{\boldmath$\theta$}^{(t+1)}\mid\mbox{\boldmath$\theta$}^{(t)})\mathop{\;\>}\nolimits\geq\mathop{\;\>}\nolimits
g(\mbox{\boldmath$\theta$}^{(t)}\mid\mbox{\boldmath$\theta$}^{(t)})\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits
f(\mbox{\boldmath$\theta$}^{(t)}),$ (3)

reflecting the definition of $\mbox{\boldmath$\theta$}^{(t+1)}$ and the
tangency and domination conditions (2). The ascent property makes the MM
algorithm remarkably stable. The validity of the descent property depends only
on increasing $g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$,
not on maximizing
$g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$. With obvious
changes, the MM algorithm also applies to minimization rather than to
maximization. To minimize a function $f(\mbox{\boldmath$\theta$})$, we
majorize it by a surrogate function
$g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$ and minimize
$g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$ to produce the
next iterate $\mbox{\boldmath$\theta$}^{(t+1)}$. The acronym should not be
confused with the maximization-maximization algorithm in the variational Bayes
context (Jeon,, 2012).

Figure 1: Log-likelihood surface of a 2-variance component model and the
surrogate functions of EM and MM minorizing the objective function at point
($\sigma_{1}^{2(t)},\sigma_{2}^{2(t)})=(18.5,0.7)$.

The MM principle (De Leeuw,, 1994; Heiser,, 1995; Kiers,, 2002; Lange et al.,,
2000; Hunter and Lange,, 2004; Wu and Lange,, 2010) finds applications in
multidimensional scaling (Borg and Groenen,, 2005), ranking of sports teams
(Hunter,, 2004), variable selection (Hunter and Li,, 2005), optimal experiment
design (Yu,, 2010), multivariate statistics (Zhou and Lange,, 2010), geometric
programming (Lange and Zhou,, 2014), and many other areas (Lange,, 2010,
Chapter 12). The celebrated EM principle (Dempster et al.,, 1977) is a special
case of the MM principle. The $Q$ function produced in the E step of an EM
algorithm minorizes the log-likelihood up to an irrelevant constant. Thus,
both EM and MM share the same advantages: simplicity, stability, graceful
adaptation to constraints, and the tendency to avoid large matrix inversion.
The more general MM perspective frees algorithm derivation from the missing
data straitjacket and invites wider applications (Wu and Lange,, 2010). Figure
1 shows the minorization functions of EM and MM for a variance components
model with $m=2$ variance components.

EM and MM algorithms often exhibit slow convergence. Fortunately, this defect
can be remedied by off-the-shelf acceleration techniques for fixed point
iterations. The recently developed squared iterative method (SQUAREM)
(Varadhan and Roland,, 2008) and the quasi-Newton acceleration method (Zhou et
al.,, 2011) are particularly attractive, given their simplicity and minimal
memory and computational costs. Our numerical experiments feature the
unadorned MM algorithm and the quasi-Newton accelerated MM (aMM) algorithm
based on one secant pair. Using more secant pairs is likely to further improve
performance.

### Convex matrix functions

For symmetric matrices we write ${\bm{A}}\preceq{\bm{B}}$ when
${\bm{B}}-{\bm{A}}$ is positive semidefinite and ${\bm{A}}\prec{\bm{B}}$ if
${\bm{B}}-{\bm{A}}$ is positive definite. A matrix-valued function $f$ is said
to be (matrix) convex if

$\displaystyle f[\lambda{\bm{A}}+(1-\lambda){\bm{B}}]$ $\displaystyle\preceq$
$\displaystyle\lambda f({\bm{A}})+(1-\lambda)f({\bm{B}})$

for all ${\bm{A}}$, ${\bm{B}}$, and $\lambda\in[0,1]$. Our derivation of the
MM variance components algorithm hinges on the convexity of the two functions
mentioned in the next lemma.

######  Lemma 1.

(a) The matrix fractional function
$f({\bm{A}},{\bm{B}})={\bm{A}}^{T}{\bm{B}}^{-1}{\bm{A}}$ is jointly convex in
the $m\times n$ matrix ${\bm{A}}$ and the $m\times m$ positive definite matrix
${\bm{B}}$. (b) The log determinant function $f({\bm{B}})=\ln\det{\bm{B}}$ is
concave on the set of positive definite matrices.

###### Proof.

The matrix fractional function is matrix convex because its epigraph

$\displaystyle\\{({\bm{A}},{\bm{B}},{\bm{C}}):{\bm{B}}\succ{\mathbf{0}},{\bm{A}}^{T}{\bm{B}}^{-1}{\bm{A}}\preceq{\bm{C}}\\}$
$\displaystyle=$
$\displaystyle\left\\{({\bm{A}},{\bm{B}},{\bm{C}}):{\bm{B}}\succ{\mathbf{0}},\begin{pmatrix}{\bm{B}}&{\bm{A}}\\\
{\bm{A}}^{T}&{\bm{C}}\end{pmatrix}\succeq{\mathbf{0}}\right\\}$

is a convex set. Here ${\bm{C}}$ varies over the set of $n\times n$ positive
semidefinite matrices. The equivalence of these two epigraph representations
is proved in Boyd and Vandenberghe, (2004, A.5.5). For the concavity of the
log determinant, see Boyd and Vandenberghe, (2004, p74). ∎

##  3 Univariate response model

Our strategy for maximizing the log-likelihood (1) is to alternate updating
the mean parameters $\beta$ and the variance components
$\mbox{\boldmath$\sigma$}^{2}$. Updating $\beta$ given
$\mbox{\boldmath$\sigma$}^{2}$ is a standard general least squares problem
with solution

$\displaystyle\mbox{\boldmath$\beta$}^{(t+1)}$ $\displaystyle=$
$\displaystyle({\bm{X}}^{T}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{X}})^{-1}{\bm{X}}^{T}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{y}}.$
(4)

Updating $\mbox{\boldmath$\sigma$}^{2}$ given $\mbox{\boldmath$\beta$}^{(t)}$
depends on two minorizations. If we assume that all of the ${\bm{V}}_{i}$ are
positive definite, then the joint convexity of the map
$({\bm{X}},{\bm{Y}})\mapsto{\bm{X}}^{T}{\bm{Y}}^{-1}{\bm{X}}$ for positive
definite ${\bm{Y}}$ implies that

$\displaystyle\mbox{\boldmath$\Omega$}^{(t)}\mbox{\boldmath$\Omega$}^{-1}\mbox{\boldmath$\Omega$}^{(t)}$
$\displaystyle=$
$\displaystyle\left(\sum_{i=1}^{m}\sigma_{i}^{2(t)}{\bm{V}}_{i}\right)\left(\sum_{i=1}^{m}\sigma_{i}^{2}{\bm{V}}_{i}\right)^{-1}\left(\sum_{i=1}^{m}\sigma_{i}^{2(t)}{\bm{V}}_{i}\right)$
$\displaystyle=$
$\displaystyle\left(\sum_{i=1}^{m}\frac{\sigma_{i}^{2(t)}}{\sum_{j}\sigma_{j}^{2(t)}}\frac{\sum_{j}\sigma_{j}^{2(t)}}{\sigma_{i}^{2(t)}}\sigma_{i}^{2(t)}{\bm{V}}_{i}\right)\cdot\left(\sum_{i=1}^{m}\frac{\sigma_{i}^{2(t)}}{\sum_{j}\sigma_{j}^{2(t)}}\frac{\sum_{j}\sigma_{j}^{2(t)}}{\sigma_{i}^{2(t)}}\sigma_{i}^{2}{\bm{V}}_{i}\right)^{-1}$
$\displaystyle\cdot\left(\sum_{i=1}^{m}\frac{\sigma_{i}^{2(t)}}{\sum_{j}\sigma_{j}^{2(t)}}\frac{\sum_{j}\sigma_{j}^{2(t)}}{\sigma_{i}^{2(t)}}\sigma_{i}^{2(t)}{\bm{V}}_{i}\right)$
$\displaystyle\preceq$
$\displaystyle\sum_{i=1}^{m}\frac{\sigma_{i}^{2(t)}}{\sum_{j}\sigma_{j}^{2(t)}}\left(\frac{\sum_{j}\sigma_{j}^{2(t)}}{\sigma_{i}^{2(t)}}\sigma_{i}^{2(t)}{\bm{V}}_{i}\right)\left(\frac{\sum_{j}\sigma_{j}^{2(t)}}{\sigma_{i}^{2(t)}}\sigma_{i}^{2}{\bm{V}}_{i}\right)^{-1}\left(\frac{\sum_{j}\sigma_{j}^{2(t)}}{\sigma_{i}^{2(t)}}\sigma_{i}^{2(t)}{\bm{V}}_{i}\right)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\frac{\sigma_{i}^{4(t)}}{\sigma_{i}^{2}}{\bm{V}}_{i}{\bm{V}}_{i}^{-1}{\bm{V}}_{i}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\frac{\sigma_{i}^{4(t)}}{\sigma_{i}^{2}}{\bm{V}}_{i}.$

When one or more of the ${\bm{V}}_{i}$ are rank deficient, we replace each
${\bm{V}}_{i}$ by ${\bm{V}}_{i,\epsilon}={\bm{V}}_{i}+\epsilon{\bm{I}}$ for
$\epsilon>0$ small. Sending $\epsilon$ to 0 in the just proved majorization

$\displaystyle\mbox{\boldmath$\Omega$}_{\epsilon}^{(t)}\mbox{\boldmath$\Omega$}_{\epsilon}^{-1}\mbox{\boldmath$\Omega$}_{\epsilon}^{(t)}$
$\displaystyle\preceq$
$\displaystyle\sum_{i=1}^{m}\frac{\sigma_{i}^{4(t)}}{\sigma_{i}^{2}}{\bm{V}}_{i,\epsilon}$

gives the desired majorization

$\displaystyle\mbox{\boldmath$\Omega$}^{(t)}\mbox{\boldmath$\Omega$}^{-1}\mbox{\boldmath$\Omega$}^{(t)}$
$\displaystyle\preceq$
$\displaystyle\sum_{i=1}^{m}\frac{\sigma_{i}^{4(t)}}{\sigma_{i}^{2}}{\bm{V}}_{i}$

in the general case. Negating both sides leads to the minorization

$\displaystyle-({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})^{T}\mbox{\boldmath$\Omega$}^{-1}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})$
$\displaystyle\succeq$
$\displaystyle-({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})^{T}\mbox{\boldmath$\Omega$}^{-(t)}\left(\sum_{i=1}^{m}\frac{\sigma_{i}^{4(t)}}{\sigma_{i}^{2}}{\bm{V}}_{i}\right)\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})$
(5)

that effectively separates the variance components
$\sigma_{1}^{2},\ldots,\sigma_{m}^{2}$ in the quadratic term of the log-
likelihood (1).

The convexity of the function ${\bm{A}}\mapsto-\log\det{\bm{A}}$ is equivalent
to the supporting hyperplane minorization

$\displaystyle-\ln\det\mbox{\boldmath$\Omega$}$ $\displaystyle\geq$
$\displaystyle-\ln\det\mbox{\boldmath$\Omega$}^{(t)}-\mathrm{tr}[\mbox{\boldmath$\Omega$}^{-(t)}(\mbox{\boldmath$\Omega$}-\mbox{\boldmath$\Omega$}^{(t)})]$
(6)

that separates $\sigma_{1}^{2},\ldots,\sigma_{m}^{2}$ in the log determinant
term of the log-likelihood (1). Combination of the minorizations (5) and (6)
gives the overall minorization

$\displaystyle
g(\mbox{\boldmath$\sigma$}^{2}\mid\mbox{\boldmath$\sigma$}^{2(t)})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}\mbox{\boldmath$\Omega$})-\frac{1}{2}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})^{T}\mbox{\boldmath$\Omega$}^{-(t)}\left(\sum_{i=1}^{m}\frac{\sigma_{i}^{4(t)}}{\sigma_{i}^{2}}{\bm{V}}_{i}\right)\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})+c^{(t)}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\left[-\frac{\sigma_{i}^{2}}{2}\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i})-\frac{1}{2}\frac{\sigma_{i}^{4(t)}}{\sigma_{i}^{2}}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})^{T}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i}\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})\right]+c^{(t)},$

where $c^{(t)}$ is an irrelevant constant. Maximization of
$g(\mbox{\boldmath$\sigma$}^{2}\mid\mbox{\boldmath$\sigma$}^{2(t)})$ with
respect to $\sigma_{i}^{2}$ yields the lovely multiplicative update

$\displaystyle\sigma_{i}^{2(t+1)}$ $\displaystyle=$
$\displaystyle\sigma_{i}^{2(t)}\sqrt{\frac{({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})^{T}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i}\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})}{\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i})}},\quad
i=1,\ldots,m.$ (8)

To preserve the uniqueness and continuity of the algorithm map, we must take
$\sigma_{i}^{2(t+1)}=0$ whenever $\sigma_{i}^{2(t)}=0$. As a sanity check on
our derivation, consider the partial derivative

$\displaystyle\frac{\partial}{\partial\sigma_{i}^{2}}L(\mbox{\boldmath$\beta$},\mbox{\boldmath$\sigma$}^{2})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-1}{\bm{V}}_{i})+\frac{1}{2}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})^{T}\mbox{\boldmath$\Omega$}^{-1}{\bm{V}}_{i}\mbox{\boldmath$\Omega$}^{-1}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}).$
(9)

Given $\sigma_{i}^{2(t)}>0$, it is clear from the update formula (8) that
$\sigma_{i}^{2(t+1)}<\sigma_{i}^{2(t)}$ when
$\frac{\partial}{\partial\sigma_{i}^{2}}L<0$. Conversely
$\sigma_{i}^{2(t+1)}>\sigma_{i}^{2(t)}$ when
$\frac{\partial}{\partial\sigma_{i}^{2}}L>0$. Algorithm 1 summarizes the MM
algorithm for MLE of the univariate response model (1).

Input : ${\bm{y}}$, ${\bm{X}}$, ${\bm{V}}_{1},\ldots,{\bm{V}}_{m}$

Output : MLE $\hat{}\mbox{\boldmath$\beta$}$,
$\hat{\sigma}_{1}^{2},\ldots,\hat{\sigma}_{m}^{2}$

1 Initialize $\sigma_{i}^{(0)}>0$, $i=1,\ldots,m$

2 repeat

3
$\mbox{\boldmath$\Omega$}^{(t)}\leftarrow\sum_{i=1}^{m}\sigma_{i}^{2(t)}{\bm{V}}_{i}$

4
$\mbox{\boldmath$\beta$}^{(t)}\leftarrow\arg\min_{\mbox{\boldmath$\beta$}}\,({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})^{T}\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})$
$\sigma_{i}^{2(t+1)}\leftarrow\sigma_{i}^{2(t)}\sqrt{\frac{({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})^{T}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i}\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})}{\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i})}},\quad
i=1,\ldots,m$

5

6until  _objective value converges_

Algorithm 1 MM algorithm for MLE of the variance components of model (1).

The update formula (8) assumes that the numerator under the square root sign
is nonnegative and the denominator is positive. The numerator requirement is a
consequence of the positive semidefiniteness of ${\bm{V}}_{i}$. The
denominator requirement can be verified through the Hadamard (elementwise)
product representation
$\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i})={\mathbf{1}}^{T}(\mbox{\boldmath$\Omega$}^{-(t)}\odot{\bm{V}}_{i}){\mathbf{1}}$.
The following lemma of Schur, (1911) is crucial. We give a self-contained
probabilistic proof.

######  Lemma 2 (Schur).

The Hadamard product of a positive definite matrix with a positive
semidefinite matrix with positive diagonal entries is positive definite.

###### Proof.

Let ${\bm{X}}=(X_{1},\ldots,X_{n})^{T}$ be a random normal vector with mean
${\mathbf{0}}$ and positive definite covariance matrix ${\bm{A}}$. Let
${\bm{Y}}=(Y_{1},\ldots,Y_{n})^{T}$ be a random normal vector independent of
${\bm{X}}$ with mean ${\mathbf{0}}$ and positive semidefinite covariance
matrix ${\bm{B}}$ having positive diagonal entries. Then
${\bm{Z}}={\bm{X}}\odot{\bm{Y}}$ has covariances
$\mathrm{E}(Z_{i}Z_{j})=\mathrm{E}(X_{i}Y_{i}X_{j}Y_{j})=\mathrm{E}(X_{i}X_{j})\mathrm{E}(Y_{i}Y_{j})=a_{ij}b_{ij}$.
It follows that $\mathrm{Cov}({\bm{Z}})={\bm{A}}\odot{\bm{B}}$. To show
${\bm{A}}\odot{\bm{B}}$ is positive definite, suppose on the contrary that
${\bm{v}}^{T}({\bm{A}}\odot{\bm{B}}){\bm{v}}=\mathrm{Var}({\bm{v}}^{T}{\bm{Z}})=0$
for some ${\bm{v}}\neq{\mathbf{0}}$. Then

$\displaystyle 0$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\mathrm{Var}({\bm{v}}^{T}{\bm{Z}})=\mathrm{E}\Big{(}\sum_{i}v_{i}X_{i}Y_{i}\Big{)}^{2}=\mathrm{E}\Big{[}\Big{(}\sum_{i}v_{i}X_{i}Y_{i}\Big{)}^{2}\mid{\bm{Y}}\Big{]}=\mathrm{E}[({\bm{v}}\odot{\bm{Y}})^{T}{\bm{A}}({\bm{v}}\odot{\bm{Y}})]$

implies ${\bm{v}}\odot{\bm{Y}}={\mathbf{0}}$ with probability 1. Since
${\bm{v}}\neq{\mathbf{0}}$, $Y_{i}=0$ with probability 1 for some $i$. This
contradicts the assumption $b_{ii}=\mathrm{Var}(Y_{i})>0$ for all $i$. ∎

We can now obtain the following characterization of the MM iterates.

######  Proposition 1.

Assume ${\bm{V}}_{i}$ has strictly positive diagonal entries. Then
$\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i})>0$ for all $t$.
Furthermore if $\sigma_{i}^{2(0)}>0$ and
$\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})\notin\text{null}({\bm{V}}_{i})$
for all $t$, then $\sigma_{i}^{2(t)}>0$ for all $t$. When ${\bm{V}}_{i}$ is
positive definite, $\sigma_{i}^{2(t)}>0$ holds if and only if
${\bm{y}}\neq{\bm{X}}\mbox{\boldmath$\beta$}^{(t)}$.

###### Proof.

The first claim follows easily from Schur’s lemma. The second claim follows by
induction. The third claim follows from the observation that
$\text{null}({\bm{V}}_{i})=\\{{\mathbf{0}}\\}$. ∎

In most applications, ${\bm{V}}_{m}={\bm{I}}$. Proposition 1 guarantees that
if $\sigma_{m}^{2(0)}>0$ and the residual vector
${\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)}$ is nonzero, then
$\sigma_{m}^{2(t)}$ remains positive and thus $\mbox{\boldmath$\Omega$}^{(t)}$
remains positive definite throughout all iterations. This fact does not
prevent any of the sequences $\sigma_{i}^{2(t)}$ from converging to 0. In this
sense, the MM algorithm acts like an interior point method, approaching the
optimum from inside the feasible region.

### Univariate response: two variance components

Input : ${\bm{y}}$, ${\bm{X}}$, ${\bm{V}}_{1},{\bm{V}}_{2}$

Output : MLE $\hat{}\mbox{\boldmath$\beta$}$,
$\hat{\sigma}_{1}^{2},\hat{\sigma}_{2}^{2}$

1 Simultaneous congruence decomposition:
$({\bm{D}},{\bm{U}})\leftarrow({\bm{V}}_{1},{\bm{V}}_{2})$

2 Transform data: $\tilde{\bm{y}}\leftarrow{\bm{U}}^{T}{\bm{y}}$,
$\tilde{\bm{X}}\leftarrow{\bm{U}}^{T}{\bm{X}}$

3 Initialize $\sigma_{1}^{(0)},\sigma_{1}^{(0)}>0$

4 repeat

5 $w_{i}^{(t)}\leftarrow(\sigma_{1}^{2(t)}d_{i}+\sigma_{2}^{2(t)})^{-1},\quad
i=1,\ldots,n$

6
$\mbox{\boldmath$\beta$}^{(t)}\leftarrow\arg\min_{\mbox{\boldmath$\beta$}}\,\sum_{i=1}^{n}w_{i}^{(t)}(\tilde{y}_{i}-\tilde{\bm{x}}_{i}^{T}\mbox{\boldmath$\beta$})^{2}$
$\sigma_{1}^{2(t+1)}\leftarrow\sigma_{1}^{2(t)}\sqrt{\frac{(\tilde{\bm{y}}-\tilde{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})^{T}(\sigma_{1}^{2(t)}{\bm{D}}+\sigma_{2}^{2(t)}{\bm{I}})^{-1}{\bm{D}}(\sigma_{1}^{2(t)}{\bm{D}}+\sigma_{2}^{2(t)}{\bm{I}})^{-1}(\tilde{\bm{y}}-\tilde{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})}{\mathrm{tr}[(\sigma_{1}^{2(t)}{\bm{D}}+\sigma_{2}^{2(t)}{\bm{I}})^{-1}{\bm{D}}]}}$

7
$\sigma_{2}^{2(t+1)}\leftarrow\sigma_{2}^{2(t)}\sqrt{\frac{(\tilde{\bm{y}}-\tilde{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})^{T}(\sigma_{1}^{2(t)}{\bm{D}}+\sigma_{2}^{2(t)}{\bm{I}})^{-2}(\tilde{\bm{y}}-\tilde{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})}{\mathrm{tr}[(\sigma_{1}^{2(t)}{\bm{D}}+\sigma_{2}^{2(t)}{\bm{I}})^{-1}]}}$

8

9until  _objective value converges_

Algorithm 2 Simplified MM algorithm for MLE of model (1) with $m=2$ variance
components and
$\mbox{\boldmath$\Omega$}=\sigma_{1}^{2}{\bm{V}}_{1}+\sigma_{2}^{2}{\bm{V}}_{2}$.

The major computational cost of Algorithm 1 is inversion of the covariance
matrix $\mbox{\boldmath$\Omega$}^{(t)}$ at each iteration. Problem specific
structures such as block diagonal matrices or a diagonal matrix plus a low-
rank matrix are often exploited to speed up matrix inversion. The special case
of $m=2$ variance components deserves attention as repeated matrix inversion
can be avoided by invoking the simultaneous congruence decomposition for two
symmetric matrices, one of which is positive definite (Rao,, 1973; Horn and
Johnson,, 1985). This decomposition is also called the generalized eigenvalue
decomposition (Golub and Van Loan,, 1996; Boyd and Vandenberghe,, 2004). If
one assumes
$\mbox{\boldmath$\Omega$}=\sigma_{1}^{2}{\bm{V}}_{1}+\sigma_{2}^{2}{\bm{V}}_{2}$
and lets $({\bm{V}}_{1},{\bm{V}}_{2})\mapsto({\bm{D}},{\bm{U}})$ be the
decomposition with ${\bm{U}}$ nonsingular,
${\bm{U}}^{T}{\bm{V}}_{1}{\bm{U}}={\bm{D}}$ diagonal, and
${\bm{U}}^{T}{\bm{V}}_{2}{\bm{U}}={\bm{I}}$, then

$\displaystyle\mbox{\boldmath$\Omega$}^{(t)}$ $\displaystyle=$
$\displaystyle{\bm{U}}^{-T}(\sigma_{1}^{2(t)}{\bm{D}}+\sigma_{2}^{2(t)}{\bm{I}}_{n}){\bm{U}}^{-1}$
$\displaystyle\mbox{\boldmath$\Omega$}^{-(t)}$ $\displaystyle=$
$\displaystyle{\bm{U}}(\sigma_{1}^{2(t)}{\bm{D}}+\sigma_{2}^{2(t)}{\bm{I}}_{n})^{-1}{\bm{U}}^{T}$
$\displaystyle\det(\mbox{\boldmath$\Omega$}^{(t)})$ $\displaystyle=$
$\displaystyle\det(\sigma_{1}^{2(t)}{\bm{D}}+\sigma_{2}^{2(t)}{\bm{I}}_{n})\det({\bm{U}}^{-T}{\bm{U}}^{-1})$
$\displaystyle=$
$\displaystyle\det(\sigma_{1}^{2(t)}{\bm{D}}+\sigma_{2}^{2(t)}{\bm{I}}_{n})\det({\bm{V}}_{2}).$

With the revised responses $\tilde{\bm{y}}={\bm{U}}^{T}{\bm{y}}$ and the
revised predictor matrix $\tilde{\bm{X}}={\bm{U}}^{T}{\bm{X}}$, the update (8)
requires only vector operations and costs $O(n)$ flops. Updating the fixed
effects is a weighted least squares problem with the transformed data
$(\tilde{\bm{y}},\tilde{\bm{X}})$ and observation weights
$w_{i}^{(t)}=(\sigma_{1}^{2(t)}d_{i}+\sigma_{2}^{2(t)})^{-1}$. Algorithm 2
summarizes the simplified MM algorithm for two variance components.

### Numerical experiments

This section compares the numerical performance of MM, quasi-Newton
accelerated MM, EM, and Fisher scoring on simulated data from a two-way ANOVA
random effects model and a genetic model. For ease of comparison, all
algorithm runs start from $\mbox{\boldmath$\sigma$}^{2(0)}={\mathbf{1}}$ and
terminate when the relative change $(L^{(t+1)}-L^{(t)})/(|L^{(t)}|+1)$ in the
log-likelihood is less than $10^{-6}$.

Two-way ANOVA: We simulated data from a two-way ANOVA random effects model

$\displaystyle
y_{ijk}=\mu+\alpha_{i}+\beta_{j}+(\alpha\beta)_{ij}+\epsilon_{ijk},$

where $\alpha_{i}\sim N(0,\sigma_{1}^{2})$, $\beta_{j}\sim
N(0,\sigma_{2}^{2})$, $(\alpha\beta)_{ij}\sim N(0,\sigma_{3}^{2})$, and
$\epsilon_{ijk}\sim N(0,\sigma_{e}^{2})$ are jointly independent. This
corresponds to $m=4$ variance components. In the simulation, we set
$\sigma_{2}^{2}=\sigma_{3}^{2}=\sigma_{e}^{2}$ and varied the ratio
$\sigma_{1}^{2}/\sigma_{e}^{2}$; the numbers of levels $a$ and $b$ in factor 1
and factor 2, respectively; and the number of observations $c$ in each
combination of factor levels. For each simulation scenario, we simulated 50
replicates. The sample size was $n=abc$ for each replicate.

Tables 1 and 2 show the average number of iterations and the average runtimes
when there are $a=b=5$ levels of each factor. Based on these results and
further results not shown for other combinations of $a$ and $b$, we draw the
following conclusions. Fisher scoring takes the fewest iterations. The MM
algorithm always takes fewer iterations than the EM algorithm. Accelerated MM
further improves the convergence rate of MM. The faster rate of convergence of
Fisher scoring is outweighed by the extra cost of evaluating and inverting the
covariance matrix. When the sample size $n=abc$ is large, Fisher scoring takes
much longer than either EM or MM.

Genetic model: We simulated a quantitative trait ${\bm{y}}$ from a genetic
model with two variance components and covariance matrix
$\mbox{\boldmath$\Omega$}=\sigma_{a}^{2}\widehat{}\mbox{\boldmath$\Phi$}+\sigma_{e}^{2}{\bm{I}}$,
where $\widehat{}\mbox{\boldmath$\Phi$}$ is a full-rank empirical kinship
matrix estimated from the genome-wide measurements of 212 individuals using
Option 29 of the Mendel software (Lange et al.,, 2013, 2005). In this example,
Fisher scoring excels at smaller $\sigma_{a}^{2}/\sigma_{e}^{2}$ ratios, while
accelerated MM is fastest at larger $\sigma_{a}^{2}/\sigma_{e}^{2}$ ratios.

In summary, the MM algorithm appears competitive even in small-scale examples.
Modern applications often involve a large number of variance components. In
this setting, the EM algorithm suffers from slow convergence and Fisher
scoring from an extremely high cost per iteration. Our genomic example in
Section 7 reinforces this point.

$\sigma^{2}_{1}/\sigma^{2}_{e}$ | Method |  $c=$ # observations per combination
---|---|---
|  | 2 | 8 | 20 | 50
0 | MM | 34.52(15.79) | 25.90(8.69) | 18.62(7.22) | 15.48(5.34)
| aMM | 16.68(5.69) | 13.48(3.47) | 11.76(3.12) | 10.88(2.43)
| EM | 123.70(63.72) | 61.58(31.36) | 38.44(18.58) | 25.66(10.31)
| FS | 6.10(1.09) | 6.74(0.99) | 6.68(0.79) | 6.36(0.72)
0.05 | MM | 27.78(13.05) | 22.82(8.96) | 19.82(6.55) | 15.48(3.97)
| aMM | 14.80(4.33) | 12.32(3.27) | 12.08(2.62) | 11.20(2.52)
| EM | 108.04(62.58) | 58.42(33.67) | 43.52(19.48) | 27.62(12.47)
| FS | 6.20(1.29) | 6.72(1.25) | 6.62(0.73) | 6.60(1.07)
0.1 | MM | 31.26(14.90) | 23.38(9.21) | 16.84(6.72) | 14.88(4.56)
| aMM | 15.96(5.65) | 12.72(3.59) | 10.36(2.51) | 10.80(2.46)
| EM | 112.12(72.70) | 62.26(28.87) | 34.86(22.61) | 24.10(11.96)
| FS | 6.10(1.25) | 6.90(0.79) | 6.48(0.86) | 6.52(0.86)
1 | MM | 29.72(15.85) | 22.72(10.86) | 17.78(8.18) | 13.94(4.73)
| aMM | 15.24(5.60) | 12.40(4.12) | 10.72(2.70) | 10.24(2.01)
| EM | 85.86(63.85) | 41.50(30.46) | 28.40(20.02) | 21.36(13.86)
| FS | 5.96(1.19) | 6.90(0.91) | 6.36(1.05) | 6.44(0.93)
10 | MM | 16.46(9.74) | 13.28(7.75) | 12.80(6.41) | 10.74(3.67)
| aMM | 11.60(3.70) | 9.36(2.78) | 9.04(3.00) | 8.68(2.54)
| EM | 24.50(32.87) | 16.18(23.06) | 15.10(16.55) | 12.36(11.13)
| FS | 6.98(0.80) | 6.96(0.70) | 6.74(0.83) | 6.76(0.52)
20 | MM | 17.34(10.70) | 14.20(6.79) | 11.58(4.46) | 10.16(4.26)
| aMM | 12.12(5.96) | 9.92(2.68) | 8.92(2.07) | 8.48(2.00)
| EM | 31.08(42.11) | 20.50(24.55) | 10.84(10.86) | 8.98(8.94)
| FS | 7.18(0.98) | 7.02(0.82) | 6.90(0.74) | 6.78(0.79)
Table 1: Average iterations until convergence for MM, quasi-Newton accelerated MM (aMM), EM, and Fisher scoring (FS) for fitting a two-way ANOVA model with $a=b=5$ levels of both factors. Standard errors are given in parentheses. $\sigma^{2}_{1}/\sigma^{2}_{e}$ | Method |  $c$ = # observations per combination
---|---|---
|  | 2 | 8 | 20 | 50
0 | MM | 30.50(81.58) | 132.50(48.51) | 739.80(272.66) | 4004.10(1317.90)
| aMM | 16.76(36.57) | 92.06(34.42) | 638.76(231.86) | 3691.90(873.70)
| EM | 70.76(43.58) | 376.82(184.27) | 1912.30(918.56) | 8276.98(3269.26)
| FS | 17.52(27.99) | 241.06(51.79) | 4039.73(6955.91) | 43315.42(76563.64)
0.05 | MM | 16.79(11.73) | 117.33(50.33) | 867.41(296.24) | 4083.60(1016.71)
| aMM | 14.23(14.21) | 80.93(31.83) | 692.54(196.19) | 3841.10(948.55)
| EM | 66.73(44.87) | 376.78(206.19) | 2291.69(1035.80) | 9054.02(3989.52)
| FS | 13.33(18.33) | 253.10(61.72) | 3198.17(379.00) | 34057.87(7132.36)
0.1 | MM | 17.01(8.97) | 122.03(55.77) | 733.37(329.08) | 3992.52(1166.67)
| aMM | 12.39(8.86) | 88.34(33.87) | 593.27(174.93) | 3745.50(922.99)
| EM | 76.65(53.37) | 389.45(179.41) | 1814.91(1152.64) | 7951.40(3810.34)
| FS | 10.24(8.57) | 257.63(45.53) | 3140.15(481.08) | 33490.67(4533.51)
1 | MM | 16.50(13.08) | 112.94(52.73) | 736.32(322.26) | 3746.87(1221.92)
| aMM | 9.86(4.10) | 80.98(39.49) | 585.98(158.65) | 3536.90(754.48)
| EM | 56.93(48.16) | 267.86(194.49) | 1465.45(986.31) | 7079.60(4430.10)
| FS | 15.75(17.26) | 262.49(44.28) | 3003.68(481.92) | 33215.82(4801.38)
10 | MM | 10.80(11.16) | 70.94(47.94) | 545.71(256.97) | 2316.96(1022.51)
| aMM | 8.64(4.17) | 62.50(26.13) | 483.63(183.47) | 2317.43(1061.57)
| EM | 21.51(31.61) | 113.76(158.82) | 803.36(816.05) | 3256.65(2624.29)
| FS | 12.32(9.81) | 261.85(37.84) | 3190.52(394.75) | 26163.81(8451.96)
20 | MM | 8.83(5.05) | 104.94(54.66) | 552.13(190.42) | 1706.71(680.84)
| aMM | 9.57(9.80) | 92.94(35.84) | 524.70(137.22) | 1750.99(489.96)
| EM | 23.13(31.17) | 175.12(198.18) | 642.39(576.82) | 2007.86(1901.66)
| FS | 12.71(11.90) | 340.81(48.29) | 3543.18(464.36) | 18796.59(2445.74)
Table 2: Average run times ($\times 10^{-3}$ seconds) of MM, quasi-Newton accelerated MM (aMM), EM, and Fisher scoring (FS) for fitting a two-way ANOVA model with $a=b=5$ levels of both factors. Standard errors are given in parentheses. $\sigma^{2}_{a}/\sigma^{2}_{e}$ | Method | Iteration | Runtime ($10^{-3}$ sec) | Objective
---|---|---|---|---
0 | MM | 88.10(29.01) | 778.24(305.37) | -374.35(9.82)
| aMM | 23.65(5.74) | 293.16(146.23) | -374.34(9.82)
| EM | 231.93(123.39) | 3509.02(1851.11) | -374.41(9.83)
| FS | 5.05(1.24) | 137.76(65.74) | -374.36(9.83)
0.05 | MM | 84.97(31.18) | 710.56(260.24) | -377.19(10.85)
| aMM | 23.05(5.45) | 272.04(67.01) | -377.18(10.85)
| EM | 220.57(124.70) | 3292.87(1865.91) | -377.25(10.85)
| FS | 5.08(1.21) | 136.47(33.18) | -377.21(10.83)
0.1 | MM | 82.45(34.39) | 673.96(268.23) | -379.62(10.54)
| aMM | 22.55(6.01) | 269.55(86.69) | -379.61(10.54)
| EM | 199.70(113.47) | 2917.71(1607.33) | -379.68(10.54)
| FS | 4.97(1.03) | 129.71(40.51) | -379.62(10.54)
1 | MM | 31.00(15.59) | 160.21(80.45) | -409.66(11.26)
| aMM | 12.55(5.38) | 90.21(43.54) | -409.66(11.26)
| EM | 51.10(28.70) | 550.55(321.89) | -409.67(11.26)
| FS | 4.60(0.59) | 80.28(25.56) | -409.66(11.26)
10 | MM | 72.67(39.23) | 374.80(209.31) | -532.57(9.11)
| aMM | 20.15(10.18) | 146.25(81.06) | -531.24(9.28)
| EM | 294.20(717.05) | 3079.82(7520.30) | -532.71(9.11)
| FS | 10.18(4.92) | 168.63(80.34) | -532.08(9.21)
20 | MM | 78.35(34.32) | 425.40(188.08) | -591.36(7.05)
| aMM | 14.80(6.53) | 117.14(71.14) | -589.13(7.15)
| EM | 362.07(764.60) | 4144.92(8862.65) | -591.62(6.82)
| FS | 10.93(4.75) | 181.48(83.96) | -590.68(7.08)
Table 3: Average performance of MM, quasi-Newton accelerated MM (aMM), EM, and
Fisher scoring (FS) for fitting a genetic model. Standard errors are given in
parentheses.

##  4 Global convergence of the MM algorithm

The KKT necessary conditions for a local maximum
$\mbox{\boldmath$\sigma$}^{2}=(\sigma_{1}^{2},\ldots,\sigma_{m}^{2})$ of the
log-likelihood (1) require each component of the score vector to satisfy

$\displaystyle\frac{\partial}{\partial\sigma_{i}^{2}}L(\mbox{\boldmath$\sigma$}^{2})$
$\displaystyle\in$ $\displaystyle\begin{cases}\\{0\\}&\sigma_{i}^{2}>0\\\
(-\infty,0]&\sigma_{i}^{2}=0\,.\end{cases}$

In this section we establish the global convergence of Algorithm 1 to a KKT
point. To reduce the notational burden, we assume that ${\bm{X}}$ is null and
omit estimation of fixed effects $\beta$. The analysis easily extends to the
MLE case. Our convergence analysis relies on characterizing the properties of
the objective function $L(\mbox{\boldmath$\sigma$}^{2})$ and the MM
algorithmic mapping $\mbox{\boldmath$\sigma$}^{2}\mapsto
M(\mbox{\boldmath$\sigma$}^{2})$ defined by equation (8). Special attention
must be paid to the boundary values $\sigma_{i}^{2}=0$. We prove convergences
for two cases, which cover most applications.

######  Assumption 1.

All ${\bm{V}}_{i}$ are positive definite.

######  Assumption 2.

${\bm{V}}_{1}$ is positive definite, each ${\bm{V}}_{i}$ is nontrivial, ${\cal
H}=\text{span}\\{{\bm{V}}_{2},\ldots,{\bm{V}}_{m}\\}$ has dimension $q<n$, and
${\bm{y}}\notin{\cal H}$.

The genetic model in Section 3 satisfies Assumption 1, while the two-way ANOVA
model satisfies Assumption 2. The key condition
${\bm{y}}\notin\text{span}\\{{\bm{V}}_{2},\ldots,{\bm{V}}_{m}\\}$ in the
second case is critical for the existence of an MLE or an REML (Demidenko and
Massam,, 1999; Grządziel and Michalski,, 2014). We will derive a sequence of
lemmas en route to the global convergence result declared in Theorem 1.

######  Lemma 3.

Under Assumption 1 or 2, the log-likelihood function (1) is coercive in the
sense that the super-level set
$S_{c}=\\{\mbox{\boldmath$\sigma$}^{2}\geq{\mathbf{0}}:L(\mbox{\boldmath$\sigma$}^{2})\geq
c\\}$ is compact for every $c$.

###### Proof.

Let us first prove the assertion when all of the covariance matrices
${\bm{V}}_{i}$ are positive definite. If we set
$r=\|\mbox{\boldmath$\sigma$}^{2}\|_{1}$ and $\alpha_{i}=r^{-1}\sigma_{i}^{2}$
for each $i$, then the log-likelihood satisfies

$\displaystyle L(\mbox{\boldmath$\sigma$}^{2})$ $\displaystyle=$
$\displaystyle-\frac{n}{2}\ln
r-\frac{1}{2}\ln\det\Big{(}\sum_{i=1}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}-\frac{1}{2r}{\bm{y}}^{T}\Big{(}\sum_{i=1}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}^{-1}{\bm{y}}.$

The functions $\ln\det\Big{(}\sum_{i=1}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}$ and
${\bm{y}}^{T}\Big{(}\sum_{i=1}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}^{-1}{\bm{y}}$
of $\alpha$ are defined and continuous on the unit simplex and hence bounded
there. The dominant term $-\frac{n}{2}\ln r$ of the loglikelihood tends to
$-\infty$ as $r$ tends to $\infty$.

To prove the assertion under Assumption 2, consider first the case
${\bm{V}}_{1}={\bm{I}}_{n}$. Setting
$\alpha_{i}=\sigma_{i}^{2}/\sigma_{1}^{2}$ for $i=2,\ldots,m$ reduces the
loglikelihood to

$\displaystyle L(\sigma_{1}^{2},\mbox{\boldmath$\alpha$})$ $\displaystyle=$
$\displaystyle-\frac{n}{2}\ln\sigma_{1}^{2}-\frac{1}{2}\ln\det\Big{(}{\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}-\frac{1}{2\sigma_{1}^{2}}{\bm{y}}^{T}\Big{(}{\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}^{-1}\\!\\!{\bm{y}}.$
(11)

The middle term on the right satisfies

$\displaystyle-\frac{1}{2}\ln\det\Big{(}{\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}$
$\displaystyle\leq$ $\displaystyle 0$

because
$\det\,({\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i})\,\geq\,\det{\bm{I}}_{n}=1.$
Now let ${\bm{U}}=({\bm{U}}_{q},{\bm{U}}_{n-q})$ be an $n\times n$ orthogonal
matrix whose left columns ${\bm{U}}_{q}$ span ${\cal H}$ and whose right
columns ${\bm{U}}_{n-q}$ span ${\cal H}^{\perp}$. The identity

$\displaystyle{\bm{U}}^{T}\Big{(}{\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}{\bm{U}}$
$\displaystyle=$
$\displaystyle\begin{pmatrix}{\bm{I}}_{q}+\sum_{i=2}^{m}\alpha_{i}{\bm{U}}_{q}^{T}{\bm{V}}_{i}{\bm{U}}_{q}&{\mathbf{0}}\\\
{\mathbf{0}}&{\bm{I}}_{n-q}\end{pmatrix}$

follows from the orthogonality relations
${\bm{U}}_{n-q}^{T}{\bm{V}}_{i}\,=\,{\bm{U}}_{n-q}^{T}{\bm{U}}_{q}\,=\,{\mathbf{0}}_{(n-q)\times
n}$. This in turn implies

$\displaystyle\Big{(}{\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}^{-1}$
$\displaystyle=$
$\displaystyle{\bm{U}}\begin{pmatrix}({\bm{I}}_{q}+\sum_{i=2}^{m}\alpha_{i}{\bm{U}}_{q}^{T}{\bm{V}}_{i}{\bm{U}}_{q})^{-1}&{\mathbf{0}}\\\
{\mathbf{0}}&{\bm{I}}_{n-q}\end{pmatrix}{\bm{U}}^{T}$ $\displaystyle\succeq$
$\displaystyle{\bm{U}}\begin{pmatrix}{\mathbf{0}}&{\mathbf{0}}\\\
{\mathbf{0}}&{\bm{I}}_{n-q}\end{pmatrix}{\bm{U}}^{T}$ $\displaystyle=$
$\displaystyle{\bm{U}}_{n-q}{\bm{U}}_{n-q}^{T}.$

Therefore the quadratic term in equation (11) is bounded below by the positive
constant

$\displaystyle{\bm{y}}^{T}\Big{(}{\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}^{-1}\\!\\!{\bm{y}}$
$\displaystyle\geq$
$\displaystyle{\bm{y}}^{T}{\bm{U}}_{n-q}{\bm{U}}_{n-q}^{T}{\bm{y}}\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits\|\mathbb{P}_{{\cal
H}^{\perp}}{\bm{y}}\|^{2}\mathop{\;\>}\nolimits>\mathop{\;\>}\nolimits 0.$

Here the assumption ${\bm{y}}\notin{\cal H}$ guarantees the projection
property $\mathbb{P}_{{\cal H}^{\perp}}{\bm{y}}\neq{\mathbf{0}}$.

Next we show that the loglikelihood tends to $-\infty$ when $\sigma_{1}^{2}$
tends to $0$ or $\infty$ or when $\|\mbox{\boldmath$\alpha$}\|_{2}$ tends to
$\infty$. The second of the two inequalities

$\displaystyle L(\sigma_{0}^{2},\mbox{\boldmath$\alpha$})$ $\displaystyle\leq$
$\displaystyle-\frac{n}{2}\ln\sigma_{1}^{2}-\frac{1}{2}\ln\det\Big{(}{\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}-\frac{1}{2\sigma_{1}^{2}}\|\mathbb{P}_{{\cal
H}^{\perp}}{\bm{y}}\|^{2}$ $\displaystyle\leq$
$\displaystyle-\frac{n}{2}\ln\sigma_{1}^{2}-\frac{1}{2\sigma_{1}^{2}}\|\mathbb{P}_{{\cal
H}^{\perp}}{\bm{y}}\|^{2}$

renders the claim about $\sigma_{1}^{2}$ obvious. To prove the claim about
$\alpha$, we make the worst case choice $\sigma_{i}^{2}=\|\mathbb{P}_{{\cal
H}^{\perp}}{\bm{y}}\|^{2}$ in the first inequality. It follows that

$\displaystyle L(\sigma_{0}^{2},\mbox{\boldmath$\alpha$})$ $\displaystyle\leq$
$\displaystyle-\frac{1}{2}\ln\det\Big{(}{\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}-\frac{n}{2}\ln\|\mathbb{P}_{{\cal
H}^{\perp}}{\bm{y}}\|^{2}-\frac{n}{2}.$

If $\alpha_{j}$ tends to $\infty$, then the inequality

$\displaystyle-\frac{1}{2}\ln\det\Big{(}{\bm{I}}_{n}+\sum_{i=2}^{m}\alpha_{i}{\bm{V}}_{i}\Big{)}$
$\displaystyle\leq$
$\displaystyle-\frac{1}{2}\ln\det\Big{(}{\bm{I}}_{n}+\alpha_{j}{\bm{V}}_{j}\Big{)}\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits-\frac{1}{2}\sum_{k=1}^{n}\ln(1+\alpha_{j}\lambda_{jk})$

holds, where the $\lambda_{jk}$ are the eigenvalues of ${\bm{V}}_{j}$. At
least one of these eigenvalues is positive because ${\bm{V}}_{j}$ is
nontrivial. It follows that $L(\sigma_{0}^{2},\mbox{\boldmath$\alpha$})$ tends
to $-\infty$ in this case as well.

For the general case where ${\bm{V}}_{1}$ is non-singular but not necessarily
${\bm{I}}_{n}$, let ${\bm{V}}^{1/2}_{1}$ be the symmetric square root of
${\bm{V}}_{1}$ and write

$\displaystyle{\bm{V}}_{1}+\sum_{i=2}^{m}\sigma_{i}^{2}{\bm{V}}_{i}={\bm{V}}^{1/2}_{1}\left({\bm{I}}+\sum_{i=2}^{m}\sigma_{i}^{2}{\bm{V}}_{1}^{-1/2}{\bm{V}}_{i}{\bm{V}}_{1}^{-1/2}\right){\bm{V}}^{1/2}_{1}.$

The above arguments still apply since each
${\bm{V}}_{1}^{-1/2}{\bm{V}}_{i}{\bm{V}}_{1}^{-1/2}$ is nontrivial and
${\bm{y}}$ belongs to the
$\text{span}\\{{\bm{V}}_{2},\ldots,{\bm{V}}_{m}\\}={\bm{S}}$ if and only if
${\bm{V}}_{1}^{-1/2}{\bm{y}}$ belongs to
${\bm{V}}^{-1/2}_{1}{\bm{S}}{\bm{V}}_{1}^{-1/2}$. ∎

######  Lemma 4.

The iterates possess the ascent property
$L(M(\mbox{\boldmath$\sigma$}^{2(t)}))\geq
L(\mbox{\boldmath$\sigma$}^{2(t)})$. Furthermore, when
$L(M(\mbox{\boldmath$\sigma$}_{*}^{2}))=L(\mbox{\boldmath$\sigma$}_{*}^{2})$,
$\mbox{\boldmath$\sigma$}_{*}^{2}$ fulfills the fixed point condition
$M(\mbox{\boldmath$\sigma$}_{*}^{2})=\mbox{\boldmath$\sigma$}_{*}^{2}$, and
each component satisfies either (i) $\sigma_{*i}^{2}=0$ or (ii)
$\sigma_{*i}^{2}>0$ and
$\frac{\partial}{\partial\sigma_{i}^{2}}L(\mbox{\boldmath$\sigma$}_{*}^{2})=0$.

###### Proof.

The ascent property is built into any MM algorithm. Suppose
$L(M(\mbox{\boldmath$\sigma$}_{*}^{2}))=L(\mbox{\boldmath$\sigma$}_{*}^{2})$
at a point $\mbox{\boldmath$\sigma$}_{*}^{2}\in\mathbb{R}_{+}^{m}$. Then
equality must hold in the string of inequalities (3). It follows that

$\displaystyle
g(M(\mbox{\boldmath$\sigma$}_{*}^{2})\mid\mbox{\boldmath$\sigma$}_{*}^{2})$
$\displaystyle=$ $\displaystyle
g(\mbox{\boldmath$\sigma$}_{*}^{2}\mid\mbox{\boldmath$\sigma$}_{*}^{2})$

and hence that
$M(\mbox{\boldmath$\sigma$}_{*}^{2})=\mbox{\boldmath$\sigma$}_{*}^{2}$. If
$\sigma_{*i}^{2}>0$, the stationarity condition

$\displaystyle\frac{\partial}{\partial\sigma_{i}^{2}}L(\mbox{\boldmath$\sigma$}_{*}^{2})$
$\displaystyle=$
$\displaystyle\frac{\partial}{\partial\sigma_{i}^{2}}g(\mbox{\boldmath$\sigma$}_{*}^{2}\mid\mbox{\boldmath$\sigma$}_{*}^{2})\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits
0$

applies. The equivalence of the two displayed partial derivatives is a
consequence of the fact that the difference
$f(\mbox{\boldmath$\sigma$}^{2})-g(\mbox{\boldmath$\sigma$}^{2}\mid\mbox{\boldmath$\sigma$}_{*}^{2})$
achieves its minimum of 0 at
$\mbox{\boldmath$\sigma$}^{2}=\mbox{\boldmath$\sigma$}_{*}^{2}$. ∎

######  Lemma 5.

The distance between successive iterates
$\|\mbox{\boldmath$\sigma$}^{2(t+1)}-\mbox{\boldmath$\sigma$}^{2(t)}\|_{2}$
converges to 0.

###### Proof.

Suppose on the contrary that
$\|\mbox{\boldmath$\sigma$}^{2(t+1)}-\mbox{\boldmath$\sigma$}^{2(t)}\|_{2}$
does not converge to 0. Then one can extract a subsequence $\\{t_{k}\\}_{k\geq
1}$ such that

$\displaystyle\|\mbox{\boldmath$\sigma$}^{2(t_{k}+1)}-\mbox{\boldmath$\sigma$}^{2(t_{k})}\|_{2}\geq\epsilon>0$
(12)

for all $k$. Let $C_{0}$ be the compact super-level set
$\\{\mbox{\boldmath$\sigma$}^{2}:L(\mbox{\boldmath$\sigma$}^{2})\geq
L(\mbox{\boldmath$\sigma$}^{2(0)})\\}$. Since the sequence
$\\{\mbox{\boldmath$\sigma$}^{2(t_{k})}\\}_{k\geq 1}$ is confined to $C_{0}$,
one can pass to a subsequence if necessary and assume that
$\mbox{\boldmath$\sigma$}^{2(t_{k})}$ converges to a limit
$\mbox{\boldmath$\sigma$}_{*}^{2}$ and that
$\mbox{\boldmath$\sigma$}^{2(t_{k}+1)}$ converges to a limit
$\mbox{\boldmath$\sigma$}_{**}^{2}$. Taking limits in the relation
$\mbox{\boldmath$\sigma$}^{2(t_{k}+1)}=M(\mbox{\boldmath$\sigma$}^{2(t_{k})})$
and invoking the continuity $M(\mbox{\boldmath$\sigma$}^{2})$ imply that
$\mbox{\boldmath$\sigma$}_{**}^{2}=M(\mbox{\boldmath$\sigma$}_{*}^{2})$.
Because the sequence $L(\mbox{\boldmath$\sigma$}^{2(t_{k})})$ is monotonically
increasing in $k$ and bounded above on $C_{0}$, it converges to a limit
$L_{*}$. Hence, the continuity of $L(\mbox{\boldmath$\sigma$}^{2})$ implies

$\displaystyle L(\mbox{\boldmath$\sigma$}_{*}^{2})$ $\displaystyle=$
$\displaystyle\mathop{\;\>}\nolimits\lim_{k}L(\mbox{\boldmath$\sigma$}^{2(t_{k})})\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits
L_{*}\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits\lim_{k}L(\mbox{\boldmath$\sigma$}^{2(t_{k}+1)})\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits
L(\mbox{\boldmath$\sigma$}_{**}^{2})\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits
L(M(\mbox{\boldmath$\sigma$}_{*}^{2})).$

Lemma 4 therefore gives
$\mbox{\boldmath$\sigma$}_{**}^{2}=M(\mbox{\boldmath$\sigma$}_{*}^{2})=\mbox{\boldmath$\sigma$}_{*}^{2}$,
contradicting the bound
$\|\mbox{\boldmath$\sigma$}_{*}^{2}-\mbox{\boldmath$\sigma$}_{**}^{2}\|_{2}\geq\epsilon$
entailed by inequality (12). ∎

######  Theorem 1.

The MM sequence $\\{\mbox{\boldmath$\sigma$}^{2(t)}\\}_{t\geq 0}$ has at least
one limit point. Every limit point is a fixed point of
$M(\mbox{\boldmath$\sigma$}^{2})$. If the set of fixed points is discrete,
then the MM sequence converges to one of them. Finally, when the iterates
converge, their limit is a KKT point.

###### Proof.

The sequence $\\{\mbox{\boldmath$\sigma$}^{2(t)}\\}_{t\geq 0}$ is contained in
the super-level compact set $C_{0}$ defined in Lemma 5 and therefore admits a
convergent subsequence $\mbox{\boldmath$\sigma$}^{2(t_{k})}$ with limit
$\mbox{\boldmath$\sigma$}^{2(\infty)}$. As argued in Lemma 5,
$L(\mbox{\boldmath$\sigma$}^{2(\infty)})=L(M(\mbox{\boldmath$\sigma$}^{2(\infty)}))$.
Lemma 4 now implies that $\mbox{\boldmath$\sigma$}^{2(\infty)}$ is a fixed
point of the algorithm map $M(\mbox{\boldmath$\sigma$}^{2})$.

According to Ostrowski’s theorem (Lange,, 2010, Proposition 8.2.1), the set of
limit points of a bounded sequence
$\\{\mbox{\boldmath$\sigma$}^{2(t)}\\}_{t\geq 0}$ is connected and compact
provided
$\|\mbox{\boldmath$\sigma$}^{2(t+1)}-\mbox{\boldmath$\sigma$}^{2(t)}\|_{2}\to
0$. If the set of fixed points is discrete, then the connected subset of limit
points reduces to a single point. Hence, the bounded sequence
$\mbox{\boldmath$\sigma$}^{2(t)}$ converges to this point. When the limit
exists, one can check that $\mbox{\boldmath$\sigma$}^{2(\infty)}$ satisfies
the KKT conditions by proving that each zero component of
$\mbox{\boldmath$\sigma$}^{2(\infty)}$ has a non-positive partial derivative.
Suppose on the contrary $\sigma_{i}^{2(\infty)}=0$ and
$\frac{\partial}{\partial\sigma_{i}^{2}}L(\mbox{\boldmath$\sigma$}^{2(\infty)})>0$.
By continuity
$\frac{\partial}{\partial\sigma_{i}^{2}}L(\mbox{\boldmath$\sigma$}^{2(t)})>0$
for all large $t$. Therefore, $\sigma_{i}^{2(t+1)}>\sigma_{i}^{2(t)}$ for all
large $t$ by the observation made after equation (9). This behavior is
inconsistent with the assumption that $\sigma_{i}^{2(t)}\to 0$. ∎

##  5 MM versus EM

Examination of Tables 2 and 3 suggests that the MM algorithm usually converges
faster than the EM algorithm. We now provide theoretical justification for
this observation. Again for notational convenience, we consider the REML case
where ${\bm{X}}$ is null. Since the EM principle is just a special instance of
the MM principle, we can compare their convergence properties in a unified
framework. Consider an MM map $M(\mbox{\boldmath$\theta$})$ for maximizing the
objective function $f(\mbox{\boldmath$\theta$})$ via the surrogate function
$g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{(t)})$. Close to the
optimal point $\mbox{\boldmath$\theta$}^{\infty}$,

$\displaystyle\mbox{\boldmath$\theta$}^{(t+1)}-\mbox{\boldmath$\theta$}^{\infty}$
$\displaystyle\approx$ $\displaystyle
dM(\mbox{\boldmath$\theta$}^{\infty})(\mbox{\boldmath$\theta$}^{(t)}-\mbox{\boldmath$\theta$}^{\infty}),$

where $dM(\mbox{\boldmath$\theta$}^{\infty})$ is the differential of the
mapping $M$ at the optimal point $\mbox{\boldmath$\theta$}^{\infty}$ of
$f(\mbox{\boldmath$\theta$})$. Hence, the local convergence rate of the
sequence $\mbox{\boldmath$\theta$}^{(t+1)}=M(\mbox{\boldmath$\theta$}^{(t)})$
coincides with the spectral radius of $dM(\mbox{\boldmath$\theta$}^{\infty})$.
Familiar calculations (McLachlan and Krishnan,, 2008; Lange,, 2010)
demonstrate that

$\displaystyle dM(\mbox{\boldmath$\theta$}^{\infty})$ $\displaystyle=$
$\displaystyle{\bm{I}}-[d^{2}g(\mbox{\boldmath$\theta$}^{\infty}\mid\mbox{\boldmath$\theta$}^{\infty})]^{-1}d^{2}f(\mbox{\boldmath$\theta$}^{\infty}).$

In other words, the local convergence rate is determined by how well the
surrogate surface
$g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{\infty})$
approximates the objective surface $f(\mbox{\boldmath$\theta$})$ near the
optimal point $\mbox{\boldmath$\theta$}^{\infty}$. In the EM literature,
$dM(\mbox{\boldmath$\theta$}^{\infty})$ is called the _rate matrix_ (Meng and
Rubin,, 1991). Fast convergence occurs when the surrogate
$g(\mbox{\boldmath$\theta$}\mid\mbox{\boldmath$\theta$}^{\infty})$ hugs the
objective $f(\mbox{\boldmath$\theta$})$ tightly around
$\mbox{\boldmath$\theta$}^{\infty}$. Figure 1 shows a case where the MM
surrogate locally dominates the EM surrogate. We demonstrate that this is no
accident.

McLachlan and Krishnan, (2008) derive the EM surrogate

$\displaystyle
g_{\text{EM}}(\mbox{\boldmath$\sigma$}^{2}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(t)})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{m}\Big{[}\mathrm{rank}({\bm{V}}_{i})\ln\sigma_{i}^{2}+\mathrm{rank}({\bm{V}}_{i})\frac{\sigma_{i}^{2(t)}}{\sigma_{i}^{2}}-\frac{\sigma_{i}^{4(t)}}{\sigma_{i}^{2}}\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i})\Big{]}$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{m}\frac{\sigma_{i}^{4(t)}}{\sigma^{2}_{i}}{\bm{y}}^{T}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{y}}$

minorizing the log-likelihood up to an irrelevant constant. Section S.3 of the
Supplementary Materials gives a detailed derivation for the more general
multivariate response case. The rank of the covariance matrix ${\bm{V}}_{i}$
appears because ${\bm{V}}_{i}$ may not be invertible. Both of the surrogates
$g_{\text{EM}}(\mbox{\boldmath$\sigma$}^{2}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(\infty)})$
and
$g_{\text{MM}}(\mbox{\boldmath$\sigma$}^{2}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(\infty)})$
are parameter separated. This implies that both second differentials
$d^{2}g_{\text{EM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(\infty)})$
and
$d^{2}g_{\text{MM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(\infty)})$
are diagonal. A small diagonal entry of either matrix indicates fast
convergence of the corresponding variance component. Our next result shows
that, under Assumption 1, on average the diagonal entries of
$d^{2}g_{\text{EM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(\infty)})$
dominate those of
$d^{2}g_{\text{MM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(\infty)})$
when $m>2$. Thus, the EM algorithm tends to converge more slowly than the MM
algorithm, and the difference is more pronounced as the number of variance
components $m$ grows.

######  Theorem 2.

Let $\mbox{\boldmath$\sigma$}^{2(\infty)}\succ{\mathbf{0}}_{m}$ be a common
limit point of the EM and MM algorithms. Then both second differentials
$d^{2}g_{\text{MM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(\infty)})$
and
$d^{2}g_{\text{EM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(\infty)})$
are diagonal with

$\displaystyle
d^{2}g_{\text{EM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\mid\mbox{\boldmath$\sigma$}^{2(\infty)})_{ii}$
$\displaystyle=$
$\displaystyle-\frac{\mathrm{rank}({\bm{V}}_{i})}{2\sigma_{i}^{4(\infty)}}$
$\displaystyle
d^{2}g_{\text{MM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\mid\mbox{\boldmath$\sigma$}^{2(\infty)})_{ii}$
$\displaystyle=$
$\displaystyle-\frac{{\bm{y}}^{T}\mbox{\boldmath$\Omega$}^{-(\infty)}{\bm{V}}_{i}\mbox{\boldmath$\Omega$}^{-(\infty)}{\bm{y}}}{\sigma_{i}^{2(\infty)}}\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits-\frac{\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(\infty)}{\bm{V}}_{i})}{\sigma_{i}^{2(\infty)}}.$

Furthermore, the average ratio

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}\frac{d^{2}g_{\text{MM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\mid\mbox{\boldmath$\sigma$}^{2(\infty)})_{ii}}{d^{2}g_{\text{EM}}(\mbox{\boldmath$\sigma$}^{2(\infty)}\mid\mbox{\boldmath$\sigma$}^{2(\infty)})_{ii}}$
$\displaystyle=$
$\displaystyle\frac{2}{mn}\sum_{i=1}^{m}\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(\infty)}\sigma_{i}^{2(\infty)}{\bm{V}}_{i})\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits\frac{2}{m}\mathop{\;\>}\nolimits<\mathop{\;\>}\nolimits
1$

for $m>2$ when all ${\bm{V}}_{i}$ have full rank $n$.

###### Proof.

See Section LABEL:sec:local-conv-rate of the Supplementary Materials. ∎

Both the EM and MM algorithms must evaluate the traces
$\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i})$ and quadratic forms
$({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})^{T}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i}\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})$
at each iteration. Since these quantities are also the building blocks of the
approximate rate matrices
$d^{2}g(\mbox{\boldmath$\sigma$}^{2(t)}\\!\mid\\!\mbox{\boldmath$\sigma$}^{2(t)})$,
one can rationally choose either the EM or MM updates based on which has
smaller diagonal entries measured by the $\ell_{1}$, $\ell_{2}$, or
$\ell_{\infty}$ norms. At negligible extra cost, this produces a hybrid
algorithm that retains the ascent property and enjoys the better of the two
convergence rates.

##  6 Extensions

Besides its competitive numerical performance, Algorithm 1 is attractive for
its simplicity and ease of generalization. In this section, we outline MM
algorithms for multivariate response models possibly with missing data, linear
mixed models, MAP estimation, penalized estimation, and generalized estimating
equations.

###  6.1 Multivariate response model

Consider the multivariate response model with $n\times d$ response matrix
${\bm{Y}}$, mean $\mathrm{E}\,{\bm{Y}}={\bm{X}}{\bm{B}}$, and covariance

$\Omega$ $\displaystyle=$
$\displaystyle\mathrm{Cov}(\mathrm{vec}{\bm{Y}})\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits\sum_{i=1}^{m}\mbox{\boldmath$\Gamma$}_{i}\otimes{\bm{V}}_{i}.$

The $p\times d$ coefficient matrix ${\bm{B}}$ collects the fixed effects, the
$\mbox{\boldmath$\Gamma$}_{i}$ are unknown $d\times d$ covariance matrices,
and the ${\bm{V}}_{i}$ are known $n\times n$ covariance matrices. If the
vector $\mathrm{vec}{\bm{Y}}$ is normally distributed, then ${\bm{Y}}$ equals
a sum of independent matrix normal distributions (Gupta and Nagar,, 1999). We
now make this assumption and pursue estimation of ${\bm{B}}$ and the
$\mbox{\boldmath$\Gamma$}_{i}$, which we collectively denote as $\Gamma$.
Under the normality assumption, the Kronecker product identity
$\mathrm{vec}({\bm{C}}{\bm{D}}{\bm{E}})=({\bm{E}}^{T}\otimes{\bm{C}})\mathrm{vec}({\bm{D}})$
yields the log-likelihood

$\displaystyle L({\bm{B}},\mbox{\boldmath$\Gamma$})$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\ln\det\mbox{\boldmath$\Omega$}-\frac{1}{2}\mathrm{vec}({\bm{Y}}-{\bm{X}}{\bm{B}})^{T}\mbox{\boldmath$\Omega$}^{-1}\mathrm{vec}({\bm{Y}}-{\bm{X}}{\bm{B}})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\ln\det\mbox{\boldmath$\Omega$}-\frac{1}{2}[\mathrm{vec}{\bm{Y}}-({\bm{I}}_{d}\otimes{\bm{X}})\mathrm{vec}{\bm{B}}]^{T}\mbox{\boldmath$\Omega$}^{-1}[\mathrm{vec}{\bm{Y}}-({\bm{I}}_{d}\otimes{\bm{X}})\mathrm{vec}{\bm{B}}].$

Updating ${\bm{B}}$ given $\mbox{\boldmath$\Gamma$}^{(t)}$ is accomplished by
solving the general least squares problem met earlier in the univariate case.
Maximization of the log-likelihood (6.1) is difficult due to the requirement
that each $\mbox{\boldmath$\Gamma$}_{i}$ be positive semidefinite. Typical
solutions involve reparameterization of the covariance matrix (Pinheiro and
Bates,, 1996). The MM algorithm derived in this section gracefully
accommodates the covariance constraints.

Updating $\Gamma$ given ${\bm{B}}^{(t)}$ requires generalizing the
minorization (5). In view of Lemma 1 and the identities
$({\bm{A}}\otimes{\bm{B}})({\bm{C}}\otimes{\bm{D}})=({\bm{A}}{\bm{C}})\otimes({\bm{B}}{\bm{D}})$
and $({\bm{A}}\otimes{\bm{B}})^{-1}={\bm{A}}^{-1}\otimes{\bm{B}}^{-1}$, we
have

$\displaystyle\mbox{\boldmath$\Omega$}^{(t)}\mbox{\boldmath$\Omega$}^{-1}\mbox{\boldmath$\Omega$}^{(t)}$
$\displaystyle=$ $\displaystyle
m\left[\frac{1}{m}\sum_{i=1}^{m}\mbox{\boldmath$\Gamma$}_{i}^{(t)}\otimes{\bm{V}}_{i}\right]\left[\frac{1}{m}\sum_{i=1}^{m}\mbox{\boldmath$\Gamma$}_{i}\otimes{\bm{V}}_{i}\right]^{-1}\left[\frac{1}{m}\sum_{i=1}^{m}\mbox{\boldmath$\Gamma$}_{i}^{(t)}\otimes{\bm{V}}_{i}\right]$
$\displaystyle\preceq$ $\displaystyle
m\frac{1}{m}\sum_{i=1}^{m}(\mbox{\boldmath$\Gamma$}_{i}^{(t)}\otimes{\bm{V}}_{i})(\mbox{\boldmath$\Gamma$}_{i}\otimes{\bm{V}}_{i})^{-1}(\mbox{\boldmath$\Gamma$}_{i}^{(t)}\otimes{\bm{V}}_{i})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{m}(\mbox{\boldmath$\Gamma$}_{i}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{-1}\mbox{\boldmath$\Gamma$}_{i}^{(t)})\otimes{\bm{V}}_{i},$

or equivalently

$\displaystyle\mbox{\boldmath$\Omega$}^{-1}$ $\displaystyle\preceq$
$\displaystyle\mbox{\boldmath$\Omega$}^{-(t)}\Big{[}\sum_{i=1}^{m}(\mbox{\boldmath$\Gamma$}_{i}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{-1}\mbox{\boldmath$\Gamma$}_{i}^{(t)})\otimes{\bm{V}}_{i}\Big{]}\mbox{\boldmath$\Omega$}^{-(t)}.$
(14)

This derivation relies on the invertibility of the matrices ${\bm{V}}_{i}$.
One can relax this assumption by substituting
${\bm{V}}_{\epsilon,i}={\bm{V}}_{i}+\epsilon{\bm{I}}_{n}$ for ${\bm{V}}_{i}$
and sending $\epsilon$ to 0.

The majorization (14) and the minorization (6) jointly yield the surrogate

$\displaystyle g(\mbox{\boldmath$\Gamma$}\mid\mbox{\boldmath$\Gamma$}^{(t)})$
$\displaystyle\\!\\!=\\!\\!$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{m}\Big{\\{}\mathrm{tr}[\mbox{\boldmath$\Omega$}^{-(t)}(\mbox{\boldmath$\Gamma$}_{i}\otimes{\bm{V}}_{i})]\\!+\\!(\mathrm{vec}\,{\bm{R}}^{(t)})^{T}[(\mbox{\boldmath$\Gamma$}_{i}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{-1}\mbox{\boldmath$\Gamma$}_{i}^{(t)})\\!\otimes\\!{\bm{V}}_{i}]\,(\mathrm{vec}\,{\bm{R}}^{(t)})\Big{\\}}\\!+\\!c^{(t)},$

where ${\bm{R}}^{(t)}$ is the $n\times d$ matrix satisfying
$\mathrm{vec}\,{\bm{R}}^{(t)}=\mbox{\boldmath$\Omega$}^{-(t)}\mathrm{vec}({\bm{Y}}-{\bm{X}}{\bm{B}}^{(t)})$
and $c^{(t)}$ is an irrelevant constant. Based on the Kronecker identities
$(\mathrm{vec}\,{\bm{A}})^{T}\mathrm{vec}\,{\bm{B}}=\mathrm{tr}({\bm{A}}^{T}{\bm{B}})$
and
$\mathrm{vec}({\bm{C}}{\bm{D}}{\bm{E}})=({\bm{E}}^{T}\otimes{\bm{C}})\mathrm{vec}({\bm{D}})$,
the surrogate can be rewritten as

$\displaystyle g(\mbox{\boldmath$\Gamma$}\mid\mbox{\boldmath$\Gamma$}^{(t)})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{m}\Big{\\{}\mathrm{tr}[\mbox{\boldmath$\Omega$}^{-(t)}(\mbox{\boldmath$\Gamma$}_{i}\otimes{\bm{V}}_{i})]+\mathrm{tr}({\bm{R}}^{(t)T}{\bm{V}}_{i}{\bm{R}}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{-1}\mbox{\boldmath$\Gamma$}_{i}^{(t)})\Big{\\}}+c^{(t)}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{m}\Big{\\{}\mathrm{tr}[\mbox{\boldmath$\Omega$}^{-(t)}(\mbox{\boldmath$\Gamma$}_{i}\otimes{\bm{V}}_{i})]+\mathrm{tr}(\mbox{\boldmath$\Gamma$}_{i}^{(t)}{\bm{R}}^{(t)T}{\bm{V}}_{i}{\bm{R}}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{-1})\Big{\\}}+c^{(t)}.$

The first trace is linear in $\mbox{\boldmath$\Gamma$}_{i}$ with the
coefficient of entry $(\mbox{\boldmath$\Gamma$}_{i})_{jk}$ equal to

$\displaystyle\mathrm{tr}(\mbox{\boldmath$\Omega$}_{jk}^{-(t)}{\bm{V}}_{i})$
$\displaystyle=$
$\displaystyle{\mathbf{1}}_{n}^{T}({\bm{V}}_{i}\odot\mbox{\boldmath$\Omega$}^{-(t)}_{jk}){\mathbf{1}}_{n},$

where $\mbox{\boldmath$\Omega$}_{jk}^{-(t)}$ is the $(j,k)$-th $n\times n$
block of $\mbox{\boldmath$\Omega$}^{-(t)}$. The matrix ${\bm{M}}_{i}$ of these
coefficients can be written as

$\displaystyle{\bm{M}}_{i}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}{\mathbf{1}}^{T}&{\mathbf{0}}^{T}&\ldots&{\mathbf{0}}^{T}\\\
{\mathbf{0}}^{T}&{\mathbf{1}}^{T}&\ldots&{\mathbf{0}}^{T}\\\
\vdots&\vdots&\ddots&\vdots\\\
{\mathbf{0}}^{T}&{\mathbf{0}}^{T}&\ldots&{\mathbf{1}}^{T}\end{pmatrix}\left[\begin{pmatrix}{\bm{V}}_{i}&\ldots&{\bm{V}}_{i}\\\
\vdots&\ddots&\vdots\\\
{\bm{V}}_{i}&\ldots&{\bm{V}}_{i}\end{pmatrix}\odot\mbox{\boldmath$\Omega$}^{(-t)}\right]\begin{pmatrix}{\mathbf{1}}&{\mathbf{0}}&\ldots&{\mathbf{0}}\\\
{\mathbf{0}}&{\mathbf{1}}&\ldots&{\mathbf{0}}\\\
\vdots&\vdots&\ddots&\vdots\\\
{\mathbf{0}}&{\mathbf{0}}&\ldots&{\mathbf{1}}\end{pmatrix}$ $\displaystyle=$
$\displaystyle({\bm{I}}_{d}\otimes{\mathbf{1}}_{n})^{T}[({\mathbf{1}}_{d}{\mathbf{1}}_{d}^{T}\otimes{\bm{V}}_{i})\odot\mbox{\boldmath$\Omega$}^{-(t)}]({\bm{I}}_{d}\otimes{\mathbf{1}}_{n}).$

The directional derivative of
$g(\mbox{\boldmath$\Gamma$}\mid\mbox{\boldmath$\Gamma$}^{(t)})$ with respect
to $\mbox{\boldmath$\Gamma$}_{i}$ in the direction
$\mbox{\boldmath$\Delta$}_{i}$ is

$\displaystyle-\frac{1}{2}\mathrm{tr}({\bm{M}}_{i}\mbox{\boldmath$\Delta$}_{i})+\frac{1}{2}\mathrm{tr}(\mbox{\boldmath$\Gamma$}_{i}^{(t)}{\bm{R}}^{(t)T}{\bm{V}}_{i}{\bm{R}}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{-1}\mbox{\boldmath$\Delta$}_{i}\mbox{\boldmath$\Gamma$}_{i}^{-1})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\mathrm{tr}({\bm{M}}_{i}\mbox{\boldmath$\Delta$}_{i})+\frac{1}{2}\mathrm{tr}(\mbox{\boldmath$\Gamma$}_{i}^{-1}\mbox{\boldmath$\Gamma$}_{i}^{(t)}{\bm{R}}^{(t)T}{\bm{V}}_{i}{\bm{R}}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{-1}\mbox{\boldmath$\Delta$}_{i}).$

Because all directional derivatives of
$g(\mbox{\boldmath$\Gamma$}\mid\mbox{\boldmath$\Gamma$}^{(t)})$ vanish at a
stationarity point, the matrix equation

$\displaystyle{\bm{M}}_{i}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\Gamma$}_{i}^{-1}\mbox{\boldmath$\Gamma$}_{i}^{(t)}{\bm{R}}^{(t)T}{\bm{V}}_{i}{\bm{R}}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{-1}$
(15)

holds. Fortunately, this equation admits an explicit solution. For positive
scalers $a$ and $b$, the solution to the equation $b=x^{-1}ax^{-1}$ is
$x=\pm\sqrt{a/b}$. The matrix analogue of this equation is the Riccati
equation ${\bm{B}}={\bm{X}}^{-1}{\bm{A}}{\bm{X}}^{-1}$, whose solution is
summarized in the next lemma.

input : ${\bm{Y}}$, ${\bm{X}}$, ${\bm{V}}_{1},\ldots,{\bm{V}}_{m}$

output : MLE $\widehat{\bm{B}}$,
$\widehat{}\mbox{\boldmath$\Gamma$}_{1},\ldots,\widehat{}\mbox{\boldmath$\Gamma$}_{m}$

1 Initialize $\mbox{\boldmath$\Gamma$}_{i}^{(0)}$ positive definite,
$i=1,\ldots,m$

2 repeat

3
$\mbox{\boldmath$\Omega$}^{(t)}\leftarrow\sum_{i=1}^{m}\mbox{\boldmath$\Gamma$}_{i}^{(t)}\otimes{\bm{V}}_{i}$

4
${\bm{B}}^{(t)}\leftarrow\arg\min_{{\bm{B}}}\,[\mathrm{vec}{\bm{Y}}-({\bm{I}}_{d}\otimes{\bm{X}})\mathrm{vec}{\bm{B}}]^{T}\mbox{\boldmath$\Omega$}^{-(t)}[\mathrm{vec}{\bm{Y}}-({\bm{I}}_{d}\otimes{\bm{X}})\mathrm{vec}{\bm{B}}]$

5
${\bm{R}}^{(t)}\leftarrow\text{reshape}(\mbox{\boldmath$\Omega$}^{-(t)}\mathrm{vec}({\bm{Y}}-{\bm{X}}{\bm{B}}^{(t)}),n,d)$

6 for  _$i=1,\ldots,m$_ do

7 Cholesky
${\bm{L}}_{i}^{(t)}{\bm{L}}_{i}^{(t)T}\leftarrow({\bm{I}}_{d}\otimes{\mathbf{1}}_{n})^{T}[({\mathbf{1}}_{d}{\mathbf{1}}_{d}^{T}\otimes{\bm{V}}_{i})\odot\mbox{\boldmath$\Omega$}^{-(t)}]({\bm{I}}_{d}\otimes{\mathbf{1}}_{n})$

8
$\mbox{\boldmath$\Gamma$}_{i}^{(t+1)}\leftarrow{\bm{L}}_{i}^{-(t)T}[{\bm{L}}_{i}^{(t)T}(\mbox{\boldmath$\Gamma$}_{i}^{(t)}{\bm{R}}^{(t)T}{\bm{V}}_{i}{\bm{R}}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{(t)}){\bm{L}}_{i}^{(t)}]^{1/2}{\bm{L}}_{i}^{-(t)}$

9 end for

10

11until  _objective value converges_

Algorithm 3 The MM algorithm for MLE of the multivariate response model (6.1).

######  Lemma 6.

Assume ${\bm{A}}$ and ${\bm{B}}$ are positive definite and ${\bm{L}}$ is the
Cholesky factor of ${\bm{B}}$. Then
${\bm{Y}}={\bm{L}}^{-T}({\bm{L}}^{T}{\bm{A}}{\bm{L}})^{1/2}{\bm{L}}^{-1}$ is
the unique positive definite solution to the matrix equation
${\bm{B}}={\bm{X}}^{-1}{\bm{A}}{\bm{X}}^{-1}$.

###### Proof.

Direct substitution shows that ${\bm{Y}}$ solves the equivalent equation
${\bm{X}}{\bm{B}}{\bm{X}}={\bm{A}}$. To show uniqueness, suppose
${\bm{Y}}^{-1}{\bm{A}}{\bm{Y}}^{-1}={\bm{B}}$ and
${\bm{Z}}^{-1}{\bm{A}}{\bm{Z}}^{-1}={\bm{B}}$. The equations

$\displaystyle({\bm{B}}^{1/2}{\bm{Y}}{\bm{B}}^{1/2})^{2}$ $\displaystyle=$
$\displaystyle{\bm{B}}^{1/2}{\bm{Y}}{\bm{B}}{\bm{Y}}{\bm{B}}^{1/2}\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits{\bm{B}}^{1/2}{\bm{A}}{\bm{B}}^{1/2}$
$\displaystyle({\bm{B}}^{1/2}{\bm{Z}}{\bm{B}}^{1/2})^{2}$ $\displaystyle=$
$\displaystyle{\bm{B}}^{1/2}{\bm{Z}}{\bm{B}}{\bm{Z}}{\bm{B}}^{1/2}\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits{\bm{B}}^{1/2}{\bm{A}}{\bm{B}}^{1/2}$

imply
${\bm{B}}^{1/2}{\bm{Y}}{\bm{B}}^{1/2}={\bm{B}}^{1/2}{\bm{Z}}{\bm{B}}^{1/2}$ by
virtue of the uniqueness of symmetric square root. Since ${\bm{B}}^{-1/2}$ is
positive definite, ${\bm{Y}}={\bm{Z}}$. ∎

The Cholesky factor ${\bm{L}}$ in Lemma 6 can be replaced by the symmetric
square root of ${\bm{B}}$. The solution, which is unique, remains the same.
The Cholesky decomposition is preferred for its cheaper computational cost and
better numerical stability.

Algorithm 3 summarizes the MM algorithm for fitting the multi-response model
(3). Each iteration invokes $m$ Cholesky decompositions and symmetric square
roots of $d\times d$ positive definite matrices. Fortunately in most
applications, $d$ is a small number. The following result guarantees the non-
singularity of the Cholesky factor throughout iterations.

######  Proposition 2.

Assume ${\bm{V}}_{i}$ has strictly positive diagonal entries. Then the
symmetric matrix
${\bm{M}}_{i}=({\bm{I}}_{d}\otimes{\mathbf{1}}_{n})^{T}[({\mathbf{1}}_{d}{\mathbf{1}}_{d}^{T}\otimes{\bm{V}}_{i})\odot\mbox{\boldmath$\Omega$}^{-(t)}]({\bm{I}}_{d}\otimes{\mathbf{1}}_{n})$
is positive definite for all $t$. Furthermore if
$\mbox{\boldmath$\Gamma$}_{i}^{(0)}\succ{\mathbf{0}}$ and no column of
${\bm{R}}^{(t)}$ lies in the null space of ${\bm{V}}_{i}$ for all $t$, then
$\mbox{\boldmath$\Gamma$}_{i}^{(t)}\succ{\mathbf{0}}$ for all $t$.

###### Proof.

If ${\bm{V}}_{i}$ has strictly positive diagonal entries, then so does
${\mathbf{1}}_{d}{\mathbf{1}}_{d}^{T}\otimes{\bm{V}}_{i}$, and the Hadamard
product
$({\mathbf{1}}_{d}{\mathbf{1}}_{d}^{T}\otimes{\bm{V}}_{i})\odot\mbox{\boldmath$\Omega$}^{-(t)}$
is positive definite by Schur’s lemma. Since the matrix
${\bm{I}}_{d}\otimes{\mathbf{1}}_{n}$ has full column rank $d$, the matrix
${\bm{M}}_{i}$ is also positive definite. Finally, if no column of
${\bm{R}}^{(t)}$ lies in the null space of ${\bm{V}}_{i}$, and
$\mbox{\boldmath$\Gamma$}^{(t)}$ is positive define, then
$\mbox{\boldmath$\Gamma$}_{i}^{(t)}{\bm{R}}^{(t)T}{\bm{V}}_{i}{\bm{R}}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{(t)}$
is positive definite. The second claim follows by induction and Lemma 6. ∎

### Multivariate response, two variance components

When there are $m=2$ variance components
$\mbox{\boldmath$\Omega$}=\mbox{\boldmath$\Gamma$}_{1}\otimes{\bm{V}}_{1}+\mbox{\boldmath$\Gamma$}_{2}\otimes{\bm{V}}_{2}$,
repeated inversion of the $nd\times nd$ covariance matrix $\Omega$ reduces to
a single $nd\times nd$ simultaneous congruence decomposition and, per
iteration, two $d\times d$ Cholesky decompositions and one $d\times d$
simultaneous congruence decomposition. The simultaneous congruence
decomposition of the matrix pair $({\bm{V}}_{1},{\bm{V}}_{2})$ involves
generalized eigenvalues ${\bm{d}}=(d_{1},\ldots,d_{n})$ and a nonsingular
matrix ${\bm{U}}$ such that
${\bm{U}}^{T}{\bm{V}}_{1}{\bm{U}}={\bm{D}}=\mathrm{diag}({\bm{d}})$ and
${\bm{U}}^{T}{\bm{V}}_{2}{\bm{U}}={\bm{I}}$. If the simultaneous congruence
decomposition of
$(\mbox{\boldmath$\Gamma$}_{1}^{(t)},\mbox{\boldmath$\Gamma$}_{2}^{(t)})$ is
$(\mbox{\boldmath$\Lambda$}^{(t)},\mbox{\boldmath$\Phi$}^{(t)})$ with
$\mbox{\boldmath$\Phi$}^{(t)T}\mbox{\boldmath$\Gamma$}_{1}^{(t)}\mbox{\boldmath$\Phi$}^{(t)}=\mbox{\boldmath$\Lambda$}^{(t)}=\mathrm{diag}(\mbox{\boldmath$\lambda$}^{(t)})$
and
$\mbox{\boldmath$\Phi$}^{(t)T}\mbox{\boldmath$\Gamma$}_{2}^{(t)}\mbox{\boldmath$\Phi$}^{(t)}={\bm{I}}_{d}$,
then

$\displaystyle\mbox{\boldmath$\Omega$}^{(t)}$ $\displaystyle=$
$\displaystyle(\mbox{\boldmath$\Phi$}^{-(t)}\otimes{\bm{U}}^{-1})^{T}(\mbox{\boldmath$\Lambda$}^{(t)}\otimes{\bm{D}}+{\bm{I}}_{d}\otimes{\bm{I}}_{n})(\mbox{\boldmath$\Phi$}^{-(t)}\otimes{\bm{U}}^{-1})$
$\displaystyle\mbox{\boldmath$\Omega$}^{-(t)}$ $\displaystyle=$
$\displaystyle(\mbox{\boldmath$\Phi$}^{(t)}\otimes{\bm{U}})(\mbox{\boldmath$\Lambda$}^{(t)}\otimes{\bm{D}}+{\bm{I}}_{d}\otimes{\bm{I}}_{n})^{-1}(\mbox{\boldmath$\Phi$}^{(t)}\otimes{\bm{U}})^{T}$
$\displaystyle\det\mbox{\boldmath$\Omega$}^{(t)}$ $\displaystyle=$
$\displaystyle\det(\mbox{\boldmath$\Lambda$}^{(t)}\otimes{\bm{D}}+{\bm{I}}_{d}\otimes{\bm{I}}_{n})\det[(\mbox{\boldmath$\Phi$}^{-(t)}\otimes{\bm{U}}^{-1})^{T}(\mbox{\boldmath$\Phi$}^{-(t)}\otimes{\bm{U}}^{-1})]$
$\displaystyle=$
$\displaystyle\det(\mbox{\boldmath$\Lambda$}^{(t)}\otimes{\bm{D}}+{\bm{I}}_{d}\otimes{\bm{I}}_{n})\det(\mbox{\boldmath$\Gamma$}_{2}^{(t)}\otimes{\bm{V}}_{2})$
$\displaystyle=$
$\displaystyle\det(\mbox{\boldmath$\Lambda$}^{(t)}\otimes{\bm{D}}+{\bm{I}}_{d}\otimes{\bm{I}}_{n})\det(\mbox{\boldmath$\Gamma$}_{2}^{(t)})^{n}\det({\bm{V}}_{2})^{d}.$

Updating the fixed effects reduces to a weighted least squares problem for the
transformed responses $\tilde{\bm{Y}}={\bm{U}}^{T}{\bm{Y}}$, transformed
predictor matrix $\tilde{\bm{X}}={\bm{U}}^{T}{\bm{X}}$, and observation
weights $(\lambda_{k}^{(t)}d_{i}+1)^{-1}$. Algorithm 4 summarizes the
simplified MM algorithm. Detailed derivations are relegated to Section
LABEL:sec:derive-mvt-2vc-algo of the Supplementary Materials.

Input : ${\bm{Y}}$, ${\bm{X}}$, ${\bm{V}}_{1}$, ${\bm{V}}_{2}$

Output : MLE $\widehat{\bm{B}}$, $\widehat{}\mbox{\boldmath$\Gamma$}_{1}$,
$\widehat{}\mbox{\boldmath$\Gamma$}_{2}$

1 Simultaneous congruence decomposition:
$({\bm{D}},{\bm{U}})\leftarrow({\bm{V}}_{1},{\bm{V}}_{2})$

2 Transform data: $\tilde{\bm{Y}}\leftarrow{\bm{U}}^{T}{\bm{Y}}$,
$\tilde{\bm{X}}\leftarrow{\bm{U}}^{T}{\bm{X}}$

3 Initialize $\mbox{\boldmath$\Gamma$}_{1}^{(0)}$,
$\mbox{\boldmath$\Gamma$}_{2}^{(0)}$ positive definite

4 repeat

5 Simultaneous congruence decomposition
$(\mbox{\boldmath$\Lambda$}^{(t)},\mbox{\boldmath$\Phi$}^{(t)})\leftarrow(\mbox{\boldmath$\Gamma$}_{1}^{(t)},\mbox{\boldmath$\Gamma$}_{2}^{(t)})$

6
${\bm{B}}^{(t)}\leftarrow\arg\min_{{\bm{B}}}\,[\mathrm{vec}(\tilde{\bm{Y}}\mbox{\boldmath$\Phi$}^{(t)})-(\mbox{\boldmath$\Phi$}^{(t)T}\otimes\tilde{\bm{X}})\mathrm{vec}{\bm{B}}]^{T}(\mbox{\boldmath$\Lambda$}^{(t)}\otimes{\bm{D}}+{\bm{I}}_{d}\otimes{\bm{I}}_{n})^{-1}[\mathrm{vec}(\tilde{\bm{Y}}\mbox{\boldmath$\Phi$}^{(t)})-(\mbox{\boldmath$\Phi$}^{(t)T}\otimes\tilde{\bm{X}})\mathrm{vec}{\bm{B}}]$

7 Cholesky
${\bm{L}}_{1}^{(t)}{\bm{L}}_{1}^{(t)T}\leftarrow\mbox{\boldmath$\Phi$}^{(t)}\mathrm{diag}\left(\mathrm{tr}\left({\bm{D}}(\lambda_{k}^{(t)}{\bm{D}}+{\bm{I}}_{n})^{-1}\right),k=1,\ldots,d\right)\mbox{\boldmath$\Phi$}^{(t)T}$

8 Cholesky
${\bm{L}}_{2}^{(t)}{\bm{L}}_{2}^{(t)T}\leftarrow\mbox{\boldmath$\Phi$}^{(t)}\mathrm{diag}\left(\mathrm{tr}\left((\lambda_{k}^{(t)}{\bm{D}}+{\bm{I}}_{n})^{-1}\right),k=1,\ldots,d\right)\mbox{\boldmath$\Phi$}^{(t)T}$

9
${\bm{N}}_{1}^{(t)}\leftarrow{\bm{D}}^{1/2}[(\tilde{\bm{Y}}-\tilde{\bm{X}}{\bm{B}}^{(t)})\mbox{\boldmath$\Phi$}^{(t)})\oslash({\bm{d}}\mbox{\boldmath$\lambda$}^{(t)T}+{\mathbf{1}}_{n}{\mathbf{1}}_{d}^{T})]\mbox{\boldmath$\Lambda$}^{(t)}\mbox{\boldmath$\Phi$}^{-(t)}$

10
${\bm{N}}_{2}^{(t)}\leftarrow[(\tilde{\bm{Y}}-\tilde{\bm{X}}{\bm{B}}^{(t)})\mbox{\boldmath$\Phi$}^{(t)})\oslash({\bm{d}}\mbox{\boldmath$\lambda$}^{(t)T}+{\mathbf{1}}_{n}{\mathbf{1}}_{d}^{T})]\mbox{\boldmath$\Phi$}^{-(t)}$

11
$\mbox{\boldmath$\Gamma$}_{i}^{(t+1)}\leftarrow{\bm{L}}_{i}^{-(t)T}({\bm{L}}_{i}^{(t)T}{\bm{M}}_{i}^{(t)T}{\bm{M}}_{i}^{(t)}{\bm{L}}_{i}^{(t)})^{1/2}{\bm{L}}_{i}^{-(t)}$,
$i=1,2$

12

13until  _objective value converges_

Algorithm 4 MM algorithm for multivariate response model
$\mbox{\boldmath$\Omega$}=\mbox{\boldmath$\Gamma$}_{1}\otimes{\bm{V}}_{1}+\mbox{\boldmath$\Gamma$}_{2}\otimes{\bm{V}}_{2}$
with two variance components matrices. Note that $\oslash$ denotes a Hadamard
quotient.

###  6.2 Multivariate response model with missing responses

In many applications the multivariate response model (6.1) involves missing
responses. For instance, in testing multiple longitudinal traits in genetics,
some trait values $y_{ij}$ may be missing due to dropped patient visits, while
their genetic covariates are complete. Missing data destroys the symmetry of
the log-likelihood (6.1) and complicates finding the MLE. Fortunately, MM
algorithm 3 easily adapts to this challenge.

The familiar EM argument (McLachlan and Krishnan,, 2008, Section 2.2) shows
that

$\displaystyle-\frac{n}{2}\ln\det\mbox{\boldmath$\Omega$}^{(t)}-\frac{1}{2}\mathrm{tr}\\{\mbox{\boldmath$\Omega$}^{-(t)}[\mathrm{vec}({\bm{Z}}^{(t)}-{\bm{X}}{\bm{B}}^{(t)})\mathrm{vec}({\bm{Z}}^{(t)}-{\bm{X}}{\bm{B}}^{(t)})^{T}+{\bm{C}}^{(t)}]\\}$
(16)

minorizes the observed log-likelihood at the current iterate
$({\bm{B}}^{(t)},\mbox{\boldmath$\Gamma$}_{1}^{(t)},\ldots,\mbox{\boldmath$\Gamma$}_{m}^{(t)})$.
Here ${\bm{Z}}^{(t)}$ is the completed response matrix given the observed
responses ${\bm{Y}}_{\text{obs}}^{(t)}$ and the current parameter values. The
complete data ${\bm{Y}}$ is assumed to be normally distributed
$N(\mathrm{vec}({\bm{X}}{\bm{B}}^{(t)}),\mbox{\boldmath$\Omega$}^{(t)})$. The
block matrix ${\bm{C}}^{(t)}$ is 0 except for a lower-right block consisting
of a Schur complement.

To maximize the surrogate (16), we invoke the familiar minorization (6) and
majorization (14) to separate the variance components
$\mbox{\boldmath$\Gamma$}_{i}$. At each iteration we impute missing entries by
their conditional means, compute their conditional variances and covariances
to supply the Schur complement, and then update the fixed effects and variance
components by the explicit updates of Algorithm 3. The required conditional
means and conditional variances can be conveniently obtained in the process of
inverting $\mbox{\boldmath$\Omega$}^{(t)}$ by the sweep operator of
computational statistics (Lange,, 2010, Section 7.3).

###  6.3 Linear mixed model (LMM)

The linear mixed model plays a central role in longitudinal data analysis. For
the sake of simplicity, consider the single-level LMM (Laird and Ware,, 1982;
Bates and Pinheiro,, 1998) for $n$ independent data clusters
$({\bm{y}}_{i},{\bm{X}}_{i},{\bm{Z}}_{i})$ with

$\displaystyle{\bm{Y}}_{i}$ $\displaystyle=$
$\displaystyle{\bm{X}}_{i}\mbox{\boldmath$\beta$}+{\bm{Z}}_{i}\mbox{\boldmath$\gamma$}_{i}+\mbox{\boldmath$\epsilon$}_{i},\quad
i=1,\ldots,n,$

where $\beta$ is a vector of fixed effects, the
$\mbox{\boldmath$\gamma$}_{i}\sim
N({\mathbf{0}},{\bm{R}}_{i}(\mbox{\boldmath$\theta$}))$ are independent random
effects, and $\mbox{\boldmath$\epsilon$}_{i}\sim
N({\mathbf{0}},\sigma^{2}{\bm{I}}_{n_{i}})$ captures random noise independent
of $\mbox{\boldmath$\gamma$}_{i}$. We assume the matrices ${\bm{Z}}_{i}$ have
full column rank. The within-cluster covariance matrices
${\bm{R}}_{i}(\mbox{\boldmath$\theta$})$ depend on a parameter vector
$\theta$; typical choices for ${\bm{R}}_{i}(\mbox{\boldmath$\theta$})$ impose
autocorrelation, compound symmetry, or unstructured correlation. It is clear
that ${\bm{Y}}_{i}$ is normal with mean ${\bm{X}}_{i}\mbox{\boldmath$\beta$}$,
covariance
$\mbox{\boldmath$\Omega$}_{i}={\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}){\bm{Z}}_{i}^{T}+\sigma^{2}{\bm{I}}_{n_{i}}$,
and log-likelihood

$\displaystyle
L_{i}(\mbox{\boldmath$\beta$},\mbox{\boldmath$\theta$},\sigma^{2})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\ln\det\mbox{\boldmath$\Omega$}_{i}-\frac{1}{2}({\bm{y}}_{i}-{\bm{X}}_{i}\mbox{\boldmath$\beta$})^{T}\mbox{\boldmath$\Omega$}_{i}^{-1}({\bm{y}}_{i}-{\bm{X}}_{i}\mbox{\boldmath$\beta$}).$

The next three facts about pseudo-inverses are used in deriving the MM
algorithm for LMM.

######  Lemma 7.

If ${\bm{A}}$ has full column rank and ${\bm{B}}$ has full row rank, then
$({\bm{A}}{\bm{B}})^{+}={\bm{B}}^{+}{\bm{A}}^{+}$.

###### Proof.

Under the hypotheses, the representations
${\bm{A}}^{+}=({\bm{A}}^{T}{\bm{A}})^{+}{\bm{A}}^{T}=({\bm{A}}^{T}{\bm{A}})^{-1}{\bm{A}}^{T}$
and ${\bm{B}}^{+}={\bm{B}}^{T}({\bm{B}}{\bm{B}}^{T})^{-1}$ are well known. The
choice
${\bm{B}}^{+}{\bm{A}}^{+}={\bm{B}}^{T}({\bm{B}}{\bm{B}}^{T})^{-1}({\bm{A}}^{T}{\bm{A}})^{-1}{\bm{A}}^{T}$
satisfies the four equations characterizing the pseudo-inverse of
${\bm{A}}{\bm{B}}$. ∎

######  Lemma 8.

If ${\bm{A}}$ and ${\bm{B}}$ are positive semidefinite matrices with the same
range, then

$\displaystyle\lim_{\epsilon\downarrow
0}\,({\bm{B}}+\epsilon{\bm{I}})({\bm{A}}+\epsilon{\bm{I}})^{-1}({\bm{B}}+\epsilon{\bm{I}})$
$\displaystyle=$ $\displaystyle{\bm{B}}{\bm{A}}^{+}{\bm{B}}.$

###### Proof.

Suppose ${\bm{A}}$ has spectral decomposition
$\sum_{i}\lambda_{i}{\bm{u}}_{i}{\bm{u}}_{i}^{T}$. The matrix
${\bm{P}}=\sum_{\lambda_{i}>0}{\bm{u}}_{i}{\bm{u}}_{i}^{T}$ projects onto the
range of ${\bm{A}}$ and therefore also projects onto the range of ${\bm{B}}$.
It follows that ${\bm{P}}{\bm{B}}={\bm{B}}$ and by symmetry that
${\bm{B}}{\bm{P}}={\bm{B}}$. This allows us to write

$\displaystyle({\bm{B}}+\epsilon{\bm{I}})({\bm{A}}+\epsilon{\bm{I}})^{-1}({\bm{B}}+\epsilon{\bm{I}})$
$\displaystyle=$
$\displaystyle{\bm{B}}{\bm{P}}({\bm{A}}+\epsilon{\bm{I}})^{-1}{\bm{P}}{\bm{B}}+\epsilon{\bm{B}}{\bm{P}}({\bm{A}}+\epsilon{\bm{I}})^{-1}+\epsilon({\bm{A}}+\epsilon{\bm{I}})^{-1}{\bm{P}}{\bm{B}}+\epsilon^{2}({\bm{A}}+\epsilon{\bm{I}})^{-1}.$

The last three of these terms vanish as $\epsilon\downarrow 0$; the first term
tends to the claimed limit. These assertions follow from the expressions

$\displaystyle{\bm{P}}({\bm{A}}+\epsilon{\bm{I}})^{-1}{\bm{P}}$
$\displaystyle=$
$\displaystyle{\bm{P}}({\bm{A}}+\epsilon{\bm{I}})^{-1}\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits({\bm{A}}+\epsilon{\bm{I}})^{-1}{\bm{P}}\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits\sum_{\lambda_{i}>0}\frac{1}{\lambda_{i}+\epsilon}{\bm{u}}_{i}{\bm{u}}_{i}^{T}$

and
$\epsilon^{2}({\bm{A}}+\epsilon{\bm{I}})^{-1}=\sum_{i}\frac{\epsilon^{2}}{\lambda_{i}+\epsilon}{\bm{u}}_{i}{\bm{u}}_{i}^{T}$.
∎

######  Lemma 9.

If ${\bm{R}}$ and ${\bm{S}}$ are positive definite matrices, and the
conformable matrix ${\bm{Z}}$ has full column rank, then the matrices
${\bm{Z}}{\bm{R}}{\bm{Z}}^{T}$ and ${\bm{Z}}{\bm{S}}{\bm{Z}}^{T}$ share a
common range.

###### Proof.

In fact, both matrices have range equal to the range of ${\bm{Z}}$. The
matrices ${\bm{Z}}$ and ${\bm{Z}}{\bm{R}}^{1/2}$ clearly have the same range.
Furthermore, the matrices ${\bm{Z}}{\bm{R}}^{1/2}$ and
${\bm{Z}}{\bm{R}}^{1/2}{\bm{R}}^{1/2}{\bm{Z}}^{T}$ also have the same range. ∎

The convexity of the map
$({\bm{X}},{\bm{Y}})\mapsto{\bm{X}}^{T}{\bm{Y}}^{-1}{\bm{X}}$ and Lemmas 7, 8,
and 9 now yield via the obvious limiting argument the majorization

$\displaystyle\mbox{\boldmath$\Omega$}^{(t)}\mbox{\boldmath$\Omega$}^{-1}\mbox{\boldmath$\Omega$}^{(t)}$
$\displaystyle\\!\\!=\\!\\!$
$\displaystyle({\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}^{(t)}){\bm{Z}}_{i}^{T}\\!+\\!\sigma^{2(t)}{\bm{I}}_{n_{i}})({\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}){\bm{Z}}_{i}^{T}\\!+\\!\sigma^{2}{\bm{I}}_{n_{i}})^{-1}({\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}^{(t)}){\bm{Z}}_{i}^{T}\\!+\\!\sigma^{2(t)}{\bm{I}}_{n_{i}})$
$\displaystyle\\!\\!\preceq\\!\\!$
$\displaystyle({\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}^{(t)}){\bm{Z}}_{i}^{T})({\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}){\bm{Z}}_{i}^{T})^{+}({\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}^{(t)}){\bm{Z}}_{i}^{T})+\frac{\sigma^{4(t)}}{\sigma^{2}}{\bm{I}}_{n_{i}}$
$\displaystyle\\!\\!=\\!\\!$
$\displaystyle[{\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}^{(t)}){\bm{Z}}_{i}^{T}{\bm{Z}}_{i}^{T+}]{\bm{R}}_{i}^{-1}(\mbox{\boldmath$\theta$})[{\bm{Z}}_{i}^{+}{\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}^{(t)}){\bm{Z}}_{i}^{T}]+\frac{\sigma^{4(t)}}{\sigma^{2}}{\bm{I}}_{n_{i}}$

In combination with the minorization (6), this gives the surrogate

$\displaystyle
g_{i}(\mbox{\boldmath$\theta$},\sigma^{2}\mid\mbox{\boldmath$\theta$}^{(t)},\sigma^{2(t)})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\mathrm{tr}({\bm{Z}}_{i}^{T}\mbox{\boldmath$\Omega$}_{i}^{-(t)}{\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}))-\frac{1}{2}{\bm{r}}_{i}^{(t)T}{\bm{R}}_{i}^{-1}(\mbox{\boldmath$\theta$}){\bm{r}}_{i}^{(t)}$
$\displaystyle-\frac{\sigma^{2}}{2}\mathrm{tr}(\mbox{\boldmath$\Omega$}_{i}^{-(t)})-\frac{\sigma^{4(t)}}{2\sigma^{2}}({\bm{y}}_{i}-{\bm{X}}_{i}\mbox{\boldmath$\beta$}^{(t)})^{T}\mbox{\boldmath$\Omega$}_{i}^{-2(t)}({\bm{y}}_{i}-{\bm{X}}_{i}\mbox{\boldmath$\beta$}^{(t)})+c^{(t)},$

for the log-likelihood $L_{i}(\mbox{\boldmath$\theta$},\sigma^{2})$, where

$\displaystyle{\bm{r}}_{i}^{(t)}$ $\displaystyle=$
$\displaystyle({\bm{Z}}_{i}^{+}{\bm{Z}}_{i}{\bm{R}}_{i}(\mbox{\boldmath$\theta$}^{(t)}){\bm{Z}}_{i}^{T})\mbox{\boldmath$\Omega$}_{i}^{-(t)}({\bm{y}}_{i}-{\bm{X}}_{i}\mbox{\boldmath$\beta$}^{(t)})\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits{\bm{R}}_{i}(\mbox{\boldmath$\theta$}^{(t)}){\bm{Z}}_{i}^{T}\mbox{\boldmath$\Omega$}_{i}^{-(t)}({\bm{y}}_{i}-{\bm{X}}_{i}\mbox{\boldmath$\beta$}^{(t)}).$

The parameters $\theta$ and $\sigma^{2}$ are nicely separated. To maximize the
overall minorization function
$\sum_{i}g_{i}(\mbox{\boldmath$\theta$},\sigma^{2}\mid\mbox{\boldmath$\theta$}^{(t)},\sigma^{2(t)})$,
we update $\sigma^{2}$ via

$\displaystyle\sigma^{2(t+1)}$ $\displaystyle=$
$\displaystyle\sigma^{2(t)}\sqrt{\frac{\sum_{i}({\bm{y}}_{i}-{\bm{X}}_{i}\mbox{\boldmath$\beta$}^{(t)})^{T}\mbox{\boldmath$\Omega$}_{i}^{-2(t)}({\bm{y}}_{i}-{\bm{X}}_{i}\mbox{\boldmath$\beta$}^{(t)})}{\sum_{i}\mathrm{tr}(\mbox{\boldmath$\Omega$}_{i}^{-(t)})}}.$

For structured models such as autocorrelation and compound symmetry, updating
$\theta$ is a low-dimensional optimization problem that can be approached
through the stationarity condition

$\displaystyle\sum_{i}\mathrm{vec}\left({\bm{Z}}_{i}^{T}\mbox{\boldmath$\Omega$}_{i}^{(t)}{\bm{Z}}_{i}-{\bm{R}}_{i}^{-1}(\mbox{\boldmath$\theta$}){\bm{r}}_{i}^{(t)}{\bm{r}}_{i}^{(t)T}{\bm{R}}_{i}^{-1}(\mbox{\boldmath$\theta$})\right)^{T}\frac{\partial}{\partial\theta_{j}}\mathrm{vec}\,{\bm{R}}_{i}(\mbox{\boldmath$\theta$})$
$\displaystyle=$ $\displaystyle 0$

for each component $\theta_{j}$. For the unstructured model with
${\bm{R}}_{i}(\mbox{\boldmath$\theta$})={\bm{R}}$ for all $i$, the
stationarity condition reads

$\displaystyle\sum_{i}{\bm{Z}}_{i}^{T}\mbox{\boldmath$\Omega$}_{i}^{(t)}{\bm{Z}}_{i}$
$\displaystyle=$
$\displaystyle{\bm{R}}^{-1}\left(\sum_{i}{\bm{r}}_{i}^{(t)}{\bm{r}}_{i}^{(t)T}\right){\bm{R}}^{-1}$

and admits an explicit solution based on Lemma 6.

Similar tactics apply to a multilevel LMM (Bates and Pinheiro,, 1998) with
responses

$\displaystyle{\bm{Y}}_{i}$ $\displaystyle=$
$\displaystyle{\bm{X}}_{i}\mbox{\boldmath$\beta$}+{\bm{Z}}_{i1}\mbox{\boldmath$\gamma$}_{i1}+\cdots{\bm{Z}}_{im}\mbox{\boldmath$\gamma$}_{im}+\mbox{\boldmath$\epsilon$}_{i}.$

Minorization separates parameters for each level (variance component).
Depending on the complexity of the covariance matrices, maximization of the
surrogate can be accomplished analytically. For the sake of brevity, details
are omitted.

###  6.4 MAP estimation

Suppose $\beta$ follows an improper flat prior, the variance components
$\sigma_{i}^{2}$ follow inverse gamma priors with shapes $\alpha_{i}>0$ and
scales $\gamma_{i}>0$, and these priors are independent. The log-posterior
density then reduces to

$\displaystyle\ln
f(\mbox{\boldmath$\beta$},\mbox{\boldmath$\sigma$}^{2}|{\bm{y}},{\bm{X}})$
(17) $\displaystyle=$
$\displaystyle-\frac{1}{2}\ln\det\mbox{\boldmath$\Omega$}-\frac{1}{2}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})^{T}\mbox{\boldmath$\Omega$}^{-1}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$})-\sum_{i=1}^{m}(\alpha_{i}+1)\ln\sigma_{i}^{2}-\sum_{i=1}^{m}\frac{\gamma_{i}}{\sigma_{i}^{2}}+c,$

where $c$ is an irrelevant constant. The MAP estimator of
$(\mbox{\boldmath$\beta$},\mbox{\boldmath$\sigma$}^{2})$ is the mode of the
posterior distribution. The update (4) of $\beta$ given
$\mbox{\boldmath$\sigma$}^{2}$ remains the same. To update
$\mbox{\boldmath$\sigma$}^{2}$ given $\beta$, apply the same minorizations (5)
and (6) to the first first two terms of equation (17). This separates
parameters and yields a convex surrogate for each $\sigma_{i}^{2}$. The
minimum of the $\sigma_{i}^{2}$ surrogate is defined by the stationarity
condition

$\displaystyle 0$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i})+\frac{\sigma_{i}^{4(t)}}{2\sigma_{i}^{4}}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})^{T}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i}\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})-\frac{\alpha_{i}+1}{\sigma_{i}^{2}}+\frac{\gamma_{i}}{\sigma_{i}^{4}}.$

Multiplying this by $\sigma_{i}^{4}$ gives a quadratic equation in
$\sigma_{i}^{2}$. The positive root should be taken to meet the nonnegativity
constraint on $\sigma_{i}^{2}$.

For the multivariate response model (6.1), we assume the variance components
$\mbox{\boldmath$\Gamma$}_{i}$ follow independent inverse Wishart
distributions with degrees of freedom $\nu_{i}>d-1$ and scale matrix
$\mbox{\boldmath$\Psi$}_{i}\succ{\mathbf{0}}$. The log density of the
posterior distribution is

$\displaystyle L({\bm{B}},\mbox{\boldmath$\Gamma$}|{\bm{X}},{\bm{Y}})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\ln\det\mbox{\boldmath$\Omega$}-\frac{1}{2}\mathrm{vec}({\bm{Y}}-{\bm{X}}{\bm{B}})^{T}\mbox{\boldmath$\Omega$}^{-1}\mathrm{vec}({\bm{Y}}-{\bm{X}}{\bm{B}})$
(18)
$\displaystyle-\frac{1}{2}\sum_{i=1}^{m}(\nu_{i}+d+1)\ln\det\mbox{\boldmath$\Gamma$}_{i}-\frac{1}{2}\sum_{i=1}^{m}\mathrm{tr}(\mbox{\boldmath$\Psi$}_{i}\mbox{\boldmath$\Gamma$}_{i}^{-1})+c,$

where $c$ is an irrelevant constant. Invoking the minorizations (6) and (14)
for the first two terms and the supporting hyperplane minorization

$\displaystyle-\ln\det\mbox{\boldmath$\Gamma$}_{i}$ $\displaystyle\geq$
$\displaystyle-\ln\det\mbox{\boldmath$\Gamma$}_{i}^{(t)}-\mathrm{tr}\\{\mbox{\boldmath$\Gamma$}_{i}^{-(t)}(\mbox{\boldmath$\Gamma$}_{i}-\mbox{\boldmath$\Gamma$}_{i}^{(t)})\\}$

for $-\ln\det\mbox{\boldmath$\Gamma$}_{i}$ gives the surrogate function

$\displaystyle g(\mbox{\boldmath$\Gamma$}|\mbox{\boldmath$\Gamma$}^{(t)})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{m}\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}(\mbox{\boldmath$\Gamma$}_{i}\otimes{\bm{V}}_{i}))-\frac{1}{2}\sum_{i=1}^{m}\mathrm{tr}\left(\mbox{\boldmath$\Gamma$}_{i}^{(t)}{\bm{R}}^{(t)T}{\bm{V}}_{i}{\bm{R}}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{-1}\right)$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{m}(\nu_{i}+d+1)\mathrm{tr}(\mbox{\boldmath$\Gamma$}_{i}^{-(t)}\mbox{\boldmath$\Gamma$}_{i})-\frac{1}{2}\sum_{i=1}^{m}\mathrm{tr}(\mbox{\boldmath$\Psi$}_{i}\mbox{\boldmath$\Gamma$}_{i}^{-1})+c^{(t)}.$

The optimal $\mbox{\boldmath$\Gamma$}_{i}$ satisfies the stationarity
condition

$\displaystyle({\bm{I}}_{d}\otimes{\mathbf{1}}_{n})^{T}[({\mathbf{1}}_{d}{\mathbf{1}}_{d}^{T}\otimes{\bm{V}}_{i})\odot\mbox{\boldmath$\Omega$}^{-(t)}]({\bm{I}}_{d}\otimes{\mathbf{1}}_{n})+(\nu_{i}+d+1)\mbox{\boldmath$\Gamma$}_{i}^{-(t)}$
$\displaystyle=$
$\displaystyle\mbox{\boldmath$\Gamma$}_{i}^{-1}(\mbox{\boldmath$\Gamma$}_{i}^{(t)}{\bm{R}}^{(t)T}{\bm{V}}_{i}{\bm{R}}^{(t)}\mbox{\boldmath$\Gamma$}_{i}^{(t)}+\mbox{\boldmath$\Psi$}_{i})\mbox{\boldmath$\Gamma$}_{i}^{-1}$

and can be found using Lemma 6.

###  6.5 Variable selection

In the statistical analysis of high-dimensional data, the imposition of
sparsity leads to better interpretation and more stable parameter estimation.
MM algorithms mesh well with penalized estimation. The simple variance
components model (1) illustrates this fact. For the selection of fixed
effects, minimizing the lasso-penalized log-likelihood

$\displaystyle-L(\mbox{\boldmath$\beta$},\mbox{\boldmath$\sigma$}^{2})+\lambda\sum_{j}|\beta_{j}|$

is often recommended (Schelldorfer et al.,, 2011). The only change to the MM
Algorithm 1 is that in estimating $\beta$, one solves a lasso penalized
general least squares problem rather than an ordinary general least squares
problem. The updates of the variance components $\sigma_{i}^{2}$ remain the
same. For selection among a large number of variance components, one can
minimize the ridge-penalized log-likelihood

$\displaystyle-L(\mbox{\boldmath$\beta$},\mbox{\boldmath$\sigma$}^{2})+\lambda\sum_{i=1}^{m}\sigma_{i}^{2}$

subject to the nonnegativity constraints $\sigma_{i}^{2}\geq 0$. Here the
standard deviations $\sigma_{i}$ are the underlying parameters. The variance
update (8) becomes

$\displaystyle\sigma_{i}^{2(t+1)}$ $\displaystyle=$
$\displaystyle\sigma_{i}^{2(t)}\sqrt{\frac{({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})^{T}\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i}\mbox{\boldmath$\Omega$}^{-(t)}({\bm{y}}-{\bm{X}}\mbox{\boldmath$\beta$}^{(t)})}{\mathrm{tr}(\mbox{\boldmath$\Omega$}^{-(t)}{\bm{V}}_{i})+2\lambda}},\quad
i=1,\ldots,m.$

The updates for the fixed effects $\beta$ are unaffected. Equation (6.5)
clearly exhibits shrinkage but no thresholding. The lasso penalized log-
likelihood

$\displaystyle-L(\mbox{\boldmath$\beta$},\mbox{\boldmath$\sigma$}^{2})+\lambda\sum_{i=1}^{m}\sigma_{i}$
(19)

subject to nonnegativity constraint $\sigma_{i}\geq 0$ achieves both ends. The
update of $\sigma_{i}$ is chosen among the positive roots of a quartic
equation and the boundary 0, whichever yields a lower objective value.

###  6.6 Beyond the linear model

One can extend the MM algorithms to binary and discrete response data with the
framework of generalized estimating equations (GEE) (Liang and Zeger,, 1986).
Again consider $n$ independent data clusters $({\bm{y}}_{i},{\bm{X}}_{i})$. In
longitudinal studies, $({\bm{y}}_{i},{\bm{X}}_{i})$ would be the responses and
clinical covariates of subject $i$ at different time points. In genetic
studies, $({\bm{y}}_{i},{\bm{X}}_{i})$ would be the trait values and
covariates of individuals within family $i$. GEE captures the within-cluster
correlation by specifying the first two moments of the conditional
distribution of ${\bm{y}}_{i}$ given ${\bm{X}}_{i}$, namely
$\mu_{ij}=E(y_{ij}|{\bm{x}}_{ij})$ and
$\sigma^{2}_{ij}=\text{Var}(y_{ij}|{\bm{x}}_{ij})$. If one assumes that
$y_{ij}$ follows an exponential family with canonical link, then

$\displaystyle\mu_{ij}(\mbox{\boldmath$\beta$})$ $\displaystyle=$
$\displaystyle\mu(\theta_{ij})\;\;\textrm{ and
}\;\;\sigma^{2}_{ij}(\mbox{\boldmath$\beta$})\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits\phi\mu^{\prime}(\theta_{ij}),\quad
i=1,\ldots,n,\;j=1,\ldots,n_{i},$

where $\mu(t)$ is a differentiable canonical link function, $\mu^{\prime}(t)$
is its first derivative,
$\theta_{ij}={\bm{x}}_{ij}^{T}\mbox{\boldmath$\beta$}$ is the linear
systematic part of $y_{ij}$ associated with the covariates, $\phi$ is an over-
dispersion parameter, and $\beta$ is the vector of fixed effects.

The GEE estimator of $\beta$ solves the equation

$\displaystyle\sum_{i=1}^{n}d\mbox{\boldmath$\mu$}_{i}(\mbox{\boldmath$\beta$})^{T}{\bm{V}}_{i}^{-1}(\mbox{\boldmath$\alpha$},\mbox{\boldmath$\beta$})[{\bm{y}}_{i}-\mbox{\boldmath$\mu$}_{i}(\mbox{\boldmath$\beta$})]$
$\displaystyle=$ $\displaystyle{\mathbf{0}}_{p},$

where ${\bm{y}}_{i}=(y_{i1},\dots,y_{in_{i}})^{T}$,
$\mbox{\boldmath$\mu$}_{i}(\mbox{\boldmath$\beta$})=[\mu_{i1}(\mbox{\boldmath$\beta$}),\ldots,\mu_{in_{i}}(\mbox{\boldmath$\beta$})]^{T}$,
${\bm{V}}_{i}=\text{cov}({\bm{y}}_{i})$ is the working covariance matrix of
the $i$-th subject, and $d\mbox{\boldmath$\mu$}_{i}(\mbox{\boldmath$\beta$})$
is the differential of $\mbox{\boldmath$\mu$}_{i}(\mbox{\boldmath$\beta$})$.
In longitudinal studies, ${\bm{V}}_{i}$ is often parameterized as
${\bm{V}}_{i}(\mbox{\boldmath$\alpha$},\mbox{\boldmath$\beta$})={\bm{A}}_{i}^{1/2}(\mbox{\boldmath$\beta$}){\bm{R}}_{i}(\mbox{\boldmath$\alpha$}){\bm{A}}_{i}^{1/2}(\mbox{\boldmath$\beta$})$,
where ${\bm{A}}^{1/2}(\mbox{\boldmath$\beta$})$ is a diagonal matrix with
standard deviations $\sigma_{ij}$ along its diagonal, and
${\bm{R}}(\mbox{\boldmath$\alpha$})$ is a correlation matrix with parameters
$\alpha$. This parameterization is too restrictive in many other applications.
For instance, in genetic studies, it is critical to dissect the variance into
different sources such as additive, dominance, and household environment
(Lange,, 2002). This suggests the variance component parameterization

$\displaystyle{\bm{R}}_{i}(\mbox{\boldmath$\sigma$}^{2})$ $\displaystyle=$
$\displaystyle\sigma_{a}^{2}2\mbox{\boldmath$\Phi$}_{i}+\sigma_{d}^{2}\mbox{\boldmath$\Delta$}_{7,i}+\sigma_{h}^{2}{\bm{H}}_{i}+\sigma_{e}^{2}{\bm{I}}_{n_{i}},\quad\sigma_{a}^{2}+\sigma_{d}^{2}+\sigma_{h}^{2}+\sigma_{e}^{2}\mathop{\;\>}\nolimits=\mathop{\;\>}\nolimits
1,$

where in the $i$-th family $\mbox{\boldmath$\Phi$}_{i}$ is the theoretical
kinship matrix, $\mbox{\boldmath$\Delta$}_{7,i}$ is the dominance variance
matrix, and ${\bm{H}}_{i}$ is the household indicator matrix. The matrices
$\mbox{\boldmath$\Phi$}_{i}$, $\mbox{\boldmath$\Delta$}_{7,i}$, and
${\bm{H}}_{i}$ are correlation matrices, and the simplex constraint ensures
${\bm{R}}_{i}$ is as well. In general, the variance component parameterization
${\bm{R}}_{i}(\mbox{\boldmath$\sigma$}^{2})=\sum_{j=1}^{m}\sigma_{j}^{2}{\bm{R}}_{ij}$
with the simplex constraint in force is reasonable. In this setting the GEE
update of $\beta$ given $\mbox{\boldmath$\sigma$}^{2}$ solves the equation

$\displaystyle\sum_{i=1}^{n}{\bm{X}}_{i}^{T}{\bm{A}}_{i}^{-1/2}(\mbox{\boldmath$\beta$}){\bm{R}}_{i}^{-1}(\mbox{\boldmath$\sigma$}^{2}){\bm{A}}^{-1/2}(\mbox{\boldmath$\beta$})[{\bm{y}}_{i}-\mbox{\boldmath$\mu$}_{i}(\mbox{\boldmath$\beta$})]$
$\displaystyle=$ $\displaystyle{\mathbf{0}}_{p}.$

This is just the classical GEE update. The difficulty lies in updating
$\mbox{\boldmath$\sigma$}^{2}$ given $\beta$. We propose minimizing the sum
$\sum_{i=1}^{n}\psi({\bm{R}}_{i})$, where $\psi(t)$ is a scalar convex loss
function. Example loss functions include the Mahalanobis criterion

$\displaystyle\sum_{i=1}^{n}\psi_{1}({\bm{R}}_{i})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}({\bm{y}}_{i}-\widehat{}\mbox{\boldmath$\mu$}_{i})^{T}\widehat{\bm{A}}_{i}^{-1/2}{\bm{R}}_{i}^{-1}(\mbox{\boldmath$\sigma$}^{2})\widehat{\bm{A}}_{i}^{-1/2}({\bm{y}}_{i}-\widehat{}\mbox{\boldmath$\mu$}_{i})$

and the sum of squared Frobenius distances

$\displaystyle\sum_{i=1}^{n}\psi_{2}({\bm{R}}_{i})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\mathrm{tr}\left[({\bm{y}}_{i}-\widehat{}\mbox{\boldmath$\mu$}_{i})({\bm{y}}_{i}-\widehat{}\mbox{\boldmath$\mu$}_{i})^{T}-\widehat{\bm{A}}_{i}^{1/2}{\bm{R}}_{i}(\mbox{\boldmath$\sigma$}^{2})\widehat{\bm{A}}_{i}^{1/2}\right]^{2}.$

The convexity of $\psi_{1}(t)$ entails a minorization similar to the
minorization (5). Minimizing the surrogate then yields the MM update

$\displaystyle\sigma_{j}^{2(t+1)}=\frac{\sum_{i=1}^{n}({\bm{y}}_{i}-\widehat{}\mbox{\boldmath$\mu$}_{i})^{T}\widehat{\bm{A}}_{i}^{-1/2}{\bm{R}}_{i}^{-1}(\mbox{\boldmath$\sigma$}^{2(t)}){\bm{R}}_{ij}{\bm{R}}_{i}^{-1}(\mbox{\boldmath$\sigma$}^{2(t)})\widehat{\bm{A}}_{i}^{-1/2}({\bm{y}}_{i}-\widehat{}\mbox{\boldmath$\mu$}_{i})}{\sum_{j=1}^{m}\sum_{i=1}^{n}({\bm{y}}_{i}-\widehat{}\mbox{\boldmath$\mu$}_{i})^{T}\widehat{\bm{A}}_{i}^{-1/2}{\bm{R}}_{i}^{-1}(\mbox{\boldmath$\sigma$}^{2(t)}){\bm{R}}_{ij}{\bm{R}}_{i}^{-1}(\mbox{\boldmath$\sigma$}^{2(t)})\widehat{\bm{A}}_{i}^{-1/2}({\bm{y}}_{i}-\widehat{}\mbox{\boldmath$\mu$}_{i})}.$

Under $\psi_{2}(t)$ the MM update boils down to projection onto the simplex.
Further exploration of these ideas probably deserves another paper and will be
omitted here for the sake of brevity.

##  7 A numerical example

Quantitative trait loci (QTL) mapping aims to identify genes associated with a
quantitative trait. Current sequencing technology measures millions of genetic
markers in study subjects. Traditional single-marker tests suffer from low
power due to the low frequency of many markers and the corrections needed for
multiple hypothesis testing. Region-based association tests are a powerful
alternative for analyzing next generation sequencing data with abundant rare
variants.

Suppose ${\bm{y}}$ is a $n\times 1$ vector of quantitative trait measurements
on $n$ people, ${\bm{X}}$ is an $n\times p$ predictor matrix (incorporating
predictors such as sex, smoking history, and principal components for ethnic
admixture), and ${\bm{G}}$ is an $n\times m$ genotype matrix of $m$ genetic
variants in a pre-defined region. The linear mixed model assumes

$\displaystyle{\bm{Y}}$ $\displaystyle=$
$\displaystyle{\bm{X}}\mbox{\boldmath$\beta$}+{\bm{G}}\mbox{\boldmath$\gamma$}+\mbox{\boldmath$\epsilon$},\quad\mbox{\boldmath$\gamma$}\sim
N({\mathbf{0}},\sigma_{g}^{2}{\bm{I}}),\quad\mbox{\boldmath$\epsilon$}\sim
N({\mathbf{0}},\sigma_{e}^{2}{\bm{I}}_{n}),$

where $\beta$ are fixed effects, $\gamma$ are random genetic effects, and
$\sigma_{g}^{2}$ and $\sigma_{e}^{2}$ are variance components for the genetic
and environmental effects, respectively. Thus, the phenotype vector ${\bm{Y}}$
has covariance
$\sigma_{g}^{2}{\bm{G}}{\bm{G}}^{T}+\sigma_{e}^{2}{\bm{I}}_{n}$, where
${\bm{G}}{\bm{G}}^{T}$ is the kernel matrix capturing the overall effect of
the $m$ variants. Current approaches test the null hypothesis
$\sigma_{g}^{2}=0$ for each region separately and then adjust for multiple
testing (Lee et al.,, 2014). Hence, we consider the joint model

$\displaystyle{\bm{y}}$ $\displaystyle=$
$\displaystyle{\bm{X}}\mbox{\boldmath$\beta$}+s_{1}^{-1/2}{\bm{G}}_{1}\mbox{\boldmath$\gamma$}_{1}+\cdots+s_{m}^{-1/2}{\bm{G}}_{m}\mbox{\boldmath$\gamma$}_{m}+\mbox{\boldmath$\epsilon$},$
$\displaystyle\mbox{\boldmath$\gamma$}_{i}$ $\displaystyle\sim$ $\displaystyle
N({\mathbf{0}},\sigma_{i}^{2}{\bm{I}}),\quad\mbox{\boldmath$\epsilon$}\sim
N({\mathbf{0}},\sigma_{e}^{2}{\bm{I}}_{n})$

and select the variance components $\sigma_{i}^{2}$ via the penalization (19).
Here $s_{i}$ is the number of variants in region $i$, and the weights
$s_{i}^{-1/2}$ put all variance components on the same scale.

We illustrate this approach using the COPDGene exome sequencing study
(http://www.copdgene.org/) (Regan et al.,, 2010). After quality control, $399$
individuals and 646,125 genetic variants remain for analysis. Genetic variants
are grouped into 16,619 genes to expose those genes associated with the
complex trait height. We include age, sex, and the top 3 principal components
in the mean effects. Because the number of genes vastly exceeds the sample
size $n=399$, we first pare the 16,619 genes down to 200 genes according to
their marginal likelihood ratio test p-values and then carry out penalized
estimation of the 200 variance components in the joint model (19). This is
similar to the sure independence screening strategy for selecting mean effects
(Fan and Lv,, 2008). Genes are ranked according to the order they appear in
the lasso solution path. Table 4 lists the top 10 genes together with their
marginal LRT p-values. Figure 2 displays the corresponding segment of the
lasso solution path. It is noteworthy that the ranking of genes by penalized
estimation differs from the ranking according to marginal p-values. The same
phenomenon occurs in selection of highly correlated mean predictors. This
penalization approach for selecting variance components warrants further
theoretical study. It is reassuring that the simple MM algorithm scales to
high-dimensional problems.

Lasso Rank | Gene | Marginal P-value | # Variants
---|---|---|---
1 | DOLPP1 | $2.35\times 10^{-6}$ | 2
2 | C9orf21 | $3.70\times 10^{-5}$ | 4
3 | PLS1 | $2.29\times 10^{-3}$ | 5
4 | ATP5D | $6.80\times 10^{-7}$ | 3
5 | ADCY4 | $1.01\times 10^{-3}$ | 11
6 | SLC22A25 | $3.95\times 10^{-3}$ | 14
7 | RCSD1 | $9.04\times 10^{-4}$ | 4
8 | PCDH7 | $1.20\times 10^{-4}$ | 7
9 | AVIL | $8.34\times 10^{-4}$ | 11
10 | AHR | $1.14\times 10^{-3}$ | 7
Table 4: Top 10 genes selected by the lasso penalized variance component model
(19) in an association study of 200 genes and the complex trait height. Figure
2: Solution path of the lasso penalized variance component model (19) in an
association study of 200 genes and the complex trait height.

##  8 Discussion

The current paper leverages the MM principle to design powerful and versatile
algorithms for variance components estimation. The MM algorithms derived are
notable for their simplicity, generality, numerical efficiency, and
theoretical guarantees. Both ordinary MLE and REML are apt to benefit. Other
extensions are possible. In nonlinear models (Bates and Watts,, 1988;
Lindstrom and Bates,, 1990), the mean response is a nonlinear function in the
fixed effects $\beta$. One can easily modify the MM algorithms to update
$\beta$ by a few rounds of Gauss-Newton iteration. The variance components
updates remain unchanged.

One can also extend our MM algorithms to elliptically symmetric densities

$\displaystyle f({\bm{y}})$ $\displaystyle=$
$\displaystyle\frac{e^{-\frac{1}{2}\kappa(\delta^{2})}}{(2\pi)^{\frac{n}{2}}(\det\mbox{\boldmath$\Omega$})^{\frac{1}{2}}}$

defined for ${\bm{y}}\in\mathbb{R}^{n}$, where
$\delta^{2}=({\bm{y}}-\mbox{\boldmath$\mu$})^{T}\mbox{\boldmath$\Omega$}^{-1}({\bm{y}}-\mbox{\boldmath$\mu$})$
denotes the Mahalanobis distance between ${\bm{y}}$ and $\mu$. Here we assume
that the function $\kappa(s)$ is strictly increasing and strictly concave.
Examples of elliptically symmetric densities include the multivariate $t$,
slash, contaminated normal, power exponential, and stable families. Previous
work (Huber and Ronchetti,, 2009; Lange and Sinsheimer,, 1993) has focused on
using the MM principle to convert parameter estimation for these robust
families into parameter estimation under the multivariate normal. One can
chain the relevant majorization
$\kappa(s)\leq\kappa(s^{(t)})+\kappa^{\prime}(s^{(t)})(s-s^{(t)})$ with our
previous minorizations and simultaneously split variance components and pass
to the more benign setting of the multivariate normal.

## Acknowledgments

The research is partially supported by NSF grant DMS-1055210 and NIH grants
R01 HG006139, R01 GM53275, and R01 GM105785. The authors thank Michael Cho,
Dandi Qiao, and Edwin Silverman for their assistance in processing and
assessing COPDGene exome sequencing data. COPDGene is supported by NIH R01
HL089897 and R01 HL089856.

## References

  * Bates and Pinheiro, (1998) Bates, D. and Pinheiro, J. (1998).  Computational methods for multilevel models.  Technical Report Technical Memorandum BL0112140-980226-01TM, Bell Labs, Lucent Technologies, Murray Hill, NJ.
  * Bates and Watts, (1988) Bates, D. M. and Watts, D. G. (1988).  Nonlinear Regression Analysis and Its Applications.  Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. John Wiley & Sons, Inc., New York.
  * Borg and Groenen, (2005) Borg, I. and Groenen, P. J. (2005).  Modern Multidimensional Scaling: Theory and Applications.  Springer Science & Business Media.
  * Boyd and Vandenberghe, (2004) Boyd, S. and Vandenberghe, L. (2004).  Convex Optimization.  Cambridge University Press, Cambridge.
  * Callanan and Harville, (1991) Callanan, T. P. and Harville, D. A. (1991).  Some new algorithms for computing restricted maximum likelihood estimates of variance components.  J. Statist. Comput. Simulation, 38(1-4):239–259.
  * De Leeuw, (1994) De Leeuw, J. (1994).  Block-relaxation algorithms in statistics.  In Information Systems and Data Analysis, pages 308–324. Springer.
  * Demidenko, (2013) Demidenko, E. (2013).  Mixed Models.  Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, second edition.  Theory and applications with R.
  * Demidenko and Massam, (1999) Demidenko, E. and Massam, H. (1999).  On the existence of the maximum likelihood estimate in variance components models.  Sankhyā Ser. A, 61(3):431–443.
  * Dempster et al., (1977) Dempster, A., Laird, N., and Rubin, D. (1977).  Maximum likelihood from incomplete data via the EM algorithm.  Journal of the Royal Statistical Soceity Series B., 39(1-38).
  * Fan and Lv, (2008) Fan, J. and Lv, J. (2008).  Sure independence screening for ultrahigh dimensional feature space (with discussion).  J. R. Statist. Soc. B, 70:849–911.
  * Fitzmaurice et al., (2011) Fitzmaurice, G. M., Laird, N. M., and Ware, J. H. (2011).  Applied Longitudinal Analysis.  Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, second edition.
  * Golub and Van Loan, (1996) Golub, G. H. and Van Loan, C. F. (1996).  Matrix Computations.  Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, third edition.
  * Grządziel and Michalski, (2014) Grządziel, M. and Michalski, A. (2014).  A note on the existence of the maximum likelihood estimate in variance components models.  Discuss. Math. Probab. Stat., 34(1-2):159–167.
  * Gupta and Nagar, (1999) Gupta, A. and Nagar, D. (1999).  Matrix Variate Distributions.  Monographs and Surveys in Pure and Applied Mathematics. Taylor & Francis.
  * Hartley and Rao, (1967) Hartley, H. O. and Rao, J. N. K. (1967).  Maximum-likelihood estimation for the mixed analysis of variance model.  Biometrika, 54:93–108.
  * Harville and Callanan, (1990) Harville, D. and Callanan, T. (1990).  Computational aspects of likelihood-based inference for variance components.  In Gianola, D. and Hammond, K., editors, Advances in Statistical Methods for Genetic Improvement of Livestock, volume 18 of Advanced Series in Agricultural Sciences, pages 136–176. Springer Berlin Heidelberg.
  * Harville, (1977) Harville, D. A. (1977).  Maximum likelihood approaches to variance component estimation and to related problems.  J. Amer. Statist. Assoc., 72(358):320–340.  With a comment by J. N. K. Rao and a reply by the author.
  * Heiser, (1995) Heiser, W. J. (1995).  Convergent computation by iterative majorization: theory and applications in multidimensional data analysis.  Recent Advances in Descriptive Multivariate Analysis, pages 157–189.
  * Horn and Johnson, (1985) Horn, R. A. and Johnson, C. R. (1985).  Matrix Analysis.  Cambridge University Press, Cambridge.
  * Huber and Ronchetti, (2009) Huber, P. J. and Ronchetti, E. M. (2009).  Robust Statistics.  Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, second edition.
  * Hunter, (2004) Hunter, D. R. (2004).  MM algorithms for generalized Bradley-Terry models.  Ann. Statist., 32(1):384–406.
  * Hunter and Lange, (2004) Hunter, D. R. and Lange, K. (2004).  A tutorial on MM algorithms.  Amer. Statist., 58(1):30–37.
  * Hunter and Li, (2005) Hunter, D. R. and Li, R. (2005).  Variable selection using MM algorithms.  Ann. Statist., 33(4):1617–1642.
  * Jeon, (2012) Jeon, M. (2012).  Estimation of Complex Generalized Linear Mixed Models for Measurement and Growth.  PhD thesis, University of California, Berkeley.
  * Khuri et al., (1998) Khuri, A. I., Mathew, T., and Sinha, B. K. (1998).  Statistical Tests for Mixed Linear Models.  Wiley Series in Probability and Statistics: Applied Probability and Statistics. John Wiley & Sons, Inc., New York.  A Wiley-Interscience Publication.
  * Kiers, (2002) Kiers, H. A. (2002).  Setting up alternating least squares and iterative majorization algorithms for solving various matrix optimization problems.  Computational Statistics & Data Analysis, 41(1):157–170.
  * Laird et al., (1987) Laird, N., Lange, N., and Stram, D. (1987).  Maximum likelihood computations with repeated measures: application of the EM algorithm.  J. Amer. Statist. Assoc., 82(397):97–105.
  * Laird and Ware, (1982) Laird, N. M. and Ware, J. H. (1982).  Random-effects models for longitudinal data.  Biometrics, 38(4):963–974.
  * Lange, (2002) Lange, K. (2002).  Mathematical and Statistical Methods for Genetic Analysis.  Statistics for Biology and Health. Springer-Verlag, New York, second edition.
  * Lange, (2010) Lange, K. (2010).  Numerical Analysis for Statisticians.  Statistics and Computing. Springer, New York, second edition.
  * Lange et al., (2000) Lange, K., Hunter, D. R., and Yang, I. (2000).  Optimization transfer using surrogate objective functions.  J. Comput. Graph. Statist., 9(1):1–59.  With discussion, and a rejoinder by Hunter and Lange.
  * Lange et al., (2013) Lange, K., Papp, J., Sinsheimer, J., Sripracha, R., Zhou, H., and Sobel, E. (2013).  Mendel: the Swiss army knife of genetic analysis programs.  Bioinformatics, 29:1568–1570.
  * Lange et al., (2005) Lange, K., Sinsheimer, J., and Sobel, E. (2005).  Association testing with Mendel.  Genetic Epidemiology, 29:36–50.
  * Lange and Sinsheimer, (1993) Lange, K. and Sinsheimer, J. S. (1993).  Normal/independent distributions and their applications in robust regression.  Journal of Computational and Graphical Statistics, 2:175–198.
  * Lange and Zhou, (2014) Lange, K. and Zhou, H. (2014).  MM algorithms for geometric and signomial programming.  Mathematical Programming Series A, 143:339–356.
  * Lee et al., (2014) Lee, S., Abecasis, G., Boehnke, M., and Lin, X. (2014).  Rare-variant association analysis: Study designs and statistical tests.  The American Journal of Human Genetics, 95(1):5 – 23.
  * Liang and Zeger, (1986) Liang, K. and Zeger, S. (1986).  Longitudinal data analysis using general linear models.  Biometrika, 73:13–22.
  * Lindstrom and Bates, (1988) Lindstrom, M. J. and Bates, D. M. (1988).  Newton-Raphson and EM algorithms for linear mixed-effects models for repeated-measures data.  J. Amer. Statist. Assoc., 83(404):1014–1022.
  * Lindstrom and Bates, (1990) Lindstrom, M. J. and Bates, D. M. (1990).  Nonlinear mixed effects models for repeated measures data.  Biometrics, 46(3):673–687.
  * McLachlan and Krishnan, (2008) McLachlan, G. J. and Krishnan, T. (2008).  The EM Algorithm and Extensions.  Wiley Series in Probability and Statistics. Wiley-Interscience [John Wiley & Sons], Hoboken, NJ, second edition.
  * Meng and Rubin, (1991) Meng, X.-L. and Rubin, D. B. (1991).  Using EM to obtain asymptotic variance-covariance matrices: the SEM algorithm.  Journal of the American Statistical Association, 86(416):899–909.
  * Pinheiro and Bates, (1996) Pinheiro, J. and Bates, D. (1996).  Unconstrained parametrizations for variance-covariance matrices.  Statistics and Computing, 6(3):289–296.
  * Rao, (1973) Rao, C. R. (1973).  Linear Statistical Inference and its Applications, 2nd ed.  John Wiley & Sons.
  * Rao and Kleffe, (1988) Rao, C. R. and Kleffe, J. (1988).  Estimation of Variance Components and Applications, volume 3 of North-Holland Series in Statistics and Probability.  North-Holland Publishing Co., Amsterdam.
  * Rao, (1997) Rao, P. S. R. S. (1997).  Variance Components Estimation, volume 78 of Monographs on Statistics and Applied Probability.  Chapman & Hall, London.  Mixed models, methodologies and applications.
  * Regan et al., (2010) Regan, E. A., Hokanson, J. E., Murphy, J. R., Make, B., Lynch, D. A., Beaty, T. H., Curran-Everett, D., Silverman, E. K., and Crapo, J. D. (2010).  Genetic epidemiology of COPD (COPDGene) study designs.  COPD, 7:32–43.
  * Schafer and Yucel, (2002) Schafer, J. L. and Yucel, R. M. (2002).  Computational strategies for multivariate linear mixed-effects models with missing values.  J. Comput. Graph. Statist., 11(2):437–457.
  * Schelldorfer et al., (2011) Schelldorfer, J., Bühlmann, P., and van de Geer, S. (2011).  Estimation for high-dimensional linear mixed-effects models using $\ell_{1}$-penalization.  Scand. J. Stat., 38(2):197–214.
  * Schur, (1911) Schur, J. (1911).  Bemerkungen zur Theorie der beschränkten Bilinearformen mit unendlich vielen Veränderlichen.  J. Reine Angew. Math., (140):1–28.
  * Searle et al., (1992) Searle, S. R., Casella, G., and McCulloch, C. E. (1992).  Variance Components.  Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. John Wiley & Sons, Inc., New York.  A Wiley-Interscience Publication.
  * Varadhan and Roland, (2008) Varadhan, R. and Roland, C. (2008).  Simple and globally convergent methods for accelerating the convergence of any EM algorithm.  Scand. J. Statist., 35(2):335–353.
  * Verbeke and Molenberghs, (2000) Verbeke, G. and Molenberghs, G. (2000).  Linear Mixed Models for Longitudinal Data.  Springer Series in Statistics. Springer-Verlag, New York.
  * Weiss, (2005) Weiss, R. E. (2005).  Modeling Longitudinal Data.  Springer Texts in Statistics. Springer, New York.
  * Wu and Lange, (2010) Wu, T. T. and Lange, K. (2010).  The MM alternative to EM.  Statistical Science, 25:492–505.
  * Yu, (2010) Yu, Y. (2010).  Monotonic convergence of a general algorithm for computing optimal designs.  Ann. Statist., 38(3):1593–1606.
  * Zhou et al., (2011) Zhou, H., Alexander, D., and Lange, K. (2011).  A quasi-Newton acceleration for high-dimensional optimization algorithms.  Statistics and Computing, 21:261–273.
  * Zhou and Lange, (2010) Zhou, H. and Lange, K. (2010).  MM algorithms for some discrete multivariate distributions.  Journal of Computational and Graphical Statistics, 19:645–665.