# Low-Rank Matrix Completion for Mobile Edge Caching in Fog-RAN via Riemannian
Optimization

Kai Yang∗, Yuanming Shi∗, and Zhi Ding†, _Fellow, IEEE_
This work is supported by Shanghai Sailing Program No. 16YF1407700. ∗School of
Information Science and Technology, ShanghaiTech University, Shanghai, China
†Dept. of ECE, University of California, Davis, California 95616, USA
E-mail: {yangkai<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

The upcoming big data era is likely to demand tremendous computation and
storage resources for communications. By pushing computation and storage to
network edges, fog radio access networks (Fog-RAN) can effectively increase
network throughput and reduce transmission latency. Furthermore, we can
exploit the benefits of cache enabled architecture in Fog-RAN to deliver
contents with low latency. Radio access units (RAUs) need content delivery
from fog servers through wireline links whereas multiple mobile devices
acquire contents from RAUs wirelessly. This work proposes a unified low-rank
matrix completion (LRMC) approach to solving the content delivery problem in
both wireline and wireless parts of Fog-RAN. To attain a low caching latency,
we present a high precision approach with Riemannian trust-region method to
solve the challenging LRMC problem by exploiting the quotient manifold
geometry of fixed-rank matrices. Numerical results show that the new approach
has a faster convergence rate, is able to achieve optimal results, and
outperforms other state-of-art algorithms.

###### Index Terms:

Edge caching, Fog-RAN, low-rank matrix completion, Riemannian optimization.

##  I Introduction

The astounding growth of smart mobile device popularity, coupled with new
types of wireless services and applications, such as Internet of Things and
mobile Cyber-Physical applications, has helped usher in the era of wireless
big data [1]. Futuristic high speed networks such as Tactile Internet [2] will
enable numerous new services and allow for new experiences, e.g., autonomous
driving, healthcare, and virtual reality. However, the rapidly growing data
and diversified services will put great strain on both storage and computation
[3] in wireless communication systems. This calls for a paradigm shift for
wireless access technologies to provide services with the stringent
requirements of ultra-low latency, high data rate, as well as massive
connectivity [4].

By leveraging mobile edge/fog computing concepts [5], the recent proposal of
fog radio access network (Fog-RAN) [6, 7] represents a disruptive wireless
network architecture to accommodate the upcoming diversified services with
ultra-low latency and high data rates. This is achieved by pushing the
computation and storage resources to network edges, thereby enabling cloud
computation within wireless networks as shown in Fig. 1.

Specifically, in Fog-RAN, each radio access unit (RAU) hosts storage and
computation entities, thereby pushing data and its processing closer to end
users. This cache enabled network architecture promises to reduce network
congestion for wireline communication scenarios [8] and enhance interference
coordination in wireless communication networks [9]. By further pushing the
content and computation resources to smart mobile devices (MDs), end-to-end
latency and network capacity can be improved, e.g., in cache-aided wireless
device-to-device (D2D) networks [10]. In particular, for wearable computing
applications, offloading computation-intensive applications to proximal smart
mobile devices or RAUs can significantly improve roundtrip delay and energy
efficiency [11].

In this paper, we mainly focus on the content-centric communication and edge
storage aspects of Fog-RAN. The fundamental goal is to leverage the network
caching capability for efficient content delivery, thereby reducing end-to-end
network latency. To understand this problem, we derive unified interference
alignment conditions for content delivery for wireline links between fog
server and RAUs, and also for wireless links between RAUs and smart mobile
devices. Surprisingly, caching problem with fixed side information is
equivalent to the index coding [12] problem, which is equivalent to many other
important problems, including topological interference management (TIM) and
network coding [13]. Unfortunately, index coding is an NP-hard problem and
only a few special cases can been solved efficiently.

In this paper, we present a low-rank matrix completion (LRMC) approach [14] to
minimize message delivery latency by exploiting side information provided by
the caching capability. The low-rank modeling framework has wide applications
in machine learning, computational big data analytics, high-dimensional
statistics. Furthermore, the LRMC approach has recently been exploited to
solve the wireless topological interference management problem [15] and index
coding problem [16]. However, these existing works mainly focus on unicast, in
which each message is desired by exactly one destination, whereas in the
caching problem each message may be desired by multiple destinations. Despite
of the NP-hardness of the resulting LRMC problem, we will reveal that the
modeling framework provides algorithmic opportunities.

Although our problem falls within the category of low-rank matrix completion,
which has stimulated a flurry of recent research activities [14, 17, 18, 19],
most prior works are not applicable in our LRMC problem because of their
special affine constraint and measurement graph structures. Specifically, the
well-known convex relaxation approach by replacing the rank function to the
nuclear norm [14] will always yield a full-rank solution in the unicast case,
which turns out to be invalid here. Furthermore, alternating minimization
approaches [17, 18] also exhibit their own limitations such as slow
convergence rate, and impractical assumptions (e.g., incoherence of the
original matrix), and available (estimated) ground-truth rank. One exception
is the conjugate gradient based Riemannian pursuit algorithm [19], which
explores the embedded manifold structures of fixed-rank matrices to design
rank estimation strategies. However, the first-order method has slow
convergence and is sensitive to initial points.

In this paper, we develop a Riemannian trust-region optimization framework
[20] to provide high accuracy solutions and faster convergence rate for
finding the minimal rank in the LRMC problem. This second-order algorithm is
more robust to initial point choices, and can obtain reliable solutions. We
first solve the fixed-rank matrices, by exploiting the quotient manifold
geometry of the search space of fixed-rank matrices [20]. We then develop an
efficient rank increase mechanism to estimate the minimum rank while
satisfying the affine constraint. This is achieved by exploiting the closure
of fixed-rank matrices. Numerical results will show that our proposed
algorithm can recover existing optimal DoF results in [12] and outperform
state-of-art algorithms.

##  II System Model and Problem Formulation

Figure 1: The network architecture of Fog-RAN, in which, the data and its
processing are pushed to the edge of networks, including RAUs and smart
devices. The benefits of caching allow for wireline communication links
between the fog server and RAUs (e.g., massive MIMO, macro base station (BS),
and femtocell access point (AP)), as well as the wireless communication links
between RAUs and small mobile devices.

Consider a cache-aided Fog-RAN as shown in Fig. 1, where both RAUs and smart
devices have caching capabilities. To evaluate the benefits of caching in
RAUs, we consider the wireline communication network between one fog data
center and multiple RAUs. To further investigate the advantages of caching in
smart devices, we consider a $K$-user interference channel with cache enabled
smart devices. We shall present that the message delivery problems are
equivalent to the index coding problem with fixed side information [21].

###  II-A The Benefits of Caching for Wireline Communications

Consider a wireline communication network between the fog data center and
cache-enabled RAUs as shown in Fig. 2(a). We assume that the file library
consists of $K$ messages $\\{W_{1},\dots,W_{K}\\}$, each with entropy of $F$
bits. Let $\mathcal{M}=\\{1,\dots,K\\}$ be the index set of the $K$
independent messages. Define the set $\mathcal{V}_{k}\subset\mathcal{M}$ as
the index set of the messages cached at the $k$-th RAU. We consider the
general multiple groupcast scenario, where each message may be desired by
multiple destination nodes [21]. Let $\mathcal{M}_{k}\subseteq\mathcal{M}$ be
the index set of messages desired by RAU $k$. We assume that the available
messages will not be further desired by the corresponding destination node,
i.e., $\mathcal{M}_{k}\cap\mathcal{V}_{k}=\emptyset$. In this paper, we are
interested in constructing the vector linear coding schemes over real field to
maximize the message delivery rate.

Specifically, the transmitted symbol sequence ${\bm{z}}\in\mathbb{R}^{N}$ over
$N$ channel uses is given by

${\bm{z}}=\sum_{m\in\mathcal{M}}{\bm{V}}_{m}{\bm{s}}_{m}.$ (1)

Here, message $W_{m}$ is split into $Q_{m}$ scalar data streams, denoted as
${\bm{s}}_{m}=[s_{i}(W_{m})]_{i=1}^{Q_{m}}\in\mathbb{R}^{Q_{m}\times 1}$, and
each of which carries one symbol from $\mathbb{R}$ and is transmitted along
the column vectors of the precoding matrix ${\bm{V}}_{m}\in\mathbb{R}^{N\times
Q_{m}}$.

Let ${\bm{U}}_{i,k}\in\mathbb{R}^{Q_{i}\times N}$ with $i\in\mathcal{M}_{k}$
be the receiver combining matrix at RAU $k$ for the desired message $W_{i}$.
We have the following decoding operation for message $W_{i}$ at RAU $k$:

$\hat{\bm{s}}_{i}=\left({\bm{U}}_{i,k}\bm{V}_{i}\right)^{-1}\bm{U}_{i,k}\left({\bm{z}}-\sum_{j\in\mathcal{V}_{k}}{\bm{V}}_{j}{\bm{s}}_{j}\right),\forall
i\in\mathcal{M}_{k}.$ (2)

To accomplish the above decoding, we impose the following interference
alignment conditions [21, 15]:

$\displaystyle{\bm{U}}_{i,k}{\bm{V}}_{j}$ $\displaystyle=$
$\displaystyle{\bm{0}},\forall i\neq
j,i\in\mathcal{M}_{k},j\notin{\mathcal{V}}_{k}$ (3)
$\displaystyle\det({\bm{U}}_{i,k}{\bm{V}}_{i})$ $\displaystyle\neq$
$\displaystyle 0,\forall i\in\mathcal{M}_{k}.$ (4)

The first condition (3) is to align and cancel the interference, and the
second condition (4) is to preserve the desired signal.

Therefore, if (3) (4) can both be satisfied over $N$ channel uses, the message
delivery rate tuple
$\mathcal{R}=(\frac{Q_{1}}{N},\frac{Q_{2}}{N},\cdots,\frac{Q_{K}}{N})$ can be
achieved. In particular, if $Q_{1}=\cdots=Q_{K}=Q_{0}$, we say that the
symmetric data rate is $\frac{Q_{0}}{N}$. Thus the smaller $N$ we can achieve,
the lower latency we can realize. Note that, with $|\mathcal{M}_{i}|=1$ and
$\mathcal{M}_{i}\cap\mathcal{M}_{j}=\emptyset,\forall i,j$, we include the
multiple unicast scenario as a special case [15]. We also observe that the
caching problem with fixed side information (i.e., cached message) is
equivalent to the index coding problem [8].

###  II-B The Benefits of Caching for Wireless Communications

Recently, caching in content-centric wireless networks is a very popular topic
[22, 23]. We consider a wireless communication modeled as the $K$-user
interference channel as shown in Fig. 2 (b). Assume that each transmitter $k$
(i.e., RAU) has the message $W_{k}$ for transmission. Let the set
$\mathcal{V}_{k}\subset\mathcal{M}$ as the index set of the messages cached at
the $k$-th receiver (i.e., smart device). Let
$\mathcal{M}_{k}\subseteq\mathcal{M}$ with
$\mathcal{M}_{k}\cap\mathcal{V}_{k}=\emptyset$ be the index set of messages
desired by receiver $k$. With message splitting as in Section II-A, over the
$N$ channel uses, the received signal $\bm{y}_{i}\in\mathbb{C}^{N}$ at
receiver $k$ is given by

${\bm{y}}_{k}=\sum_{i=1}^{K}{\bm{H}}^{[ki]}{\bm{V}}_{i}{\bm{s}}_{i}+{\bm{n}}_{k},\forall
k,$ (5)

where ${\bm{n}}_{k}\sim\mathcal{CN}({\bm{0}},{\bm{I}}_{N})$ is the additive
isotropic white Gaussian noise, ${\bm{s}}_{i}\in\mathbb{C}^{Q_{i}\times 1}$
represents message $W_{i}$, ${\bm{H}}^{[ki]}={H}_{ki}{\bm{I}}_{N}$ is an
$N\times N$ diagonal matrix with $H_{ki}\in\mathbb{C}$ as the constant channel
coefficient between transmitter $i$ and receiver $k$ over $N$ channel uses in
the considered block.

Let ${\bm{U}}_{i,k}\in\mathbb{C}^{Q_{i}\times N}$ with $i\in\mathcal{M}_{k}$
be the receiver combining matrix at receiver $k$ for the desired message
$W_{i}$. We assume that the channel state information $H_{ki}$ with
$i\in\mathcal{V}_{k}$ is available at receiver $k$. Therefore, we can first
eliminate the undesired messages available in $\mathcal{V}_{k}$, resulting the
interference space as $\sum_{j\notin\mathcal{V}_{k},i\neq
j}{\bm{H}}^{[kj]}{\bm{V}}_{j}$ for the message $W_{i}$. In the regime of
asymptotically high SNR, to accomplish decoding, we impose the constraint that
the desired signal space ${\bm{H}}^{[ki]}{\bm{V}}_{i}$ is complementary to the
interference space, resulting the following interference alignment condition
(see [12, 15])

$\displaystyle{\bm{U}}_{i,k}\bm{H}^{[kj]}{\bm{V}}_{j}$ $\displaystyle=$
$\displaystyle{\bm{0}},\forall i\neq
j,i\in\mathcal{M}_{k},j\notin{\mathcal{V}}_{k}$ (6)
$\displaystyle\det({\bm{U}}_{i,k}\bm{H}^{[ki]}{\bm{V}}_{i})$
$\displaystyle\neq$ $\displaystyle 0,\forall i\in\mathcal{M}_{k}.$ (7)

Since ${\bm{H}}^{[ki]}=H_{ki}{\bf{I}}_{N}$ for the constant channel over the
$N$ channel uses, conditions (6) and (7) are equivalent to conditions (3) and
(4), respectively.

Therefore, if (6) (7) are both satisfied over $N$ channel uses, the parallel
interference-free channels can be obtained. Therefore, the DoF tuples
$\mathcal{D}=(\frac{Q_{1}}{N},\frac{Q_{2}}{N},\cdots,\frac{Q_{K}}{N})$ for the
message delivery over wireless interference channel. In particular, if
$Q_{1}=\cdots=Q_{K}=Q_{0}$, we say that the the symmetric DoF
$\frac{Q_{0}}{N}$ can be achieved for each message. As DoF provides a first-
order characterization for capacity, we shall minimize channel uses $N$ to
improve the message delivery data rate.

##  III Low-Rank Matrix Completion Approach for Mobile Edge Caching

In this section, we present a low-rank matrix completion modeling framework to
find the minimal number of channel uses such that the interference alignment
conditions (3) and (4) can be satisfied. Specifically, let
$\bm{X}_{ij,k}={\bm{U}}_{i,k}{\bm{V}}_{j}$. Without lose of generality, we set
${\bm{U}}_{i,k}{\bm{V}}_{i}=\bm{I}_{Q_{i}\times Q_{i}},\forall
i\in\mathcal{M}_{k}$. We thus can restrict $\bm{X}_{ij,k}$’s to the real field
without losing any performance in terms of achievable DoFs [15]. Therefore,
for $\bm{X}_{ij,k}\in\mathbb{R}^{Q_{i}\times Q_{j}}$, we have

$\displaystyle{\bm{X}}_{ij,k}=\left\\{\begin{array}[]{ll}\bm{0}&~{}~{}\forall
i\neq j,i\in\mathcal{M}_{k},j\notin{\mathcal{V}}_{k}\\\ \bm{I}&~{}~{}\forall
i=j,i\in\mathcal{M}_{k}\\\ \star&~{}~{}\textrm{otherwise},\end{array}\right.$
(11)

where $\star$ represents arbitrary $Q_{i}\times Q_{j}$ vectors.

Let $\bm{X}=[\bm{X}_{ij,k}]\in\mathbb{R}^{M\times Q}$ with
$M=\sum_{i=1}^{K}|\mathcal{M}_{i}|$ and $Q=\sum_{i=1}^{K}Q_{i}$. Based on
(11), the interference alignment conditions (3) and (4) can be rewritten as

$\displaystyle\mathcal{P}_{\Omega}(\bm{X})=\bm{J},$ (12)

where $\mathcal{P}_{\Omega}:\mathbb{R}^{M\times
Q}\rightarrow\mathbb{R}^{M\times Q}$ is the orthogonal projection operator
onto the subspace of matrices which vanish outside $\Omega$ such that the
$(i,j)$-th component of $\mathcal{P}_{\Omega}(\bm{X})$ equals to $X_{ij}$ if
$(i,j)\in\Omega$ and zero otherwise. Here, the constant $M\times Q$ matrix
$\bm{J}=[\bm{J}_{mj}]$ is given by

$\displaystyle{\bm{J}}_{mj}=\left\\{\begin{aligned}
\bm{I}&~{}~{}{\textrm{if}}~{}\sum_{i=1}^{k-1}|\mathcal{M}_{i}|+1\leq
m\leq\sum_{i=1}^{k}|\mathcal{M}_{i}|,j\in\mathcal{M}_{k}\\\
\bm{0}&~{}~{}\textrm{otherwise},\end{aligned}\right.~{}~{}$ (13)

where $\bm{I}$ is the $Q_{j}\times Q_{j}$ identity matrix. The set $\Omega$ is
constructed as
$\Omega=\\{\mathcal{S}_{m}\times\mathcal{S}_{j}|\sum_{i=1}^{k-1}|\mathcal{M}_{i}|+1\leq
m\leq\sum_{i=1}^{k}|\mathcal{M}_{i}|,~{}{\textrm{with}}~{}j\notin\mathcal{V}_{k}~{}\textrm{or}~{}j\in\mathcal{M}_{k}\\}$,
where
$\mathcal{S}_{i}=\\{1+\sum_{k=1}^{i-1}Q_{k},\dots,\sum_{k=1}^{i}Q_{k}\\}$.

###  III-A Low-Rank Matrix Completion Modeling Framework

Observe that the rank of matrix $\bm{X}\in\mathbb{R}^{M\times Q}$ equals $N$,
we propose to solve the following LRMC problem

$\displaystyle\mathop{\textrm{minimize }}_{\bm{X}}$
$\displaystyle\textrm{rank}(\bm{X})$ subject to
$\displaystyle\mathcal{P}_{\Omega}(\bm{X})=\bm{J}$ (14)

to maximize the message delivery rate while satisfying interference alignment
conditions (3) (4), thereby reducing latency.

An illustrative example of a caching problem with $5$ messages and $3$
destination nodes is given in Fig. 2.

Figure 2: Cache enabled Fog-RAN. (a) A wireline communication network between
fog data center and cache-enabled RAUs. (b) Wireless cache-enabled networks
between RAUs and cache-enabled smart devices. (c) The corresponding index
coding problem. Red lines represent the side information. For example,
receiver$~{}1$ has $\mathcal{V}_{1}=\\{W_{2},W_{5}\\}$ as side information.
(d) Associated low-rank matrix completion model. In this example no message
splitting is taken into consideration, therefore $Q_{i}=1$. Gray cells
represent elements to be arbitrary values. For example, $(3,1)$-th element of
$\bm{X}$ must be $0$ since $W_{1}$ is not in the side information of
receiver$~{}3$.

In this setting, we have $\mathcal{V}_{1}=\\{2,5\\}$,
$\mathcal{V}_{2}=\\{1,5\\}$, $\mathcal{V}_{3}=\\{2,4\\}$,
$\mathcal{V}_{4}=\\{2,3\\}$, $\mathcal{V}_{5}=\\{1,3,4\\}$. Then all elements
of the set $\Omega$ are $\\{(1,1),(1,3),(1,4),(2,2),(2,3),(2,4)$,
$(3,1),(3,3),(3,5),(4,1),(4,4),(4,5),(5,2),(5,5)\\}$.

###  III-B Problem Analysis

Although LRMC problem is NP-hard generally, there are several popular methods
to yield a feasible solution. The most common approach is to relax the rank
function to nuclear norm $\bm{X}_{*}$ [14]. But obviously in the unicast case,
because $\bm{X}_{*}\geq|\textrm{Tr}(\bm{X})|=Q$ always holds, we would always
end up with a full rank solution. Another common way to solve the LRMC problem
applies alternating minimization. In this algorithm, under a fixed rank $r$ we
can factorize $\bm{X}$ as $\bm{X}=\bm{U}\bm{V}^{T}$. Next we optimize
$\bm{U}\in\mathbb{R}^{M\times r}$ and $\bm{V}\in\mathbb{R}^{Q\times r}$
alternatively. Unfortunately, this method has a very slow convergence rate and
the rank assumption is impractical since the rank is unknown for our problem.
Additionally, the algorithm calls for the incoherence property of matrix
$\bm{X}$, which is also invalid. Interestingly, the authors of [19] presented
a Riemannian pursuit algorithm based on conjugate gradient method, termed as
EmbG. EmbG takes a rank increase strategy by exploiting the embedded manifold
structures of fixed rank matrices. As a first-order algorithm, the convergence
of EmbG is slow, and is sensitive to initial points. Therefore, we hope to
find a method with higher accuracy, faster convergence rate and more robust to
initial points.

##  IV Riemannian Trust-Region Method for Low-Rank Matrix Completion on
Quotient Manifolds

In this section, we shall present a Riemannian trust-region framework which is
a second-order algorithm for faster convergence rate and higher resolution. We
first recast the LRMC problem as a set of fixed-rank subproblems and present
the corresponding Riemannian trust region approach. Then we develop our rank-
increasing strategy by exploiting the closure of fixed-rank matrices to ensure
that the objective function monotonically decreases.

###  IV-A Fixed-Rank Riemannian Trust-Region for Matrix Completion

By defining the cost function with Frobenius norm

$\displaystyle
f(\bm{X})=\frac{1}{2}\|\mathcal{P}_{\Omega}(\bm{X})-{\bm{J}}\|_{F}^{2},$ (15)

problem (III-A) can be reformulated as a sequence of subproblems $\mathscr{P}$
given the fixed rank $r$:

$\displaystyle\mathscr{P}:\quad\mathop{\textrm{minimize}}_{\bm{X}\in\mathcal{M}_{r}}\frac{1}{2}\|\mathcal{P}_{\Omega}(\bm{X})-{\bm{J}}\|_{F}^{2}.$
(16)

With a little abuse of notation here,
$\mathcal{M}_{r}:=\\{\bm{X}\in\mathbb{R}^{M\times
Q}:~{}\textrm{rank}(\bm{X})=r\\}$ is a smooth ($C^{\infty}$) manifold.

The matrix $\bm{X}$ can be represented by a SVD-type factorization of well
defined manifold
$\bm{X}=\bm{U}\bm{\Sigma}\bm{V}^{T},~{}\bm{U}\in\textrm{St}(r,M),\bm{\Sigma}\in\textrm{GL}(r),\bm{V}\in\textrm{St}(r,Q)$.
Here $\textrm{St}(r,M)=\\{\bm{Y}\in\mathbb{R}^{M\times
r}:\bm{Y}^{T}\bm{Y}=\bm{I}_{r}\\}$ is the compact Stiefel manifold of
orthonormal $M\times r$ matrices, and
$\textrm{GL}(r)=\\{\bm{Y}\in\mathbb{R}^{r\times r}:\textrm{rank}(\bm{Y})=r\\}$
is the non-compact Stiefel manifold of all $r\times r$ nonsingular matrices.
Thus $\bm{X}=(\bm{U},\bm{\Sigma},\bm{V})$ and computational space is
$\mathcal{M}_{r}:=\textrm{St}(r,M)\times\textrm{GL}(r)\times\textrm{St}(r,Q)$,
on which our algorithm is built. Then the quotient space is
$\mathcal{M}_{r}/\sim:=\mathcal{M}_{r}/(\mathcal{O}(r)\times\mathcal{O}(r))$,
$\mathcal{O}(r)=\\{\bm{Y}\in\mathbb{R}^{r\times
r}:\bm{Y}^{T}\bm{Y}=\bm{I}_{r}\\}$ denote the $r$-order orthogonal group.

Our algorithm is based on trust region method on Riemannian manifold. The
trust-region subproblem in the quotient manifold $\mathcal{M}_{r}/\sim$ is
horizontally lifted to the horizontal space $\mathcal{H}_{\bm{X}}$, and
formulated as

$\displaystyle\mathop{\rm{minimize}}_{{\bm{\xi}}_{\bm{X}}\in\mathcal{H}_{\bm{X}}\mathcal{M}_{r}}$
$\displaystyle{m({\bm{\xi}}_{\bm{X}})}$ $\displaystyle{\rm{subject~{}to}}$
$\displaystyle
g_{\bm{X}}({\bm{\xi}}_{\bm{X}},{\bm{\xi}}_{\bm{X}})\leq\delta^{2}.$ (17)

The trust-region radius is $\delta$ and the cost function is given by

$\displaystyle m({\bm{\xi}}_{\bm{X}})=$ $\displaystyle
f(\bm{X})+g_{\bm{X}}({\bm{\xi}}_{\bm{X}},\textrm{grad}_{\bm{X}}f)+$ (18)
$\displaystyle\frac{1}{2}g_{\bm{X}}({\bm{\xi}}_{\bm{X}},\textrm{Hess}_{\bm{X}}f\left[\bm{\xi}_{\bm{X}}\right]).$

To develop an effective algorithm, we need to determine

  * •

Riemannian metric: $g_{\bm{X}}$

  * •

Riemannian gradient: $\textrm{grad}_{\bm{X}}f$

  * •

Riemannian Hessian: $\textrm{Hess}_{\bm{X}}f[\bm{\xi}_{\bm{X}}]$

$\textrm{grad}_{\bm{X}}f$ and $\textrm{Hess}_{\bm{X}}f$ are the generalization
of Euclidean gradient and Hessian to the Riemannian manifold [20].

1) Riemannian metric: The Riemannian metric must have identical matrix
representation along the equivalent class $[\bm{X}]$. By taking second-order
derivative of $f(\bm{X})$, namely $\bm{\mathcal{L}}_{\bm{X}\bm{X}}(\bm{X})$,
we can induce its following Riemannian metric from the approximation [15]

$\displaystyle g_{\bm{X}}(\bm{\xi}_{\bm{X}},\bm{\zeta}_{\bm{X}})\approx$
$\displaystyle\langle\bm{\xi}_{U},\bm{\zeta}_{U}\bm{\Sigma}\bm{\Sigma}^{T}\rangle+\langle\bm{\xi}_{\Sigma},\bm{\zeta}_{\Sigma}\rangle$
(19)
$\displaystyle+\langle\bm{\xi}_{V},\bm{\zeta}_{V}\bm{\Sigma}^{T}\bm{\Sigma}\rangle,$

in which
$\bm{\xi}_{\bm{X}}=(\bm{\xi}_{U},\bm{\xi}_{\Sigma},\bm{\xi}_{V}),\bm{\zeta}_{\bm{X}}=(\bm{\zeta}_{U},\bm{\zeta}_{\Sigma},\bm{\zeta}_{V})\in
T_{\bm{X}}\mathcal{M}_{r}$ and $\bm{X}\in(\bm{U},\bm{\Sigma},\bm{V})$.

2) Riemannian gradient: Euclidean gradient is $\nabla
f(\bm{X})=\mathcal{P}_{\Omega}(\bm{X})-{\bm{J}}$. Riemannian gradient
$\textrm{grad}f(\bm{X})$ gives an orthogonal projection of $\nabla f(\bm{X})$
to the tangent space, i.e.

$\displaystyle\textrm{grad}f(\bm{X})=\Pi_{T_{\bm{X}}\mathcal{M}_{r}}(\nabla
f(\bm{X}))$ (20)
$\displaystyle\Pi_{T_{\bm{X}}\mathcal{M}_{r}}(\bm{G}):\bm{G}\rightarrowtail
P_{U}\bm{G}P_{V}+{P_{U}^{\perp}}\bm{G}P_{V}+{P_{U}}\bm{G}{P_{V}^{\perp}},$
(21)

in which $P_{U}=\bm{U}\bm{U}^{T}$ and $P_{U}^{\perp}=\bm{I}-\bm{U}\bm{U}^{T}$
for $\bm{X}=\bm{U}\textrm{diag}(\bm{\sigma})\bm{V}^{T}$ and
$\bm{\sigma}\in\mathbb{R}^{r}$.

3) Riemannian Hessian: Much like Riemannian gradient, Riemannian Hessian has
the form of

$\textrm{Hess}_{\bm{X}}f[\bm{\xi}_{\bm{X}}]=\Pi_{\mathcal{H}_{\bm{X}}\mathcal{M}_{r}}(\nabla_{\bm{\xi}_{\bm{X}}}\textrm{grad}f(\bm{X})).$
(22)

Let us define the Riemannian connection
$\nabla_{\bm{\eta}_{\bm{X}}}\bm{\xi}_{\bm{X}}$, and substitute the Riemannian
gradient and Riemannian connection into (22). From the Koszul formula, its
connection to Euclidean directional derivative
$D\bm{\xi}_{\bm{X}}[\bm{\eta}_{\bm{X}}]$

$\nabla_{\bm{\eta}_{\bm{X}}}\bm{\xi}_{\bm{X}}=D\bm{\xi}_{\bm{X}}[\bm{\eta}_{\bm{X}}]+(\bm{\theta}_{U},\bm{\theta}_{\Sigma},\bm{\theta}_{V}),$
(23)

where $D\bm{\xi}_{\bm{X}}$ denotes the Euclidean directional derivative,
$\bm{\theta}=(\bm{\theta}_{U},\bm{\theta}_{\Sigma},\bm{\theta}_{V})$ is
related to the solutions of Lyapunov equations.

With Riemannian metric, Riemannian gradient and Riemannian Hessian defined, a
fixed-rank Riemannian trust-region algorithm has been well developed and can
be implemented in Manopt package [24].

###  IV-B Rank Increase Algorithm

To solve the rank minimization problem (III-A) with the fixed rank trust-
region algorithm, we develop a rank-increasing strategy based on the closure
of the set of fixed-rank matrices, $\mathcal{M}_{\leq
r}=\\{\bm{X}\in\mathbb{R}^{M\times Q}:\textrm{rank}(\bm{X})\leq r\\}$.

1 Input: $\Omega,{\bm{J}}$, accuracy $\epsilon$.

2 Initialize: $\bm{X}_{0}^{[1]}\in\mathbb{R}^{m\times n}$,
$(m,n)=\textrm{size}({\bm{J}})$.

3 while  _rank-deficient and not converged_ do

4 Compute a solution to the trust-region subproblem $\mathscr{P}$ with initial
point $\bm{X}_{0}^{[r]}$ for fixed rank $r$ with Riemannian optimization.

5 Compute the value of cost function $f(\bm{X})$.

6 Update rank $r$ to $r+1$ with rank increase strategy (25) in Section IV-B.

7 end while

Output: $\bm{X}^{[r]}$ and rank $r$.

Algorithm 1 Riemannian Trust-Region Algorithm for LRMC Problem (III-A) TABLE I: Recover the Optimal DoF Results in Cases with/without Message Splitting in [12] via LRMC Network Topology in [12] | Fig. 4 | Fig. 6 | Fig. 12(a) | Fig. 12(b) | Fig. 9(a) | Fig. 10(a)
---|---|---|---|---|---|---
DoF in [12] | 1/2 | 1/2 | 1/4 | 1/3 | 3/7 | 2/5
Number of Data Streams ($Q_{i}$) | 1 | 1 | 1 | 1 | 3 | 2
Optimal Rank by LRMC Approach ($n$) | 2 | 2 | 4 | 3 | 7 | 5

Consider the $(r+1)$-th step of iteration where we compute $\bm{X}_{r+1}$ from
$\bm{X}_{r}$. To escape from $\mathcal{M}_{r}$ and embed $\bm{X}_{r}$ to
$\mathcal{M}_{r+1}$, we can give a good initial point by using the linear-
search method

$\bm{X}_{r+1}=\mathcal{R}_{{r+1}}(\bm{X}_{r}+\alpha_{r}\bm{\Xi}_{r}),$ (24)

where we choose the negative Euclidean gradient in the tangent cone
$T_{\bm{X}_{r}}\mathcal{M}_{r+1}$ as our search direction $\bm{\Xi}_{r}$ at
$\bm{X}_{r}$ [25] with step-size $\alpha_{r}$ and $\mathcal{R}_{r+1}$ denotes
the retraction to computation space $\mathcal{M}_{r+1}$. Then
$\bm{\Xi}=\arg\min_{\bm{\Xi}\in T_{\bm{X}_{r}}\mathcal{M}_{\leq
r+1}}\|-\nabla_{\bm{X}_{r}}f-\bm{\Xi}\|_{F}=-\textrm{grad}_{\bm{X}_{r}}f+\bm{\Xi}_{r}^{(1)}$.
$\bm{\Xi}_{r}^{(1)}$ is the orthogonal projection on the tangent space
$T_{\bm{X}_{r}}\mathcal{M}_{r}$. Hence, our rank updating strategy can be
formulated as

$\bm{X}_{r+1}=\mathcal{R}_{r+1}(\bm{X}_{r}+\alpha_{r}(\bm{\Xi}_{r}^{(1)}-\textrm{grad}_{\bm{X}_{r}}f)),$
(25)

which keeps the cost function $f$ decreasing monotonically.

Thus far, through solving fixed-rank subproblem and rank increase strategy, a
complete Riemannian trust-region algorithm can be used to solve the caching
problem in Fog-RAN.

##  V Numerical Results

In this section we first test the convergence rate of our trust-region
approach. We then show that our algorithm can achieve optimality in cases that
[21] refers to. Finally, we run the whole simulation on $20$ messages and
destination nodes with $3$ data streams for varying cache size varies from no
caching to just one message needs to be transmitted. The following two known
algorithms are compared:

  * •

EmbG: This algorithm [19] is developed on the embedded manifold via fixed-rank
optimization [26] with the Riemannian pursuit rank increase strategy [19].

  * •

LMaFit: This algorithm introduces the alternating minimization scheme to solve
problem $\mathscr{P}$ [17].

###  V-A Convergence Rate

Consider a caching problem where the cache size is $10$ and each message is
exactly desired by one destination node, the number of messages and
destination nodes are both $30$ and each message is split to $5$ data streams.
We investigate the convergence rate of all three algorithms given a fixed rank
$r=40$. Fig. 3 shows that the trust-region method outperforms the two
competing schemes in convergence rate and can achieve a higher precision
solution as expected. The convergence rate of LMaFit is lower than other two
manifold approach.

Figure 3: Convergence rate when fixing the rank of $\bm{X}$ as $40$.

###  V-B Achieve the Existing Optimal DoF Results

In this part of experiments, we consider equivalent model for topological
interference management in [12] where optimal DoF is given. To check whether
or not the Riemannian trust-region method achieves the optimal DoF, both of
the unicast cases with and without message splitting are tested. Table I shows
that the LRMC algorithm can numerically achieve all of the optimal DoF results
in [12].

###  V-C Data Rates for Different Cache Size

In order to evaluate the performance of the Riemannian trust-region method in
caching networks with different cache sizes, we test a simple unicast model
with 20 messages and users with 3 data streams for each message. The data rate
we can achieve is $3/n$ where $n$ is the solution of algorithms. The cache
size $m$ varies from $0$ to $19$ since $m=20$ represents that no message is
needed. The average data rates results are shown in Fig. 4. In our tests, the
side information is generated randomly and 50 experiments for each cache size
are performed. When the cost falls below $\epsilon=10^{-7}$, rank increase
strategy stops.

Figure 4: Average data rates varying from cache sizes. Each message has 3 data
streams. The cache size is normalized and tolerance of the cost is $10^{-7}$.

From the experimental result, we observe that Riemannian trust-region approach
can achieve a higher data rates than two other algorithms. And EmbG
outperforms alternating projection method (LMaFit).

##  VI Conclusions

In this paper, we propose a unified low-rank matrix completion approach for
content delivering of caching problem in Fog-RAN. We connected the caching
problem with LRMC problem, and presented a Riemannian trust-region algorithms,
solving the fixed-rank subproblem and giving an rank increase strategy that
can guarantee monotonic decrease of objective function. Numerical results show
that our approach can achieve optimal solution in existing cases for which the
optimal value is known and outperforms algorithms such as EmbG and LMaFit. One
interesting future work is the determination of conditions under which the
Riemannian trust-region algorithm can guarantee optimality.

## References

  * [1] S. Bi, R. Zhang, Z. Ding, and S. Cui, “Wireless communications in the era of big data,” IEEE Commun. Mag., vol. 53, pp. 190–199, October 2015.
  * [2] M. Simsek, A. Aijaz, M. Dohler, J. Sachs, and G. Fettweis, “5G-enabled tactile internet,” IEEE J. Sel. Areas Commun., vol. 34, pp. 460–473, March 2016.
  * [3] Y. Shi, J. Zhang, B. O’Donoghue, and K. B. Letaief, “Large-scale convex optimization for dense wireless cooperative networks,” IEEE Trans. Signal Process., vol. 63, pp. 4729–4743, Sept 2015.
  * [4] J. G. Andrews, S. Buzzi, W. Choi, S. V. Hanly, A. Lozano, A. C. K. Soong, and J. C. Zhang, “What will 5G be?,” IEEE J. Sel. Areas Commun., vol. 32, pp. 1065–1082, June 2014.
  * [5] F. Bonomi, R. Milito, P. Natarajan, and J. Zhu, “Fog computing: A platform for internet of things and analytics,” in Big Data and Internet of Things: A Roadmap for Smart Environments, pp. 169–186, Springer, 2014.
  * [6] Y. Shi, J. Zhang, K. B. Letaief, B. Bai, and W. Chen, “Large-scale convex optimization for ultra-dense Cloud-RAN,” IEEE Commun. Mag., vol. 22, pp. 84–91, June 2015.
  * [7] M. Peng, S. Yan, K. Zhang, and C. Wang, “Fog-computing-based radio access networks: issues and challenges,” IEEE Network, vol. 30, pp. 46–53, July 2016.
  * [8] M. A. Maddah-Ali and U. Niesen, “Fundamental limits of caching,” IEEE Trans. Inf. Theory, vol. 60, pp. 2856–2867, May 2014.
  * [9] A. Liu and V. K. N. Lau, “Exploiting base station caching in MIMO cellular networks: Opportunistic cooperation for video streaming,” IEEE Trans. Signal Process., vol. 63, pp. 57–69, Jan 2015.
  * [10] M. Ji, G. Caire, and A. F. Molisch, “Fundamental limits of caching in wireless D2D networks,” IEEE Trans. Inf. Theory, vol. 62, pp. 849–869, Feb 2016\.
  * [11] S. Barbarossa, S. Sardellitti, and P. D. Lorenzo, “Communicating while computing: Distributed mobile cloud computing over 5G heterogeneous networks,” IEEE Signal Process. Mag., vol. 31, pp. 45–55, Nov 2014.
  * [12] S. A. Jafar, “Topological interference management through index coding,” IEEE Trans. Inf. Theory, vol. 60, pp. 529–568, Jan 2014.
  * [13] M. Effros, S. E. Rouayheb, and M. Langberg, “An equivalence between network coding and index coding,” IEEE Trans. Inf. Theory, vol. 61, pp. 2478–2487, May 2015.
  * [14] E. J. Candès and B. Recht, “Exact matrix completion via convex optimization,” Found. Comput. Math., vol. 9, no. 6, pp. 717–772, 2009\.
  * [15] Y. Shi, J. Zhang, and K. B. Letaief, “Low-rank matrix completion for topological interference management by Riemannian pursuit,” IEEE Trans. Wireless Commun., vol. 15, pp. 4703–4717, July 2016.
  * [16] H. Esfahanizadeh, F. Lahouti, and B. Hassibi, “A matrix completion approach to linear index coding problem,” in IEEE Information Theory Workshop (ITW), pp. 531–535, IEEE, 2014.
  * [17] Z. Wen, W. Yin, and Y. Zhang, “Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm,” Mathematical Programming Computation, vol. 4, no. 4, pp. 333–361, 2012.
  * [18] P. Jain, P. Netrapalli, and S. Sanghavi, “Low-rank matrix completion using alternating minimization,” in ACM Symp. Theory Comput., pp. 665–674, ACM, 2013.
  * [19] M. Tan, I. W. Tsang, L. Wang, B. Vandereycken, and S. J. Pan, “Riemannian pursuit for big matrix recovery,” in Proc. Int. Conf. Mach. Learn. (ICML), pp. 1539–1547, 2014.
  * [20] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization algorithms on matrix manifolds.  Princeton University Press, 2009.
  * [21] H. Maleki, V. R. Cadambe, and S. A. Jafar, “Index coding–an interference alignment perspective,” IEEE Trans. Inf. Theory, vol. 60, pp. 5402–5432, Sept 2014.
  * [22] R. Wang, X. Peng, J. Zhang, and K. B. Letaief, “Mobility-aware caching for content-centric wireless networks: modeling and methodology,” IEEE Commun. Mag., vol. 54, pp. 77–83, August 2016.
  * [23] M. Tao, E. Chen, H. Zhou, and W. Yu, “Content-centric sparse multicast beamforming for cache-enabled Cloud-RAN,” IEEE Trans. Wireless Commun., to appear, 2016.
  * [24] N. Boumal, B. Mishra, P.-A. Absil, and R. Sepulchre, “Manopt, a matlab toolbox for optimization on manifolds,” J. Mach. Learn. Res., vol. 15, no. 1, pp. 1455–1459, 2014.
  * [25] R. Schneider and A. Uschmajew, “Convergence results for projected line-search methods on varieties of low-rank matrices via lojasiewicz inequality,” SIAM Journal on Optimization, vol. 25, no. 1, pp. 622–646, 2015.
  * [26] B. Vandereycken, “Low-rank matrix completion by riemannian optimization,” SIAM Journal on Optimization, vol. 23, no. 2, pp. 1214–1236, 2013.


#  confintROB: An R Package for Confidence Intervals in Robust Linear Mixed
Models

Fabio Mason
University of Geneva &Manuel Koller
&Eva Cantoni
University of Geneva &Paolo Ghisletta
University of Geneva<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Fabio Mason, Manuel Koller, Eva Cantoni, Paolo Ghisletta confintROB Package:
Confindence Intervals in robust linear mixed models confintROB Package
Statistical inference is a major scientific endeavor for many researchers. In
terms of inferential methods implemented to mixed-effects models, significant
progress has been made in the R software. However, these advances primarily
concern classical estimators (ML, REML) and mainly focus on fixed effects. In
the confintROB package, we have implemented various bootstrap methods for
computing confidence intervals (CIs) not only for fixed effects but also for
variance components. These methods can be implemented with the widely used
lmer() function from the lme4 package, as well as with the rlmer() function
from the robustlmm package and the varComprob() function from the
robustvarComp package. These functions implement robust estimation methods
suitable for data with outliers. The confintROB package implements the Wald
method for fixed effects, whereas for both fixed effects and variance
components, two bootstrap methods are implemented: the parametric bootstrap
and the wild bootstrap. Moreover, the confintROB package can obtain both the
percentile and the bias-corrected accelerated versions of CIs. robust
estimators; bootstrap; confidence intervals; linear mixed models; lme4;
robustlmm; robustvarComp, R robust estimators; bootstrap; confidence
intervals; linear mixed models; lme4; robustlmm; robustvarComp, R Fabio Mason
Faculty of Psychology and Educational Sciences
University of Geneva
Boulevard du Pont d’Arve 40, 1211 Geneva, Switzerland
E-mail:

Paolo Ghisletta
Faculty of Psychology and Educational Sciences
University of Geneva
Boulevard du Pont d’Arve 40, 1211 Geneva, Switzerland
E-mail:

Eva Cantoni
Research Center for Statistics and Geneva School of Economics and Management
University of Geneva
Boulevard du Pont d’Arve 40, 1211 Geneva, Switzerland
E-mail:

Manuel Koller
Zürich, Switzerland
E-mail:

##  1 Introduction

The linear mixed model (LMM) is frequently applied to correlated multilevel
data. The package lme4 (Bates _et al._ , 2015) has notably contributed to the
dissemination of LMM in various research fields. Indeed, it allows the
analysis of data of variable degrees of complexity (e.g., hierarchical, cross-
effects or mixed), with or without missing data, allows to compute confidence
intervals (CIs) for the different types of parameter via the confint.merMod()
function, and performs testing (and compute corresponding $p$-values) with the
lmerTest package (Kuznetsova _et al._ , 2017).

Recently, the field of robustness in the context of LMM has undergone
significant developments (e.g. Pinheiro _et al._ , 2001; Copt and Victoria-
Feser, 2006; Chervoneva and Vishnyakov, 2014; Koller, 2013; Agostinelli and
Yohai, 2016; Saraceno _et al._ , 2023; Lucadamo _et al._ , 2021; Zhang _et
al._ , 2020; Burger and Schall, 2018; Mason _et al._ , 2021, 2024). However,
only a few of these proposals are implemented in R: the Robust Scoring
Equations (RSE) estimator of Koller (2013) and Koller and Stahel (2023) is
implemented in the package robustlmm (Koller, 2016), the composite-$\tau$
estimator (cTAU) of Agostinelli and Yohai (2016) and the S-estimator (S) of
Copt and Victoria-Feser (2006) are implemented in the package robustvarComp
(Agostinelli and Yohai, 2016) and the multivariate $t$ maximum likelihood
estimator (tML) of Pinheiro _et al._ (2001) is implemented in the heavy
package (Osorio, 2019). For all these estimators, the Wald test (heavyLme()
and varComprob() return $p$-values in the summary, whereas rlmer() returns
standard errors) is the only statistical inference option, and with the
exception of the estimators implemented in robustvarComp, this test only
applies to the fixed effects. Unfortunately, the robust estimators are not
compatible with the lmerTest package or the confint.merMod() function. This
gap in existing software inspired us to create the confintROB package, which
allows computing confidence intervals (CIs) based on the parametric and the
wild bootstrap for both fixed effects and variance components of models
estimated with the lmer() function from lme4, the rlmer() function from
robustlmm, and with the varComprob() function from robustvarComp. Our package
was inspired by the confint.merMod() and confint() functions for the arguments
and for the output, with the difference that it implements two bootstrap
methods also applicable to the robust estimators implemented in the robustlmm
package and robustvarComp package, as well as the Wald type confidence
intervals.

This article is intended as a tutorial for confintROB(), a function to compute
CI with robust estimators implemented in the R environment (Singmann _et al._
, 2015) in the homonymous package. The function can currently be applied to 5
classes of estimators and allows one to choose between 5 choices of CI types,
including 4 bootstrap methods. The package is available on CRAN at
https://cran.r-project.org/package=confintROB under the GPL-2 license.

In the next Section (2), we introduce the technicalities germane to the model
and the estimators, and provide details on the wild and the parametric
bootstrap schemes in their _percentile_ and _Bias-corrected and accelerated
(BCa)_ versions. In Section 3, we introduce a real data application with the
function rlmer(). Limitation and future extensions are discussed in Section 4.

##  2 Technicalities

###  2.1 The Linear Mixed Model (LMM) and its (robust) estimation

The LMM can be written as:

$\begin{array}[]{cc}\boldsymbol{y}_{i}=X_{i}\boldsymbol{\gamma}+Z_{i}\boldsymbol{b}_{i}+\boldsymbol{\varepsilon}_{i},\end{array}$
(1)

where $\boldsymbol{y}_{i}$ is a vector of length $J_{i}$ containing the
responses of cluster $i$ (e.g. participants in a longitudinal study),
$\boldsymbol{\gamma}$ is a vector of coefficients for the p fixed effects,
$X_{i}$ is the ($J_{i}$ $\times$ $p$) design matrix for fixed effects,
$\boldsymbol{b}_{i}$ is the vector of random effects of length q, independent
of the errors $\boldsymbol{\varepsilon}_{i}$ and $Z_{i}$ is the ($J_{i}$
$\times$ $q$) design matrix for the random effects. The vector
$\boldsymbol{\varepsilon}_{i}=(\varepsilon_{i1},\ldots,\varepsilon_{iJ_{i}})^{T}$
contains the $J_{i}$ individual error terms for cluster $i$. Generally, it is
assumed that random effects and error terms are normally distributed around
zero (Laird and Ware, 1982):

$\boldsymbol{b}_{i}\sim\mathcal{N}(\boldsymbol{0},\Sigma)\ \ \mbox{ and }\ \
\boldsymbol{\varepsilon}_{i}\sim\mathcal{N}(\boldsymbol{0},\sigma_{\varepsilon}^{2}I),$
(2)

where $\Sigma=\Sigma(\boldsymbol{\theta})$ is the $q\times q$ (parameterized)
covariance matrix of the random effects and $\sigma_{\varepsilon}^{2}I$ is the
$(J_{i}\times J_{i})$ (diagonal) covariance matrix of the error terms. That
is, the model assumes that the random effects $\boldsymbol{b}_{i}$ and errors
$\boldsymbol{\varepsilon}_{i}$ stem from zero-centered normal distributions.
The entire set of parameters to be estimated is
$(\gamma,\sigma_{\varepsilon}^{2},\theta)^{T}$. We consider two general
classes of estimators: _classical_ and _robust_. The classical estimators are
based on the likelihood of observing the data, assuming the model to be
correct in the population. The full maximum likelihood (ML) approach and the
restricted maximum likelihood (REML) approach assume that both random effects
and error terms are normally distributed (see Equation (2)). Because classical
estimators may lead to biased estimates and wrong inferential conclusions with
contaminated data (i.e., data that contain outliers, defined as observations
with extremely low probability of occurrence assuming normality of random
effects and errors terms as in Equation (2); cf. Mason _et al._ (2021, 2024)),
robust estimators have been developed (e.g. Copt and Victoria-Feser, 2006;
Agostinelli and Yohai, 2016; Koller, 2016). These estimators are thus
constructed to reduce bias in estimation of parameters in the presence of
outliers.

####  2.1.1 Classical Estimation

The classical estimator of the LMM is maximum likelihood (ML), which assumes
multivariate normality (e.g. Richardson and Welsh, 1995; Welsh and Richardson,
1997). The (log-)likelihood function expresses the likelihood of the data as a
function of the model parameters, assuming the model is correct. The estimates
are the parameter values that maximize this function, for both fixed effects
$\boldsymbol{\gamma}$ and variance components $\sigma_{\varepsilon}^{2}$ and
$\boldsymbol{\theta}$.

The ML can be implemented in two versions, either full (ML), or restricted
(REML). The main difference is that the latter estimates the variance
components contingent on the fixed effects, while the former estimates both
sets of parameters simultaneously. The function lmer() from the package lme4
(Bates _et al._ , 2015) in R implements both versions.

####  2.1.2 Robust Estimation

Several robust methods have been proposed for the LMM. The confintROB package
is available with two functions that implement a total of 8 robust estimators.
The rlmer() function implements the RSE estimator and its fast approximative
version (fast RSE), whereas the varComprob() function implements the $\tau$-,
S- and MM-estimators and their _composite_ versions (see Table 1). All
estimators presume a so-called central model (see, e.g., Koller, 2013), where
a central model $G$ is assumed for the majority of data (e.g., a normal
distribution) and an alternative unknown model $H$ is assumed for the
outliers, which are hence assumed from a different population. Typically, the
substantive focus of the analysis is on the parameters of $G$. All 8 robust
estimators also provide robustness weights that are applied to the data during
the estimation procedure. The weights reveal important information about the
possible outliers, because they indicate the extent to which the data are
supposed to stem from the $G$ (weights close to $1$) rather than the $H$
(weights close to $0$) model. Whereas some robust methods estimate weights at
the single observation level, to accommodate within-subject outliers (i.e.,
outlying points within the vector of an otherwise non-outlying subject),
others only estimate subject-specific weights, which are applied to all
observations of a given subject and are thus ideal for identifying between-
subject outliers (i.e., subjects whose entire vector of data are atypical for
$G$). Other methods yet estimate simultaneously both types of weights to
accommodate both within- and between-subject outliers.

In this paper we discuss the RSE estimator implemented in the rlmer()
function. This function allows to handle data with complex structure: nested
or non-nested observations, balanced or unbalanced designs, with or without
missing data, correlated or non-correlated random effects, etc.

Type | Function | Package | Estimator |  Argument
---|---|---|---|---
Classical |  lmer() |  lme4 | ML |  REML = FALSE
|  |  | REML |  REML = TRUE
Robust |  rlmer() |  robustlmm | RSE |  method = "DAStau"
|  |  | fast RSE |  method = "DASvar"
|  varComprob() |  robustvarComp | $\tau$ |  method = "Tau"
|  |  | composite-$\tau$ |  method = "compositeTau"
|  |  | S |  method = "S"
|  |  | composite-S |  method = "compositeS"
|  |  | MM |  method = "MM"
|  |  | composite-MM |  method = "compositeMM"
Table 1:  Overview of available functions and estimators with confintROB().

###  2.2 Statistical Inference in the LMM

In the frequentist approach to statistical inference, assessing a null
hypothesis is usually performed by either performing an appropriate
statistical test about a parameter of interest and comparing its $p$-value to
a given statistical significance level, or by building a CI for that parameter
and checking whether it covers the value 0.

####  2.2.1 Tests

For LMM, Kuznetsova _et al._ (2017) developed the lmerTest package that
implements F-, and $t$-tests for the fixed effects parameters with degrees of
freedom calculated following the Satterthwaite (Satterthwaite, 1946) or the
Kenward-Roger (Kenward and Roger, 1997) approximations. For the fixed effect
parameters, these tests are largely considered as the benchmark for classical
LMM estimators (Kuznetsova _et al._ , 2017). The package also implements the
likelihood ratio test (LRT) for the random effects. These methods are
applicable only to the classical LMM estimators and objects of class ‘lmerMod’
obtained via the lmer() function.

With robust estimators, standard errors are computed for the fixed effects
with rlmer() and varComprob(), and for random effects with the latter only.
With these quantities, one can then construct $t$-test statistics and compute
associated $p$-values according to the reference distribution (Student $t$).

####  2.2.2 Confidence Intervals

Confidence Intervals (CIs) are an alternative approach to perform statistical
inference. The lme4 package implements three kinds of CI for objects of class
‘lmerMod’ with the
confint.merMod() function: the _profile_ method based on the LRT test, which
is available only with ML estimates, the _Wald_ method based on the standard
error estimates, and available only for fixed effects, and the _bootstrap_
method, which is often recommended to obtain CIs for fixed effects and
variance components for LMM parameters (e.g., Koller, 2016; Modugno and
Giannerini, 2013).

To obtain bootstrap CIs, first, the empirical distribution of the parameter of
interest, let’s say $\beta$, is constructed, based on the multiple estimations
$\hat{\beta}_{1}^{*},\ldots,\hat{\beta}_{B}^{*}$ of the parameter in the $B$
resamples. With the percentile version, implemented in confint.merMod(), the
lower and upper bounds are obtained from the $\frac{\alpha}{2}\times 100$ and
$(1-\frac{\alpha}{2})\times 100$ empirical percentiles of the bootstrap
distribution. When the bootstrap distribution is skewed and/or biased,
Tibshirani and Efron (1993) proposed a refined version, called bias-corrected
and accelerated (BCa), which is however computationally very intensive. To
mitigate potential bias, the BCa estimates the parameter $\widehat{z}_{0}$
that corresponds to the proportion of estimates from the resamples smaller
than the estimate of interest from the original sample. To correct for
skewness, the “accelerated” parameter $\widehat{a}$ is estimated using the
Jacknife procedure (e.g. Efron, 1979), resulting in

$\widehat{a}=\frac{\sum_{i=1}^{n}\left(\widehat{\beta}_{(\cdot)}-\widehat{\beta}_{(i)}\right)^{3}}{6\left[\sum_{i=1}^{n}\left(\widehat{\beta}_{(\cdot)}-\widehat{\beta}_{(i)}\right)^{2}\right]^{3/2}},$

where $n$ is the number of subjects, $\widehat{\beta}_{(i)}$ is the estimate
of $\beta$ on the sample without subject $i$ and $\widehat{\beta}_{(\cdot)}$
is the mean of the $n$ $\widehat{\beta}_{(i)}$ obtained. The BCa interval is
then constructed as follows:

$[\widehat{\beta}^{*}_{\alpha 1};\widehat{\beta}^{*}_{\alpha 2}],$

where

${\alpha
1}=\Phi\left(\widehat{z}_{0}+\frac{\widehat{z}_{0}+z_{\alpha/2}}{1-\widehat{a}(\widehat{z}_{0}+z_{\alpha/2})}\right)$

and

${\alpha
2}=\Phi\left(\widehat{z}_{0}+\frac{\widehat{z}_{0}+z_{1-\alpha/2}}{1-\widehat{a}(\widehat{z}_{0}+z_{1-\alpha/2})}\right),$

where $z_{\alpha/2}$ and $z_{1-\alpha/2}$ are the $\alpha/2$ and
$(1-\alpha/2)$ quantiles of a Gaussian distribution respectively, $\alpha$ is
the level of confidence, and $\Phi(\cdot)$ is the cumulative distribution
function of a Gaussian distribution.

The confintROB() function implements two kinds of bootstrap (parametric and
wild) in both versions (percentile and BCa), plus the Wald-type confidence
interval, for a total of 5 options available with objects from lmer(), rlmer()
and varComprob(). All methods are summarized in Table 2.

##### Parametric Bootstrap

First, the procedure estimates fixed effects ($\boldsymbol{\gamma}$) and
variance component parameters $(\sigma_{\varepsilon}^{2},\boldsymbol{\theta})$
from model (1) and (2) on the original data
$(\boldsymbol{y}_{i},X_{i},Z_{i};i=1,\ldots,n)$, to produce
$\hat{\boldsymbol{\gamma}}$, $\hat{\sigma}_{\varepsilon}^{2}$, and
$\Sigma(\hat{\boldsymbol{\theta}})$. Then, for $b=1,\ldots,B$, the procedure :

  1. 1.

for $i=1,\ldots,n,$ generates
$\boldsymbol{\varepsilon}_{i}^{*}\sim\mathcal{N}(\hat{\sigma}_{\varepsilon}^{2})$
and $\boldsymbol{b}_{i}^{*}\sim\mathcal{N}(\Sigma(\hat{\boldsymbol{\theta}}))$
to build

$\boldsymbol{y}_{i}^{*}=X_{i}\widehat{\boldsymbol{\gamma}}+Z_{i}\boldsymbol{b}_{i}^{*}+\boldsymbol{\varepsilon}_{i}^{*}$

  2. 2.

estimates $\boldsymbol{\gamma}$, $\sigma_{\varepsilon}^{2}$, and
$\boldsymbol{\theta}$ from model (1) and (2) on the resample
$(\boldsymbol{y}_{i}^{*},X_{i},Z_{i};i=1,\ldots,n)$ to produce
$\hat{\boldsymbol{\gamma}}^{*}$, $(\hat{\sigma}_{\varepsilon}^{2})^{*}$, and
$\Sigma(\hat{\boldsymbol{\theta}}^{*})$.

The estimator used at step 2 can either be the same as, or differ from, that
used on the original sample. Because usually $B=5000$ and because robust LMM
estimation can be computationally very heavy111 Over 500 large simulated
samples (40 subjects measured 80 times following the design described in
Modugno and Giannerini (2013)) the mean computing time for a CI is 3.034 days
with the first method and 0.9 days with the second method on the most powerful
computing facility at the University of Geneva (4300 processors, 50 GPUs, up
to 1.5TB RAM per calculation server, 1.2Po of dedicated storage)., we choose
to apply ML at step 2 also when a robust estimator was used at step 1 (cf.
Mason _et al._ (2024)).

##### Wild Bootstrap

The wild bootstrap method replaces step 1 of the parametric bootstrap
computing “residuals” $\tilde{\boldsymbol{\upsilon}}_{i}$ as follows (Modugno
and Giannerini, 2013):

$\tilde{\boldsymbol{\upsilon}}_{i}=\mbox{diag}(I-H_{i})^{-1/2}\circ(\boldsymbol{y}_{i}-X_{i}\widehat{\boldsymbol{\gamma}}),$

with $H_{i}=X_{i}(X^{T}X)^{-1}X_{i}^{T}$, where the operator “$\circ$” denotes
the (element-wise) Hadamard product.

For $i=1,\ldots,n$, the procedure samples independently $w^{*}_{i}$ from the
following standardized distribution (Mammen, 1993):

$w^{*}_{i}=\left\\{\begin{array}[]{rl}-\dfrac{\sqrt{5}-1}{2}&\mbox{with
probability}\ p=(\sqrt{5}+1)/(2\sqrt{5})\\\ \dfrac{\sqrt{5}+1}{2}&\mbox{with
probability}\ q=1-p\end{array}\right.$

and builds individual responses

$\boldsymbol{y}^{*}_{i}=X_{i}\widehat{\boldsymbol{\gamma}}+\tilde{\boldsymbol{\upsilon}_{i}}w^{*}_{i}.$

The wild bootstrap can be applied to all LMM estimators considered here.

Modugno and Giannerini (2013) compared several methods to generate the
resamples, including the parametric bootstrap and the semi-parametric wild
bootstrap. The assumptions are stronger for the former method than for the
latter (i.e., fixed covariates, correct specification of both fixed effects
and variance components, and homoscedasticity and normality for random
effects). The authors found that the semi-parametric wild method was superior
in coverage in large samples and with heteroscedastic random effects.

Method | Type | Estimates |  Function specification
---|---|---|---
Wald | - | Fixed effects |  confintROB(object = model.output, method = "Wald")
Percentile | wild | All |  confintROB(object = model.output, boot.type = "wild")
| parametric | All |  confintROB(object = model.output, boot.type = "parametric")
BCa | wild | All |  confintROB(object = model.output, clusterID = "id", method = "BCa", boot.type = "wild")
| parametric | All |  confintROB(object = model.output, clusterID = "id", method = "BCa", boot.type = "parametric")
Table 2:  Overview of available methods to compute CI with confintROB().

###  2.3 Design of the confintROB function

The structure of the confintROB() function is inspired by the confint.merMod()
function from the lme4 package, which computes CIs from objects of class
‘lmerMod’. Similarly to confint.merMod(), the value of a confintROB object is
a matrix containing the names of parameters in the first column, the lower and
upper bounds of CIs in the second and third column, respectively, labelled
with the corresponding percentages. The argument parm can be used to select a
subset of parameters of interest using a vector of string characters
indicating their names, or numerical values indicating the corresponding rows
in the full matrix. The first argument is object, the fitted model with either
lmer() (as confint.merMod()), rlmer() or varComprob(). The type of CI is
chosen with the argument method, which can take the values "Wald", "boot"
(default), or "BCa", which is an addition to confint.merMod(). If "boot" or
"BCa" are selected, then argument boot.type needs to be specified, choosing
from "wild" (default) or "parametric". The argument clusterID is needed only
with method="BCa" to specify the cluster variable with a string character
corresponding to the variable name in the dataset. In the example of Section
3, participants are the clusters, and the variable name that identifies the
cluster is "id". The numerical argument level is required to define the level
of confidence and must be positive and smaller to 1, with default value .95.
For bootstrap CIs (i.e., with method = "boot" or method = "BCa"), the number
of bootstrap samples can be changed (the default value is 5000) with the
argument nsim. There is an optional argument verify.saved that allows checking
for a correspondence between the current results and those obtained, for
example, in a previous version of a package in dependency. It requires an
object produced by the confintROB() function.

##  3 Illustration

Here, we analyze the medication dataset, often used for didactic purposes
(Singer and Willett, 2003), and originally discussed in Tomarken _et al._
(1997). The dataset is included in the package and contains 1242 observations
with 5 variables:

R> library(confintROB)R> data(medication)R> str(medication)

’data.frame’:1242 obs. of 5 variables: $ obs : int 1 2 3 4 5 6 7 8 9 10 ... $
id : int 1 1 1 1 1 1 1 1 1 1 ... $ treat: int 1 1 1 1 1 1 1 1 1 1 ... $ time :
num 0 0.333 0.667 1 1.333 ... $ pos : num 107 100 100 100 100 ... where obs
denotes the observation number, id identifies the 64 subjects of the sample ,
treat is a dichotomous variable to define control (treat=0 for 27 subjects)
and treatment (treat=1 for 37 subjects), time takes 21 values (from 0 to 6.67
increasing by 0.33 at each time point), pos is the positive mood score,
varying from 100 to 500 (mean = 167.79). Although the data set does not
contain any NAs, it is unbalanced because 33 subjects were assessed fewer than
21 times, which amounts to 7.59% of data missing.

The group-specific representation is presented in Figure 1. Here, we see each
subject’s measurements represented by lineplots. We can see that individual
trajectories are extremely variable between subjects (between-subject
variability) and unstable in time (within-subject variability). If we try to
represent the trajectories by straight lines, the two groups seem to differ
with respect to their slopes, but not on their starting levels. To emphasize
this, we plotted with a solid (treat = 0) and a dotted (treat = 1) line the
group-implied average straight-line trajectories. The shade of the
individuals’ lineplots provides information about the weight(s) assigned to
each individual by the estimator. The lower the weight(s) of an individual,
the darker their line, thus the greater the chance that the individual be
considered a between-subject outlier (Koller, 2016). Two individual
trajectories are particularly darker, because those subjects’ data appear
extreme with respect to their group, probably due to their extremely high
slope scores (Koller, 2016). Their weights are 0.41 and 0.12, for the
individual in the treat = 0 and treat = 1 group, respectively. A second
individual in the treat = 1 group has a weight of 0.90, rendering their line
less distinguishable than those of the two individuals mentioned before.
According to Mason _et al._ (2024), these are the outliers that have the
greatest impact on the estimates of the time by group interaction fixed
parameter, which is the parameter of main interest in many group comparison
studies, because it captures the different effect of time across the two
groups (e.g., differential treatment effect). The inferential results of this
particular parameter are also greatly affected by this type of outlier. The
presence of random slope between-subject outliers can greatly affect the
rejection rates regarding the fixed interaction effect Mason _et al._ (2024).

Figure 1: The individual positive mood trajectories (pos) as a function of
waves of assessment (time) and treatment (treat = 0 on the left and treat = 1
on the right panel). The regression lines are based on the Maximum Likelihood
(ML) estimates of the model in Equation (1). The solid black line is the
predicted average trajectory for individuals of treat = 0 and the dotted line
is the predicted average trajectory for individuals of treat = 1. The
lineplots represent the individual trajectories (The clearer the lines, the
closer the individual estimated weights (by RSE) to 1.).

A first step towards assessing whether outliers influenced estimation consists
in comparing parameter estimates from a classical estimator, here ML, with
those from a robust alternative, here RSE. In this example, we are
particularly motivated by the effect of time on pos possibly moderated by
treat. Below, we show the code to model this hypothesis taking into account
that participants may differ in their intercept and slope values.

R> library(lmerTest)R> model.ML <\- lmer(pos ~ treat * time + (time|id), data
= medication,\+ REML = "FALSE") The object model.ML contains the parameter
estimates obtained with the lmer() function based on the ML estimator because
of the argument REML=FALSE (for more details about this function see Bates _et
al._ (2015)). To obtain the equivalent model with RSE estimates, we use the
function rlmer().

R> library(robustlmm)R> model.RSE <\- rlmer(pos ~ treat * time + (time|id),
data = medication,\+ method = "DAStau", init = inits.RSE) An easy way to
compare estimates is to produce a table using the compare() function from
package robustlmm calling objects model.ML and model.RSE:

R> compare(model.ML, model.RSE)

| model.ML | model.RSE
---|---|---
Coefficients (Std. Error) |  |
(Intercept) | 167.46 ( 9.33) | 163.831 ( 9.77)
treat | -3.11 (12.33) | 0.232 (12.90)
time | -2.42 ( 1.73) | -2.551 ( 1.40)
treat:time | 5.54 ( 2.28) | 4.206 ( 1.85)
Variance components |  |
(Intercept) $|$ id | 45.95 | 48.51
time $|$ id | 7.98 | 6.38
Correlations |  |
(Intercept) $\times$ time $|$ id | -0.332 | -0.431
$\sigma$ | 35.1 | 27.8
deviance | 12680 |
rho.e |  | smoothed Huber (k = 1.345, s = 10)
rho.$\sigma$.e |  | smoothed Huber, Proposal 2 (k = 1.345, s = 10)
rho.b_1 |  | smoothed Huber (k = 5.14, s = 10)
rho.$\sigma$.b_1 |  | smoothed Huber (k = 5.14, s = 10)
Table 3: Comparison of classical and robust fits. Rendering of output of
xtable(compare(model.ML,model.RSE)).

Table 3 shows the output of xtable in rendered form. The _Coefficients_
section contains the estimates of the fixed parameters, the estimated variance
of the random effects are in the section _Variance components_ , their
estimated correlation in the _Correlations_ section, and the standard
deviation of error terms is called $\sigma$. The last part of the output gives
the deviance for model.ML and the functions used in the algorithm to obtain
the robust estimates (see Koller (2016) for more details). The estimates of
the variance components and the fixed-effect coefficients for the intercept,
treat and time are at times very similar (e.g., fixed and random effects of
Intercept) and at times different (e.g., fixed and random effects of time)
across the two fits. A meaningful difference in estimates is the fixed effect
of the interaction, which is slightly smaller with RSE. This is coherent with
the fact that the very high slope values of two outliers in group treat=1 tend
to increase the difference in the average slope between the two groups. These
results can be complemented by CIs, which are also useful for inferential
decisions. We therefore introduce the confintROB() function that applies the
percentile wild bootstrap to the model estimated by ML.

R> set.seed(3)R> wild.ML <\- confintROB(object = model.ML, boot.type =
"wild",\+ verify.saved = wild.ML) The function confintROB() can be applied to
any object of class ‘lmerMod’, ‘lmerModLmerTest’, ‘rlmerMod’ or ‘varComprob’.
In this example, model.ML is an object of class ‘lmerMod’ (see Table 2). The
percentile bootstrap CI being the default argument, we specify that we want
the wild bootstrap with the boot.type argument. We use the default values for
argument nsim and the argument level (respectively, 5000 bootstrap samples and
.95 confidence level).

R> wild.ML

2.5 % 97.5 %(Intercept) 150.1803810 185.45316421treat -26.9182945
20.44265435time -5.3947424 0.67222061treat:time 1.3825904 9.89946119Sigma id
(Intercept) 37.1847684 53.02209619Sigma id time 5.2272994 10.76968251Sigma id
(Intercept) time -0.5461602 -0.07403813Sigma Residual 28.7413919
41.48149925attr(,"fullResults") Full results of confintROB, a list with
components: "Percentile", "bootstrap_estimates" Results are saved in the
object wild.ML. With the print() function, the matrix of selected CIs is
returned. In this example, the CIs for all the parameters are provided (the
default). If the interest is on a subset of the parameters, the argument parm
can be used. For instance, to obtain the CI for the variable treat only, one
would specify parm = 2 or parm = "treat". To obtain the CIs for treat and
Sigma id time, one would specify parm = c(2, 6) or parm = c("treat", "Sigma id
time"). The first column of the output matrix always contains the names of the
parameters. When the name begins with Sigma, the corresponding CIs are for the
standard deviations of residuals (Sigma Residual), random effect, or
correlation between random effects. For random effects, the string Sigma is
followed by the name of the cluster (here, id (Intercept) indicates that it is
the random effect of the intercept, when random slopes are supposed, the name
of the variable follows id (here, time) and correlation is indicated with both
names (in this example, (Intercept) time). The remaining effects are fixed.
The second and third columns represent, respectively, the lower and upper
bound of the CIs. The string attr(,"fullResults") indicates that the object
wild.ML also contains an attribute of class ‘fullResults’, which is a list
containing the complete matrix with the percentile version of the CIs
("Percentile") and all the parameter estimates of the simulated samples
("bootstrap_estimates"). We now apply the wild bootstrap to RSE estimates by
replacing the model.ML object by model.RSE.

R> set.seed(3)R> wild.RSE <\- confintROB(object = model.RSE,\+ boot.type =
"wild",\+ verify.saved = wild.RSE)R> wild.RSE

2.5 % 97.5 %(Intercept) 146.44579838 181.89124806treat -23.80293651
23.55395617time -5.50695767 0.57216356treat:time 0.02290454 8.63920558Sigma id
(Intercept) 37.14917091 53.38396261Sigma id time 5.22385401 10.89180572Sigma
id (Intercept) time -0.54931053 -0.07081393Sigma Residual 28.69255451
41.42750808attr(,"fullResults") Full results of confintROB, a list with
components: "Percentile", "bootstrap_estimates" Again, estimates are extremely
similar to the ones obtained with ML, with a slight negative shift in the CIs
for the fixed-effect parameters of the Intercept, treat, and the interaction
between treat and time.

Both classical and robust estimators agree in concluding that there is no
effect of treat, time, and a moderately strong, negative correlation between
random effects in the original sample. However, whereas the classical
estimator clearly concludes that there is a time by group interaction effect,
the conclusion drawn from the robust estimator is not as clearcut, with a
lower bound much closer to 0. In substantive terms, this difference in results
is of utmost importance, as it may potentially invalidate a treatment.

Obviously, in this first example with real data it is impossible to know the
true values of the parameters in the population, and some may suspect that the
robust estimator is simply more conservative or less powerful than the
classical estimator. Therefore, we introduce a second example (medsim)
inspired by the first one, but based on simulated data, for which we know the
true parameter values in model (1) and (2):
$\boldsymbol{\gamma}=(167.46,~{}-3.11,~{}-2.42,~{}4.00)$,
$\boldsymbol{\theta}=(2111.54,~{}63.74,~{}-121.63)^{T}$ and
$\sigma_{\varepsilon}^{2}=1229.93$.

R> str(medsim)

’data.frame’:420 obs. of 5 variables: $ obs : int 1 2 3 4 5 6 7 8 9 10 ... $
id : int 1 1 1 1 1 1 1 2 2 2 ... $ time : int 0 3 6 9 12 15 18 0 3 6 ... $
treat: int 0 0 0 0 0 0 0 0 0 0 ... $ pos : num 205.8 192.5 126 51.8 59.6 ...
The sample size is similar to that of the real data example ($N=60$) but with
balanced and complete data and with fewer repeated assessments per subjects
($J_{i}=J=7$), to reduce the file size. Time is equal to (0,3,6,9,12,15 and
18). Two participants from to the control group were originally in the
treatment group: we intentionally changed their respective label to create
outliers (see the dark trajectories in the left panel of the Figure 2). Unlike
the “outliers” in (treat = 1) highlighted in Figure 1, which tended to
accentuate the difference in trajectories between participants in the two
groups and thereby increase the probability of obtaining a significant effect
for the interaction parameter $\gamma_{3}$, these outliers should increase the
confusion between the two groups and, in this case, decrease the probability
of rejecting the related null hypothesis (similarly to the outlier in (treat =
0) highlighted in Figure 1).

Figure 2: The individual positive mood trajectories (pos) as a function of
waves of assessment (time) and treatment (treat = 0 on the left and treat = 1
on the right panel). The regression lines are based on the Maximum Likelihood
(ML) estimates of the model in Equation (1). The solid black line is the
predicted average trajectory for individuals of treat = 0 and the dotted line
is the predicted average trajectory for individuals of treat = 1. The
lineplots represent the individual trajectories (The clearer the lines, the
closer the individual estimated weights (by RSE) are to 1.).

After estimating the same model on the simulated data by creating objects
modelSim.ML and modelSim.RSE respectively, the estimation of the parameters
can once again be directly compared using the compare() function:

R> modelSim.ML <\- lmer(pos ~ treat * time + (time|id), data = medsim, REML =
"FALSE")

R> modelSim.RSE <\- rlmer(pos ~ treat * time + (time|id), data = medsim,\+
method = "DAStau", init = inits.SimRSE)

R> compare(modelSim.ML, modelSim.RSE)

| modelSim.ML | modelSim.RSE
---|---|---
Coefficients (Std. Error) |  |
(Intercept) | 179.50 (10.18) | 178.55 (10.64)
treat | -20.54 (14.40) | -21.07 (15.04)
time | -2.90 ( 1.14) | -3.00 ( 1.21)
treat:time | 3.03 ( 1.61) | 3.19 ( 1.71)
Variance components |  |
(Intercept) $|$ id | 51.13 | 52.83
time $|$ id | 5.90 | 6.17
Correlations |  |
(Intercept) $\times$ time $|$ id | -0.234 | -0.248
$\sigma$ | 32.6 | 32.8
deviance | 4428 |
rho.e |  | smoothed Huber (k = 1.345, s = 10)
rho.$\sigma$.e |  | smoothed Huber, Proposal 2 (k = 1.345, s = 10)
rho.b_1 |  | smoothed Huber (k = 5.14, s = 10)
rho.$\sigma$.b_1 |  | smoothed Huber (k = 5.14, s = 10)
Table 4: Comparison of classical and robust fits on the simulated dataset.
Rendering of output, xtable(compare(modelSim.ML, modelSim.RSE)).

Again, the parameter estimates for fixed effects, their standard errors, and
variance components are very similar for both estimators. By applying the wild
bootstrap using the confintROB() function to these two new objects, we once
again obtain a CI for each parameter of the model.

R> set.seed(3)R> wildSim.ML <\- confintROB(object = modelSim.ML,\+ boot.type =
"wild",\+ verify.saved = wildSim.ML)R> wildSim.ML

2.5 % 97.5 %(Intercept) 158.5044035 201.29909601treat -48.8550232
6.93072963time -5.2582788 -0.51693479treat:time -0.1187373 6.17566692Sigma id
(Intercept) 38.4990062 62.30562438Sigma id time 4.6069444 7.04286448Sigma id
(Intercept) time -0.4960982 0.05952634Sigma Residual 28.0164266
37.43709604attr(,"fullResults") Full results of confintROB, a list with
components: "Percentile", "bootstrap_estimates"

R> set.seed(3)R> wildSim.RSE <\- confintROB(object = modelSim.RSE,\+ boot.type
= "wild",\+ verify.saved = wildSim.RSE)R> wildSim.RSE

2.5 % 97.5 %(Intercept) 157.54765680 199.99603194treat -49.50558142
6.54352734time -5.33278702 -0.59838226treat:time 0.04245368 6.33599786Sigma id
(Intercept) 38.49685977 62.11219643Sigma id time 4.60446985 7.04626841Sigma id
(Intercept) time -0.49917653 0.06116436Sigma Residual 27.94773579
37.41651808attr(,"fullResults") Full results of confintROB, a list with
components: "Percentile", "bootstrap_estimates" Once again, the confidence
intervals (CIs) are very similar for each parameter, and again, the
interaction effect is marginally significant only with the robust estimator,
thus leading to an opposite substantive finding concerning the treatment. This
result suggests that it is not a matter of a more conservative or less
powerful method, but rather the nature of the outliers that determines the
outcome of the inferential analysis.

##  4 Summary and discussion

The main purpose of the current package is to provide a complementary tool to
the robust estimators of LMM obtained with the packages robustlmm and
robustvarComp, so as to facilitate comparisons with classical estimators, and
to allow for inferential decisions based on CIs for both fixed and variance
components. The output and options are constructed to be as similar as
possible to the function confint() from package stats and confint.merMod()
from package lme4 so that users of those packages can readily become
accustomed to results from the confintROB package.

In the future, it would be worthwhile to investigate also non-parametric
tests. While parametric tests are more powerful when the assumptions are met,
these assumptions can be quite demanding, and in cases where they are not met,
parametric tests are generally less effective (Anderson, 1961). In the context
of the LMM, the cluster bootstrap is available for classical estimator in R
through the ClusterBootstrap package (Deen and de Rooij, 2020). Unfortunately,
its performance in terms of coverage rates was found to be inferior to the
parametric and the wild bootstrap CIs (Modugno and Giannerini, 2013).
Moreover, with contaminated data, the proportion of outliers can be very high
in the bootstrap samples, thus biasing parameter estimation (Salibian-Barrera
and Yohai, 2006; Salibián-Barrera _et al._ , 2008).

Permutation tests are also non-parametric alternatives that, like the wild
bootstrap, have the advantage, compared to the cluster bootstrap, of
maintaining the same proportion of outliers in the resamples as in the
original sample. In permutation tests, the same observations are retained but
with permuted characteristic, such as the group labels in a two-sample
independent t-test or the predictor-response pairs in simple linear
regressions, as if the link between the predictor and the response variable
was null, to obtain an empirical distribution of the parameter (or statistic)
of interest under the null hypothesis. Despite the complexity of the repeated
measures data structure, permutation methods have been developed for classical
LMMs and made accessible, for example, through the lmm.perm function in the
minque package in R (Wu and Wu, 2019). These tests are well-suited for group
comparisons Shields-Cutler _et al._ (2018). Under the null hypothesis,
individuals from both groups come from the same population and are therefore
independent and identically distributed. It suffices to permute the group
labels. Several permutation methods for LMMs are based on the likelihood of
the model and are thus not applicable to robust estimators (e.g., Krzciuk,
2014; Lee and Braun, 2012; Rao _et al._ , 2019). In their first study, Krzciuk
(2014) also combined the Wald statistic with the permutation method. The Wald
statistic was computed in the original sample and in the permuted samples. The
$p$-value was defined as the proportion of Wald statistic values obtained from
the permuted samples that were greater than the one obtained from the original
sample. However, the results regarding the variances were significantly poorer
than the two other methods based on likelihood models, so this method was
discarded. Nevertheless, Krzciuk’s method appears to be a good alternative to
tests of fixed effects with robust estimators, especially considering the
excellent performance of the original test with this class of estimators.

We hope that with the confintROB package, LMM users who are wary of outliers
can have the best of both worlds: a robust estimation procedure for their LMM
estimates (via the robustlmm or robustvarComp packages) coupled with state-of-
the-art confidence intervals, as currently implemented in the stats and the
lme4 packages for classical estimators. We believe that being able to compare
parameter estimates as well as their CIs across classical and robust
estimators is a valuable tool when testing research hypotheses that could
easily go astray because of a few, odd, and yet influential observations.

## Computational details

The results in this paper were obtained using R 4.3.1 with the robustlmm 3.2-3
and lme4 1.1-34 packages. R itself and all packages used are available from
the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/.

## Acknowledgments

We thank Drs. Lucia Modugno and Simone Giannerini for helping us with their R
script to compute wild bootstraps.

## References

  * Agostinelli and Yohai (2016) Agostinelli C, Yohai VJ (2016).  “Composite Robust Estimators for Linear Mixed Models.”  _Journal of the American Statistical Association_ , 111(516), 1764–1774.  10.1080/01621459.2015.1115358.
  * Anderson (1961) Anderson NH (1961).  “Scales and statistics: parametric and nonparametric.”  _Psychological bulletin_ , 58(4), 305.  10.1037/h0042576.
  * Bates _et al._ (2015) Bates D, Mächler M, Bolker B, Walker S (2015).  “Fitting Linear Mixed-Effects Models Using lme4.”  _Journal of Statistical Software_ , 67(1), 1–48.  10.18637/jss.v067.i01.
  * Burger and Schall (2018) Burger DA, Schall R (2018).  “Robust fit of Bayesian mixed effects regression models with application to colony forming unit count in tuberculosis research.”  _Statistics in medicine_ , 37(4), 544–556.  10.1002/sim.7529.
  * Chervoneva and Vishnyakov (2014) Chervoneva I, Vishnyakov M (2014).  “Generalized S-estimators for linear mixed effects models.”  _Statistica Sinica_ , 24(3), 1257–1276.  ISSN 10170405.  10.5705/ss.2012.285.
  * Copt and Victoria-Feser (2006) Copt S, Victoria-Feser MP (2006).  “High-breakdown inference for mixed linear models.”  _Journal of the American Statistical Association_ , 101(473), 292–300.  10.1198/016214505000000772.
  * Deen and de Rooij (2020) Deen M, de Rooij M (2020).  “ClusterBootstrap: An R package for the analysis of hierarchical data using generalized linear models with the cluster bootstrap.”  _Behavior Research Methods_ , 52(2), 572–590.  10.3758/s13428-019-01252-y.
  * Efron (1979) Efron B (1979).  “Bootstrap methods: Another look at the jacknife.”  _The Annals of Statistics_ , 7(1), 1–26.  10.1214/aos/1176344552.
  * Kenward and Roger (1997) Kenward MG, Roger JH (1997).  “Small Sample Inference for Fixed Effects from Restricted Maximum Likelihood.”  _Biometrics_ , 53(3), 983–997.  ISSN 0006-341X.  10.2307/2533558.
  * Koller (2013) Koller M (2013).  _Robust Estimation of Linear Mixed Models_.  Ph.D. thesis, Eidgenössische Technische Hochschule Zurich.  10.3929/ethz-a-007632241.
  * Koller (2016) Koller M (2016).  “robustlmm: An R Package for Robust Estimation of Linear Mixed-Effects Models.”  _Journal of Statistical Software_ , 75(6), 1–24.  ISSN 1548-7660.  10.18637/jss.v075.i06.
  * Koller and Stahel (2023) Koller M, Stahel WA (2023).  “Robust Estimation of General Linear Mixed Effects Models.”  In PM Yi, PK Nordhausen (eds.), _Robust and Multivariate Statistical Methods_. Springer Nature Switzerland AG.
  * Krzciuk (2014) Krzciuk M (2014).  “On some tests of fixed effects for linear mixed models.”  _Studia Ekonomiczne_ , 189, 49–57.
  * Kuznetsova _et al._ (2017) Kuznetsova A, Brockhoff PB, Christensen RHB (2017).  “lmerTest package: tests in linear mixed effects models.”  _Journal of Statistical Software_ , 82(13).  10.18637/jss.v082.i13.
  * Laird and Ware (1982) Laird NM, Ware JH (1982).  “Random-effects models for longitudinal data.”  _Biometrics_ , pp. 963–974.  10.2307/2529876.
  * Lee and Braun (2012) Lee OE, Braun TM (2012).  “Permutation tests for random effects in linear mixed models.”  _Biometrics_ , 68(2), 486–493.  10.1111/j.1541-0420.2011.01675.x.
  * Lucadamo _et al._ (2021) Lucadamo A, Greco L, Amenta P, Crisci A (2021).  “Snipped robust inference in mixed linear models.”  _Preface XIX 1 Plenary Sessions_ , p. 695.
  * Mammen (1993) Mammen E (1993).  “Bootstrap and Wild Bootstrap for High Dimensional Linear Models.”  _The Annals of Statistics_ , 21(1), 255–285.  10.1214/aos/1176349025.
  * Mason _et al._ (2021) Mason F, Cantoni E, Ghisletta P (2021).  “Parametric and Semi-Parametric Bootstrap-Based Confidence Intervals for Robust Linear Mixed Models.”  _Methodology_ , 17(4), 271–295.  10.5964/meth.6607.
  * Mason _et al._ (2024) Mason F, Cantoni E, Ghisletta P (2024).  “Linear mixed models and latent growth curve models for group comparison studies contaminated by outliers.”  _Psychological Methods_ , Advance online publication.  10.1037/met0000643.
  * Modugno and Giannerini (2013) Modugno L, Giannerini S (2013).  “The Wild Bootstrap for Multilevel Models.”  _Communications in Statistics - Theory and Methods_ , 44(22), 4812–4825.  ISSN 0361-0926, 1532-415X.  10.1080/03610926.2013.802807.
  * Osorio (2019) Osorio F (2019).  _heavy: Robust estimation using heavy-tailed distributions_.
  * Pinheiro _et al._ (2001) Pinheiro J, Liu C, Wu YN (2001).  “Efficient Algorithms for Robust Estimation in Linear Mixed-Effects Models Using the Multivariate t Distribution.”  _Journal of Computational and Graphical Statistics_ , 10(2), 249–276.  ISSN 1061-8600, 1537-2715.  10.1198/10618600152628059.
  * Rao _et al._ (2019) Rao K, Drikvandi R, Saville B (2019).  “Permutation and Bayesian tests for testing random effects in linear mixed-effects models.”  _Statistics in medicine_ , 38(25), 5034–5047.
  * Richardson and Welsh (1995) Richardson AM, Welsh AH (1995).  “Robust restricted maximum likelihood in mixed linear models.”  _Biometrics_ , pp. 1429–1439.  10.2307/2533273.
  * Salibián-Barrera _et al._ (2008) Salibián-Barrera M, Van Aelst S, Willems G (2008).  “Fast and robust bootstrap.”  _Statistical Methods and Applications_ , 17(1), 41–71.  10.1007/s10260-007-0048-6.
  * Salibian-Barrera and Yohai (2006) Salibian-Barrera M, Yohai VJ (2006).  “A fast algorithm for S-regression estimates.”  _Journal of computational and Graphical Statistics_ , 15(2), 414–427.  10.1198/106186006X113629.
  * Saraceno _et al._ (2023) Saraceno G, Ghosh A, Basu A, Agostinelli C (2023).  “Robust estimation of fixed effect parameters and variances of linear mixed models: the minimum density power divergence approach.”  _AStA Advances in Statistical Analysis_ , pp. 1–31.
  * Satterthwaite (1946) Satterthwaite FE (1946).  “An approximate distribution of estimates of variance components.”  _Biometrics bulletin_ , 2(6), 110–114.  10.2307/3002019.
  * Shields-Cutler _et al._ (2018) Shields-Cutler RR, Al-Ghalith GA, Yassour M, Knights D (2018).  “Splinectomer enables group comparisons in longitudinal microbiome studies.”  _Frontiers in microbiology_ , 9, 785.  10.3389/fmicb.2018.00785.
  * Singer and Willett (2003) Singer JD, Willett JB (2003).  _Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence_.  Oxford University Press.  ISBN 978-0-19-515296-8.
  * Singmann _et al._ (2015) Singmann H, Bolker B, Westfall J, Aust F, Ben-Shachar MS (2015).  “afex: Analysis of factorial experiments.”  _R package version 0.28–1_.
  * Tibshirani and Efron (1993) Tibshirani RJ, Efron B (1993).  “An introduction to the bootstrap.”  _Monographs on statistics and applied probability_ , 57, 1–436.  10.1007/978-1-4899-4541-9_1.
  * Tomarken _et al._ (1997) Tomarken A, Shelton R, Elkins L, Anderson T (1997).  “Sleep deprivation and anti-depressant medication: unique effects on positive and negative affect.”  In _American Psychological Society Meeting, Washington, DC_.
  * Welsh and Richardson (1997) Welsh A, Richardson A (1997).  “13 Approaches to the robust estimation of mixed models.”  _Handbook of statistics_ , 15, 343–384.  10.1016/S0169-7161(97)15015-5.
  * Wu and Wu (2019) Wu J, Wu MJ (2019).  “Package ‘minque’.”  _Available online: cran. rproject. org/web/packages/minque/minque. pdf_.
  * Zhang _et al._ (2020) Zhang Y, Qin G, Zhu Z, Fu B (2020).  “Robust estimation of models for longitudinal data with dropouts and outliers.”  _Journal of Applied Statistics_ , pp. 1–24.  10.1080/02664763.2020.1845623.


* Johnson et al. (2018b) —. 2018b, ApJ, 865, 104.  https://arxiv.org/abs/1808.08966
  * Jorstad et al. (2022) Jorstad, S., et al. 2022, in prep.
  * Kim et al. (2020) Kim, J.-Y., Krichbaum, T. P., Broderick, A. E., et al. 2020, A&A, 640, A69, doi: 10.1051/0004-6361/202037493
  * Kostić et al. (2009) Kostić, U., Čadež, A., Calvani, M., & Gomboc, A. 2009, A&A, 496, 307, doi: 10.1051/0004-6361/200811059
  * Krichbaum et al. (1997) Krichbaum, T. P., Graham, D. A., Greve, A., et al. 1997, A&A, 323, L17
  * Krichbaum et al. (1998) Krichbaum, T. P., Graham, D. A., Witzel, A., et al. 1998, A&A, 335, L106
  * Lee et al. (2014) Lee, S.-S., Petrov, L., Byun, D.-Y., et al. 2014, AJ, 147, 77, doi: 10.1088/0004-6256/147/4/77
  * Li et al. (2013) Li, Z., Morris, M. R., & Baganoff, F. K. 2013, ApJ, 779, 154, doi: 10.1088/0004-637X/779/2/154
  * Liu et al. (2016) Liu, H. B., Wright, M. C. H., Zhao, J.-H., et al. 2016, A&A, 593, A107, doi: 10.1051/0004-6361/201628731
  * Liu et al. (2004) Liu, S., Petrosian, V., & Melia, F. 2004, ApJ, 611, L101, doi: 10.1086/423985
  * Lo & Claussen (1983) Lo, K. Y., & Claussen, M. J. 1983, Nature, 306, 647, doi: 10.1038/306647a0
  * Lo et al. (1998) Lo, K. Y., Shen, Z.-Q., Zhao, J.-H., & Ho, P. T. P. 1998, ApJ, 508, L61, doi: 10.1086/311726
  * Lu et al. (2018) Lu, R.-S., Krichbaum, T. P., Roy, A. L., et al. 2018, ApJ, 859, 60, doi: 10.3847/1538-4357/aabe2e
  * Macquart & Bower (2006) Macquart, J.-P., & Bower, G. C. 2006, ApJ, 641, 302, doi: 10.1086/500317
  * Marrone et al. (2006a) Marrone, D. P., Moran, J. M., Zhao, J.-H., & Rao, R. 2006a, ApJ, 640, 308, doi: 10.1086/500106
  * Marrone et al. (2006b) Marrone, D. P., Moran, J. M., Zhao, J.-H., & Rao, R. 2006b, in Journal of Physics Conference Series, Vol. 54, Journal of Physics Conference Series, 354–362, doi: 10.1088/1742-6596/54/1/056
  * Marrone et al. (2007) —. 2007, ApJ, 654, L57, doi: 10.1086/510850
  * Marrone et al. (2008) Marrone, D. P., Baganoff, F. K., Morris, M. R., et al. 2008, ApJ, 682, 373, doi: 10.1086/588806
  * Martí-Vidal et al. (2016) Martí-Vidal, I., Roy, A., Conway, J., & Zensus, A. J. 2016, A&A, 587, A143, doi: 10.1051/0004-6361/201526063
  * Matthews et al. (2018) Matthews, L. D., Crew, G. B., Doeleman, S. S., et al. 2018, PASP, 130, 015002, doi: 10.1088/1538-3873/aa9c3d
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Michail et al. (2021) Michail, J. M., Wardle, M., Yusef-Zadeh, F., & Kunneriath, D. 2021, ApJ, 923, 54, doi: 10.3847/1538-4357/ac2d2c
  * Miyazaki et al. (2004) Miyazaki, A., Tsutsumi, T., & Tsuboi, M. 2004, ApJ, 611, L97, doi: 10.1086/424004
  * Mori et al. (2013) Mori, K., Gotthelf, E. V., Zhang, S., et al. 2013, ApJ, 770, L23, doi: 10.1088/2041-8205/770/2/L23
  * Morris et al. (2012) Morris, M. R., Meyer, L., & Ghez, A. M. 2012, Research in Astronomy and Astrophysics, 12, 995, doi: 10.1088/1674-4527/12/8/007
  * Narayan (1992) Narayan, R. 1992, Royal Society of London Philosophical Transactions Series A, 341, 151, doi: 10.1098/rsta.1992.0090
  * Neilsen et al. (2013a) Neilsen, J., Nowak, M. A., Gammie, C., et al. 2013a, ApJ, 774, 42, doi: 10.1088/0004-637X/774/1/42
  * Neilsen et al. (2013b) —. 2013b, ApJ, 774, 42, doi: 10.1088/0004-637X/774/1/42
  * Neilsen et al. (2013c) —. 2013c, ApJ, 774, 42, doi: 10.1088/0004-637X/774/1/42
  * Neilsen et al. (2015) Neilsen, J., Markoff, S., Nowak, M. A., et al. 2015, ApJ, 799, 199, doi: 10.1088/0004-637X/799/2/199
  * Niinuma et al. (2015) Niinuma, K., Lee, S.-S., Kino, M., & Sohn, B. W. 2015, Publication of Korean Astronomical Society, 30, 637, doi: 10.5303/PKAS.2015.30.2.637
  * Nowak et al. (2012a) Nowak, M. A., Neilsen, J., Markoff, S. B., et al. 2012a, ApJ, 759, 95, doi: 10.1088/0004-637X/759/2/95
  * Nowak et al. (2012b) —. 2012b, ApJ, 759, 95, doi: 10.1088/0004-637X/759/2/95
  * Padin et al. (1990) Padin, S., Woody, D. P., Hodges, M. W., et al. 1990, ApJ, 360, L11, doi: 10.1086/185800
  * Peißker et al. (2021) Peißker, F., Zajaček, M., Eckart, A., et al. 2021, ApJ, 923, 69, doi: 10.3847/1538-4357/ac23df
  * Ponti et al. (2017) Ponti, G., George, E., Scaringi, S., et al. 2017, MNRAS, 468, 2447, doi: 10.1093/mnras/stx596
  * Porquet et al. (2008) Porquet, D., Grosso, N., Predehl, P., et al. 2008, A&A, 488, 549, doi: 10.1051/0004-6361:200809986
  * Primiani et al. (2016) Primiani, R. A., Young, K. H., Young, A., et al. 2016, Journal of Astronomical Instrumentation, 5, 1641006, doi: 10.1142/S2251171716410063
  * Psaltis et al. (2018) Psaltis, D., Johnson, M., Narayan, R., et al. 2018, ArXiv e-prints.  https://arxiv.org/abs/1805.01242
  * Quataert (2002) Quataert, E. 2002, ApJ, 575, 855, doi: 10.1086/341425
  * Rauch et al. (2016) Rauch, C., Ros, E., Krichbaum, T. P., et al. 2016, A&A, 587, A37, doi: 10.1051/0004-6361/201527286
  * Rea et al. (2013) Rea, N., Esposito, P., Pons, J. A., et al. 2013, ApJ, 775, L34, doi: 10.1088/2041-8205/775/2/L34
  * Rea et al. (2020) Rea, N., Coti Zelati, F., Viganò, D., et al. 2020, ApJ, 894, 159, doi: 10.3847/1538-4357/ab8387
  * Reid & Brunthaler (2004) Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872, doi: 10.1086/424960
  * Rickett (1990) Rickett, B. J. 1990, ARA&A, 28, 561, doi: 10.1146/annurev.aa.28.090190.003021
  * Roelofs et al. (2017) Roelofs, F., Johnson, M. D., Shiokawa, H., Doeleman, S. S., & Falcke, H. 2017, ApJ, 847, 55, doi: 10.3847/1538-4357/aa8455
  * Scargle (1998) Scargle, J. D. 1998, The Astrophysical Journal, 504, 405, doi: 10.1086/306064
  * Scargle et al. (2013) Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ApJ, 764, 167, doi: 10.1088/0004-637X/764/2/167
  * Schödel et al. (2011) Schödel, R., Morris, M. R., Muzic, K., et al. 2011, A&A, 532, A83, doi: 10.1051/0004-6361/201116994
  * Schwab & Cotton (1983) Schwab, F. R., & Cotton, W. D. 1983, AJ, 88, 688, doi: 10.1086/113360
  * Serabyn et al. (1997) Serabyn, E., Carlstrom, J., Lay, O., et al. 1997, ApJ, 490, L77, doi: 10.1086/311010
  * Shcherbakov et al. (2012) Shcherbakov, R. V., Penna, R. F., & McKinney, J. C. 2012, ApJ, 755, 133, doi: 10.1088/0004-637X/755/2/133
  * Shen et al. (2005) Shen, Z.-Q., Lo, K. Y., Liang, M.-C., Ho, P. T. P., & Zhao, J.-H. 2005, Nature, 438, 62, doi: 10.1038/nature04205
  * Smirnov (2011a) Smirnov, O. M. 2011a, A&A, 527, A106, doi: 10.1051/0004-6361/201016082
  * Smirnov (2011b) —. 2011b, A&A, 527, A107, doi: 10.1051/0004-6361/201116434
  * Smirnov (2011c) —. 2011c, A&A, 527, A108, doi: 10.1051/0004-6361/201116435
  * Smirnov (2011d) —. 2011d, A&A, 531, A159, doi: 10.1051/0004-6361/201116764
  * Steel et al. (2019) Steel, S., Wielgus, M., Blackburn, L., Issaoun, S., & Johnson, M. 2019, EHT Memo Series, 2019-CE-03
  * Stone et al. (2016) Stone, J. M., Marrone, D. P., Dowell, C. D., et al. 2016, ApJ, 825, 32, doi: 10.3847/0004-637X/825/1/32
  * Su & Finkbeiner (2012) Su, M., & Finkbeiner, D. P. 2012, ApJ, 753, 61, doi: 10.1088/0004-637X/753/1/61
  * Telesco et al. (1996) Telesco, C. M., Davidson, J. A., & Werner, M. W. 1996, ApJ, 456, 541, doi: 10.1086/176678
  * Thompson et al. (2017) Thompson, A. R., Moran, J. M., & Swenson, Jr., G. W. 2017, Interferometry and Synthesis in Radio Astronomy, 3rd Edition (Springer International Publishing), doi: 10.1007/978-3-319-44431-4
  * van Bemmel et al. (2019) van Bemmel, I., Small, D., Kettenis, M., et al. 2019, PoS(EVN2018)079, arXiv:1904.11747, doi: 10.22323/1.344.0079
  * van den Eijnden et al. (2021) van den Eijnden, J., A, A., B, B., & C, C. 2021, in prep., 00, 0, doi: TBD
  * van Langevelde et al. (1992) van Langevelde, H. J., Frail, D. A., Cordes, J. M., & Diamond, P. J. 1992, ApJ, 396, 686, doi: 10.1086/171750
  * Čadež et al. (2008) Čadež, A., Calvani, M., & Kostić, U. 2008, A&A, 487, 527, doi: 10.1051/0004-6361:200809483
  * Verner et al. (1996) Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487, doi: 10.1086/177435
  * Vertatschitsch et al. (2015) Vertatschitsch, L., Primiani, R., Young, A., et al. 2015, PASP, 127, 1226, doi: 10.1086/684513
  * von Fellenberg et al. (2018) von Fellenberg, S. D., Gillessen, S., Graciá-Carpio, J., et al. 2018, ApJ, 862, 129, doi: 10.3847/1538-4357/aacd4b
  * Wajima et al. (2016) Wajima, K., Hagiwara, Y., An, T., et al. 2016, Astronomical Society of the Pacific Conference Series, Vol. 502, The East-Asian VLBI Network, ed. L. Qain & D. Li, 81
  * Wang et al. (2013) Wang, Q. D., Nowak, M. A., Markoff, S. B., et al. 2013, Science, 341, 981, doi: 10.1126/science.1240755
  * Weisskopf et al. (2002) Weisskopf, M. C., Brinkman, B., Canizares, C., et al. 2002, PASP, 114, 1, doi: 10.1086/338108
  * Whitney et al. (2004) Whitney, A. R., Cappallo, R., Aldrich, W., et al. 2004, Radio Science, 39, RS1007, doi: 10.1029/2002RS002820
  * Whitney et al. (2013) Whitney, A. R., Beaudoin, C. J., Cappallo, R. J., et al. 2013, PASP, 125, 196, doi: 10.1086/669718
  * Wielgus et al. (2019) Wielgus, M., Blackburn, L., Issaoun, S., et al. 2019, EHT data set validation and characterization of errors, Tech. Rep. 2019-CE-02, EHT Memo Series
  * Wielgus et al. (2022) Wielgus, M., Marchili, N., Martí-Vidal, I., et al. 2022, in prep., 00, 0, doi: TBD
  * Williams et al. (2017) Williams, P. K. G., Clavel, M., Newton, E., & Ryzhkov, D. 2017, pwkit: Astronomical utilities in Python.  http://ascl.net/1704.001
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914, doi: 10.1086/317016
  * Witzel et al. (2012) Witzel, G., Eckart, A., Bremer, M., et al. 2012, ApJS, 203, 18, doi: 10.1088/0067-0049/203/2/18
  * Witzel et al. (2018a) Witzel, G., Martinez, G., Hora, J., et al. 2018a, ApJ, 863, 15, doi: 10.3847/1538-4357/aace62
  * Witzel et al. (2018b) —. 2018b, ApJ, 863, 15, doi: 10.3847/1538-4357/aace62
  * Witzel et al. (2021) Witzel, G., Martinez, G., Willner, S. P., et al. 2021, ApJ, 917, 73, doi: 10.3847/1538-4357/ac0891
  * Wright et al. (1987) Wright, M. C. H., Genzel, R., Güsten, R., & Jaffe, D. T. 1987, in American Institute of Physics Conference Series, Vol. 155, The Galactic Center, ed. D. C. Backer, 133–137, doi: 10.1063/1.36411
  * Young et al. (2016) Young, A., Primiani, R., Weintroub, J., et al. 2016, in Phased Array Systems and Technology (ARRAY, E1, doi: 10.1109/ARRAY.2016.7832576
  * Yuan et al. (2003) Yuan, F., Quataert, E., & Narayan, R. 2003, ApJ, 598, 301, doi: 10.1086/378716
  * Yusef-Zadeh et al. (2015) Yusef-Zadeh, F., Bushouse, H., Schödel, R., et al. 2015, ApJ, 809, 10, doi: 10.1088/0004-637X/809/1/10
  * Yusef-Zadeh et al. (1990) Yusef-Zadeh, F., Morris, M., & Ekers, R. D. 1990, Nature, 348, 45, doi: 10.1038/348045a0
  * Yusef-Zadeh et al. (2011) Yusef-Zadeh, F., Wardle, M., Miller-Jones, J. C. A., et al. 2011, ApJ, 729, 44, doi: 10.1088/0004-637X/729/1/44
  * Yusef-Zadeh et al. (2012) Yusef-Zadeh, F., Arendt, R., Bushouse, H., et al. 2012, ApJ, 758, L11, doi: 10.1088/2041-8205/758/1/L11
  * Zhang et al. (2017a) Zhang, S., Baganoff, F. K., Ponti, G., et al. 2017a, ApJ, 843, 96, doi: 10.3847/1538-4357/aa74e8
  * Zhang et al. (2017b) —. 2017b, ApJ, 843, 96, doi: 10.3847/1538-4357/aa74e8
  * Zhao et al. (2003) Zhao, J.-H., Young, K. H., Herrnstein, R. M., et al. 2003, ApJ, 586, L29, doi: 10.1086/374581
  * Zheng (2015) Zheng, W. 2015, in IAU General Assembly, Vol. 29, 2255896
  * Zhu et al. (2019) Zhu, Z., Li, Z., Morris, M. R., Zhang, S., & Liu, S. 2019, ApJ, 875, 44, doi: 10.3847/1538-4357/ab0e05
  * Zubovas et al. (2012) Zubovas, K., Nayakshin, S., & Markoff, S. 2012, Monthly Notices of the Royal Astronomical Society, 421, 1315, doi: 10.1111/j.1365-2966.2011.20389.x
  * Zylka & Mezger (1988) Zylka, R., & Mezger, P. G. 1988, A&A, 190, L25
  * Zylka et al. (1992) Zylka, R., Mezger, P. G., & Lesch, H. 1992, A&A, 261, 119
  * Zylka et al. (1995) Zylka, R., Mezger, P. G., Ward-Thompson, D., Duschl, W. J., & Lesch, H. 1995, A&A, 297, 83.  https://arxiv.org/abs/astro-ph/9410086

  *[EHT]: Event Horizon Telescope
  *[Sgr A∗]: Sagittarius A*
  *[VLBI]: very-long-baseline interferometry
  *[ALMA]: Atacama Large Millimeter/submillimeter Array
  *[APEX]: Atacama Pathfinder Experiment
  *[LMT]: Large Millimeter Telescope Alfonso Serrano
  *[JCMT]: James Clerk Maxwell Telescope
  *[SMA]: Submillimeter Array
  *[PV]: IRAM 30 m Telescope on Pico Veleta
  *[SMT]: Submillimeter Telescope
  *[SPT]: South Pole Telescope
  *[$S/N$]: signal-to-noise ratio
  *[HOPS]: Haystack Observatory Postprocessing System
  *[CASA]: Common Astronomy Software Applications
  *[SEFD]: system equivalent flux density
  *[SEFDs]: system equivalent flux density
  *[µas]: microarcseconds