# S-estimation in Linear Models with Structured Covariance Matrices *

Hendrik, P. Lopuhaä, Valérie Gares, Anne Ruiz-Gazen

► To cite this version:

Hendrik, P. Lopuhaä, Valérie Gares, Anne Ruiz-Gazen. S-estimation in Linear Models with Structured Covariance Matrices *. Annals of Statistics, 2023, 51 (6), pp.2415-2439. 10.1214/23-aos2334. hal-04397997

HAL Id: hal-04397997

https://hal.science/hal-04397997v1

Submitted on 16 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# S-estimation in Linear Models with Structured Covariance Matrices *

Hendrik Paul Lopuhaä¹, Valerie Gares², and Anne Ruiz-Gazen³

¹Delft University of Technology

²Institut National des Sciences Appliqués de Rennes

³Toulouse School of Economics

August 4, 2022

## Abstract

We provide a unified approach to S-estimation in balanced linear models with structured covariance matrices. Of main interest are S-estimators for linear mixed effects models, but our approach also includes S-estimators in several other standard multivariate models, such as multiple regression, multivariate regression, and multivariate location and scatter. We provide sufficient conditions for the existence of S-functionals and S-estimators, establish asymptotic properties such as consistency and asymptotic normality, and derive their robustness properties in terms of breakdown point and influence function. All the results are obtained for general identifiable covariance structures and are established under mild conditions on the distribution of the observations, which goes far beyond models with elliptically contoured densities. Some of our results are new and others are more general than existing ones in the literature. In this way this manuscript completes and improves results on S-estimation in a wide variety of multivariate models. We illustrate our results by means of a simulation study and an application to data from a trial on the treatment of lead-exposed children.

## 1 Introduction

Linear models are widely used and provide a versatile approach for analyzing correlated responses, such as longitudinal data, growth data or repeated measurements. In such models, each subject $i$, $i = 1, ..., n$, is observed at $k_i$ occasions, and the vector of responses $\mathbf{y}_i$ is assumed to arise from the model

$$ \mathbf{y}_i = \mathbf{X}_i\beta + \mathbf{u}_i, $$

where $\mathbf{X}_i$ is the design matrix for the $i$th subject and $\mathbf{u}_i$ is a vector whose covariance matrix can be used to model the correlation between the responses. One possibility is the linear mixed effects model, in which the random effects together with the measurement error yields a specific covariance structure depending on a vector $\theta$ consisting of some unknown covariance parameters. Other covariance structures may arise, for example if the $\mathbf{u}_i$ are the outcome of a time series, see e.g., [14] or [10], for different possible covariance structures.

Maximum likelihood estimation of $\beta$ and $\theta$ has been studied, e.g., in [12, 25, 15], see also [10, 7]. To be resistant against outliers, robust methods have been investigated for linear mixed effects models, e.g., in [22, 5, 4, 13, 1, 3]. This mostly concerns S-estimators, originally introduced in the multiple regression context by Rousseeuw and Yohai [27] and extended to multivariate location and scatter in [6, 17], to multivariate regression in [29], and to linear mixed effects models in [5, 13, 3].

<footnote>*This work has been partly supported by the French Agence Nationale de la Recherche through the Investments for the Future (Investissements d'Avenir) program, grant ANR-17-EURE-0010.</footnote>
S-estimators are well known smooth versions of the minimum volume ellipsoid estimator [26] that are highly resistant against outliers. As such, S-estimators have gained popularity as robust estimators, but they may also serve as initial estimators to further improve the efficiency. However, the theory about these estimators is far from complete, even in balanced models where the number of observed responses is the same for all subjects.

In view of this, we provide a unified approach to S-estimation in balanced linear models with structured covariance matrices, and postpone a unified approach for unbalanced models to a future paper. The balanced setup is already quite flexible and includes several specific multivariate statistical models. Of main interest are S-estimators for linear mixed effects models, but our approach also includes S-estimators in several other standard multivariate models, such as multiple regression, multivariate regression, and multivariate location and scatter. We provide sufficient conditions for the existence of S-functionals and S-estimators, establish their asymptotic properties, such as consistency and asymptotic normality, and derive their robustness properties in terms of breakdown point and influence function. All results are obtained for a large class of identifiable covariance structures, and are established under very mild conditions on the distribution of the observations, which goes far beyond models with elliptically contoured densities. In this way, some of our results are new and others are more general than existing ones in the literature.

Existence of S-estimators and S-functionals is established under mild conditions. Although existence of the estimators seems a basic requirement, such results are missing for instance for multivariate regression in [30] and for linear mixed effects models in [5, 3]. We obtain robustness properties for S-estimators, such as breakdown point and influence function, under mild conditions on collections of observations and under mild conditions on the distribution of the observations. High breakdown and a bounded influence function seem basic requirements for a robust method, but both properties are not available for linear mixed effects models [5, 3]. For multivariate regression [30], the influence function is only determined at distributions with an elliptical contoured density. Finally, we establish consistency and asymptotic normality for S-estimators under mild conditions on the distribution of the observations. A rigorous derivation is missing for multivariate regression [30], or is only available for observations from a normal distribution [27, 3].

We apply our asymptotic results, such as influence function and asymptotic normality, to the special case for which the distribution of the observations corresponds to an elliptically contoured density. In this way we retrieve earlier results found in [27, 17, 30]. Somewhat surprisingly, the asymptotic variances of our S-estimators for linear mixed effects models in which the response has an elliptically contoured density, differ from the ones found in [5]. We investigate this difference by means of a simulation study.

The paper is organized as follows. In Section 2, we explain the model in detail and provide some examples of standard multivariate models that are included in our setup. In Section 3 we define the S-estimator and S-functional and in Section 4 we give conditions under which they exist. In Section 5 we establish continuity of the S-functional, which is then used to obtain consistency of the S-estimator. Section 6 deals with the breakdown point. Section 7 provides the preparation for Sections 8 and 9 in which we obtain the influence function and establish asymptotic normality. Finally, in Section 10, we illustrate our results by means of a simulation and investigate the performance of our estimators by means of an application to data from a trial on the treatment of lead-exposed children. All proofs and some technical lemmas are put in an Appendix at the end of the paper. Other long and technical proofs are available as supplemental material [16].

## 2 Balanced models with structured covariances

We consider independent observations $(y_1, X_1), \dots, (y_n, X_n)$, for which we assume the following model

$$y_i = X_i\beta + u_i, \quad i = 1, \dots, n, \qquad (2.1)$$

where $y_i \in \mathbb{R}^k$ contains repeated measurements for the $i$-th subject, $\beta \in \mathbb{R}^q$ is an unknown parameter vector, $X_i \in \mathbb{R}^{k \times q}$ is a known design matrix, and $u_i \in \mathbb{R}^k$ are unobservable independent mean zero random vectors with covariance matrix $V \in \text{PDS}(k)$, the class of positive definite
symmetric $k \times k$ matrices. The model is balanced in the sense that all $\mathbf{y}_i$ have the same dimension. Furthermore, we consider a structured covariance matrix, that is, the matrix $\mathbf{V} = \mathbf{V}(\boldsymbol{\theta})$ is a known function of unknown covariance parameters combined in a vector $\boldsymbol{\theta} \in \mathbb{R}^l$. We first discuss some examples that are covered by this setup.

**Example 1.** An important case of interest is the (balanced) linear mixed effects model

$$ \mathbf{y}_i = \mathbf{X}_i\boldsymbol{\beta} + \mathbf{Z}\gamma_i + \boldsymbol{\epsilon}_i, \quad i = 1, \dots, n. \qquad (2.2) $$

This model arises from $\mathbf{u}_i = \mathbf{Z}\gamma_i + \boldsymbol{\epsilon}_i$, for $i = 1, \dots, n$, where $\mathbf{Z} \in \mathbb{R}^{k \times g}$ is known and $\gamma_i \in \mathbb{R}^g$ and $\boldsymbol{\epsilon}_i \in \mathbb{R}^k$ are independent mean zero random variables, with unknown covariance matrices $\mathbf{G}$ and $\mathbf{R}$, respectively. In this case $\mathbf{V}(\boldsymbol{\theta}) = \mathbf{Z}\mathbf{G}\mathbf{Z}^T + \mathbf{R}$ and $\boldsymbol{\theta} = (\text{vech}(\mathbf{G})^T, \text{vech}(\mathbf{R})^T)^T$, where

$$ \text{vech}(\mathbf{A}) = (a_{11}, \dots, a_{k1}, a_{22}, \dots, a_{kk}) \qquad (2.3) $$

is the unique $k(k+1)/2$-vector that stacks the columns of the lower triangle elements of a symmetric matrix $\mathbf{A}$. In full generality, the model is usually overparametrized and one may run into identifiability problems. A more feasible example is obtained by taking $\mathbf{R} = \sigma_0^2\mathbf{I}_k$, $\mathbf{Z} = [\mathbf{Z}_1 \cdots \mathbf{Z}_r]$ and $\gamma_i = (\gamma_{i1}, \dots, \gamma_{ir})^T$, where the $\mathbf{Z}_j$'s are known $k \times g_j$ design matrices and the $\gamma_{ij} \in \mathbb{R}^{g_j}$ are independent mean zero random variables with covariance matrix $\sigma_j^2\mathbf{I}_{g_j}$, for $j = 1, \dots, r$. This leads to

$$ \mathbf{y}_i = \mathbf{X}_i\boldsymbol{\beta} + \sum_{j=1}^{r} \mathbf{Z}_j\gamma_{ij} + \boldsymbol{\epsilon}_i, \quad i = 1, \dots, n, \qquad (2.4) $$

with $\mathbf{V}(\boldsymbol{\theta}) = \sum_{j=1}^r \sigma_j^2 \mathbf{Z}_j \mathbf{Z}_j^T + \sigma_0^2 \mathbf{I}_k$ and $\boldsymbol{\theta} = (\sigma_0^2, \sigma_1^2, \dots, \sigma_r^2)$.

**Example 2.** An example with an unstructured covariance is the multivariate linear regression model

$$ \mathbf{y}_i = \mathbf{B}^T\mathbf{x}_i + \mathbf{u}_i, \quad i = 1, \dots, n, \qquad (2.5) $$

where $\mathbf{B} \in \mathbb{R}^{q \times k}$ is a matrix of unknown parameters, $\mathbf{x}_i \in \mathbb{R}^q$ is known, and $\mathbf{u}_i$, for $i = 1, \dots, n$, are independent mean zero random variables with covariance matrix $\mathbf{V}(\boldsymbol{\theta}) = \mathbf{C} \in \text{PDS}(k)$. In this case, the vector of unknown covariance parameters is given by

$$ \boldsymbol{\theta} = \text{vech}(\mathbf{C}) = (c_{11}, \dots, c_{1k}, c_{22}, \dots, c_{kk})^T \in \mathbb{R}^{\frac{1}{2}k(k+1)}. \qquad (2.6) $$

The model can be obtained as a special case of (2.1), by taking $\mathbf{X}_i = \mathbf{x}_i^T \otimes \mathbf{I}_k$ and $\boldsymbol{\beta} = \text{vec}(\mathbf{B}^T)$, where $\otimes$ denotes the Kronecker product and $\text{vec}(\cdot)$ is the $k^2$-vector that stacks the columns of a matrix. Clearly, the linear multiple regression model is a special case with $k = 1$.

**Example 3.** Model (2.1) also includes examples, for which $\mathbf{u}_1, \dots, \mathbf{u}_n$ are generated from a time series. One example, is the case where $\mathbf{u}_i$ has a covariance matrix with elements

$$ v_{st} = \sigma^2 \rho^{|s-t|}, \quad s, t = 1, \dots, k. \qquad (2.7) $$

This arises when the $\mathbf{u}_i$'s are generated by an autoregressive process of order one. The vector of unknown covariance parameters is $\boldsymbol{\theta} = (\sigma^2, \rho) \in (0, \infty) \times [-1, 1]$. A general stationary process leads to

$$ v_{st} = \theta_{|s-t|+1}, \quad s, t = 1, \dots, k, \qquad (2.8) $$

in which case $\boldsymbol{\theta} = (\theta_1, \dots, \theta_k)^T \in \mathbb{R}^k$, where $\theta_{|s-t|+1}$ represents the autocovariance over lag $|s-t|$.

**Example 4.** Also the multivariate location-scale model can be obtained as a special case of (2.1), by taking $\mathbf{X}_i = \mathbf{I}_k$, the $k \times k$ identity matrix. In this case, $\boldsymbol{\beta} \in \mathbb{R}^k$ is the unknown location parameter and $\mathbf{V}(\boldsymbol{\theta})$ is the unstructured covariance matrix as in Example 2, with $\boldsymbol{\theta}$ as in (2.6).

Throughout the manuscript we will assume that the parameter $\boldsymbol{\theta}$ is identifiable in the sense that,

$$ \mathbf{V}(\boldsymbol{\theta}_1) = \mathbf{V}(\boldsymbol{\theta}_2) \quad \Rightarrow \quad \boldsymbol{\theta}_1 = \boldsymbol{\theta}_2. \qquad (2.9) $$

This is true for all models in Examples 2, 3 and 4. This may not be true in general for the linear mixed effects model in Example 1 with unknown $\text{vech}(\mathbf{G})$ and $\text{vech}(\mathbf{R})$. For linear mixed effects models in (2.4), identifiability of $\boldsymbol{\theta} = (\sigma_0^2, \sigma_1^2, \dots, \sigma_r^2)$ holds for particular choices of the design matrices $\mathbf{Z}_1, \dots, \mathbf{Z}_r$.
# 3 Definitions

We start by representing our observations as points in $\mathbb{R}^k \times \mathbb{R}^{kq}$ in the following way. For $r = 1, \dots, k$, let $\mathbf{x}_i^T$ denote the $r$-th row of the $k \times q$ matrix $\mathbf{X}$, so that $\mathbf{x}_r \in \mathbb{R}^q$. We represent the pair $\mathbf{s} = (\mathbf{y}, \mathbf{X})$ as an element in $\mathbb{R}^k \times \mathbb{R}^{kq}$ defined by $\mathbf{s}^T = (\mathbf{y}^T, \mathbf{x}_1^T, \dots, \mathbf{x}_k^T)$. In this way our observations can be represented as $\mathbf{s}_1, \dots, \mathbf{s}_n$, with $\mathbf{s}_i = (\mathbf{y}_i, \mathbf{X}_i) \in \mathbb{R}^k \times \mathbb{R}^{kq}$.

## 3.1 S-estimator

S-estimators are defined by means of a function $\rho: \mathbb{R} \to [0, \infty)$ that satisfies the following properties

(R1) $\rho$ is symmetric around zero with $\rho(0) = 0$ and $\rho$ is continuous at zero;

(R2) There exists a finite constant $c_0 > 0$, such that $\rho$ is non-decreasing on $[0, c_0]$ and constant on $[c_0, \infty)$; put $a_0 = \sup \rho$.

The S-estimator $\xi_n = (\beta_n, \theta_n)$ is defined as the solution to the following minimization problem

$$ \begin{array}{ll} \min_{\beta, \theta} & \det(\mathbf{V}(\boldsymbol{\theta})) \\ \text{subject to} & \frac{1}{n} \sum_{i=1}^{n} \rho \left( \sqrt{(\mathbf{y}_i - \mathbf{X}_i \boldsymbol{\beta})^T \mathbf{V}(\boldsymbol{\theta})^{-1} (\mathbf{y}_i - \mathbf{X}_i \boldsymbol{\beta})} \right) \le b_0, \end{array} \quad (3.1) $$

where the minimum is taken over all $\beta \in \mathbb{R}^q$ and $\theta \in \mathbb{R}^l$, such that $\mathbf{V}(\boldsymbol{\theta}) \in \text{PDS}(k)$, with $\rho$ satisfying (R1)-(R2).

The S-estimator defined by (3.1) for the setup in (2.1) includes several specific cases that have been considered in the literature. The original regression S-estimator introduced by Rousseeuw and Yohai [27] is obtained as a special case by taking $\mathbf{X}_i = \mathbf{x}_i^T$ a $1 \times q$ vector and $\mathbf{V}(\boldsymbol{\theta}) = \sigma^2 > 0$. S-estimators for multivariate location and scale, as considered in Davies [6] and Lopuhaä [17] can be obtained by taking $\mathbf{X}_i$ and $\mathbf{V}(\boldsymbol{\theta})$ as in Example 4. For the multivariate regression model in Example 2, S-estimators have been considered by Van Aelst and Willems [30]. Copt and Victoria-Feser [5] and Chervoneva and Vishnyakov [3] consider S-estimators for the parameters in the linear mixed effects model (2.4).

The constant $0 < b_0 < a_0$ in (3.1) can be chosen in agreement with an assumed underlying distribution. For the multivariate regression model in [30], it is assumed that $\mathbf{y}_i | \mathbf{X}_i$ has an elliptically contoured density of the form

$$ f_{\mu, \Sigma}(\mathbf{y}) = \det(\Sigma)^{-1/2} h ((\mathbf{y} - \boldsymbol{\mu})^T \Sigma^{-1} (\mathbf{y} - \boldsymbol{\mu})), \quad (3.2) $$

with $\boldsymbol{\mu} = \mathbf{X}_i\boldsymbol{\beta}$ and $\boldsymbol{\Sigma} = \mathbf{V}(\boldsymbol{\theta})$ and $h: [0, \infty) \to [0, \infty)$. For the linear mixed effects model in [5], it is assumed that $\mathbf{y}_i | \mathbf{X}_i$ has a multivariate normal distribution, which is a special case of (3.2) with $h(t) = (2\pi)^{-k/2} \exp(-t/2)$. When the underlying distribution corresponds to a density of the form (3.2), then a natural choice is $b_0 = \mathbb{E}_{\mathbf{0}, \mathbf{I}} \rho(\|\mathbf{z}\|)$, where $\mathbf{z}$ has density (3.2) with $(\boldsymbol{\mu}, \boldsymbol{\Sigma}) = (\mathbf{0}, \mathbf{I}_k)$. Finally, it should be emphasized that the ratio $b_0/a_0$ determines the breakdown point of the S-estimator (see Theorem 4), as well as its limiting variance (see Corollary 6). By choosing the constant $c_0$ in (R2) one then has to make a trade-off between robustness and efficiency.

Note that at this point we do not assume smoothness of $\rho$ or strict monotonicity on $[0, c_0]$. This means that (R1)-(R2) allow the function $\rho(d) = 1 - \mathbb{1}_{[-c_0, c_0]}(d)$, which corresponds to the minimum volume ellipsoid estimator in location-scale models (see [26]) and to the least median of squares estimator in linear regression models (see [28]). Indeed, with $\rho(d) = 1 - \mathbb{1}_{[-c_0, c_0]}(d)$, the S-estimator $(\beta_n, \theta_n)$ corresponds to the smallest cylinder

$$ C(\boldsymbol{\beta}, \boldsymbol{\theta}, c_0) = \{ (\mathbf{y}, \mathbf{X}) \in \mathbb{R}^k \times \mathbb{R}^{kq} : (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}(\boldsymbol{\theta})^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \le c_0^2 \} \quad (3.3) $$

that contains at least $n - nb_0$ points.
**Remark 3.1.** Clearly, the definition of the S-estimator in (3.1) has great similarities with the S-estimator for multivariate location and covariance (see [6] and [17]), defined as the solution $(\mathbf{t}_n, \mathbf{C}_n)$ to the minimization problem

$$ \begin{array}{ll} \min_{\mathbf{t}, \mathbf{C}} & \det(\mathbf{C}) \\ \text{subject to} & \frac{1}{n} \sum_{i=1}^{n} \rho \left( \sqrt{(\mathbf{y}_i - \mathbf{t})^T \mathbf{C}^{-1} (\mathbf{y}_i - \mathbf{t})} \right) \le b_0, \end{array} \qquad (3.4) $$

where the minimum is taken over all $\mathbf{t} \in \mathbb{R}^k$ and $\mathbf{C} \in \text{PDS}(k)$. Even more so, if all $\mathbf{X}_i$ are assumed to be equal to the same design matrix $\mathbf{X}$ of full rank, as was done in [5, 4]. However, there is a subtle, but important difference between minimization problems (3.4) and (3.1). The important difference is that in (3.4) we minimize over all positive definite symmetric $k \times k$ matrices $\mathbf{C}$, whereas in (3.1), we only minimize over positive definite symmetric $k \times k$ matrices $\mathbf{V}(\boldsymbol{\theta})$, which can arise as the image of the mapping $\boldsymbol{\theta} \mapsto \mathbf{V}(\boldsymbol{\theta})$. The latter collection is a subset of the other:

$$ \{ \mathbf{V}(\boldsymbol{\theta}) \in \text{PDS}(k) : \boldsymbol{\theta} \in \mathbb{R}^l \} \subset \text{PDS}(k), $$

and will typically be a strictly smaller subset. This means that the properties of $\mathbf{V}(\boldsymbol{\theta}_n)$ and $\mathbf{C}_n$ are related, but the properties of $\mathbf{V}(\boldsymbol{\theta}_n)$ cannot simply be derived from properties of $\mathbf{C}_n$, not even in the case where all $\mathbf{X}_i$ are equal to the same $\mathbf{X}$. In fact, this will lead to limiting covariances that differ from the ones found in [5], see Corollary 6.

## 3.2 S-functional

The concept of S-functional is needed to investigate local robustness properties of the corresponding S-estimator, such as the influence function (see Section 8). Let $\mathbf{s} = (\mathbf{y}, \mathbf{X})$ have a probability distribution $P$ on $\mathbb{R}^k \times \mathbb{R}^{kq}$. The S-functional $\xi(P) = (\beta(P), \theta(P))$ is defined as the solution to the following minimization problem:

$$ \begin{array}{ll} \min_{\boldsymbol{\beta}, \boldsymbol{\theta}} & \det(\mathbf{V}(\boldsymbol{\theta})) \\ \text{subject to} & \displaystyle \int \rho \left( \sqrt{(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}(\boldsymbol{\theta})^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})} \right) dP(\mathbf{y}, \mathbf{X}) \le b_0, \end{array} \qquad (3.5) $$

where the minimum is taken over all $\boldsymbol{\beta} \in \mathbb{R}^q$ and $\boldsymbol{\theta} \in \mathbb{R}^l$, such that $\mathbf{V}(\boldsymbol{\theta}) \in \text{PDS}(k)$, with $\rho$ satisfying (R1)-(R2).

As a special case, we obtain the S-estimator $\xi_n = (\beta_n, \theta_n)$ by taking $P = \mathbb{P}_n$, the empirical measure corresponding to the observations $(\mathbf{y}_1, \mathbf{X}_1), \dots, (\mathbf{y}_n, \mathbf{X}_n)$. In view of this connection, existence and consistency of solutions to (3.1) will follow from general results on the existence and the continuity of solutions to (3.5).

The definition of the S-functionals for the multivariate location-scale model given in Lopuhaä [17] and for the multivariate regression model given by Van Aelst and Willems [30] can be obtained as special cases of (3.5), by choosing $\mathbf{X}, \boldsymbol{\beta}$ and $\mathbf{V}(\boldsymbol{\theta})$ as in Examples 4 and 2, respectively. Copt and Victoria-Feser [5] do not pay attention to S-functionals or the influence function in the linear mixed effects model (2.4). However, S-functionals for linear mixed effects models can be also be obtained as a special case of (3.5), by choosing $\mathbf{X}, \boldsymbol{\beta}$ and $\mathbf{V}(\boldsymbol{\theta})$ as in Example 1.

## 4 Existence

We will first establish existence of the S-functional $\xi(P)$ defined by (3.5), under particular conditions on the probability measure $P$. As a consequence, this will also yield the existence of
the S-estimator, defined by (3.1). Recall that $(\mathbf{y}_1, \mathbf{X}_1), \dots, (\mathbf{y}_n, \mathbf{X}_n)$ are represented as points in $\mathbb{R}^k \times \mathbb{R}^{kq}$. Note however, that for linear models with intercept the first column of each $\mathbf{X}_i$ consists of 1's. This means that the points $(\mathbf{y}_i, \mathbf{X}_i)$ are concentrated in a lower dimensional subset of $\mathbb{R}^k \times \mathbb{R}^{kq}$. A similar situation occurs when all $\mathbf{X}_i$ are equal to the same design matrix, such as in [5]. In view of this, define $\mathcal{X} \subset \mathbb{R}^{kq}$ as the subset with the lowest dimension $p = \dim(\mathcal{X}) \le kq$ satisfying

$$P(\mathbf{X} \in \mathcal{X}) = 1. \qquad (4.1)$$

Hence, $P$ is then concentrated on the subset $\mathbb{R}^k \times \mathcal{X}$ of $\mathbb{R}^k \times \mathbb{R}^{kq}$, which is of dimension $k+p$, which may be of smaller than $k+qk$.

The first condition we require, expresses the fact that $P$ cannot have too much mass at infinity, in relation to the ratio $r = b_0/a_0$.

(C1$_{\epsilon}$) There exists a compact set $K_{\epsilon} \subset \mathbb{R}^k \times \mathcal{X}$, such that $P(K_{\epsilon}) \ge r + \epsilon$.

The second condition requires that $P$ cannot have too much mass at arbitrarily thin strips in $\mathbb{R}^k \times \mathcal{X}$. For $\alpha \in \mathbb{R}^{k+qk}$, such that $\|\alpha\| = 1$, $\ell \in \mathbb{R}$ and $\delta \ge 0$, we define a strip $H(\alpha, \ell, \delta)$ as follows:

$$H(\alpha, \ell, \delta) = \{\mathbf{s} \in \mathbb{R}^k \times \mathbb{R}^{kq} : \ell - \delta/2 \le \alpha^T \mathbf{s} \le \ell + \delta/2\}. \qquad (4.2)$$

Defined in this way, a strip is the area between two parallel hyperplanes which are symmetric around the hyperplane $H(\alpha, \ell, 0) = \{\mathbf{s} \in \mathbb{R}^k \times \mathbb{R}^{kq} : \alpha^T \mathbf{s} = \ell\}$. Since the distance between two parallel hyperplanes $\alpha^T \mathbf{s} = \ell_1$ and $\alpha^T \mathbf{s} = \ell_2$ is $|\ell_1 - \ell_2|$, the strip $H(\alpha, \ell, \delta)$ defined as in (4.2) has width $\delta$. We require the following condition

(C2$_{\epsilon}$) The value

$$\delta_{\epsilon} = \inf \{ \delta : P(H(\alpha, \ell, \delta)) \ge \epsilon, \alpha \in \mathbb{R}^{k+qk}, \|\alpha\| = 1, \ell \in \mathbb{R}, \delta \ge 0 \}$$

is strictly positive.

According to (4.1), in (C2$_{\epsilon}$) one only needs to consider strips in $\mathbb{R}^k \times \mathcal{X}$.

Both conditions are satisfied for any $0 < \epsilon \le 1 - r$ by any probability measure $P$ that is absolutely continuous. Clearly, condition (C1$_{\epsilon}$) holds for any $0 \le \epsilon \le 1 - r$ for the empirical measure $P_n$ corresponding to a collection of $n$ points $S_n = \{\mathbf{s}_1, \dots, \mathbf{s}_n\} \subset \mathbb{R}^k \times \mathcal{X}$. Condition (C2$_{\epsilon}$) for $\epsilon = (k+p+1)/n$ is also satisfied by the empirical measure $P_n$, when the collection $S_n$ is in general position, i.e., no subset $J \subset S_n$ of $k+p+1$ points is contained in the same hyperplane in $\mathbb{R}^k \times \mathcal{X}$. Conditions (C1$_{\epsilon}$) and (C2$_{\epsilon}$) together, are similar to condition (C$_{\epsilon}$) in [17]. The reason that (C1$_{\epsilon}$) slightly deviates from [17], is to handle the presence of $\mathbf{X}$ in minimization problem (3.5).

**Remark 4.1.** Note that condition (C2$_{\epsilon}$) is equivalent with

$$\omega_{\epsilon} = \inf_{P(J) \ge \epsilon} \inf_{\|\alpha\|=1} \inf_{\ell \in \mathbb{R}} \sup_{\mathbf{s} \in J} |\alpha^T \mathbf{s} - \ell| > 0, \qquad (4.3)$$

where the infima are taken over all subsets $J \subset \mathbb{R}^k \times \mathcal{X}$ with $P(J) \ge \epsilon$, all vectors $\alpha \in \mathbb{R}^{k+qk}$, with $\|\alpha\| = 1$, and levels $\ell \in \mathbb{R}$. Details can be found in [16].

To establish existence of the S-functional, we follow the reasoning in [17]. The idea is to argue that one can restrict oneself to a compact set for finding solutions to (3.5). When the object function in (3.5) is continuous, this immediately yields existence of a solution of (3.5). To this end, we assume the following condition.

(V1) The mapping $\theta \mapsto \mathbf{V}(\theta)$ is continuous.

The lemma below is fundamental for the existence of the S-functional. It requires that the identity is in $\mathcal{V} = \{\mathbf{V}(\theta) \in \text{PDS}(k) : \theta \in \mathbb{R}^l\}$ and that $\mathcal{V}$ is closed under multiplication with a positive scalar.
(V2) There exists a $\boldsymbol{\theta} \in \mathbb{R}^l$, such that $\mathbf{V}(\boldsymbol{\theta}) = \mathbf{I}_k$. For any $\mathbf{V}(\boldsymbol{\theta}) \in \mathcal{V}$ and any $\alpha > 0$, it holds that $\alpha \mathbf{V}(\boldsymbol{\theta}) = \mathbf{V}(\boldsymbol{\theta}')$, for some $\boldsymbol{\theta}' \in \mathbb{R}^l$.

Conditions (V1)-(V2) are not very restrictive. For example, all models in Examples 1 to 4 satisfy these conditions.

For any $k \times k$ matrix $\mathbf{A}$, let $\lambda_k(\mathbf{A}) \le \dots \le \lambda_1(\mathbf{A})$ denote the eigenvalues of $\mathbf{A}$. We then have the following key lemma for the existence of S-functionals. The lemma is similar to Lemma 1 in [17] and its proof can be found in [16].

**Lemma 1.** Let $(\boldsymbol{\beta}, \boldsymbol{\theta}) \in \mathbb{R}^q \times \mathbb{R}^l$, $0 < m_0 < \infty$, $0 < c < \infty$, and $0 < \epsilon < 1$, and suppose that the mapping $\boldsymbol{\theta} \mapsto \mathbf{V}(\boldsymbol{\theta})$ satisfies (V2). Then the following properties hold.

(i) If $P$ satisfies $(C2_\epsilon)$ and $P(C(\boldsymbol{\beta}, \boldsymbol{\theta}, c)) \ge \epsilon$, then $\lambda_k(\mathbf{V}(\boldsymbol{\theta})) \ge a_1 > 0$, where $a_1$ only depends on $c$ and the width $\delta_\epsilon$ from condition (C2$_\epsilon$).

(ii) Suppose $\int \rho(\|\mathbf{y}\|/m_0) dP(\mathbf{s}) \le b_0$. Then for any solution $(\boldsymbol{\beta}, \boldsymbol{\theta})$ of (3.5), which is such that $\lambda_k(\mathbf{V}(\boldsymbol{\theta})) \ge a_1 > 0$, it holds that $\lambda_1(\mathbf{V}(\boldsymbol{\theta})) \le a_2 < \infty$, where $a_2$ only depends on $a_1$ and $m_0$.

(iii) Let $P$ satisfy (C2$_\epsilon$) and suppose that $P(C(\boldsymbol{\beta}, \boldsymbol{\theta}, c)) \ge a > 0$. Suppose there exists a compact set $K \in \mathbb{R}^k \times \mathcal{X}$, such that $P(K) \ge 1 - a + \epsilon$. If $\lambda_1(\mathbf{V}(\boldsymbol{\theta})) \le a_2 < \infty$, then $\|\boldsymbol{\beta}\| \le M < \infty$, where $M$ only depends on $c$, $a_2$, the set $K$, and a constant $\gamma_\epsilon > 0$ that can be deduced from condition (C2$_\epsilon$).

Lemma 1 will ensure that there exists a compact set that contains all pairs $(\boldsymbol{\beta}, \mathbf{V}(\boldsymbol{\theta}))$ that correspond to possible solutions $(\boldsymbol{\beta}, \boldsymbol{\theta})$ of (3.5). To establish that possible solutions $(\boldsymbol{\beta}, \boldsymbol{\theta})$ of (3.5) are in a compact set, we need that the pre-image $\{\boldsymbol{\theta} \in \mathbb{R}^l : \mathbf{V}(\boldsymbol{\theta}) \in K\}$ of a compact set $K \subset \mathbb{R}^{k \times k}$ is again compact. Recall that subsets of $\mathbb{R}^l$ are compact if and only if they are closed and bounded, and note that the pre-image of a continuous mapping of a closed set is closed. Hence, in view of condition (V1), it suffices to require the following condition.

(V3) The mapping $\boldsymbol{\theta} \mapsto \mathbf{V}(\boldsymbol{\theta})$ is such that the pre-image of a bounded set is bounded.

Condition (V3) is satisfied by all models in Examples 1 to 4, including the linear mixed effects model of Example 1, as long as the matrix $\mathbf{Z}$ is of full rank. We then have the following theorem.

**Theorem 1.** Consider minimization problem (3.5) with $\rho$ satisfying (R1)-(R2). Suppose that $P$ satisfies (C1$_\epsilon$) and (C2$_\epsilon$), for some $0 < \epsilon \le 1 - r$, where $r = b_0/a_0$, and suppose that $\mathbf{V}$ satisfies (V1)-(V3). Then there exists at least one solution to (3.5).

The theorem has a direct corollary for the existence of the S-estimator, when dealing with a collections of points. Let $S_n = \{\mathbf{s}_1, \dots, \mathbf{s}_n\}$, with $\mathbf{s}_i = (\mathbf{y}_i, \mathbf{X}_i)$ be a collection of $n$ points in $\mathbb{R}^k \times \mathcal{X}$. Define

$$ \kappa(S_n) = \text{maximal number of points of } S_n \text{ lying on the same hyperplane in } \mathbb{R}^k \times \mathcal{X}. \quad (4.4) $$

For example, if the distribution $P$ is absolutely continuous, then $\kappa(S_n) \le k + p$ with probability one. We then have the following corollary.

**Corollary 1.** Consider minimization problem (3.1) with $\rho$ satisfying (R1)-(R2), for a collection $S_n = \{\mathbf{s}_1, \dots, \mathbf{s}_n\} \subset \mathbb{R}^k \times \mathcal{X}$, with $\mathbf{s}_i = (\mathbf{y}_i, \mathbf{X}_i)$, for $i = 1, \dots, n$. Suppose that $\mathbf{V}$ satisfies (V1)-(V3). If $\kappa(S_n) + 1 \le n(1-r)$, where $r = b_0/a_0$, then there exists at least one solution $\xi_n = (\boldsymbol{\beta}_n, \boldsymbol{\theta}_n)$ to the minimization problem (3.1).

Copt and Victoria-Feser [5] consider S-estimators for the linear mixed effects model (2.4). Despite their Proposition 1 about the asymptotic behavior of solutions to their S-minimization problem [5, equation (7)], the actual existence of such a solution is not established. However, this now follows from our Corollary 1. In their case, $\mathbf{V}(\boldsymbol{\theta})$ satisfies conditions (V1) and (V2). It can be seen, that if all matrices $\mathbf{Z}_j$, for $j = 1, \dots, r$, are of full rank, then $\mathbf{V}(\boldsymbol{\theta})$ also satisfies (V3). The translated bi-weight $\rho$-function proposed in [5] satisfies (R1)-(R2). Finally, under their assumption
that $\mathbf{X}_i = \mathbf{X}$ is the same and $\mathbf{y}_i | \mathbf{X} \sim N_k(\mathbf{X}\beta, \mathbf{V}(\theta))$, it follows that $\kappa(\mathcal{S}_n) \le k$. It then follows from Corollary 1 that with $b_0 \le a_0(n - k - 1)/n$, at least one solution to their S-minimization problem exists.

For the multivariate regression model from Example 2, Van Aelst & Willems [30] do not explicitly prove existence of the S-estimator. Since in their case, $\mathbf{V}(\theta) = \mathbf{C} \in \text{PDS}(k)$ satisfies (V1)-(V3) and the conditions imposed in [30] on the $\rho$-function satisfy (R1)-(R2), the existence of their S-estimator now also follows from Corollary 1, when $b_0$ is chosen suitably.

Existence of S-estimators is obtained from existence of S-functionals at the empirical measure $\mathbb{P}_n$. The following corollary shows that existence can be established in general, for probability measures that are close to $P$. It requires the following condition on $P$.

(C3) Let $\mathcal{C}$ be the class of all measurable convex subsets of $\mathbb{R}^k \times \mathbb{R}^{kq}$. Every $C \in \mathcal{C}$ is a $P$-continuity set, i.e., $P(\partial C) = 0$, where $\partial C$ denotes the boundary of $C$.

**Corollary 2.** Suppose that $\rho$ satisfies (R1)-(R2) and $\mathbf{V}$ satisfies (V1)-(V3). Let $P_t, t \ge 0$ be a sequence of probability measures on $\mathbb{R}^k \times \mathbb{R}^{kq}$ that converges weakly to $P$, as $t \to \infty$. Suppose that $P$ satisfies (C3), as well as $(C1_{\epsilon'})$ and $(C2_{\epsilon})$, for some $0 < \epsilon < \epsilon' \le 1 - r = b_0/a_0$. Then, for $t$ sufficiently large, the minimization problem (3.5) with probability measure $P_t$ has at least one solution $\xi(P_t)$.

Condition (C3) is needed to apply (A.2). Clearly, this condition is satisfied if $P$ is absolutely continuous.

# 5 Continuity and consistency

Consider a sequence $P_t, t \ge 0$, of probability measures on $\mathbb{R}^k \times \mathbb{R}^{kq}$ that converges weakly to $P$, as $t \to \infty$. By continuity of the S-functional $\xi(P)$ we mean that $\xi(P_t) \to \xi(P)$, as $t \to \infty$. An example of such a sequence is the sequence of empirical measures $\mathbb{P}_n, n = 1, 2, \dots$, that converges weakly to $P$, almost surely. Continuity of the S-functional for this sequence would then mean that the S-estimator $\xi_n$ is consistent, i.e., $\xi_n = \xi(\mathbb{P}_n) \to \xi(P)$, almost surely.

We require an additional condition for the function $\rho$.

(R3) $\rho$ is continuous and strictly increasing on $[0, c_0]$.

For $\mathbf{s} = (\mathbf{y}, \mathbf{X})$ and $\xi = (\beta, \theta)$, define the Mahalanobis distances by

$$d^2(\mathbf{s}, \xi) = d^2(\mathbf{s}, \beta, \theta) = (\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}(\theta)^{-1} (\mathbf{y} - \mathbf{X}\beta). \quad (5.1)$$

We then have the following theorem for the S-functional $\xi(P) = (\beta(P), \theta(P))$.

**Theorem 2.** Let $P_t, t \ge 0$ be a sequence of probability measures on $\mathbb{R}^k \times \mathbb{R}^{kq}$ that converges weakly to $P$, as $t \to \infty$, and let $\xi(P_t)$ be a solution to minimization problem (3.5) with probability measure $P_t$. Suppose that $\rho$ satisfies (R1)-(R3) and $\mathbf{V}$ satisfies (V1)-(V3). Suppose that $P$ satisfies (C3), as well as $(C1_{\epsilon'})$ and $(C2_{\epsilon})$, for some $0 < \epsilon < \epsilon' \le 1 - r = b_0/a_0$. If the solution $\xi(P)$ of (3.5) is unique, then for any sequence of solutions $\xi(P_t), t \ge 0$, it holds

$$\lim_{t \to \infty} \xi(P_t) = \xi(P).$$

Theorem 2 is an extension of Theorem 3.1 in [17] on the continuity of S-functionals for multivariate location and scale. Continuity of S-functionals for multiple regression has been investigated in [9].

Continuity of the S-functional will be used to derive the influence function of the S-estimator in Section 8. Another nice consequence of the continuity of the S-functional is, that one can directly obtain consistency of the S-estimator. Consider the S-estimator $\xi_n$ defined by minimization problem (3.1). Recall that $\xi_n = \xi(\mathbb{P}_n)$, so that we can use Theorem 2 to establish consistency of the S-estimator.
**Corollary 3.** Let $\xi_n$ be a solution to minimization problem (3.1). Suppose $\rho$ satisfies (R1)-(R3) and $\mathbf{V}$ satisfies (V1)-(V3). Suppose that $P$ satisfies (C3) as well as (C1$_{\epsilon}$') and (C2$_{\epsilon}$), for some $0 < \epsilon < \epsilon' \le 1 - r = b_0/a_0$. If the solution $\xi(P)$ of (3.5) is unique, then

$$ \lim_{n \to \infty} \xi_n = \xi(P), $$

with probability one.

Theorem 2 and Corollary 3 require that $\xi(P)$ is the unique solution to minimization problem (3.5). An example of a distribution $P$ for which $\xi(P)$ is unique, is when $P$ is such that $\mathbf{y} | \mathbf{X}$ has an elliptically contoured density (3.2). This situation is very similar to that of multivariate location-scale S-estimators, for which Davies [6, Theorem 1] shows that the corresponding S-minimization problem (3.5) has a unique solution. The next theorem is a direct consequence of that result. Its proof can be found in [16].

**Theorem 3.** Suppose that $\rho: \mathbb{R} \to [0, \infty)$ satisfies (R1)-(R2) and suppose that the probability distribution $P$ of $(\mathbf{y}, \mathbf{X})$ is such that $\mathbf{y} | \mathbf{X}$ has an elliptically contoured density $f_{\mu, \Sigma}$ from (3.2), with $\mu = \mathbf{X}\beta_0$ and $\Sigma = \mathbf{V}(\theta_0)$. Suppose that $h$ in (3.2) is non-increasing and such that the functions $-\rho$ and $h$ have at least one common point of decrease $d_0 > 0$, i.e.,

$$ \rho(s) < \rho(d_0) < \rho(t) \quad \text{and} \quad h(s) > h(d_0) > h(t) $$

for all $s, t \ge 0$, such that $s < d_0 < t$. If $\mathbf{X}^T\mathbf{X}$ is non-singular with probability one, then the minimization problem

$$ \begin{array}{ll} \min_{\beta, \theta} & \det(\mathbf{V}(\theta)) \\ \text{subject to} & \\ & \int \rho \left( \sqrt{(\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}(\theta)^{-1} (\mathbf{y} - \mathbf{X}\beta)} \right) f_{\mu, \Sigma}(\mathbf{y}) \mathrm{d}\mathbf{y} \le b_0, \end{array} \qquad (5.2) $$

where the minimum is taken over all $\beta \in \mathbb{R}^q$ and $\theta \in \mathbb{R}^l$, such that $\mathbf{V}(\theta) \in \text{PDS}(k)$, has the unique solution $(\beta, \theta) = (\beta_0, \theta_0)$ with probability one.

Minimization problem (5.2) seems to be slightly different from the one in (3.5). However, note that when $P$ is such that $\mathbf{X}$ is equal to a single value with probability one, both minimization problems are identical. This situation was considered, e.g., in [5].

An elliptically contoured density for $\mathbf{y}_i | \mathbf{X}_i$ in the context of S-estimators for specific cases of the model (2.1) has been assumed in [6] for the multivariate location-scale model of Example 4, in [30] for the multivariate regression model of Example 2, and in [5] for the linear mixed effects model (2.4). More precisely, in [5] it is assumed that $\mathbf{X}_i = \mathbf{X}$ and that $\mathbf{y}_i | \mathbf{X}$ has a multivariate normal distribution. In that case, the function $h$ in (3.2) satisfies all the conditions of Theorem 3.

# 6 Global robustness: the breakdown point

Consider a collection of points $S_n = \{\mathbf{s}_i = (\mathbf{y}_i, \mathbf{X}_i), i = 1, \dots, n\} \subset \mathbb{R}^k \times \mathcal{X}$. To emphasize the dependence on the collection $S_n$, denote by $\xi_n(S_n) = (\beta_n(S_n), \theta_n(S_n))$, the S-estimator, as defined in (3.1). To investigate the global robustness of S-estimators, we compute that finite-sample (replacement) breakdown point. For a given collection $S_n$ the finite-sample breakdown point (see Donoho and Huber [8]) of a regression S-estimator $\beta_n$ is defined as the smallest proportion of points from $S_n$ that one needs to replace in order to carry the estimator over all bounds. More precisely,

$$ \epsilon_n^*(\beta_n, S_n) = \min_{1 \le m \le n} \left\{ \frac{m}{n} : \sup_{S'_m} \|\beta_n(S_n) - \beta_n(S'_m)\| = \infty \right\}, \qquad (6.1) $$
where the minimum runs over all possible collections $S'_m$ that can be obtained from $S_n$ by replacing $m$ points of $S_n$ by arbitrary points in $\mathbb{R}^k \times \mathcal{X}$.

The estimator $\theta_n$ determines the covariance estimator $V_n = V(\theta_n)$. For this reason it seems natural to let the breakdown point of $\theta_n$ correspond to the breakdown of a covariance estimator. We define the finite sample (replacement) breakdown point of the S-estimator $\theta_n$ at a collection $S_n$, as

$$ \epsilon_n^*(\theta_n, S_n) = \min_{1 \le m \le n} \left\{ \frac{m}{n} : \sup_{S'_m} \text{dist}(V(\theta_n(S_n))), V(\theta_n(S'_m)) = \infty \right\}, \quad (6.2) $$

with $\text{dist}(\cdot, \cdot)$ defined as $\text{dist}(\mathbf{A}, \mathbf{B}) = \max\{|\lambda_1(\mathbf{A}) - \lambda_1(\mathbf{B})|, |\lambda_k(\mathbf{A})^{-1} - \lambda_k(\mathbf{B})^{-1}|\}$, where the minimum runs over all possible collections $S'_m$ that can be obtained from $S_n$ by replacing $m$ points of $S_n$ by arbitrary points in $\mathbb{R}^k \times \mathcal{X}$. So the breakdown point of $\theta_n$ is the smallest proportion of points from $S_n$ that one needs to replace in order to make the largest eigenvalue of $V(\theta(S'_m))$ arbitrarily large (explosion), or to make the smallest eigenvalue of $V(\theta(S'_m))$ arbitrarily small (implosion).

Good global robustness is illustrated by a high breakdown point. The breakdown point of the S-estimators is given the theorem below. It extends the results for S-estimators of multivariate location and scale, see [6] and [19], and S-estimators for multivariate regression, see [30]. For S-estimators in the linear mixed effects model considered in [5], the breakdown point has not been established. This will now follow as a special case from the next theorem. Its proof can be found in [16].

**Theorem 4.** Consider the minimization problem (3.1) with $\rho$ satisfying (R1)-(R2). Suppose that $V$ satisfies (V1)-(V3). Let $S_n \subset \mathbb{R}^k \times \mathcal{X}$ be a collection of $n$ points $\mathbf{s}_i = (\mathbf{y}_i, \mathbf{X}_i)$, $i = 1, \dots, n$. Let $r = b_0/a_0$ and suppose that $0 < r \le (n - \kappa(S_n))/(2n)$, where $\kappa(S_n)$ is defined by (4.4). Then for any solution $(\beta_n, \theta_n)$ of minimization problem (3.5),

$$ \frac{\lfloor (n+1)/2 \rfloor}{n} \ge \epsilon_n^*(\beta_n, S_n) \ge \frac{\lfloor nr \rfloor}{n}, \\ \epsilon_n^*(\theta_n, S_n) = \frac{\lfloor nr \rfloor}{n}. $$

The largest possible value of the breakdown point occurs when $r = (n - \kappa(S_n))/(2n)$, in which case $\lfloor nr \rfloor/n = \lfloor (n - \kappa(S_n))/2 \rfloor/n = \lfloor (n - \kappa(S_n) + 1)/2 \rfloor/n$. When the collection $S_n$ is in general position, then $\kappa(S_n) = k + p$. In that case the breakdown point of both estimators is at least equal to $\lfloor (n - k - p + 1)/2 \rfloor/n$. When all $\mathbf{X}_i$ are equal to the same $\mathbf{X}$, in [5, 4], one has $p = 0$ and $\kappa(S_n) = k$. In that case, the breakdown point of $\theta_n$ is equal to $\lfloor (n - k + 1)/2 \rfloor/n$. This coincides with the maximal breakdown point for affine equivariant estimators for $k \times k$ covariance matrices (see [6, Theorem 6]).

**Remark 6.1.** Van Aelst & Willems [30] also take into account the case $r > (n - \kappa(S_n))/(2n)$. For this case, by replacing $\lfloor n - nr \rfloor - \kappa(S_n)$ points, a specific solution to the S-minimization problem is constructed that breaks down. However, since there may be multiple solutions to the S-minimization problem, this does not necessarily mean that all solutions break down. In the proof of our Theorem 4, for the case $r \le (n - \kappa(S_n))/(2n)$, we show that all solutions to (3.1) do not break down, when replacing at most $\lfloor nr \rfloor - 1$ points, and that the covariance part of all solutions do break down, when replacing $\lfloor nr \rfloor$ points. For the case $r > (n - \kappa(S_n))/(2n)$, we can show that all solutions to (3.1) do not break down, when replacing at most $\lfloor n - nr \rfloor - \kappa(S_n) - 1$ points.

# 7 Score equations

Up to this point, properties of S-functionals and S-estimators have been derived from the minimization problems (3.1) and (3.5). To obtain the influence function and to establish the limiting distribution of S-estimators, we use the score equations that can be found by differentiation of the Lagrangian corresponding to the constrained minimization problems. To this end, we require the following additional condition on the function $\rho$,
(R4) $\rho$ is continuously differentiable and $u(s) = \rho'(s)/s$ is continuous,
and the following condition on the mapping $\theta \mapsto \mathbf{V}(\theta)$,

(V4) $\mathbf{V}(\theta)$ is continuously differentiable.

Obviously, condition (V4) implies the former condition (V1).

## 7.1 General covariance structures

Let $\xi_P = (\beta_P, \theta_P)$ be a solution to minimization problem (3.5). If we denote the corresponding Lagrange multiplier by $\lambda_P$, then the pair $(\xi_P, \lambda_P)$ is a zero of all partial derivatives $\partial L_P/\partial \beta$, $\partial L_P/\partial \theta$, and $\partial L_P/\partial \lambda$, where $L_P$ is the Lagrangian given by

$$L_P(\xi, \lambda) = \log \det(\mathbf{V}(\theta)) - \lambda \left\{ \int \rho \left( \sqrt{(\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}(\theta)^{-1} (\mathbf{y} - \mathbf{X}\beta)} \right) dP(\mathbf{y}, \mathbf{X}) - b_0 \right\}.$$

If $\mathbb{E}_P \|\mathbf{X}\| < \infty$, then under conditions (R4) and (V4), one may interchange the order of integration and differentiation in $\partial L_P/\partial \beta$ and $\partial L_P/\partial \theta$, on a neighborhood of $\xi_P$. It follows that besides the constraint in (3.5), the pair $(\xi_P, \lambda_P)$ satisfies

$$\int u(d) \mathbf{X}^T \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\beta) dP(\mathbf{s}) = \mathbf{0} \\ \text{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) + \frac{\lambda}{2} \int u(d) (\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\beta) dP(\mathbf{s}) = 0, \quad (7.1)$$

for $j = 1, \dots, l$, where $u(s) = \rho'(s)/s$ and $d = d(\mathbf{s}, \xi)$ is defined by (5.1), and where we abbreviate $\mathbf{V}(\theta)$ by $\mathbf{V}$. To solve $\lambda_P$ from the second set of equations, we multiply the $j$-th equation by $\theta_j$ and then sum over $j = 1, \dots, l$. This leads to

$$\text{tr} \left( \mathbf{V}^{-1} \sum_{j=1}^{l} \theta_j \frac{\partial \mathbf{V}}{\partial \theta_j} \right) + \frac{\lambda}{2} \int u(d) (\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}^{-1} \left( \sum_{j=1}^{l} \theta_j \frac{\partial \mathbf{V}}{\partial \theta_j} \right) \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\beta) dP(\mathbf{s}) = 0,$$

which is solved by

$$\lambda_P = \frac{-2\text{tr}\left(\mathbf{V}^{-1}\sum_{j=1}^{l}\theta_j(\partial\mathbf{V}/\partial\theta_j)\right)}{\int u(d)(\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}^{-1}\left(\sum_{j=1}^{l}\theta_j(\partial\mathbf{V}/\partial\theta_j)\right) \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\beta) dP(\mathbf{s})}.$$

When we insert this back into the second equation in (7.1), we find

$$\text{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) \int u(d) (\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}^{-1} \left( \sum_{t=1}^{l} \theta_t \frac{\partial \mathbf{V}}{\partial \theta_t} \right) \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\beta) dP(\mathbf{s}) \\ - \text{tr} \left( \mathbf{V}^{-1} \sum_{t=1}^{l} \theta_t \frac{\partial \mathbf{V}}{\partial \theta_t} \right) \int u(d) (\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\beta) dP(\mathbf{s}) = 0,$$

or briefly

$$\int u(d)(\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\beta) dP(\mathbf{s}) = 0, \quad j = 1, \dots, l, \quad (7.2)$$

where

$$\mathbf{H}_j = \text{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) \left( \sum_{t=1}^{l} \theta_t \frac{\partial \mathbf{V}}{\partial \theta_t} \right) - \text{tr} \left( \mathbf{V}^{-1} \sum_{t=1}^{l} \theta_t \frac{\partial \mathbf{V}}{\partial \theta_t} \right) \frac{\partial \mathbf{V}}{\partial \theta_j}. \quad (7.3)$$
Because $\sum_{j=1}^l \theta_j \mathbf{H}_j = \mathbf{0}$, the system of equations (7.2) is linearly dependent. Similar to [17] we subtract the S-constraint from each equation. For each $j = 1, \dots, l$, we subtract the term

$$ \mathrm{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) (\rho(d) - b_0) $$

from the left hand side of equation (7.2). We then find that any solution $\xi_P$ of (3.5) satisfies the following equation

$$ \int \Psi(\mathbf{s}, \boldsymbol{\xi}) \, dP(\mathbf{s}) = \mathbf{0}, \qquad (7.4) $$

where $\Psi = (\Psi_{\beta}, \Psi_{\theta})$, with $\Psi_{\theta} = (\Psi_{\theta,1}, \dots, \Psi_{\theta,l})$, where

$$ \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi}) = u(d) \mathbf{X}^T \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \qquad (7.5) $$

$$ \Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi}) = u(d)(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) - \mathrm{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) (\rho(d) - b_0), $$

for $j = 1, \dots, l$, where $\mathbf{H}_j$ and $d = d(\mathbf{s}, \boldsymbol{\xi})$ are defined in (7.3) and (5.1), respectively, and where we abbreviate $\mathbf{V}(\boldsymbol{\theta})$ by $\mathbf{V}$.

The regression score equation for $\Psi_{\beta}$ with the empirical measure $\mathbb{P}_n$ for $P$ in (7.4) coincides with the one for the regression S-estimator in the linear mixed effects model (2.4) considered in [5] (see their equation (10)). The empirical regression score equation also coincides with the one for the regression S-estimator in the multivariate regression model of Example 2 considered in [30] (see equation (2.2) in [29]). Similarly, the empirical score equation for $\Psi_{\beta}$ coincides with the one for the location S-estimator of Example 4 considered in [17].

For general covariance structures the empirical covariance score equation for $\Psi_{\theta}$ does not compare directly to existing equations in the literature. However, as we will see in the next subsection, similar comparisons are available for models with a linear covariance structure.

## 7.2 Linear covariance structures

In the previous section, we solved $\lambda$ from (7.1) and subtracted the S-constraint, leading to score equation (7.4) with $\Psi$ given in (7.5). The fact that this was done in a specific way has the following reason. In cases where $\mathbf{V}(\boldsymbol{\theta})$ is linear, say

$$ \mathbf{V}(\boldsymbol{\theta}) = \sum_{j=1}^{l} \theta_j \mathbf{L}_j, \qquad (7.6) $$

the function $\Psi_{\theta}$ simplifies a lot and can also be related to the covariance psi-function in [17]. Typical models of interest that have a covariance matrix of this type are the mixed linear effects model from Example 1 and the multivariate regression model from Example 2. But also the multivariate location-scale model from Example 4 and the time series model (2.8) from Example 3 have linear covariance structures.

When $\mathbf{V}$ is of the form (7.6), then $\partial \mathbf{V}/\partial \theta_j = \mathbf{L}_j$ and $\sum_{j=1}^l \theta_j (\partial \mathbf{V}/\partial \theta_j) = \mathbf{V}$. In this case, (7.3) simplifies to $\mathbf{H}_j = \mathrm{tr} (\mathbf{V}^{-1} \mathbf{L}_j) \mathbf{V} - k \mathbf{L}_j$, and $\Psi_{\theta,j}$ in (7.5) becomes

$$ \Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi}) = \mathrm{tr} (\mathbf{V}^{-1} \mathbf{L}_j) v(d) - ku(d)(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{L}_j \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}), $$

where $u(s)$ is defined in (R4) and

$$ v(s) = u(s)s^2 - \rho(s) + b_0. \qquad (7.7) $$

Using that $\mathrm{tr}(\mathbf{A}^T \mathbf{B}) = \mathrm{vec}(\mathbf{A})^T \mathrm{vec}(\mathbf{B})$, this can be written as

$$ \Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi}) = -\mathrm{vec}(ku(d)(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T - v(d)\mathbf{V})^T \mathrm{vec}(\mathbf{V}^{-1}\mathbf{L}_j\mathbf{V}^{-1}). $$
On the right hand side we recognize $ku(d)(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T - v(d)\mathbf{V}$, being the covariance psi-function that also appears in (2.8) in [17]. For our purposes we define

$$ \Psi_{\mathbf{V}}(\mathbf{s}, \boldsymbol{\xi}) = ku(d(\mathbf{s}, \boldsymbol{\xi}))(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T - v(d(\mathbf{s}, \boldsymbol{\xi}))\mathbf{V}. \quad (7.8) $$

The functions $\Psi_{\theta,j}$, for $j = 1, \dots, l$, can be combined in one expression for the vector valued function $\Psi_{\theta}$ as follows. First note that

$$ \mathrm{vec}(\mathbf{V}^{-1}\mathbf{L}_j\mathbf{V}^{-1}) = (\mathbf{V}^{-1} \otimes \mathbf{V}^{-1}) \mathrm{vec}(\mathbf{L}_j) $$

for $j = 1, \dots, l$. Define the $k^2 \times l$ matrix

$$ \mathbf{L} = \left[ \mathrm{vec}(\mathbf{L}_1) \quad \dots \quad \mathrm{vec}(\mathbf{L}_l) \right]. \quad (7.9) $$

Then, the column vector $\Psi_{\theta} = (\Psi_{\theta,1}, \dots, \Psi_{\theta,l})$ can be written as

$$ \Psi_{\theta}(\mathbf{s}, \boldsymbol{\xi}) = -\mathbf{L}^T (\mathbf{V}^{-1} \otimes \mathbf{V}^{-1}) \mathrm{vec}(\Psi_{\mathbf{V}}(\mathbf{s}, \boldsymbol{\xi})), $$

where $\Psi_{\mathbf{V}}$ is defined in (7.8) and $\mathbf{L}$ in (7.9). Note that the dependence on $\mathbf{s} = (\mathbf{y}, \mathbf{X})$ in $\Psi_{\theta}$ is only through the function $\Psi_{\mathbf{V}}$. We conclude that in the case of a linear covariance structure, any solution $\boldsymbol{\xi}_P$ of (3.5) satisfies (7.4), where $\Psi = (\Psi_{\beta}, \Psi_{\theta})$, with

$$ \begin{aligned} \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi}) &= u(d)\mathbf{X}^T\mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ \Psi_{\theta}(\mathbf{s}, \boldsymbol{\xi}) &= -\mathbf{L}^T(\mathbf{V}^{-1} \otimes \mathbf{V}^{-1})\mathrm{vec}(\Psi_{\mathbf{V}}(\mathbf{s}, \boldsymbol{\xi})) \end{aligned} \quad (7.10) $$

where $d = d(\mathbf{s}, \boldsymbol{\xi})$ is defined in (5.1), and where we abbreviate $\mathbf{V}(\boldsymbol{\theta})$ by $\mathbf{V}$.

For the multivariate regression model in Example 2, one has $\mathbf{V}(\boldsymbol{\theta}) = \mathbf{C}$, where $\boldsymbol{\theta} = \mathrm{vech}(\mathbf{C})$. The matrix $\mathbf{L} = \partial\mathrm{vec}(\mathbf{V})/\partial\boldsymbol{\theta}^T$ is then equal to the so-called duplication matrix $\mathcal{D}_k$, which is the unique $k^2 \times k(k+1)/2$ matrix, with the properties $\mathcal{D}_k\mathrm{vech}(\mathbf{C}) = \mathrm{vec}(\mathbf{C})$ and $(\mathcal{D}_k\mathcal{D}_k)^{-1}\mathcal{D}_k^T\mathrm{vec}(\mathbf{C}) = \mathrm{vech}(\mathbf{C})$ (e.g., see [20, Ch. 3, Sec. 8]). Because $\mathbf{V}$ has full rank, it follows that equation (7.4) holds for $\Psi = (\Psi_{\beta}, \Psi_{\mathbf{V}})$. The resulting score equations for the empirical measure $\mathbb{P}_n$ corresponding to observations $(\mathbf{y}_i, \mathbf{X}_i)$, for $i = 1, \dots, n$, are then equivalent with the ones found in [30].

For the linear mixed effects model (2.4), the covariance matrix $\mathbf{V}(\boldsymbol{\theta})$ has a linear structure with the vector $\boldsymbol{\theta} = (\sigma_0^2, \dots, \sigma_r^2)$ of unknown covariance parameters. The matrix $\mathbf{L}$ is then a $k^2 \times (r+1)$ matrix and will typically be of rank $r+1 < k^2$. As a consequence, in this case one cannot further simplify equation (7.4), by removing the factor $\mathbf{L}^T (\mathbf{V}^{-1} \otimes \mathbf{V}^{-1})$ from the function $\Psi_{\theta}$. The score equation for $\Psi_{\beta}$ resulting from the empirical measure $\mathbb{P}_n$ corresponding to observations $(\mathbf{y}_i, \mathbf{X}_i)$, for $i = 1, \dots, n$, is the same as the one obtained in [5]. The corresponding score equation for $\Psi_{\theta}$ differs slightly from the one in [5], because the authors do not subtract a term with $\rho(d) - b_0$ to remove the linear dependency of the equations (7.2).

# 8 Local robustness: the influence function

For $0 < h < 1$ and $\mathbf{s} = (\mathbf{y}, \mathbf{X}) \in \mathbb{R}^k \times \mathbb{R}^{kq}$ fixed, define the perturbed probability measure

$$ P_{h,s} = (1-h)P + h\delta_s, $$

where $\delta_s$ denotes the Dirac measure at $\mathbf{s} \in \mathbb{R}^k \times \mathbb{R}^{kq}$. The influence function of the functional $\boldsymbol{\xi}(\cdot)$ at probability measure $P$, is defined as

$$ \mathrm{IF}(\mathbf{s}; \boldsymbol{\xi}, P) = \lim_{h \downarrow 0} \frac{\boldsymbol{\xi}((1-h)P + h\delta_s) - \boldsymbol{\xi}(P)}{h}, \quad (8.1) $$

if this limit exists. In contrast to the global robustness measured by the breakdown point, the influence function measures the local robustness. It describes the effect of an infinitesimal contamination at a single point $\mathbf{s}$ on the functional (see Hampel [11]). Good local robustness is therefore illustrated by a bounded influence function.
## 8.1 The general case

The theorem below gives the influence function for the S-functional $\xi$. It extends the result for S-functionals of multivariate location and scale [17]. Under the assumption that the limit in (8.1) exists and $P$ has an elliptical contoured density (3.2), Van Aelst and Willems [30] relate the influence function for S-functionals of multivariate regression to that of S-functionals of multivariate location and scale. For the linear mixed effects model considered in [5], the influence function has not been established. The influence function for these functionals now follows as a special case from the theorem below.

We will show that the limit in (8.1) exists and derive its expression at general $P$. Since the value of $\theta$ determines the covariance matrix $V(\theta)$, we also include the influence function of the covariance functional. Consider the S-functional at $P_{h,s_0}$. From the Portmanteau theorem [2, Theorem 2.1] it can easily be seen that $P_{h,s_0} \to P$, weakly, as $h \downarrow 0$. Therefore, under the conditions of Corollary 2 and Theorem 2, it follows that there exist solutions $\xi(P_{h,s_0})$ and $\xi(P)$ to minimization problems (3.5) at $P_{h,s_0}$ and $P$, respectively, and that $\xi(P_{h,s_0}) \to \xi(P)$, as $h \downarrow 0$.

**Theorem 5.** Let $\xi(P_{h,s_0})$ and $\xi(P)$ be solutions to minimization problems (3.5) at $P_{h,s_0}$ and $P$, respectively, and suppose that $\xi(P_{h,s_0}) \to \xi(P)$, as $h \downarrow 0$. Suppose that $\rho$ satisfies (R4) and $V$ satisfies (V4). Let $\Psi$ be defined in (7.5) and suppose that

$$ \Lambda(\xi) = \int \Psi(s, \xi) dP(s), \qquad (8.2) $$

is continuously differentiable with a non-singular derivative $D(P)$ at $\xi(P)$. Then for $s_0 \in \mathbb{R}^k \times \mathbb{R}^{kq}$,

$$ \mathrm{IF}(s_0; \xi, P) = -D(P)^{-1} \Psi(s_0, \xi(P)). $$

For the covariance functional $C(P) = V(\theta(P))$, it holds that

$$ \mathrm{IF}(s_0; \mathrm{vec}(C), P) = \left( \frac{\partial \mathrm{vec}(V(\theta(P)))}{\partial \theta^T} \right) \mathrm{IF}(s_0; \theta, P). $$

To investigate the local robustness of S-estimators, we derive the following bound on the influence function for $\xi(P)$.

**Corollary 4.** Suppose that $\rho$ satisfies (R2) and (R4), and $V$ satisfies (V4). Then there exist $0 < C_1 < \infty$ and $0 < C_2 < \infty$, only depending on $P$, such that for $s = (y, X)$ it holds that $\|\mathrm{IF}(s, \xi(P))\| \le C_1 + C_2\|X\|$.

Its proof can be found in [16].

## 8.2 Elliptically contoured densities

When $P$ is such that $y | X$ has an elliptically contoured density (3.2) and $V(\theta)$ is linear, we can obtain a more detailed expression for the influence function, This requires the following condition on the function $\rho$,

(R5) $\rho$ is twice continuously differentiable,

and the following condition on the mapping $\theta \mapsto V(\theta)$,

(V5) $V(\theta)$ is twice continuously differentiable.

Conditions (R5) and (V5) are needed to establish that $\Lambda$, as defined in (8.2), is continuously differentiable. Clearly, condition (V5) implies former conditions (V4) and (V1).

Suppose that $P$ is such that $y | X$ has an elliptically contoured density $f_{\mu, \Sigma}$ from (3.2), with $\mu \in \mathbb{R}^k$ and $\Sigma \in \mathrm{PDS}(k)$. When the S-functional is affine equivariant, it suffices to determine the influence function for the case $(\mu, \Sigma) = (0, I_k)$. However, this does not hold in general for the
S-functionals in our setting. The reason is that, for a $k \times k$ non-singular matrix $\mathbf{A}$ and $\boldsymbol{\theta} \in \mathbb{R}^l$, the matrix $\mathbf{AV}(\boldsymbol{\theta})\mathbf{A}^T$ may not be of the form $\mathbf{V}(\boldsymbol{\theta}')$, for some $\boldsymbol{\theta}' \in \mathbb{R}^l$. Examples are the (linear) covariance structure that corresponds to the linear mixed effects model (2.4) considered in [5] or the models discussed in Example 3.

Nevertheless, note that for the general case with $\boldsymbol{\mu} \in \mathbb{R}^k$ and $\boldsymbol{\Sigma} \in \text{PDS}(k)$, we can still use the fact that, conditionally on $\mathbf{X}$, the distribution of $\mathbf{y}$ is the same as that of $\boldsymbol{\Sigma}^{1/2}\mathbf{z} + \boldsymbol{\mu}$, where $\mathbf{z}$ has a spherical density $f_{0,\mathbf{I}_k}$. As a consequence, we can still obtain the following result, which enables one to determine the influence functions of the functionals $\boldsymbol{\beta}(P)$ and $\boldsymbol{\theta}(P)$ separately.

If $P$ itself is also absolutely continuous, then it satisfies (C3), as well as (C1$_{\epsilon}$') and (C2$_{\epsilon}$), for any $0 < \epsilon' < \epsilon \le 1-r$. When $\rho$ and $\mathbf{V}$ satisfy (R1)-(R3) and (V1)-(V3), it follows from Theorem 1 and Corollary 2 that $\boldsymbol{\xi}(P)$ and $\boldsymbol{\xi}(P_{h,s})$ exist, for $h$ sufficiently small. If $h$ in (3.2) is non-increasing and not constant on $[0, c_0^2]$, then $\boldsymbol{\xi}(P)$ is unique, according to Theorem 3, so that $\boldsymbol{\xi}(P_{h,s}) \to \boldsymbol{\xi}(P)$, as $h \downarrow 0$. Hence, in order to apply Theorem 5, it remains to show that $\Lambda$ in (8.2) is continuously differentiable with a non-singular derivative at $\boldsymbol{\xi}(P)$. As a first step we obtain that the derivative of $\Lambda$ is a block matrix.

**Lemma 2.** Suppose that $P$ is such that $\mathbf{y} | \mathbf{X}$ has an elliptically contoured density $f_{\boldsymbol{\mu},\boldsymbol{\Sigma}}$ from (3.2) and $\mathbb{E}\|\mathbf{X}\|^2 < \infty$. Suppose that $\boldsymbol{\xi}(P)$ is a solution to the corresponding minimization problem (3.5), such that $(\mathbf{X}\boldsymbol{\beta}(P), \mathbf{V}(\boldsymbol{\theta}(P))) = (\boldsymbol{\mu}, \boldsymbol{\Sigma})$. Suppose that $\rho$ satisfies (R2), (R4)-(R5) and that $\mathbf{V}$ satisfies (V5) and has a linear structure (7.6). Let $\Lambda$ be defined in (8.2) with $\Psi$ defined in (7.10). Then

$$ \frac{\partial \Lambda(\boldsymbol{\xi}(P))}{\partial \boldsymbol{\xi}} = \begin{pmatrix} \frac{\partial \Lambda_{\boldsymbol{\beta}}(\boldsymbol{\xi}(P))}{\partial \boldsymbol{\beta}} & \mathbf{0} \\ \mathbf{0} & \frac{\partial \Lambda_{\boldsymbol{\theta}}(\boldsymbol{\xi}(P))}{\partial \boldsymbol{\theta}} \end{pmatrix}, $$

where

$$ \frac{\partial \Lambda_{\boldsymbol{\beta}}(\boldsymbol{\xi}(P))}{\partial \boldsymbol{\beta}} = -\alpha \mathbb{E} \left[ \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right], \quad (8.3) $$

with

$$ \alpha = \mathbb{E}_{0, \mathbf{I}_k} \left[ \left( 1 - \frac{1}{k} \right) \frac{\rho'(\|\mathbf{z}\|)}{\|\mathbf{z}\|} + \frac{1}{k} \rho''(\|\mathbf{z}\|) \right], \quad (8.4) $$

and

$$ \frac{\partial \Lambda_{\boldsymbol{\theta}}(\boldsymbol{\xi}(P))}{\partial \boldsymbol{\theta}} = \gamma_1 \mathbf{L}^T (\boldsymbol{\Sigma}^{-1} \otimes \boldsymbol{\Sigma}^{-1}) \mathbf{L} - \gamma_2 \mathbf{L}^T \text{vec}(\boldsymbol{\Sigma}^{-1}) \text{vec}(\boldsymbol{\Sigma}^{-1})^T \mathbf{L}. $$

where $\mathbf{L} = \partial\text{vec}(\mathbf{V}(\boldsymbol{\theta}(P)))/\partial\boldsymbol{\theta}^T$ is the $k^2 \times l$ matrix given in (7.9) and

$$ \begin{align} \gamma_1 &= \frac{\mathbb{E}_{0, \mathbf{I}_k} [\rho''(\|\mathbf{z}\|) \|\mathbf{z}\|^2 + (k+1)\rho'(\|\mathbf{z}\|) \|\mathbf{z}\|]}{k+2} \tag{8.5} \\ \gamma_2 &= \frac{\mathbb{E}_{0, \mathbf{I}_k} [2\rho''(\|\mathbf{z}\|) \|\mathbf{z}\|^2 + k\rho'(\|\mathbf{z}\|) \|\mathbf{z}\|]}{2k(k+2)}, \notag \end{align} $$

The proof is tedious, but straightforward, and can be found in [16].

**Remark 8.1.** The proof of Lemma 2 uses the fact that

$$ \frac{\partial \Lambda(\boldsymbol{\xi})}{\partial \boldsymbol{\xi}} = \int \frac{\partial \Psi(\mathbf{s}, \boldsymbol{\xi})}{\partial \boldsymbol{\xi}} dP(\mathbf{s}). $$

for all $\boldsymbol{\xi}$ in a neighborhood of $\boldsymbol{\xi}(P)$. This holds for general $P$ and any covariance structure $\mathbf{V}(\boldsymbol{\theta})$ that satisfies (V2)-(V3) and (V5), see Lemma B.3 in [16]. Furthermore, Lemma 2 is obtained for a linear covariance structure. However, with some additional technicalities, this result can also be
shown to hold for $\Psi$ defined in (7.5) corresponding to general covariance structures. For general covariance structures one still obtains (8.3), and that

$$ \frac{\partial \Lambda_{\theta,j}(\xi(P))}{\partial \theta_s} = -\alpha_1 \mathrm{tr} \left( \Sigma^{-1} \frac{\partial \mathbf{V}(\theta(P))}{\partial \theta_s} \Sigma^{-1} \mathbf{H}_j \right) \\ + \alpha_2 \mathrm{tr} \left( \Sigma^{-1} \frac{\partial \mathbf{V}(\theta(P))}{\partial \theta_s} \right) \mathrm{tr} \left( \Sigma^{-1} \frac{\partial \mathbf{V}(\theta(P))}{\partial \theta_j} \right), $$

for $j, s = 1, \dots, l$, and where $\alpha_1 = \gamma_1/k$ and $\alpha_2 = \gamma_1/k - \gamma_2$, with $\gamma_1, \gamma_2$ from (8.5), and where $\mathbf{H}_j$ is defined in (7.3).

The next corollary gives expressions for the influence functions of the functionals $\beta(P)$ and $\theta(P)$ separately, at a distribution $P$ that is such that $\mathbf{y} | \mathbf{X}$ has an elliptically contoured density. The proof is tedious, but straightforward, and can be found in [16].

**Corollary 5.** Suppose that $P$ is such that $\mathbf{y} | \mathbf{X}$ has an elliptically contoured density $f_{\mu, \Sigma}$ from (3.2), such that $(\mathbf{X}\beta(P), \mathbf{V}(\theta(P))) = (\mu, \Sigma)$. Let $\xi(P_{h,s_0})$ and $\xi(P)$ be a solution to minimization problem (3.5) at $P_{h,s_0}$ and $P$, respectively, and suppose that $\xi(P_{h,s_0}) \to \xi(P)$, as $h \downarrow 0$. Suppose that $\mathbb{E}\|\mathbf{X}\|^2 < \infty$ and suppose that $\rho$ satisfies (R2)-(R5) and that $\mathbf{V}$ satisfies (V5), and has a linear structure (7.6). Let $\alpha, \gamma_1$, and $\gamma_2$ be defined in (8.4) and (8.5), and suppose that $\mathbb{E}_{0, \mathbf{I}_k}[\rho''(\|\mathbf{z}\|)] > 0$. If $\mathbf{X}$ has full rank with probability one, then

$$ \mathrm{IF}(\mathbf{s}_0, \beta, P) = \frac{u(d_0)}{\alpha} \left( \mathbb{E} \left[ \mathbf{X}^T \Sigma^{-1} \mathbf{X} \right] \right)^{-1} \mathbf{X}_0^T \Sigma^{-1} (\mathbf{y}_0 - \mathbf{X}_0 \beta(P)) $$

where $d_0^2 = (\mathbf{y}_0 - \mathbf{X}_0\beta(P))^T \Sigma^{-1}(\mathbf{y}_0 - \mathbf{X}_0\beta(P))$ and $u(s) = \rho'(s)/s$. If $\gamma_1 > 0$ and the $k^2 \times l$ matrix $\mathbf{L}$, as defined in (7.9), has full rank, then $\mathrm{IF}(\mathbf{s}_0, \theta, P)$ is given by

$$ \frac{ku(d_0)}{\gamma_1} \left( \mathbf{L}^T (\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1}) \mathbf{L} \right)^{-1} \mathbf{L}^T \mathrm{vec} \left( \mathbf{\Sigma}^{-1} (\mathbf{y}_0 - \mathbf{X}_0 \beta(P)) (\mathbf{y}_0 - \mathbf{X}_0 \beta(P))^T \mathbf{\Sigma}^{-1} \right) \\ + \left( -\frac{u(d_0)d_0^2}{\gamma_1} + \frac{\rho(d_0) - b_0}{\gamma_1 - k\gamma_2} \right) \theta(P). $$

Note that since $\mathbf{L}\theta(P) = \mathrm{vec}(\mathbf{V}(\theta(P))) = \mathrm{vec}(\mathbf{\Sigma})$, we can immediately obtain the influence function for the covariance functional $\mathbf{C}(P) = \mathbf{V}(\theta(P))$. From Theorem 5 it immediately follows that $\mathrm{IF}(\mathbf{s}_0, \mathrm{vec}(\mathbf{C}), P)$ is given by

$$ \frac{ku(d_0)}{\gamma_1} \mathbf{L} \left( \mathbf{L}^T (\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1}) \mathbf{L} \right)^{-1} \mathbf{L}^T \mathrm{vec} \left( \mathbf{\Sigma}^{-1} (\mathbf{y}_0 - \mathbf{X}_0 \beta)(\mathbf{y}_0 - \mathbf{X}_0 \beta)^T \mathbf{\Sigma}^{-1} \right) \\ + \left( -\frac{u(d_0)d_0^2}{\gamma_1} + \frac{\rho(d_0) - b_0}{\gamma_1 - k\gamma_2} \right) \mathrm{vec}(\mathbf{\Sigma}). $$

Since the functions $u(s)s = \rho'(s)$, $u(s)s^2 = \rho'(s)s$, and $\rho(s)$ are bounded, it follows that $\mathrm{IF}(\mathbf{s}, \theta, P)$ and $\mathrm{IF}(\mathbf{s}, \mathrm{vec}(\mathbf{C}), P)$ are bounded uniformly in both $\mathbf{y}$ and $\mathbf{X}$, whereas $\mathrm{IF}(\mathbf{s}, \beta, P)$ is bounded uniformly in $\mathbf{y}$, but not in $\mathbf{X}$. This illustrates the phenomenon in linear regression that leverage points can have a high effect on the regression S-estimator.

For the S-estimators in the linear mixed effects model (2.4) with normal errors considered in [5], the influence function is not available. The expression can now be obtained from Corollary 5. The expression for $\mathrm{IF}(\mathbf{s}, \beta, P)$ in Corollary 5 coincides with the one found for the multivariate regression S-functional in [30], where $\alpha > 0$ is the same constant as the one in the expression of the influence function for the location S-functional in [17]. Furthermore, for the multivariate regression model, one has $\theta = \mathrm{vech}(\mathbf{C})$ and the matrix $\mathbf{L}$ is equal to the duplication matrix $\mathcal{D}_k$. From the properties of $\mathcal{D}_k$, the expressions for the influence functions simplify. One finds in this case that

$$ \mathrm{IF}(\mathbf{s}, \theta, P) = \frac{ku(d)}{\gamma_1} \mathrm{vech} \left( (\mathbf{y} - \mathbf{X}\beta(P))(\mathbf{y} - \mathbf{X}\beta(P))^T \right) + \left( -\frac{u(d)d^2}{\gamma_1} + \frac{\rho(d) - b_0}{\gamma_1 - k\gamma_2} \right) \theta(P) $$
and the influence function of the covariance functional $C(P) = V(\theta(P))$ itself is given by

$$ \mathrm{IF}(\mathbf{s}, \mathbf{C}, P) = \frac{ku(d)}{\gamma_1}(\mathbf{y} - \mathbf{X}\beta(P))(\mathbf{y} - \mathbf{X}\beta(P))^T + \left(-\frac{u(d)d^2}{\gamma_1} + \frac{\rho(d) - b_0}{\gamma_1 - k\gamma_2}\right) \mathbf{\Sigma}. $$

This coincides with the expressions found for the covariance S-functionals in [30] and in [17].

# 9 Asymptotic normality

To establish asymptotic normality of the S-estimators, we use the score equations obtained from differentiation of the Lagrangian corresponding to the minimization problem (3.1). In the same way as before, we obtain score equation (7.4), with $P$ equal to the empirical measure $\mathbb{P}_n$ corresponding to observations $\mathbf{s}_1, \dots, \mathbf{s}_n$, with $\mathbf{s}_i = (\mathbf{y}_i, \mathbf{X}_i) \in \mathbb{R}^k \times \mathbb{R}^{kq}$. From (7.4), we see that any solution $\xi_n = \xi(\mathbb{P}_n)$ to the S-minimization problem (3.1) must satisfy

$$ \int \Psi(\mathbf{s}, \xi_n) d\mathbb{P}_n(\mathbf{s}) = 0, \qquad (9.1) $$

where $\Psi = (\Psi_\beta, \Psi_\theta)$ is defined in (7.5).

## 9.1 General case

Writing $\xi_P = \xi(P)$, we decompose (9.1) as follows

$$ \begin{aligned} 0 = \int \Psi(\mathbf{s}, \xi_n) dP(\mathbf{s}) &+ \int \Psi(\mathbf{s}, \xi_P) d(\mathbb{P}_n - P)(\mathbf{s}) \\ &+ \int (\Psi(\mathbf{s}, \xi_n) - \Psi(\mathbf{s}, \xi_P)) d(\mathbb{P}_n - P)(\mathbf{s}). \end{aligned} \qquad (9.2) $$

The essential step in establishing asymptotic normality of $\xi_n$ is to show that the third term on the right hand side of (9.2) is of the order $o_P(n^{-1/2})$. To this end we will apply results from empirical process theory as developed in Pollard [23]. This leads to the following theorem.

**Theorem 6.** Suppose that $\rho$ satisfies (R1)-(R2) and (R4), such that $u(s)$ is of bounded variation, and suppose that $V$ satisfies (V4). Let $\xi_n$ and $\xi(P)$ be solutions to minimization problems (3.1) and (3.5), and suppose that $\xi_n \to \xi(P)$ in probability. Suppose that $\Lambda$, as defined in (8.2) with $\Psi$ defined in (7.5), is continuously differentiable with a non-singular derivative $D(P)$ at $\xi(P)$ and suppose that $E\|\mathbf{X}\|^2 < \infty$. Then $\sqrt{n}(\xi_n - \xi(P))$ is asymptotically normal with mean zero and covariance matrix

$$ D(P)^{-1}E[\Psi(\mathbf{s}, \xi(P))\Psi(\mathbf{s}, \xi(P))^T] D(P)^{-1}. $$

Theorem 6 is similar to Theorem 4.1 in [17]. Note that Theorem 6 confirms the well known heuristic that relates the limiting covariance of $\sqrt{n}(\xi_n - \xi(P))$ to the influence function of the functional $\xi(\cdot)$ given in Theorem 5,

$$ D(P)^{-1}E[\Psi(\mathbf{s}, \xi(P))\Psi(\mathbf{s}, \xi(P))^T] D(P)^{-1} = E[\mathrm{IF}(\mathbf{s}, \xi, P)\mathrm{IF}(\mathbf{s}, \xi, P)^T]. \qquad (9.3) $$

Van Aelst and Willems [30] consider the limiting behavior of S-estimators in the multivariate regression model of Example 2, but only under $P$ for which $\mathbf{y} | \mathbf{X}$ has an elliptical contoured density. Copt and Victoria-Feser [5] consider asymptotic normality for S-estimators in the linear mixed effects model (2.4) with a constant design matrix $\mathbf{X}_i = \mathbf{X}$ and only consider $P$ for which $\mathbf{y} | \mathbf{X}$ has an multivariate normal distribution.
## 9.2 Elliptically contoured densities

Consider the special case that $P$ is such that $y \mid X$ has an elliptically contoured density $f_{\mu,\Sigma}$ from (3.2), with $\mu \in \mathbb{R}^k$ and $\Sigma \in \text{PDS}(k)$. As before, in determining the limiting normal distribution of the individual S-estimators, we cannot use affine equivariance and restrict ourselves to the case $(0, I_k)$. Instead, we use some of the results obtained in Section 8.2 to establish the limiting normal distributions of the S-estimators $\beta_n = \beta(\mathbb{P}_n)$, $\theta_n = \theta(\mathbb{P}_n)$, and $C_n = V(\theta(\mathbb{P}_n))$.

**Corollary 6.** Suppose that $P$ is such that $y \mid X$ has an elliptically contoured density $f_{\mu,\Sigma}$ from (3.2), such that $(X\beta(P), V(\theta(P))) = (\mu, \Sigma)$. Let $\xi_n$ and $\xi(P)$ be solutions to minimization problems (3.1) and (3.5), and suppose that $\xi_n \to \xi(P)$ in probability. Suppose that $E\|X\|^2 < \infty$ and suppose that $\rho$ satisfies (R2)-(R5), such that $u(s)$ is of bounded variation. Suppose that $V$ satisfies (V5), and has a linear structure (7.6). Let $\alpha, \gamma_1$, and $\gamma_2$ be defined in (8.4) and (8.5), and suppose that $E_{0, I_k}[\rho''(\|z\|)] > 0$. If $X$ has full rank with probability one, then $\sqrt{n}(\beta_n - \beta(P))$ is asymptotically normal with mean zero and covariance matrix

$$ \frac{\mathbb{E}_{0, I_k}[\rho'(\|z\|)^2]}{k\alpha^2} (\mathbb{E}[X^T\Sigma^{-1}X])^{-1}. $$

If $\gamma_1 > 0$ and the $k^2 \times l$ matrix $L$, as defined in (7.9), has full rank, then $\sqrt{n}(\theta_n - \theta(P))$ is asymptotically normal with mean zero and covariance matrix

$$ 2\sigma_1 (L^T (\Sigma^{-1} \otimes \Sigma^{-1}) L)^{-1} + \sigma_2 \theta(P) \theta(P)^T, $$

where

$$ \begin{aligned} \sigma_1 &= \frac{k(k+2)\mathbb{E}_{0, I_k}[u(\|z\|)^2 \|z\|^4]}{(\mathbb{E}_{0, I_k}[\rho''(\|z\|) \|z\|^2] + (k+1)\rho'(\|z\|) \|z\|)^2} \\ \sigma_2 &= -\frac{2}{k}\sigma_1 + \frac{4\mathbb{E}_{0, I_k}[(\rho(\|z\|) - b_0)^2]}{(\mathbb{E}_{0, I_k}[\rho'(\|z\|)^2])^2} \end{aligned} $$

Due to the linearity of $V$, we can immediately establish asymptotic normality of the covariance estimator $C_n = V(\theta_n)$. From Corollary 6 it follows that

$$ \sqrt{n}(\text{vec}(C_n) - \text{vec}(\Sigma)) = \sqrt{n}(\mathbf{L}\theta_n - \mathbf{L}\theta(P)) = \mathbf{L}\sqrt{n}(\theta_n - \theta(P)). $$

It follows that the limiting covariance of $\sqrt{n}(\text{vec}(V(\theta_n)) - \text{vec}(\Sigma))$ is given by

$$ 2\sigma_1 \mathbf{L} (\mathbf{L}^T (\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1}) \mathbf{L})^{-1} \mathbf{L}^T + \sigma_2 \text{vec}(\mathbf{\Sigma}) \text{vec}(\mathbf{\Sigma})^T. $$

**Corollary 6** is a direct consequence of Theorem 6. Its proof, in particular the derivations of the expressions for the limiting covariances, can be found in [16]. Note that the constants $\mathbb{E}_{0, I_k}[\rho'(\|z\|)^2] / (k\alpha^2)$, $\sigma_1$ and $\sigma_2$, are the same as the ones found in [17] for the location and covariance S-estimators, respectively. In fact, Corollary 6 is an extension of Corollary 5.1 in [17] for S-estimators in the multivariate location-scale model of Example 4.

Asymptotic normality of S-estimators in the multivariate regression model of Example 2 follows from Corollary 6. These estimators have been considered in [30], but asymptotic normality has not been established. Under the assumption that the heuristic (9.3) holds, asymptotic relative efficiencies are computed on the basis of this heuristic. Indeed, now that Corollary 6 has been established, one may check that (9.3) holds.

Finally, note that the limiting covariances of $\sqrt{n}(\beta_n - \beta(P))$ and $\sqrt{n}(\theta_n - \theta(P))$ in Corollary 6 differ from the ones found in [5] for the linear mixed effects model (2.4) with $X_i = X$, for $i = 1, \dots, n$. The results in [5] are obtained by re-parameterizing $X\beta = \mu$ and interpreting the model as a multivariate location-scale model. Then building on the results in [17] for S-estimators of multivariate location-scale, the limiting covariances in [5] are found by application of the delta method. However, in view of Remark 3.1 this does not seem to be a correct approach.
**Remark 9.1.** Although our expressions for the limiting covariances in Corollary 6 differ from the ones found in Proposition 1 in [5], somewhat surprisingly, they yield the same matrices for the example discussed in Section 5.1 in [5]. However, this is a consequence of the specific structure of the design matrices X and Z in this example. One can easily find other design matrices for which the limiting covariances in Corollary 6 yield different matrices as the ones found in [5]. Moreover, the corresponding confidence regions based on the expressions in Corollary 6 can be substantially smaller than the ones based on the expressions found in [5]. See the simulation in Section 10.

# 10 Simulation and data example

We compare the asymptotic results of the S-estimators with their finite sample behavior by means of a simulation. Moreover we investigate the differences between the expressions found in Corollary 6 and the ones in Copt and Victoria-Feser [5]. To this end we will study the behavior of the estimators for samples generated from a model that is close to the one in [5]:

$$y_i = X\beta + \gamma_i Z + \epsilon_i, \quad i = 1, \dots, n, \qquad (10.1)$$

a linear mixed effects model with $y_i$ in dimension $k = 4$ and all subjects with the same design matrix X for the fixed effects $\beta = (\beta_1, \beta_2)^T$. Following the setup in [5], the matrix X is built as follows. The first column of X is taken to be a vector 1 consisting of ones of length four. The four x-values in the second column are generated from a standard normal, and then X is rescaled to a new matrix $X = [1 \ x]$, such that $X^T X = 4I_2$. For our simulation we used

$$X = \begin{pmatrix} 1 & -0.9504967 & & \\ 1 & -0.5428346 & & \\ 1 & 1.6650521 & & \\ 1 & -0.1717207 & & \end{pmatrix}.$$

The random effects $\gamma_i$ are independent $N(0, \sigma_\gamma^2)$ distributed random variables, which are independent from the measurement error $\epsilon_i \sim N(0, \sigma_\epsilon^2 R)$. This leads to a structured covariance $\Sigma = \sigma_\gamma^2 ZZ^T + \sigma_\epsilon^2 R$, with covariance parameter vector $\theta = (\theta_1, \theta_2)^T$, where $\theta_1 = \sigma_\gamma^2$ and $\theta_2 = \sigma_\epsilon^2$. Following the setup in [5], we set $\beta_1 = \beta_2 = 1$ and $\theta_1 = \theta_2 = 1$.

In [5], the authors took $Z = (1, 1, 1, 1)^T$ and $R = I_4$. With these choices the expression

$$\mathrm{Var}_{\mathrm{CVF}}(\beta_n) = \frac{\mathbb{E}_{0, I_k} [\rho' (\|z\|^2)]}{k\alpha^2} (X^T X)^{-1} (X^T \Sigma X (X^T X)^{-1}} \quad (10.2)$$

found in [5] for the limiting covariance matrix of $\sqrt{n}(\beta_n - \beta)$ (see (14) in [5]), is equal to our expression

$$\mathrm{Var}_{\mathrm{LGRG}}(\beta_n) = \frac{\mathbb{E}_{0, I_k} [\rho' (\|z\|^2)]}{k\alpha^2} (X^T \Sigma^{-1} X)^{-1}, \qquad (10.3)$$

found in Corollary 6, and similarly for the limiting covariance matrix of $\sqrt{n}(\theta_n - \theta)$. However, this is just the consequence of the extreme simple choices for X, Z and R. Already, if we keep X as it is, and only take a slight variation of either Z or R, one finds severe differences between (10.2) and (10.3), and similarly for the expression of the limiting covariance matrix of $\sqrt{n}(\theta_n - \theta)$.

We considered the following two alternatives

1. take $Z = (1, 2, 3, 4)^T$ and leave X and R = I_4 as they are;

2. take R = (1, 4, 9, 16)^T and leave X and Z = (1, 1, 1, 1)^T as they are.

We generated 10 000 samples of size $n = 100$ according to model (10.1) and computed the value of S-estimators $\beta_n$ and $\theta_n$ by means of Tukey's bi-weight

$$\rho_B(s; c) = \begin{cases} s^2/2 - s^4/(2c^2) + s^6/(6c^4), & |s| \le c \\ c^2/6 & |s| > c. \end{cases} \qquad (10.4)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Empirical marginal and joined distributions together with limiting marginal and joined distributions of $\sqrt{n}(\beta_n - \beta)$ (first row) and $\sqrt{n}(\theta_n - \theta)$ (second row).</figcaption></figure>

and $b_0 = E_{0, \mathbf{I}_k}[\rho_B(\|z\|; c_0)]$ with the cut-off value $c_0$ chosen such that $b_0/a_0 = 0.5$. According to Theorem 4, this corresponds to (asymptotic) breakdown point 50%.

Figure 1 displays the limiting marginal and joined distributions of $\sqrt{n}(\beta_n - \beta)$ in the first row, where we generated the samples with alternative 1. The histograms and scatterplot correspond to the 10 000 different values of $\sqrt{n}(\beta_n - \beta)$. The dashed curves correspond to the densities and 95% contourlines of the theoretical limiting marginal and joined normal distributions using the covariance matrix in (10.2). The solid curves correspond to the marginal and joined normal distributions using the covariance matrix in (10.3). The empirical contourlines based on the sample mean and sample covariance of the 10 000 estimates are plotted in dotted lines, but they almost indistinguishable from the solid contourlines. We find

$$ \mathrm{Var}_{\mathrm{CVF}}(\beta_n) = \begin{pmatrix} 8.13 & 1.78 \\ 1.78 & 0.72 \end{pmatrix} \quad \text{and} \quad \mathrm{Var}_{\mathrm{LGRG}}(\beta_n) = \begin{pmatrix} 1.97 & 0.38 \\ 0.38 & 0.40 \end{pmatrix}. $$

Clearly, the histograms of the repeated estimates for $\beta_1$ and $\beta_2$ match the graphs of the (marginal) normal densities with the variances given by $\mathrm{Var}_{\mathrm{LGRG}}(\beta_n)$, and the scatterplot matches with the 95% contourline corresponding to $\mathrm{Var}_{\mathrm{LGRG}}(\beta_n)$. Note that the differences with $\mathrm{Var}_{\mathrm{CVF}}(\beta_n)$ are quite severe. For example, this yields that the length of the confidence interval for $\beta_1$ based on $\mathrm{Var}_{\mathrm{CVF}}(\beta_n)$ will be two times larger than the one based on $\mathrm{Var}_{\mathrm{LGRG}}(\beta_n)$.

The second row in Figure 1 displays the limiting distributions of $\sqrt{n}(\theta_n - \theta)$, where we generated the samples with alternative 2. In [5], the limiting covariance matrix was given by (see (15) in [5]) $(\mathbf{L}^T\mathbf{L})^{-1}\mathbf{L}^T\mathbf{V}_{\Sigma}\mathbf{L}(\mathbf{L}^T\mathbf{L})^{-1}$, where $\mathbf{V}_{\Sigma} = \sigma_1(\mathbf{I}_{k^2} + \mathbf{K}_{k,k})(\mathbf{\Sigma} \otimes \mathbf{\Sigma}) + \sigma_2\mathrm{vec}(\mathbf{\Sigma})\mathrm{vec}(\mathbf{\Sigma})^T$, see Corollary 5.1 in [17]. Because

$$ (\mathbf{L}^T\mathbf{L})^{-1}\mathbf{L}^T(\mathbf{I}_{k^2} + \mathbf{K}_{k,k}) = 2(\mathbf{L}^T\mathbf{L})^{-1}\mathbf{L}^T, \\ (\mathbf{L}^T\mathbf{L})^{-1}\mathbf{L}^T\mathrm{vec}(\mathbf{\Sigma}) = \mathbf{\theta}(P), $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Left picture: standardized residuals for the S-estimates (horizontal axis) and the ML estimates (vertical axis). Right picture: observations for the subjects in the treatment group.</figcaption></figure>

the expression given in [5] becomes

$$ \mathrm{Var}_{\mathrm{CVF}}(\theta_n) = 2\sigma_1(\mathbf{L}^T\mathbf{L})^{-1}\mathbf{L}^T(\mathbf{\Sigma} \otimes \mathbf{\Sigma})\mathbf{L}(\mathbf{L}^T\mathbf{L})^{-1} + \sigma_2\theta(P)\theta(P)^T. $$

This differs from our Corollary 6, which gives

$$ \mathrm{Var}_{\mathrm{LGRG}}(\theta_n) = 2\sigma_1(\mathbf{L}^T(\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1})\mathbf{L})^{-1} + \sigma_2\theta(P)\theta(P)^T. $$

For the choices of **X**, **Z** and **R** in [5], both expressions are equal. However, for the alternative choice for **R** made in alternative 2, one finds

$$ \mathrm{Var}_{\mathrm{CVF}}(\theta_n) = \begin{pmatrix} 20.63 & -1.22 \\ -1.22 & 1.77 \end{pmatrix} \quad \text{and} \quad \mathrm{Var}_{\mathrm{LGRG}}(\theta_n) = \begin{pmatrix} 8.57 & -0.82 \\ -0.82 & 0.80 \end{pmatrix}. $$

Again the differences are quite large. For example, as a consequence the length of the confidence interval for $\theta_1$ based on $\mathrm{Var}_{\mathrm{CVF}}(\theta_n)$ will be 1.5 times larger than the one based on $\mathrm{Var}_{\mathrm{LGRG}}(\theta_n)$.

Finally, we illustrate the performance of S-estimators by an application to data from a trial on the treatment of lead-exposed children. This dataset is discussed in [10] and consists of four repeated measurements of blood lead levels obtained at baseline (or week 0), week 1, week 4, and week 6 on 100 children who were randomly assigned to chelation treatment with succimer (a chelation agent) or placebo. On the basis of a graphical display of the mean response over time, it is suggested in [10] that a quadratic trend over time seems suitable. We fitted the following model

$$ y_{ij} = \beta_0 + \beta_1\delta_i + (\beta_3 + \beta_4\delta_i)t_j + (\beta_5 + \beta_6\delta_i)t_j^2 + \gamma_{1i} + \gamma_{2i}t_j + \gamma_{3i}t_j^2 + \epsilon_{ij}, $$

for $i = 1, \dots, 100$ and $j = 1, \dots, 4$, where $(t_1, \dots, t_4) = (0, 1, 4, 6)$ refer to the different weeks, $y_{ij}$ is the blood lead level (mcg/dL) of subject $i$ obtained at time $t_j$, and $\delta_i = 0$ if the $i$-th subject is in the placebo group and $\delta_i = 1$, otherwise. The random effects $\gamma_i = (\gamma_{1i}, \gamma_{2i}, \gamma_{3i})$, $i = 1, \dots, 100$, are assumed to be independent mean zero normal random vectors with a diagonal covariance matrix consisting of variances $\sigma_{\gamma_1}^2$, $\sigma_{\gamma_2}^2$ and $\sigma_{\gamma_3}^2$, respectively. The measurement errors $\epsilon_i = (\epsilon_{i1}, \dots, \epsilon_{i4})$, $i = 1, \dots, 100$, are assumed to be independent mean zero random vectors with covariance matrix $\sigma_{\epsilon}^2\mathbf{I}_4$, also being independent of the random effects. In this way we are fitting a balanced linear
mixed effects model with unknown parameters $\beta = (\beta_1, \dots, \beta_6)$ and $\theta = (\sigma_{\gamma_1}^2, \sigma_{\gamma_2}^2, \sigma_{\gamma_3}^2, \sigma_{\epsilon}^2)$, and a linear covariance structure.

We estimated $(\beta, \theta)$ by means of maximum likelihood and by means of the S-estimator corresponding to Tukey's bi-weight defined in (10.4). The tuning-constant was chosen to be $c = 4.097$, which corresponds to asymptotic breakdown point 0.5. For each estimate $(\hat{\beta}, \hat{\theta})$, we determined the estimate $\mathbf{V}(\hat{\theta})$ for the structured covariance and the standardized residuals for each subject

$$ \mathrm{RES}_i = \sqrt{(\mathbf{y}_i - \mathbf{X}_i\hat{\beta})^T \mathbf{V}(\hat{\theta})^{-1} (\mathbf{y}_i - \mathbf{X}_i\hat{\beta})} $$

The residuals for both estimation procedures are visible in the left picture of Figure 2, with the residuals determined from the S-estimate on the horizontal axis and the ones determined from the ML estimate on the vertical axis. Both estimates identify subject 40 as an outlier, but only the robust S-estimate also clearly identifies observation 98 as outlier. The extreme large observation in week 6 seems to be the reason that observation 40 is identified as outlier by both methods. See the right picture in Figure 2. Observation 98 also seems to deviate from the overall quadratic trend, by having a suspicious low observation in week 6. The corresponding S-residual clearly sticks out from the other S-residuals, whereas this is much less so for the corresponding ML residual.

# A Proofs and technical lemmas

## Proof of Theorem 1

*Proof.* Let $(\beta, \theta) \in \mathbb{R}^q \times \mathbb{R}^l$ satisfy the S-constraint in (3.5). Then from (R1)-(R2) it follows that

$$ P(C(\beta, \mathbf{V}(\theta), c_0)) \ge 1 - \frac{1}{a_0} \int \rho \left( \sqrt{(\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}(\theta)^{-1} (\mathbf{y} - \mathbf{X}\beta)} \right) dP(s) \ge 1 - r. \quad (A.1) $$

Since $1 - r \ge \epsilon$, Lemma 1(i) then implies that $\lambda_k(\mathbf{V}(\theta)) \ge a_1 > 0$. Because

$$ \lim_{m \to \infty} \int \rho (\|\mathbf{y}\|/m) dP(\mathbf{y}, \mathbf{X}) = 0, $$

we can find $m_0 > 0$, such that $\int \rho(\|\mathbf{y}\|/m_0) dP(\mathbf{y}, \mathbf{X}) \le b_0$. Lemma 1(ii) then yields that $\lambda_1(\mathbf{V}(\theta)) \le a_2 < \infty$. Application of Lemma 1(iii), with $a = 1 - r$ together with $(C1_\epsilon)$, implies that $\|\beta\| \le M < \infty$. It follows that $\beta$ is in a compact subset of $\mathbb{R}^q$ and $\mathbf{V}(\theta)$ is in a compact set $K \subset \mathbb{R}^{k \times k}$.

According to (2.9), the mapping $\theta \mapsto \mathbf{V}(\theta)$ is one-to-one, so that we can restrict $\theta$ to the pre-image $\mathbf{V}^{-1}(K)$. Then with conditions (V1) and (V3) it follows that also $\mathbf{V}^{-1}(K)$ is compact in $\mathbb{R}^l$. We conclude that for solving minimization problem (3.5), we can restrict ourselves to a compact set $K' \subset \mathbb{R}^q \times \mathbb{R}^l$. As det$(\mathbf{V}(\theta))$ is a continuous function of $(\beta, \theta)$, due to condition (V1), it must attain a minimum on $K'$. $\square$

## Proof of Corollary 1

*Proof.* Let $\mathbb{P}_n$ be the empirical measure corresponding to the collection $\mathcal{S}_n$. Then $\mathbb{P}_n$ satisfies $(C1_\epsilon)$ for any $0 < \epsilon \le 1 - r$ and satisfies $(C2_\epsilon)$, for $\epsilon = (\kappa(\mathcal{S}_n) + 1)/n$. Clearly $0 < \epsilon \le 1 - r$, where $r = b_0/a_0$, so according to Theorem 1 there exists at least one solution to (3.5) with $P = \mathbb{P}_n$. This means that there exists at least one solution to (3.1). $\square$

## Proof of Corollary 2

*Proof.* First note there exists $0 < \eta < \epsilon' - \epsilon$. According to Ranga Rao [24, Theorem 4.2] we have

$$ \sup_{C \in \mathcal{C}} |P_t(C) - P(C)| \to 0, \quad \text{as } t \to \infty. \quad (A.2) $$
Because strips $H(\alpha, \ell, \delta) \in \mathbb{C}$, property (A.2) implies that every strip with $P_t(H(\alpha, \ell, \delta)) \ge \epsilon + \eta$, for $t$ sufficiently large, must also satisfy $P(H(\alpha, \ell, \delta)) \ge \epsilon$. This means that

$$ \inf \{\delta : P_t(H(\alpha, \ell, \delta)) \ge \epsilon + \eta\} \ge \inf \{\delta : P(H(\alpha, \ell, \delta)) \ge \epsilon\} > 0. $$

It follows that, for $t$ sufficiently large, $P_t$ satisfies condition (C2$_{\epsilon+\eta}$). Next, consider the compact set $K$ from (C1$_{\epsilon}$'). Without loss of generality we may assume that it belongs to $\mathbb{C}$. Therefore, as $P(K) \ge r + \epsilon'$, for $t$ sufficiently large $P_t(K) \ge r + \epsilon + \eta$. It follows that, for $t$ sufficiently large, $P_t$ satisfies condition (C1$_{\epsilon+\eta}$). Since $\epsilon+\eta < 1-r$, according to Theorem 1 at least one solution $\xi(P_t)$ exists, for $t$ sufficiently large. $\square$

## Proof of Theorem 2

*Proof.* First note that there exists $0 < \eta < \epsilon' - \epsilon$. Denote $\xi(P_t) = \xi_t = (\beta_t, \theta_t)$. Similar to (A.1) we find that $P_t(C(\beta_t, \theta_t, c_0)) \ge 1-r$. Therefore, as $C(\beta_t, \theta_t, c_0) \in \mathbb{C}$ and $1-r > 1-r-\eta$, it follows from (A.2) that

$$ P(C(\beta_t, \theta_t, c_0)) \ge P_t(C(\beta_t, \theta_t, c_0)) - \sup_{C \in \mathcal{C}} |P_t(C) - P(C)| \ge 1 - r - \eta, \quad (A.3) $$

for $t$ sufficiently large. Since $1-r-\eta > \epsilon$, this means that, according to Lemma 1(i), there exists $a_1 > 0$ only depending on $c_0$ and $P$, such that for $t$ sufficiently large,

$$ \lambda_k(\mathbf{V}(\theta_t)) \ge a_1 > 0. $$

Denote $\xi(P) = \xi_0 = (\beta_0, \theta_0)$ and let $d_0(s) = d(s, \beta_0, \theta_0)$. Then according to Lemma B.1 in [16], for any $\sigma > -1$,

$$ \int \rho \left( \frac{d_0(s)}{1+\sigma} \right) dP_t(s) \to \int \rho \left( \frac{d_0(s)}{1+\sigma} \right) dP(s), $$

as $t \to \infty$. As the limit is strictly decreasing at $\sigma = 0$, and $\xi_0$ satisfies the constraint in (3.5), we find that for all $\sigma > 0$,

$$ \int \rho \left( \frac{d_0(s)}{1+\sigma} \right) dP_t(s) \le b_0, $$

for $t$ sufficiently large. Hence, similar to the proof of Lemma 1(ii) we conclude that for any possible solution $\xi_t = (\beta_t, \theta_t)$, it must hold that

$$ \det(\mathbf{V}(\theta_t)) \le (1+\sigma)^{2k}\det(\mathbf{V}(\theta_0)), $$

for $t$ sufficiently large. As $\sigma > 0$ can be taken arbitrarily small, we conclude that

$$ \limsup_{t \to \infty} \det(\mathbf{V}(\theta_t)) \le \det(\mathbf{V}(\theta_0)), \quad (A.4) $$

and we find that $\lambda_1(\mathbf{V}(\theta_t)) \le \det(\mathbf{V}(\theta_0))/a_1^{k-1} < \infty$, for $t$ large sufficiently large. Finally, let $K$ be the compact set from (C1$_{\epsilon}$'), so that $P(K) \ge r + \epsilon' > r + \epsilon + \eta$. Then, according to (A.3), it follows from Lemma 1(iii) with $a = 1-r-\eta$, that there exists $0 < M < \infty$ such that $\|\beta_t\| \le M$, for $t$ sufficiently large. This means that there exists a compact set $K'$, such that for $t$ sufficiently large the sequence $\{(\beta_t, \mathbf{V}(\theta_t))\} \subset K'$. Then, similar to the second part of the proof of Theorem 1, the conditions on the mapping $\theta \mapsto \mathbf{V}(\theta)$ yield that there exists a compact set $K'' \subset \mathbb{R}^{q+l}$, such that for $t$ sufficiently large, the sequence $\{\xi_t\} \subset K''$.

Consider a convergent subsequence $\{\xi_{t_j}\}$ with $\xi_{t_j} \to \xi_L$. With Lemma B.1 in [16] and the fact that $\xi_{t_j}$ satisfies the S-constraint in (3.5) at $P = P_{t_j}$, we find

$$ \int \rho(d(s, \xi_L)) dP(s) = \lim_{j \to \infty} \int \rho(d(s, \xi_{t_j})) dP_{t_j}(s) \le b_0. $$

Hence, $\xi_L$ satisfies the S-constraint in (3.5), which has solution $\xi_0$. This means that $\det(\mathbf{V}(\theta_L)) \ge \det(\mathbf{V}(\theta_0))$. But then from (A.4), it follows $\det(\mathbf{V}(\theta_L)) = \det(\mathbf{V}(\theta_0))$. Uniqueness of $\xi_0$ together with identifiability (2.9) then implies that $\xi_L = \xi_0$. Because $\{\xi_t\}$ eventually stays in a compact set, this means that we must have $\lim_{t \to \infty} \xi_t = \xi_0$. $\square$
### Proof of Corollary 3

*Proof.* We apply Theorem 2 to the sequence $\mathbb{P}_n$, $n = 1, 2, \dots$, of probability measures, where $\mathbb{P}_n$ is the empirical measure corresponding to $(\mathbf{y}_1, \mathbf{X}_1), \dots, (\mathbf{y}_n, \mathbf{X}_n)$. According to the Portmanteau Theorem (e.g., see Theorem 2.1 in [2]), $\mathbb{P}_n$ converges weakly to $P$, with probability one. The corollary then follows from Theorem 2. □

### Proof of Theorem 5

*Proof.* Denote $\xi_{h,s_0} = \xi(P_{h,s_0})$. This solution satisfies the score equation (7.4) for the regression S-functional at $P_{h,s_0}$, that is

$$ \int \Psi(\mathbf{s}, \xi_{h,s_0}) dP_{h,s_0}(\mathbf{s}) = 0. $$

We decompose as follows

$$ \begin{aligned} \mathbf{0} &= \int \Psi(\mathbf{s}, \xi_{h,s_0}) dP_{h,s_0}(\mathbf{s}) \\ &= (1-h) \int \Psi(\mathbf{s}, \xi_{h,s_0}) dP(\mathbf{s}) + h \Psi(\mathbf{s}_0, \xi_{h,s_0}) \\ &= (1-h) \Lambda(\xi_{h,s_0}) + h \left( \Psi(\mathbf{s}_0, \xi_{h,s_0}) - \Psi(\mathbf{s}_0, \xi(P)) \right) + h \Psi(\mathbf{s}_0, \xi(P)). \end{aligned} $$

We first determine the order of $\xi_{h,s_0} - \xi(P)$, as $h \downarrow 0$. Because $\xi \mapsto \Psi(\mathbf{s}_0, \xi)$ is continuous, it follows that

$$ \Psi(\mathbf{s}_0, \xi_{h,s_0}) = \Psi(\mathbf{s}_0, \xi(P)) + o(1), \quad \text{as } h \downarrow 0. $$

Furthermore, because $\xi \mapsto \Lambda(\xi)$ is continuously differentiable at $\xi(P)$, we have that

$$ \Lambda(\xi_{h,s_0}) = \Lambda(\xi(P)) + \mathbf{D}(P)(\xi_{h,s_0} - \xi(P)) + o(\|\xi_{h,s_0} - \xi(P)\|). $$

Since $\xi(P)$ is the S-functional at $P$, it is a zero of the corresponding score equation, i.e., $\Lambda(\xi(P)) = 0$. It follows that

$$ \mathbf{0} = (1-h)\mathbf{D}(P)(\xi_{h,s_0} - \xi(P)) + o(\|\xi_{h,s_0} - \xi(P)\|) + o(h) + h\Psi(\mathbf{s}_0, \xi(P)). $$

Because $\mathbf{D}(P)$ is non-singular and $\Psi(\mathbf{s}_0, \xi(P))$ is fixed, this implies $\xi_{h,s_0} - \xi(P) = O(h)$. After inserting this in the previous equality, it follows that

$$ \begin{aligned} \mathbf{0} &= (1-h)\mathbf{D}(P)(\xi_{h,s_0} - \xi(P)) + h\Psi(\mathbf{s}_0, \xi(P)) + o(h) \\ &= \mathbf{D}(P)(\xi_{h,s_0} - \xi(P)) + h\Psi(\mathbf{s}_0, \xi(P)) + o(h). \end{aligned} $$

We conclude

$$ \frac{\xi_{h,s_0} - \xi(P)}{h} = -\mathbf{D}(P)^{-1}\Psi(\mathbf{s}_0, \xi(P)) + o(1), \quad \text{as } h \downarrow 0. $$

This means that the limit of the left hand side exists and

$$ \mathrm{IF}(\mathbf{s}_0; \xi, P) = \lim_{h \downarrow 0} \frac{\xi((1-h)P + h\delta_{\mathbf{s}_0}) - \xi(P)}{h} = -\mathbf{D}(P)^{-1}\Psi(\mathbf{s}_0, \xi(P)). $$

Next, consider the covariance functional $\mathbf{C}(P) = \mathbf{V}(\theta(P))$. By definition

$$ \mathrm{IF}(\mathbf{s}_0; \mathrm{vec}(\mathbf{C}), P) = \lim_{h \downarrow 0} \frac{\mathrm{vec}(\mathbf{C}(P_{h,s_0})) - \mathrm{vec}(\mathbf{C}(P))}{h}. $$

Due to (V4), by applying the chain rule (e.g., see [20, Theorem 12, page 108]), we find

$$ \begin{aligned} \mathrm{vec}(\mathbf{C}(P_{h,s_0})) - \mathrm{vec}(\mathbf{C}(P)) &= \frac{\partial \mathrm{vec}(\mathbf{V}(\theta(P)))}{\partial \theta^T} (\theta(P_{h,s_0}) - \theta(P)) + o(\|\theta(P_{h,s_0}) - \theta(P)\|) \\ &= \left( \frac{\partial \mathrm{vec}(\mathbf{V}(\theta(P)))}{\partial \theta^T} \right) (\theta(P_{h,s_0}) - \theta(P)) + o(h). \end{aligned} $$

Dividing by $h$ and letting $h \downarrow 0$, finishes the proof. □
## Proof of Theorem 6

*Proof.* Write $\xi_P = \xi(P)$. Then from (9.2) and Lemma B.8 in [16], it follows that

$$0 = \Lambda(\xi_n) + \int \Psi(s, \xi_P) d(\mathbb{P}_n - P)(s) + o_P(n^{-1/2}). \quad (A.5)$$

From Theorem 2, we know that $\xi_n \to \xi_P$ with probability one. For the first term on the right hand side we have that

$$\Lambda(\xi_n) = \mathbf{D}(P)(\xi_n - \xi_P) + o_P(\|\xi_n - \xi_P\|) = \mathbf{D}(P)(\xi_n - \xi_P) + o_P(\|\xi_n - \xi_P\|),$$

using that $\Lambda(\xi_P) = 0$, due to the fact that $\xi_P$ is the solution of (3.5), see also (7.4). From Lemma B.2 in [16] we have that $\|\Psi_\beta\| \le C_1 \|\mathbf{X}\|$ and $\|\Psi_\theta\| \le C_2$, for universal constants $0 < C_1, C_2 < \infty$. Hence, from the conditions of the theorem, it follows that $\mathbb{E}\|\Psi(s, \xi(P))\|^2 < \infty$. This means that for the second term on the right hand side of (A.5),

$$\int \Psi(s, \xi_P) d(\mathbb{P}_n - P)(s) = \frac{1}{n} \sum_{i=1}^{n} (\Psi(s_i, \xi_P) - \mathbb{E}\Psi(s, \xi_P)) = O_P(n^{-1/2}),$$

according to the central limit theorem. It follows that

$$0 = \mathbf{D}(P)(\xi_n - \xi_P) + o_P(\|\xi_n - \xi_P\|) + O_P(n^{-1/2}),$$

so that $\|\xi_n - \xi_P\| = O_P(n^{-1/2})$. If we insert this is in (A.5), we obtain

$$0 = \mathbf{D}(P)(\xi_n - \xi_P) + \frac{1}{n} \sum_{i=1}^{n} (\Psi(s_i, \xi_P) - \mathbb{E}\Psi(s, \xi_P)) + o_P(n^{-1/2}),$$

from which it follows that

$$\sqrt{n}(\xi_n - \xi_P) = -\mathbf{D}(P)^{-1} \sqrt{n} \left( \frac{1}{n} \sum_{i=1}^{n} (\Psi(s_i, \xi_P) - \mathbb{E}\Psi(s, \xi_P)) \right) + o_P(1).$$

After application of the central limit theorem, we conclude that $\sqrt{n}(\xi_n - \xi_P)$ is asymptotically normal with mean zero and covariance matrix

$$\mathbf{D}(P)^{-1} \mathbb{E} [\Psi(s, \xi_P) \Psi(s, \xi_P)^T] \mathbf{D}(P)^{-1},$$

where we use that $\mathbb{E}\Psi(s, \xi_P) = \Lambda(\xi_P) = 0$. This proves the theorem. $\square$

## References

[1] C. Agostinelli and V. J. Yohai. Composite robust estimators for linear mixed models. *Journal of the American Statistical Association*, 111(516):1764–1774, 2016.

[2] P. Billingsley. *Convergence of probability measures*. John Wiley & Sons, Inc., New York-London-Sydney, 1968.

[3] I. Chervoneva and M. Vishnyakov. Generalized s-estimators for linear mixed effects models. *Statistica Sinica*, 24(3):1257–1276, 2014.

[4] S. Copt and S. Heritier. Robust alternatives to the f-test in mixed linear models based on mm-estimates. *Biometrics*, 63(4):1045–1052, 2007.

[5] S. Copt and M.-P. Victoria-Feser. High-breakdown inference for mixed linear models. *Journal of the American Statistical Association*, 101(473):292–300, 2006.
[6] P. L. Davies. Asymptotic behaviour of S-estimates of multivariate location parameters and dispersion matrices. *Ann. Statist.*, 15(3):1269–1292, 1987.

[7] E. Demidenko. *Mixed models*. Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, second edition, 2013. Theory and applications with R.

[8] D. Donoho and P. J. Huber. The notion of breakdown point. In *A Festschrift for Erich L. Lehmann*, Wadsworth Statist./Probab. Ser., pages 157–184. Wadsworth, Belmont, CA, 1983.

[9] M. V. Fasano, R. A. Maronna, M. Sued, and V. J. Yohai. Continuity and differentiability of regression M functionals. *Bernoulli*, 18(4):1284–1309, 2012.

[10] G. M. Fitzmaurice, N. M. Laird, and J. H. Ware. *Applied longitudinal analysis*. Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, second edition, 2011.

[11] F. R. Hampel. The influence curve and its role in robust estimation. *J. Amer. Statist. Assoc.*, 69:383–393, 1974.

[12] H. O. Hartley and J. N. K. Rao. Maximum-likelihood estimation for the mixed analysis of variance model. *Biometrika*, 54:93–108, 1967.

[13] S. Heritier, E. Cantoni, S. Copt, and M.-P. Victoria-Feser. *Robust methods in biostatistics*. Wiley Series in Probability and Statistics. John Wiley & Sons, Ltd., Chichester, 2009.

[14] R. I. Jennrich and M. D. Schluchter. Unbalanced repeated-measures models with structured covariance matrices. *Biometrics*, 42(4):805–820, 1986.

[15] N. M. Laird and J. H. Ware. Random-effects models for longitudinal data. *Biometrics*, 38(4):963–974, 1982.

[16] H. Lopuhaä, V. Gares, and A. Ruiz-Gazen. Supplement to “S-estimation in linear models with structured covariance matrices”. 2022.

[17] H. P. Lopuhaä. On the relation between S-estimators and M-estimators of multivariate location and covariance. *Ann. Statist.*, 17(4):1662–1683, 1989.

[18] H. P. Lopuhaä. Asymptotic expansion of S-estimators of location and covariance. *Statist. Neerlandica*, 51(2):220–237, 1997.

[19] H. P. Lopuhaä and P. J. Rousseeuw. Breakdown points of affine equivariant estimators of multivariate location and covariance matrices. *Ann. Statist.*, 19(1):229–248, 1991.

[20] J. R. Magnus and H. Neudecker. *Matrix differential calculus with applications in statistics and econometrics*. Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. John Wiley & Sons, Ltd., Chichester, 1988.

[21] D. Nolan and D. Pollard. *U*-processes: rates of convergence. *Ann. Statist.*, 15(2):780–799, 1987.

[22] J. C. Pinheiro, C. Liu, and Y. N. Wu. Efficient algorithms for robust estimation in linear mixed-effects models using the multivariate *t* distribution. *J. Comput. Graph. Statist.*, 10(2):249–276, 2001.

[23] D. Pollard. *Convergence of stochastic processes*. Springer Series in Statistics. Springer-Verlag, New York, 1984.

[24] R. Ranga Rao. Relations between weak and uniform convergence of measures with applications. *Ann. Math. Statist.*, 33:659–680, 1962.

[25] C. R. Rao. Estimation of variance and covariance components in linear models. *J. Amer. Statist. Assoc.*, 67:112–115, 1972.
[26] P. Rousseeuw. Multivariate estimation with high breakdown point. In *Mathematical statistics and applications, Vol. B (Bad Tatzmannsdorf, 1983)*, pages 283–297. Reidel, Dordrecht, 1985.

[27] P. Rousseeuw and V. Yohai. Robust regression by means of S-estimators. In *Robust and nonlinear time series analysis (Heidelberg, 1983)*, volume 26 of *Lect. Notes Stat.*, pages 256–272. Springer, New York, 1984.

[28] P. J. Rousseeuw. Least median of squares regression. *J. Amer. Statist. Assoc.*, 79(388):871–880, 1984.

[29] S. Van Aelst and G. Willems. Multivariate regression S-estimators for robust estimation and inference, 2004, preprint received by personal communication.

[30] S. Van Aelst and G. Willems. Multivariate regression S-estimators for robust estimation and inference. *Statist. Sinica*, 15(4):981–1001, 2005.
# B Supplemental Material

For later use we first define two important matrix norms and mention some useful properties. For $m \times n$ real-valued matrices $\mathbf{A}$, we define the Euclidean norm or Frobenius norm as

$$ \|\mathbf{A}\| = \sqrt{\sum_{i=1}^{m} \sum_{j=1}^{n} a_{ij}^2}. $$

and the spectral norm by

$$ \|\mathbf{A}\|_2 = \sup_{\|\mathbf{u}\|=1} \|\mathbf{A}\mathbf{u}\|. $$

Recall that for real-valued $\mathbf{A}$, the largest eigenvalue is defined by

$$ \lambda_1(\mathbf{A}) = \sup_{\|\mathbf{u}\|=1} \mathbf{u}^T \mathbf{A} \mathbf{u}. \qquad (B.1) $$

This means that $\|\mathbf{A}\|_2^2 = \lambda_1(\mathbf{A}^T \mathbf{A})$. Other useful properties are

$$ \|\mathbf{A}\|_2 \le \|\mathbf{A}\| \le \sqrt{\min(m, n)} \|\mathbf{A}\|_2 \qquad (B.2) $$

and

$$ \|\mathbf{A}\|^2 = \mathrm{tr}(\mathbf{A}^T \mathbf{A}) = \mathrm{tr}(\mathbf{A} \mathbf{A}^T), \qquad (B.3) $$

and

$$ \|\mathbf{u}\mathbf{v}^T\| = \|\mathbf{u}\|\|\mathbf{v}\| \qquad (B.4) $$

for $\mathbf{u} \in \mathbb{R}^m$ and $\mathbf{v} \in \mathbb{R}^n$. When $\mathbf{A}$ is symmetric then $\lambda_1(\mathbf{A}^T \mathbf{A}) = \lambda_1(\mathbf{A}^2) = \lambda_1(\mathbf{A})^2$. In that case

$$ |\lambda_1(\mathbf{A})| = \|\mathbf{A}\|_2 \le \|\mathbf{A}\| \le \sqrt{\min(m, n)} \|\mathbf{A}\|_2 = \sqrt{\min(m, n)} |\lambda_1(\mathbf{A})|. \qquad (B.5) $$

Finally, note that both matrix norms are submultiplicative, that is

$$ \|\mathbf{A}\mathbf{B}\|_2 \le \|\mathbf{A}\|_2 \|\mathbf{B}\|_2, \qquad (B.6) $$

and

$$ \|\mathbf{A}\mathbf{B}\| \le \|\mathbf{A}\|\|\mathbf{B}\|. \qquad (B.7) $$

## B.1 Proofs of Section 4

### Proof of Lemma 1

*Proof.* Note that cylinder $C(\beta, \theta, c)$ is contained in some strip $H(\alpha, \ell, 2c\sqrt{\lambda_k(\mathbf{V}(\theta))})$. It then follows from (C2$_{\epsilon}$) that $\lambda_k(\mathbf{V}(\theta)) \ge \delta_{\epsilon}^2/4c^2$. This proves (i). Let $\theta_0$ be such that the pair $(\beta, \mathbf{V}(\theta_0)) = (0, m_0^2 \mathbf{I}_k)$ satisfies the S-constraint in (3.5), which is possible due to condition (V2). It then follows that for any solution $(\beta, \mathbf{V}(\theta))$ of (3.5), one must have $\det(\mathbf{V}(\theta)) \le m_0^{2k}$. The lower bound on the smallest eigenvalue then implies $\lambda_1(\mathbf{V}(\theta)) \le m_0^{2k}/a_1^{k-1} < \infty$, which proves (ii). Note that

$$ \gamma_{\epsilon} = \inf_{P(J) \ge \epsilon} \inf_{\|\gamma\|=1} \sup_{\mathbf{s} \in J} \|\mathbf{X}\gamma\| > 0, \qquad (B.8) $$

where the infima are taken over all subsets $J \subset \mathbb{R}^k \times \mathcal{X}$ with $P(J) \ge \epsilon$ and all unit vectors $\gamma \in \mathbb{R}^q$. This can be seen as follows. Take $\alpha^T = (\mathbf{0}^T, \gamma^T, \dots, \gamma^T)/k$, where $\mathbf{0}$ is a $k$-vector of zeros, so $\alpha \in \mathbb{R}^k \times \mathbb{R}^{kq}$ and $\|\alpha\| = 1$. Then, with $\ell = 0$, we have $\alpha^T \mathbf{s} - \ell = (\gamma^T \mathbf{x}_1 + \dots + \gamma^T \mathbf{x}_k)/k$. Note that $\|\mathbf{X}\gamma\|^2 = \sum_{j=1}^k (\mathbf{x}_j^T \gamma)^2$. Therefore, if $\gamma_{\epsilon} = 0$, then $\gamma^T \mathbf{x}_j = 0$, for all $j = 1, \dots, k$, which means that $\alpha^T \mathbf{s} - \ell = 0$. This would be in contradiction with (4.3), which is equivalent to (C2$_{\epsilon}$), according to Remark 4.1. Now, note that

$$ P(C(\beta, \theta, c) \cap K) \ge P(C(\beta, \theta, c)) - P(K^c) \ge a - 1 + 1 - a + \epsilon = \epsilon. $$
Hence, according to (B.8), there exists $s_0 = (y_0, X_0) \in C(\beta, \theta, c) \cap K$, such that $\|X_0\gamma\| \ge \gamma_\epsilon > 0$, for all $\gamma$ with $\|\gamma\| = 1$. Because $s_0 \in C(\beta, \theta, c)$, it holds

$$ \|y_0 - X_0\beta\|^2 \le (y_0 - X_0\beta)^T V(\theta)^{-1} (y_0 - X_0\beta) \lambda_1(V(\theta)) \le c\lambda_1(V(\theta)) \le ca_2, $$

due to part (ii). Because $s_0 = (y_0, X_0) \in K$, this means that

$$ \|X_0\beta\| \le \|y_0 - X_0\beta\| + \|y_0\| \le \sqrt{ca_2} + \sup_{(y, X) \in K} \|y\|. $$

Because $s_0 \in C(\beta, \theta, c)$, together with (B.8) we conclude that

$$ \|\beta\| = \|X_0\beta\| \times \frac{\|\beta\|}{\|X_0\beta\|} \le \frac{1}{\gamma_\epsilon} \|X_0\beta\| \le \frac{1}{\gamma_\epsilon} \left( \sqrt{ca_2} + \sup_{(y, X) \in K} \|y\| \right) < \infty. $$

This proves part (iii).

□

### Proof of Remark 4.1

*Proof.* Suppose that $(C2_\epsilon)$ holds and suppose that $\omega_\epsilon = 0$. Then there exists a sequence $(J_n, \alpha_n, \ell_n)$, with $J_n \subset \mathbb{R}^k \times \mathcal{X}$, $\alpha_n \in \mathbb{R}^{k+\ell_n}$, $\|\alpha_n\| = 1$, and $\ell_n \in \mathbb{R}$, such that $P(J_n) \ge \epsilon$, for all $n = 1, 2, \dots$, and

$$ \sup_{s \in J_n} |\alpha_n^T s - \ell_n| \to 0. $$

This means that for $n$ sufficiently large $\delta_n = 2\sup_{s \in J_n} |\alpha_n^T s - \ell_n| < \delta_\epsilon$. Then, there exists a strip $H(\alpha_n, \ell_n, \delta_n)$ containing $J_n$, such that $\delta_n < \delta_\epsilon$ and $P(H(\alpha_n, \ell_n, \delta_n) \cap (\mathbb{R}^k \times \mathcal{X})) \ge P(J_n) \ge \epsilon$. This would be in contradiction with the definition of $\delta_\epsilon$ in (C2$_\epsilon$). On the other hand, suppose that (4.3) holds and suppose that $\delta_\epsilon = 0$. That means that we can find a sequence $(\alpha_n, \ell_n, \delta_n)$, with $\alpha_n \in \mathbb{R}^k \times \mathbb{R}^{\ell_n}$, $\|\alpha_n\| = 1$, and $\ell_n \in \mathbb{R}$, such that $\delta_n \downarrow 0$ and $P(H(\alpha_n, \ell_n, \delta_n) \cap (\mathbb{R}^k \times \mathcal{X})) \ge \epsilon$, for all $n = 1, 2, \dots$. This means that for $n$ sufficiently large $\delta_n < 2\omega_\epsilon$ and then

$$ \sup_{s \in H(\alpha_n, \ell_n, \delta_n) \cap \mathbb{R}^k \times \mathcal{X}} |\alpha_n^T s - \ell_n| = \delta_n/2 < \omega_\epsilon, $$

which would contradict definition (4.3).

□

**Lemma B.1.** Let $P_t$, $t \ge 0$ be a sequence of probability measures on $\mathbb{R}^k \times \mathbb{R}^{\ell_n}$ that converges weakly to $P$, as $t \to \infty$. Let $\xi_t = (\beta_t, \theta_t)$, $t \ge 0$, be a sequence in $\mathbb{R}^q \times \mathbb{R}^l$, such that $\xi_t \to \xi_L$, as $t \to \infty$. If $g(s, \xi) = \rho(d(s, \xi)/\alpha)$, for some $\alpha > 0$ fixed, where $\rho$ satisfies (R2)-(R3) and $V$ satisfies (V1), then

$$ \lim_{t \to \infty} \int g(s, \xi_t) dP_t(s) = \int g(s, \xi_L) dP(s). $$

*Proof.* Let $g_t(s) = g(s, \xi_t)$ and $g_L(s) = g(s, \xi_L)$. Then for every sequence $\{s_t\}$, such that $s_t \to s$, we have

$$ \lim_{t \to \infty} g_t(s_t) = g_L(s). $$

Now, apply Theorem 5.5 from [2]. Let $\Gamma: [0, \infty) \to [0, \infty)$ be the function

$$ \Gamma(u) = u1_{[0, a_0]}(u) + a_0 1_{(a_0, \infty)}(u), $$

which is bounded and uniformly continuous. Then as a consequence of $P_t \to P$ weakly, we have

$$ \lim_{t \to \infty} \int g(s, \xi_t) dP_t(s) = \lim_{t \to \infty} \int \Gamma(g_t(s)) dP_t(s) = \int \Gamma(g_L(s)) dP(s) = \int g(s, \xi_L) dP(s). $$

□
## B.2 Proofs of Section 5

### Proof of Theorem 3

*Proof.* Davies [6] defines location-scale S-estimators by means of a function $\kappa: [0, \infty) \to [0, 1]$. It relates to our $\rho$-function as $\rho(d) = a_0(1 - \kappa(d^2))$. The S-minimization problem considered in [6] can be formulated in our notation as follows

$$ \begin{array}{ll} \min_{\alpha, \mathbf{A}} & \det(\mathbf{A}) \\ \text{subject to} & \displaystyle\int \rho \left( \sqrt{(\mathbf{y} - \alpha)^T \mathbf{A}^{-1} (\mathbf{y} - \alpha)} \right) f_{\mu, \Sigma}(\mathbf{y}) \, d\mathbf{y} \le b_0, \end{array} \qquad (B.9) $$

where the minimum is taken over all $\alpha \in \mathbb{R}^k$ and $\mathbf{A} \in \text{PDS}(k)$. The conditions (R1)-(R2) imply the conditions on $\kappa(s) = 1 - \rho(\sqrt{s})/a_0$ imposed in [6], and $\kappa$ and $h$ have a common point of decrease. It then follows from Theorem 1 in [6] that (B.9) has a unique solution

$$ (\alpha^*, \mathbf{A}^*) = (\mu, \Sigma) = (\mathbf{X}\beta_0, \mathbf{V}(\theta_0)). $$

Since this solution is unique, candidate solutions to (B.9) must be of the form $(\alpha, \mathbf{A}) = (\mathbf{X}\beta, \mathbf{V}(\theta))$, for some $\beta \in \mathbb{R}^q$ and $\theta \in \mathbb{R}^l$. It follows that minimization problem (B.9) is equivalent to minimization problem (5.2). As a consequence, minimization problem (5.2) has a unique solution $(\beta^*, \theta^*)$ for which $\mathbf{X}\beta^* = \mathbf{X}\beta_0$ and $\mathbf{V}(\theta^*) = \mathbf{V}(\theta_0)$. Since $\mathbf{X}^T\mathbf{X}$ is non-singular we can multiply $\mathbf{X}\beta^* = \mathbf{X}\beta_0$ from the left by $(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$. It then follows that $\beta^* = \beta_0$. Finally, from (2.9) we find that $\theta^* = \theta_0$. This proves the theorem. $\square$

## B.3 Proofs for Section 6

### Proof of Theorem 4

*Proof.* Without loss of generality we may assume that $c_0 = 1$ and $\sup \rho = 1$, so that $r = b_0$. The first step is to show that for both estimators $\epsilon_n^* \ge [nr]/n$. To this end, consider a collection $S'_m$ obtained from the original collection $S_n$ by replacing at most $m = [nr] - 1$ number of points in $\mathbb{R}^k \times \mathcal{X}$. We must show that at least one solution $(\beta_n(S'_m), \theta_n(S'_m))$ to the S-minimization problem (3.1) exists for the corrupted collection $S'_m$, and that all possible solutions do not break down.

Denote a possible solution to the S-minimization problem (3.1) corresponding to the corrupted collection $S'_m$, by

$$ \xi'_m = (\beta'_m, \theta'_m) = (\beta_n(S'_m), \theta_n(S'_m)), $$

and consider the corresponding cylinder $C(\beta'_m, \theta'_m, 1)$ defined by (3.3). We apply Lemma 1 to the empirical measure $P'_m$ corresponding to the corrupted collection $S'_m$ of $n$ points. Because $\xi'_m = (\beta'_m, \theta'_m)$ must satisfy the S-constraint in (3.1), one can argue as in (A.1),

$$ \begin{aligned} P'_m(C(\beta'_m, \theta'_m, 1)) &= \frac{1}{n} \sum_{\mathbf{s}_i \in S'_m} 1 \{\mathbf{s}_i \in C(\beta'_m, \theta'_m, 1)\} \\ &\ge 1 - \frac{1}{n} \sum_{\mathbf{s}_i \in S'_m} \rho(d(\mathbf{s}_i, \xi'_m)) \ge 1 - b_0, \end{aligned} $$

where $d(\mathbf{s}_i, \xi'_m)$ is defined in (5.1). It follows that the cylinder $C(\beta'_m, \theta'_m, 1)$ must contain at least $[n - nb_0] = [n - nr]$ number of points from the corrupted collection $S'_m$. Furthermore, since $r \le (n - \kappa(S_n))/(2n)$, for any such subset of $S'_m$ it holds that it contains

$$ [n - nr] - m = n - \lfloor nr \rfloor - \lfloor nr \rfloor + 1 \ge \kappa(S_n) + 1 \qquad (B.10) $$
points of the original collection $S_n$. It follows that the measure $\mathbb{P}'_m$ satisfies condition (C2$\epsilon$), for $\epsilon = (\kappa(S_n) + 1)/n$ and with the value $\delta_\epsilon > 0$ only depending on the original collection $S_n$. Moreover, we also have that $\mathbb{P}'_m(C(\beta'_m, \theta'_m, 1)) \ge \epsilon$, for $\epsilon = (\kappa(S_n) + 1)/n > 0$. According to Lemma 1(i), it then follows that $\lambda_k(\mathbf{V}(\theta'_m)) \ge a_1 > 0$, where $a_1$ only depends on the original collection $S_n$.

Because $nb_0 - m = nr - \lfloor nr \rfloor + 1 > 0$ and

$$ \lim_{R \to \infty} \sum_{(\mathbf{y}_i, \mathbf{X}_i) \in S_n} \rho \left( \frac{\|\mathbf{y}_i\|}{R} \right) = 0, $$

we can find an $R_0 > 0$, only depending on the original collection $S_n$, such that

$$ \sum_{(\mathbf{y}_i, \mathbf{X}_i) \in S_n} \rho \left( \frac{\|\mathbf{y}_i\|}{R_0} \right) \le nb_0 - m. $$

The collection $S'_m$ contains $n-m$ points of the original collection $S_n$. Consider the smallest $M > 0$, such that

$$ \sum_{(\mathbf{y}_i, \mathbf{X}_i) \in S'_m \cap S_n} \rho \left( \frac{\|\mathbf{y}_i\|}{M} \right) \le nb_0 - m. $$

Because $S'_m \cap S_n$ has less points than $S_n$, it holds that $M \le R_0$. It follows that

$$ \begin{aligned} \int \rho \left( \frac{\|\mathbf{y}\|}{R_0} \right) d\mathbb{P}'_m(\mathbf{s}) &= \frac{1}{n} \sum_{(\mathbf{y}_i, \mathbf{X}_i) \in S'_m} \rho \left( \frac{\|\mathbf{y}_i\|}{R_0} \right) \\ &\le \frac{1}{n} \sum_{(\mathbf{y}_i, \mathbf{X}_i) \in S'_m} \rho \left( \frac{\|\mathbf{y}_i\|}{M} \right) \\ &\le \frac{1}{n} \left( \sum_{(\mathbf{y}_i, \mathbf{X}_i) \in S'_m \cap S_n} \rho \left( \frac{\|\mathbf{y}_i\|}{M} \right) + m \right) \le b_0. \end{aligned} $$

According to Lemma 1(ii), it then follows that $\lambda_1(\mathbf{V}(\theta'_m)) \le a_2 < \infty$. where $a_2$ only depends on $a_1$ and the collection $S_n$.

To show that the estimate $\beta'_m$ stays bounded, recall that the cylinder $C(\beta'_m, \theta'_m, 1)$ contains a subset $J_0$ of $\kappa(S_n) + 1$ points from the original collection $S_n$, according to (B.10). By definition, $\kappa(S_n) + 1$ original points cannot be on the same hyperplane, so that

$$ \gamma_n = \inf_{J \subset S_n} \inf_{\|\gamma\|=1} \max_{\mathbf{s} \in J} \|\mathbf{X}\gamma\| > 0. $$

where the first infimum runs over all subsets $J \subset S_n$ of $\kappa(S_n) + 1$ points. By definition of $\gamma_n$, there exists an original point $\mathbf{s}_0 \in J_0 \subset S_n \cap C(\beta'_m, \theta'_m, 1)$, such that

$$ \|\beta'_m\| = \|\mathbf{X}_0\beta'_m\| \times \frac{\|\beta'_m\|}{\|\mathbf{X}_0\beta'_m\|} \le \frac{1}{\gamma_n} \|\mathbf{X}_0\beta'_m\|. $$

Because $\mathbf{s}_0 \in C(\beta'_m, \theta'_m, 1)$, similar to the proof of Lemma 1(iii), it follows that $\|\mathbf{y}_0 - \mathbf{X}_0\beta'_m\|^2 \le a_2$, and because $\mathbf{s}_0 \in S_n$, we have that

$$ \|\mathbf{X}_0\beta'_m\| \le \sqrt{a_2} + \max_{(\mathbf{y}_i, \mathbf{X}_i) \in S_n} \|\mathbf{y}_i\| < \infty. $$

We conclude that there exists a compact set $K_n$, only depending on the original collection $S_n$, that contains the pair $(\beta'_m, \mathbf{V}(\theta'_m))$. Similar to the reasoning in the proof of Theorem 1, it follows that at least one solution $(\beta_n(S'_m), \theta_n(S'_m))$ to the S-minimization problem (3.1) exists for the collection $S'_m$, and that all possible solutions $\beta_n(S'_m)$ and $\theta_n(S'_m)$ do not break down.
We continue by showing $\epsilon_n^*(\theta_n) \le \lceil nr \rceil/n$. Replace $m = \lceil nr \rceil$ points of $S_n$ to obtain a corrupted collection $S'_m$ of $n$ points. Suppose that a solution $\xi'_m = (\beta'_m, \theta'_m) = (\beta_n(S'_m), \theta_n(S'_m))$ exists to the S-minimization problem (3.1) corresponding to the corrupted collection $S'_m$. We must show that the estimate $\theta'_m$ breaks down. Note that the estimates $\beta'_m$ and $\theta'_m$ satisfy the S-constraint in (3.1) for the corrupted collection,

$$ \sum_{s_i \in S'_m} \rho(d(s_i, \xi'_m)) \le nr. \qquad (B.11) $$

If all $m = \lceil nr \rceil$ replaced points are outside the cylinder $C(\beta'_m, \theta'_m, 1)$, then

$$ \sum_{s_i \in S'_m} \rho(d(s_i, \xi'_m)) = \sum_{s_i \in S'_m \cap S_n} \rho(d(s_i, \xi'_m)) + \lceil nr \rceil > nr, $$

when $nr \notin \mathbb{N}$. When $nr \in \mathbb{N}$, then by assumption $n - m = n - \lceil nr \rceil \ge \kappa(S_n) + 1$. Hence, there is at least one point $s_i \in S_n$, for which $d(s_i, \xi'_m) > 0$. Because $\rho$ is strictly increasing on $[0, 1]$, this implies

$$ \sum_{s_i \in S'_m} \rho(d(s_i, \xi'_m)) > nr. $$

We conclude that at least one replaced point $s'_i = (y'_i, X'_i) \in C(\beta'_m, \theta'_m, 1)$. Similarly, if all original points in $S'_m$ are outside $C(\beta'_m, \theta'_m, 1)$, then

$$ \sum_{s_i \in S'_m} \rho(d(s_i, \xi'_m)) \ge n - \lceil nr \rceil \ge n - \lceil n - nr \rceil + \kappa(S_n) > nr, $$

which is in contradiction with (B.11). Therefore, the cylinder $C(\beta'_m, \theta'_m, 1)$ must contain a point $s_0 = (y_0, X_0)$ from the original collection $S_n$ as well as one replaced point $s'_i = (y'_i, X'_i)$.

Note that for each point $s = (y, X) \in C(\beta'_m, \theta'_m, 1)$ it holds that

$$ \frac{|\mathbf{q}_1^T(y - \mathbf{X}\beta'_m)|^2}{\lambda_1(\mathbf{V}(\theta'_m))} + \dots + \frac{|\mathbf{q}_k^T(y - \mathbf{X}\beta'_m)|^2}{\lambda_k(\mathbf{V}(\theta'_m))} \le 1, \qquad (B.12) $$

where $\lambda_j(\mathbf{V}(\theta'_m)) > 0$, for $j = 1, \dots, k$, are the eigenvalues of $\mathbf{V}(\theta'_m)$ and $\mathbf{q}_1, \dots, \mathbf{q}_k$ are the corresponding orthonormal eigenvectors. Now, replace all $m$ points by $s'_i = (y'_i, X'_i) = (\mathbf{z}, \mathbf{0})$, where $\mathbf{0}$ is a $k \times q$ matrix of zeros and $\mathbf{z} = t \sum_{j=1}^k \mathbf{q}_j$, so that $\mathbf{q}_j^T \mathbf{z} = t$, for each $j = 1, \dots, k$. By sending $t \to \infty$ and the fact that at least one replaced point $(\mathbf{z}, \mathbf{0})$ satisfies (B.12), it follows that $\lambda_j(\mathbf{V}(\theta'_m)) \to \infty$, for each $j = 1, \dots, k$. This means $\theta'_m$ breaks down.

The upper bound for $\epsilon_n^*(\beta_n, S_n)$ follows from the fact that $\beta_n$ is regression equivariant. Similar to Theorem 2 in [19] it can be shown that the maximal breakdown point of regression equivariant estimators is $\lfloor (n+1)/2 \rfloor/n$. This proves the theorem. $\square$

## B.4 Proofs for Section 8

**Lemma B.2.** Suppose that $\rho$ satisfies (R2), (R4) and $\mathbf{V}$ satisfies (V4). Let $\Psi = (\Psi_\beta, \Psi_\theta)$, as defined in (7.5). Then there exist $0 < C_1 < \infty$ and $0 < C_2 < \infty$, only depending on $P$, such that $\|\Psi_\beta(\mathbf{s}, \xi(P))\| \le C_1 \|\mathbf{X}\|$ and $\|\Psi_\theta(\mathbf{s}, \xi(P))\| \le C_2$.

*Proof.* Consider the expression of $\Psi_\beta$ in (7.5). Consecutively, we apply (B.3), (B.1), (B.5), (B.6), and (B.2). This gives

$$ \begin{aligned} \|\mathbf{X}^T \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\beta)\|^2 &= (\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}^{-1} \mathbf{X} \mathbf{X}^T \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\beta) \\ &\le \lambda_1(\mathbf{V}^{-1/2} \mathbf{X} \mathbf{X}^T \mathbf{V}^{-1/2})(\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\beta) \\ &= d^2 \|\mathbf{V}^{-1/2} \mathbf{X} \mathbf{X}^T \mathbf{V}^{-1/2}\|_2 \\ &\le d^2 \|\mathbf{V}^{-1/2}\|_2^2 \|\mathbf{X} \mathbf{X}^T\|_2 \\ &\le d^2 \|\mathbf{X}\|^2 \lambda_1(\mathbf{V}^{-1}), \end{aligned} \qquad (B.13) $$
where $d = d(\mathbf{s}, \boldsymbol{\xi})$ as defined by (5.1), and where we abbreviate $\mathbf{V}(\boldsymbol{\theta})$ by $\mathbf{V}$. This means that

$$ \|\Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi}(P))\| \le |u(d)d| \times \|\mathbf{X}\| \times \sqrt{\lambda_1(\mathbf{V}(\boldsymbol{\theta}(P)))^{-1}}. $$

From (R2) and (R4) it follows that $u(s)s = \rho'(s)$ is bounded. This means that there exists a universal constant $0 < C_1 < \infty$, such that

$$ \|\Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi}(P))\| \le C_1 \|\mathbf{X}\|. $$

For $\Psi_{\boldsymbol{\theta}} = (\Psi_{\boldsymbol{\theta},1}, \dots, \Psi_{\boldsymbol{\theta},l})$, we have

$$ \Psi_{\boldsymbol{\theta},j}(\mathbf{s}, \boldsymbol{\xi}) = u(d)(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) - \text{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \right) (\rho(d) - b_0), $$

for $j = 1, \dots, l$, where we write $\mathbf{V}$ instead of $\mathbf{V}(\boldsymbol{\theta})$. Recall that

$$ \mathbf{H}_j = \text{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \right) \left( \sum_{t=1}^{l} \theta_t \frac{\partial \mathbf{V}}{\partial \theta_t} \right) - \text{tr} \left( \mathbf{V}^{-1} \sum_{t=1}^{l} \theta_t \frac{\partial \mathbf{V}}{\partial \theta_t} \right) \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j}. $$

To bound $\mathbf{H}_j$, we first obtain a bound on

$$ (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}). \quad (B.14) $$

Note that $\partial \mathbf{V}/\partial \boldsymbol{\theta}_j$ is symmetric, but not necessarily positive definite. When (B.14) is positive, then application of (B.1) gives

$$ \begin{aligned} 0 < (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) &\le d^2 \lambda_1 \left( \mathbf{V}^{-1/2} \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \mathbf{V}^{-1/2} \right) \\ &\le d^2 \left\| \mathbf{V}^{-1/2} \right\|^2 \left\| \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \right\| \\ &\le d^2 \left\| \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \right\| \lambda_1(\mathbf{V}^{-1}), \end{aligned} $$

according to (B.7) and (B.5). When (B.14) is negative, then similarly

$$ \begin{aligned} 0 < (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \left(-\frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j}\right) \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) &\le d^2 \lambda_1 \left( \mathbf{V}^{-1/2} \left(-\frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j}\right) \mathbf{V}^{-1/2} \right) \\ &\le d^2 \left\| \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \right\| \lambda_1(\mathbf{V}^{-1}). \end{aligned} $$

It follows that

$$ \left| (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \right| \le d^2 \left\| \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \right\| \lambda_1(\mathbf{V}^{-1}). \quad (B.15) $$

Furthermore, according to (V4), the mapping $\boldsymbol{\theta} \mapsto \mathbf{V}(\boldsymbol{\theta})$ is continuously differentiable. This means that there exists a universal constant $0 < M_1 < \infty$, such that

$$ \max_{1 \le j \le l} \left\| \frac{\partial \mathbf{V}(\boldsymbol{\theta}(P))}{\partial \boldsymbol{\theta}_j} \right\| \le M_1. \quad (B.16) $$

Because $\partial \mathbf{V}/\partial \boldsymbol{\theta}_j$ is symmetric, according to (B.5) and (B.6)

$$ \left| \lambda_1 \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \right) \right| \le \| \mathbf{V}^{-1} \|_2 \left\| \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \right\|_2 \le \sqrt{k} \lambda_1 (\mathbf{V}^{-1}) \left\| \frac{\partial \mathbf{V}}{\partial \boldsymbol{\theta}_j} \right\|. $$
Together with (B.16), we find

$$ \left| \mathrm{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) \right| \leq k \left| \lambda_1 \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) \right| \leq k^{3/2} M_1 \lambda_1 (\mathbf{V}^{-1}), \quad (B.17) $$

where we abbreviate $\mathbf{V}(\theta(P))$ by $\mathbf{V}$. It then follows there exists a constant $0 < M_2 < \infty$, only depending on $P$, such that at $\theta(P)$,

$$ \max_{1 \le j \le l} |(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})| \le d^2 M_2. \quad (B.18) $$

From (R2) and (R4) it follows that $u(s)s^2 = \rho'(s)s$ and $\rho(s) - b_0$ are bounded. Together with (B.17), it follows that there exists a universal constant $0 < C_2 < \infty$, such that

$$ \|\Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi}(P))\| \le C_2, $$

for all $j = 1, \dots, l$. This finishes the proof. $\square$

### Proof of Corollary 4

*Proof.* Take $\mathbf{s} = (\mathbf{y}, \mathbf{X})$ fixed and consider $\mathrm{IF}(\mathbf{s}; \boldsymbol{\xi}, P)$. Since $\mathbf{D}(P)$ does not depend on $\mathbf{s}$, from Theorem 5 and Lemma B.2, it follows immediately that $\mathrm{IF}(\mathbf{s}; \boldsymbol{\xi}, P)$ remains bounded in $\mathbf{y}$, but not necessarily in $\mathbf{X}$. $\square$

**Lemma B.3.** Consider $\Lambda$ as defined by (8.2) with $\Psi$ defined in (7.5). Suppose that $\rho$ satisfies (R2) and (R5) and $\mathbf{V}$ satisfies (V5). Furthermore, suppose that $\mathbb{E}\|\mathbf{X}\|^2 < \infty$. Let $\boldsymbol{\xi}(P)$ be a solution to (3.5) and let $N$ be an open neighborhood of $\boldsymbol{\xi}(P)$. Then, $\Lambda$ is continuous differentiable at $\boldsymbol{\xi}(P)$ and for all $\boldsymbol{\xi} \in N$,

$$ \frac{\partial \Lambda(\boldsymbol{\xi})}{\partial \boldsymbol{\xi}} = \int \frac{\partial \Psi(\mathbf{s}, \boldsymbol{\xi})}{\partial \boldsymbol{\xi}} dP(\mathbf{s}). $$

*Proof.* Write $\partial\Lambda/\partial\boldsymbol{\xi}$ as the block matrix

$$ \frac{\partial \Lambda(\boldsymbol{\xi})}{\partial \boldsymbol{\xi}} = \begin{pmatrix} \frac{\partial \Lambda_{\beta}(\boldsymbol{\xi})}{\partial \boldsymbol{\beta}} & \frac{\partial \Lambda_{\beta}(\boldsymbol{\xi})}{\partial \boldsymbol{\theta}} \\ \frac{\partial \Lambda_{\theta}(\boldsymbol{\xi})}{\partial \boldsymbol{\beta}} & \frac{\partial \Lambda_{\theta}(\boldsymbol{\xi})}{\partial \boldsymbol{\theta}} \end{pmatrix}, \quad (B.19) $$

where

$$ \Lambda_{\beta}(\boldsymbol{\xi}) = \int \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi}) dP(\mathbf{s}), $$

$$ \Lambda_{\theta}(\boldsymbol{\xi}) = \int \Psi_{\theta}(\mathbf{s}, \boldsymbol{\xi}) dP(\mathbf{s}). $$

We prove the lemma for each block separately. Consider $\partial\Lambda_{\beta}/\partial\boldsymbol{\beta}$. We have

$$ \frac{\partial \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi})}{\partial \boldsymbol{\beta}} = - \frac{u'(d)}{d} \mathbf{X}^T \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{X} - u(d) \mathbf{X}^T \mathbf{V}^{-1} \mathbf{X}, \quad (B.20) $$

where $d = d(\mathbf{s}, \boldsymbol{\xi})$ is defined by (5.1) and where we abbreviate $\mathbf{V}(\boldsymbol{\theta})$ by $\mathbf{V}$. First note that, according to (R5), $\boldsymbol{\xi} \mapsto u'(d(\mathbf{s}, \boldsymbol{\xi}))/d(\mathbf{s}, \boldsymbol{\xi})$ is continuous at $\boldsymbol{\xi}(P)$, for each $\mathbf{s}$ fixed such that $d(\mathbf{s}, \boldsymbol{\xi}(P)) \neq 0$. Together with (V1), this means that for such $\mathbf{s}$ fixed, $\boldsymbol{\xi} \mapsto \partial\Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi})/\partial\boldsymbol{\beta}$ is continuous at $\boldsymbol{\xi}(P)$. For the first term on the right hand side of (B.20), we apply (B.4) and (B.13). This gives

$$ \|\mathbf{X}^T \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{X}\| = \|\mathbf{X}^T \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\|^2 \le d^2 \|\mathbf{X}\|^2 \lambda_1 (\mathbf{V}^{-1}). $$
Similarly, for the second term on the right hand side of (B.20), after application of (B.7) and (B.5), we get

$$ \|\mathbf{X}^T \mathbf{V}^{-1} \mathbf{X}\| \le \|\mathbf{X}\|^2 \|\mathbf{V}^{-1}\| \le \sqrt{q} \|\mathbf{X}\|^2 \lambda_1(\mathbf{V}^{-1}). $$

Since $\lambda_1(\mathbf{V}^{-1})$ is bounded uniformly on the neighborhood $N$ of $\xi(P)$ and because $u(s)$ and $u'(s)s = \rho''(s) - u(s)$ are bounded, due to (R2), it follows that there exists a constant $0 < C_1 < \infty$, only depending on $P$, such that

$$ \left\| \frac{\partial \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi})}{\partial \beta^T} \right\| \le C_1 \|\mathbf{X}\|^2. $$

Since $\mathbb{E}\|\mathbf{X}\|^2 < \infty$, it follows by dominated convergence that for $\boldsymbol{\xi}$ in the neighborhood $N$ of $\xi(P)$, it holds that

$$ \frac{\partial \Lambda_{\beta}(\boldsymbol{\xi})}{\partial \boldsymbol{\beta}} = \int \frac{\partial \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi})}{\partial \boldsymbol{\beta}} dP(\mathbf{s}), \qquad (B.21) $$

and that $\partial \Lambda_{\beta}/\partial \beta$ is continuous at $\xi(P)$.

Next consider $\partial \Psi_{\beta}/\partial \theta$. For each $j = 1, \dots, l$ fixed, we have

$$ \frac{\partial \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi})}{\partial \theta_j} = \frac{u'(d)}{2d} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \cdot \mathbf{X}^T \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ + u(d) \cdot \mathbf{X}^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}). \qquad (B.22) $$

First note that, similar to (B.20), $\boldsymbol{\xi} \mapsto \partial \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi})/\partial \theta_j$ is continuous at $\xi(P)$, for each $\mathbf{s}$ fixed such that $d(\mathbf{s}, \xi(P)) \ne 0$. Consider the first term on the right hand side of (B.22). From (B.13), we have

$$ \|\mathbf{X}^T \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\| \le d \|\mathbf{X}\| \sqrt{\lambda_1(\mathbf{V}^{-1})}. $$

Moreover, similar to the reasoning in (B.15), we find

$$ \left| (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \right| \le d^2 \left\| \frac{\partial \mathbf{V}}{\partial \theta_j} \right\| \lambda_1(\mathbf{V}^{-1}). \qquad (B.23) $$

For the second term on the right hand side of (B.22), similar to the reasoning in (B.13), we have

$$ \begin{aligned} \left\| \mathbf{X}^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \right\|^2 &= (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} \mathbf{X} \mathbf{X}^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ &\le d^2 \lambda_1 \left( \mathbf{V}^{-1/2} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} \mathbf{X} \mathbf{X}^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1/2} \right) \\ &\le d^2 \left\| \mathbf{V}^{-1/2} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1} \mathbf{X} \mathbf{X}^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \mathbf{V}^{-1/2} \right\|_2^2 \\ &\le d^2 \|\mathbf{X}\|^2 \left\| \frac{\partial \mathbf{V}}{\partial \theta_j} \right\|^2 \lambda_1(\mathbf{V}^{-1})^3. \end{aligned} $$

According to (V4), the mapping $\mathbf{V}(\theta)$ is continuously differentiable. This means that $\|\partial \mathbf{V}/\partial \theta_j\|$ is bounded on the neighborhood $N$ of $\xi(P)$. Since $\lambda_1(\mathbf{V}^{-1})$ is bounded uniformly on $N$ and because $u(s)s = \rho'(s)$ and $u'(s)s^2 = \rho''(s)s - \rho'(s)$ are bounded, it follows that there exists a constant $0 < C_2 < \infty$, only depending on $P$, such that for $j = 1, \dots, l$,

$$ \left\| \frac{\partial \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi})}{\partial \theta_j} \right\| \le C_2 \|\mathbf{X}\|^2. $$

As before, it follows by dominated convergence that for $\boldsymbol{\xi}$ in the neighborhood $N$ of $\xi(P)$, it holds that

$$ \frac{\partial \Lambda_{\beta}(\boldsymbol{\xi})}{\partial \boldsymbol{\theta}} = \int \frac{\partial \Psi_{\beta}(\mathbf{s}, \boldsymbol{\xi})}{\partial \boldsymbol{\theta}} dP(\mathbf{s}), \qquad (B.24) $$
and that $\partial\Lambda_{\beta}/\partial\theta$ is continuous at $\xi(P)$.

Next consider $\partial\Psi_{\theta,j}/\partial\beta$, for $j = 1, \dots, l$. We have

$$ \frac{\partial \Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi})}{\partial \boldsymbol{\beta}} = \frac{u'(d)}{d} \mathbf{X}^T \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \cdot (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ - u(d) \cdot 2 \mathbf{X}^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ + \operatorname{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) u(d) \mathbf{X}^T \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}), \qquad (B.25) $$

where $\mathbf{H}_j$ is defined in (7.3). As before, $\xi \mapsto \partial\Psi_{\theta_j}(\mathbf{s}, \boldsymbol{\xi})/\partial\beta$ is continuous at $\xi(P)$, for each $\mathbf{s}$ fixed such that $d(\mathbf{s}, \boldsymbol{\xi}(P)) \neq 0$. Consider the first term on the right hand side of (B.25). From (B.18),

$$ \left| (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \right| \le d^2 M_2. \qquad (B.26) $$

Because $u'(s)s^2$ is bounded, together with (B.13), the norm of the first term on the right hand side of (B.25) is bounded by a constant times $\|\mathbf{X}\|\lambda_1(\mathbf{V}^{-1})^{1/2}$. Similar to (B.13), for the second term on the right hand side of (B.25),

$$ \begin{aligned} \| \mathbf{X}^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \|^2 & \le d^2 \left\| \mathbf{V}^{-1/2} \mathbf{H}_j \mathbf{V}^{-1} \right\|_2^2 \left\| \mathbf{X}\mathbf{X}^T \right\|_2^2 \\ & \le k d^2 \|\mathbf{X}\|^2 \|\mathbf{H}_j\|^2 \lambda_1 (\mathbf{V}^{-1})^3, \end{aligned} $$

and for the third term on the right hand side of (B.25), we can use (B.13) and (B.17). As before, since $u'(s)s = \rho''(s) - u(s)$ and $u(s)s^2 = \rho'(s)s$ are bounded, it follows that there exists a constant $0 < C_3 < \infty$, only depending on $P$, such that for $j = 1, \dots, l$,

$$ \left\| \frac{\partial \Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi})}{\partial \boldsymbol{\beta}} \right\| \le C_3 \|\mathbf{X}\|. $$

Since $\mathbb{E}\|\mathbf{X}\| < \infty$, it follows by dominated convergence that for $\xi$ in the neighborhood $N$ of $\xi(P)$, it holds that

$$ \frac{\partial \Lambda_{\theta}(\boldsymbol{\xi})}{\partial \boldsymbol{\beta}} = \int \frac{\partial \Psi_{\theta}(\mathbf{s}, \boldsymbol{\xi})}{\partial \boldsymbol{\beta}} dP(\mathbf{s}), \qquad (B.27) $$

and that $\partial\Lambda_{\theta}/\partial\beta$ is continuous at $\xi(P)$.

Finally, consider $\partial\Psi_{\theta,j}/\partial\theta_t$, for $j, t = 1, \dots, l$. We find

$$ \begin{aligned} \frac{\partial \Psi_{\theta,j}}{\partial \theta_t} = & -\frac{u'(d)}{2d} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ & \cdot (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ & - u(d) (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} \cdot \mathbf{H}_j \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ & + u(d) (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \frac{\partial \mathbf{H}_j}{\partial \theta_t} \cdot \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ & - u(d) (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \mathbf{H}_j \cdot \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ & + \operatorname{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} \cdot \frac{\partial \mathbf{V}}{\partial \theta_j} \right) (\rho(d) - b_0) \\ & - \operatorname{tr} \left( \mathbf{V}^{-1} \cdot \left( \frac{\partial^2 \mathbf{V}}{\partial \theta_j \partial \theta_t} \right) \right) (\rho(d) - b_0) \\ & + \operatorname{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \cdot \frac{u(d)}{2} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \right). \end{aligned} \qquad (B.28) $$

As before, together with (V5), $\xi \mapsto \partial\Psi_{\theta_j}(\mathbf{s}, \boldsymbol{\xi})/\partial\theta_t$ is continuous at $\xi(P)$, for each $\mathbf{s}$ fixed such that $d(\mathbf{s}, \boldsymbol{\xi}(P)) \neq 0$. From (B.15) and (B.18), it follows that the first term on the right hand side
of (B.28) is bounded by

$$ \frac{|u'(d)d^3|}{2} \left\| \frac{\partial \mathbf{V}}{\partial \theta_t} \right\| \lambda_1(\mathbf{V}^{-1}) M_2. $$

Because $u'(s)s^3 = \rho''(s)s^2 - \rho'(s)s$ is bounded, together with (B.16), we conclude this is bounded on the neighborhood $N$ of $\xi(P)$. Similar to (B.13), the second term on the right hand side of (B.28) is bounded by

$$ |u(d)|d^2 \left\| \mathbf{V}^{-1/2} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1/2} \right\|_2^2 \le |u(d)|d^2 \left\| \frac{\partial \mathbf{V}}{\partial \theta_t} \right\|^2 \|\mathbf{H}_j\|^2 \lambda_1 (\mathbf{V}^{-1})^4. $$

Since $u(s)s^2 = \rho'(s)s$ is bounded, together with (B.17), we again find that this is bounded on the neighborhood $N$ of $\xi(P)$. The same holds for the fourth term on the right hand side of (B.28).

We continue with the third term on the right hand side of (B.28). Similar to (B.13), this is bounded by

$$ |u(d)|d^2 \left\| \mathbf{V}^{-1/2} \frac{\partial \mathbf{H}_j}{\partial \theta_t} \mathbf{V}^{-1/2} \right\|_2^2 \le |u(d)|d^2 \lambda_1 (\mathbf{V}^{-1})^2 \left\| \frac{\partial \mathbf{H}_j}{\partial \theta_t} \right\|. $$

We have that

$$ \begin{aligned} \frac{\partial \mathbf{H}_j}{\partial \theta_t} = & \left\{ \mathrm{tr} \left( -\mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) + \mathrm{tr} \left( \mathbf{V}^{-1} \frac{\partial^2 \mathbf{V}}{\partial \theta_j \partial \theta_t} \right) \right\} \left( \sum_{s=1}^l \theta_s \frac{\partial \mathbf{V}}{\partial \theta_s} \right) \\ & + \mathrm{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) \left\{ \left( \frac{\partial \mathbf{V}}{\partial \theta_t} \right) + \left( \sum_{s=1}^l \theta_s \frac{\partial^2 \mathbf{V}}{\partial \theta_s \partial \theta_t} \right) \right\} \\ & - \left\{ \mathrm{tr} \left( -\mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} \sum_{s=1}^l \theta_s \frac{\partial \mathbf{V}}{\partial \theta_s} \right) + \mathrm{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \right) + \mathrm{tr} \left( \mathbf{V}^{-1} \sum_{s=1}^l \theta_s \frac{\partial^2 \mathbf{V}}{\partial \theta_s \partial \theta_t} \right) \right\} \frac{\partial \mathbf{V}}{\partial \theta_j} \\ & - \mathrm{tr} \left( \mathbf{V}^{-1} \sum_{s=1}^l \theta_s \frac{\partial \mathbf{V}}{\partial \theta_s} \right) \left( \frac{\partial^2 \mathbf{V}}{\partial \theta_j \partial \theta_t} \right). \end{aligned} $$

With (B.5), (B.6), and (B.16), we find

$$ \begin{aligned} \left| \mathrm{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_s} \right) \right| & \le k \left\| \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_s} \right\|_2 \\ & \le k \left\| \mathbf{V}^{-1} \right\|_2^2 \left\| \frac{\partial \mathbf{V}}{\partial \theta_t} \right\| \left\| \frac{\partial \mathbf{V}}{\partial \theta_s} \right\| \\ & \le k \lambda_1 (\mathbf{V}^{-1}) M_1^2, \end{aligned} \quad (B.29) $$

which is uniformly bounded on the neighborhood $N$ of $\xi(P)$, and similar to (B.17) we find

$$ \left| \mathrm{tr} \left( \mathbf{V}^{-1} \frac{\partial^2 \mathbf{V}}{\partial \theta_j \partial \theta_t} \right) \right| \le k^{3/2} \lambda_1 (\mathbf{V}^{-1}) \left\| \frac{\partial^2 \mathbf{V}}{\partial \theta_j \partial \theta_t} \right\|. \quad (B.30) $$

Because, according to (V5), the mapping $\theta \mapsto \mathbf{V}(\theta)$ is twice continuously differentiable, it follows that $\|\partial^2 \mathbf{V}/\partial \theta_j \partial \theta_t\|$ is uniformly bounded on the neighborhood $N$ of $\xi(P)$. Together with the fact that with (B.16),

$$ \left\| \sum_{s=1}^{l} \theta_s \frac{\partial \mathbf{V}}{\partial \theta_s} \right\| \le \sum_{s=1}^{l} \|\theta_s\| \left\| \frac{\partial \mathbf{V}}{\partial \theta_s} \right\| \le M_1 \sum_{s=1}^{l} \|\theta_s\|, $$

it follows that the first term of $\partial \mathbf{H}_j/\partial \theta_t$ is bounded on the neighborhood $N$ of $\xi(P)$. The traces in the other terms can be handled in the same way, which yields that $\|\partial \mathbf{H}_j/\partial \theta_t\|$ is bounded on the neighborhood $N$ of $\xi(P)$. Because $u(s)s^2 = \rho'(s)s$ is bounded, it follows that the third term on the right hand side of (B.28) is bounded.
Next, consider the fifth term on the right hand side of (B.28). From (B.29),

$$ \left| \operatorname{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_t} \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) \right| \leq k \lambda_1 (\mathbf{V}^{-1}) M_1^2, $$

which is uniformly bounded on the neighborhood N of $\xi(P)$. Because $\rho(s)$ is bounded, it follows that the fifth term on the right hand side of (B.28) is bounded. For the sixth term on the right hand side of (B.28), from (B.30) we have

$$ \left| \operatorname{tr} \left( \mathbf{V}^{-1} \frac{\partial^2 \mathbf{V}}{\partial \theta_j \theta_t} \right) \right| \leq k^{3/2} \lambda_1 (\mathbf{V}^{-1}) \left\| \frac{\partial^2 \mathbf{V}}{\partial \theta_j \theta_t} \right\|. $$

Because $\mathbf{V}(\theta)$ is twice continuously differentiable and $\rho(s)$ is bounded, we conclude that the sixth term on the right hand side of (B.28) is bounded. Finally, from (B.17) and (B.23) together with the fact that $u(s)s^2 = \rho'(s)s$ is bounded, it also follows that the last term on the right hand side of (B.28) is bounded. By putting everything together, it follows that there exists a constant $0 < C_4 < \infty$, only depending on P, such that for $j, t = 1, \dots, l$,

$$ \left\| \frac{\partial \Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi})}{\partial \theta_t} \right\| \leq C_4. $$

It follows by dominated convergence that for $\boldsymbol{\xi}$ in the neighborhood N of $\xi(P)$, it holds

$$ \frac{\partial \Lambda_{\theta}(\boldsymbol{\xi})}{\partial \boldsymbol{\theta}} = \int \frac{\partial \Psi_{\theta}(\mathbf{s}, \boldsymbol{\xi})}{\partial \boldsymbol{\theta}} dP(\mathbf{s}), \qquad (B.31) $$

and that $\partial \Lambda_{\theta}/\partial \theta$ is continuous at $\xi(P)$. This finishes the proof. $\square$

For convenience we state the following result from [17] about spherically contoured densities, see Lemma 5.1 in [17]. This lemma uses the commutation matrix $\mathbf{K}_{k,k}$, which is the $k^2 \times k^2$ block matrix with the $(i,j)$-block being equal to the $k \times k$ matrix $\mathbf{\Delta}_{ji}$ consisting of zero's except a 1 at entry $(j,i)$. A useful property (e.g., see [20, Section 3.7]) is that for any $k \times k$ matrix $\mathbf{A}$, it holds that

$$ \mathbf{K}_{k,k} \operatorname{vec}(\mathbf{A}) = \operatorname{vec}(\mathbf{A}^T). \qquad (B.32) $$

**Lemma B.4.** Suppose that $\mathbf{z}$ has a k-variate elliptical contoured density defined in (3.2), with parameters $\boldsymbol{\mu} = \mathbf{0}$ and $\boldsymbol{\Sigma} = \mathbf{I}_k$. Then $\mathbf{u} = \mathbf{z}/\|\mathbf{z}\|$ is independent of $\|\mathbf{z}\|$, has mean zero and covariance matrix $(1/k)\mathbf{I}_k$. Furthermore, $\mathbb{E}_{\mathbf{0},\mathbf{I}_k} \mathbf{u}\mathbf{u}^T \mathbf{u} = \mathbf{0}$ and

$$ \mathbb{E}_{\mathbf{0},\mathbf{I}_k} \operatorname{vec}(\mathbf{u}\mathbf{u}^T) \operatorname{vec}(\mathbf{u}\mathbf{u}^T)^T = \sigma_1(\mathbf{I}_{k^2} + \mathbf{K}_{k,k}) + \sigma_2 \operatorname{vec}(\mathbf{I}_k) \operatorname{vec}(\mathbf{I}_k)^T, $$

$$ \text{where } \sigma_1 = \sigma_2 = (k(k + 2))^{-1}. $$

**Proof of Lemma 2**

*Proof.* Write $\partial \Lambda/\partial \boldsymbol{\xi}$ as in (B.19). We determine each block separately and apply Lemma B.3. Let $\boldsymbol{\xi}_P = (\boldsymbol{\beta}_P, \boldsymbol{\theta}_P) = (\boldsymbol{\beta}(P), \boldsymbol{\theta}(P))$. When we also write $\mathbf{V}_P$ instead of $\mathbf{V}(\boldsymbol{\theta}(P))$, then according to Lemma B.3 we have

$$ \begin{aligned} \frac{\partial \Lambda_{\boldsymbol{\beta}}(\boldsymbol{\xi}_P)}{\partial \boldsymbol{\beta}} &= \int \frac{\partial \Psi_{\boldsymbol{\beta}}(\mathbf{s}, \boldsymbol{\xi}_P)}{\partial \boldsymbol{\beta}} dP(\mathbf{s}) \\ &= -\mathbb{E} \left[ \frac{u'(d)}{d} \mathbf{X}^T \mathbf{V}_P^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}_P) (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}_P)^T \mathbf{V}_P^{-1} \mathbf{X} + u(d) \cdot \mathbf{X}^T \mathbf{V}_P^{-1} \mathbf{X} \right] \\ &= -\mathbb{E} \left[ \mathbb{E} \left[ \frac{u'(d)}{d} \mathbf{X}^T \boldsymbol{\Sigma}^{-1} (\mathbf{y} - \boldsymbol{\mu})(\mathbf{y} - \boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} \mathbf{X} + u(d) \cdot \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \middle| \mathbf{X} \right] \right], \end{aligned} $$
where $d^2 = (\mathbf{y} - \mathbf{X}\beta_P)^T \mathbf{V}_P^{-1} (\mathbf{y} - \mathbf{X}\beta_P)$. The inner expectation on the right hand side is the conditional expectation of $\mathbf{y} | \mathbf{X}$, which has the same distribution as $\mathbf{\Sigma}^{1/2}\mathbf{z} + \mu$, where $\mathbf{z}$ has a spherical density $f_{0,\mathbf{I}_k}$. This implies that the inner expectation on the right hand side is equal to

$$ \mathbf{X}^T \mathbf{\Sigma}^{-1/2} \mathbb{E}_{0,\mathbf{I}_k} \left[ \frac{u'(\|\mathbf{z}\|)}{\|\mathbf{z}\|} \mathbf{z}\mathbf{z}^T + u(\|\mathbf{z}\|) \mathbf{I}_k \right] \mathbf{\Sigma}^{-1/2} \mathbf{X}. $$

From Lemma B.4, we find

$$
\begin{align*}
\mathbb{E}_{0,\mathbf{I}_k} \left[ \frac{u'(\|\mathbf{z}\|)}{\|\mathbf{z}\|} \mathbf{z}\mathbf{z}^T + u(\|\mathbf{z}\|) \mathbf{I}_k \right] &= \mathbb{E}_{0,\mathbf{I}_k} \left[ u'(\|\mathbf{z}\|) \|\mathbf{z}\| \frac{\mathbf{z}\mathbf{z}^T}{\|\mathbf{z}\|^2} + u(\|\mathbf{z}\|) \mathbf{I}_k \right] \\
&= \mathbb{E}_{0,\mathbf{I}_k} \left[ u'(\|\mathbf{z}\|) \|\mathbf{z}\| \right] \mathbb{E}_{0,\mathbf{I}_k} \left[ \mathbf{u}\mathbf{u}^T \right] + \mathbb{E}_{0,\mathbf{I}_k} \left[ u(\|\mathbf{z}\|) \mathbf{I}_k \right] \\
&= \mathbb{E}_{0,\mathbf{I}_k} \left[ u'(\|\mathbf{z}\|) \|\mathbf{z}\| \right] \frac{1}{k} \mathbf{I}_k + \mathbb{E}_{0,\mathbf{I}_k} \left[ u(\|\mathbf{z}\|) \mathbf{I}_k \right] \\
&= \alpha \mathbf{I}_k,
\end{align*}
$$

where $\mathbf{u} = \mathbf{z}/\|\mathbf{z}\|$ and

$$ \alpha = \mathbb{E}_{0,\mathbf{I}_k} \left[ \frac{1}{k} u'(\|\mathbf{z}\|) \|\mathbf{z}\| + u(\|\mathbf{z}\|) \right] = \mathbb{E}_{0,\mathbf{I}_k} \left[ \left(1 - \frac{1}{k}\right) \frac{\rho'(\|\mathbf{z}\|)}{\|\mathbf{z}\|} + \frac{1}{k} \rho''(\|\mathbf{z}\|) \right]. $$

It follows that

$$ \frac{\partial \Lambda_{\beta}(\xi_P)}{\partial \beta} = -\alpha \mathbb{E} \left[ \mathbf{X}^T \mathbf{\Sigma}^{-1} \mathbf{X} \right]. $$

Next, for $j = 1, \dots, l$, consider

$$
\begin{align}
\frac{\partial \Lambda_{\beta}(\xi_P)}{\partial \theta_j} &= \int \frac{\partial \Psi_{\beta}(\mathbf{s}, \xi_P)}{\partial \theta_j} dP(\mathbf{s}) \nonumber \\
&= \mathbb{E} \left[ \mathbb{E} \left[ \frac{u'(d)}{2d} (\mathbf{y} - \mathbf{X}\beta_P)^T \mathbf{V}_P^{-1} \frac{\partial \mathbf{V}_P}{\partial \theta_j} \mathbf{V}_P^{-1} (\mathbf{y} - \mathbf{X}\beta_P) \cdot \mathbf{X}^T \mathbf{V}_P^{-1} (\mathbf{y} - \mathbf{X}\beta_P) \middle| \mathbf{X} \right] \right] \tag{B.33} \\
&\quad + \mathbb{E} \left[ \mathbb{E} \left[ u(d) \cdot \mathbf{X}^T \mathbf{V}_P^{-1} \frac{\partial \mathbf{V}_P}{\partial \theta_j} \mathbf{V}_P^{-1} (\mathbf{y} - \mathbf{X}\beta_P) \middle| \mathbf{X} \right] \right]. \nonumber
\end{align}
$$

According to Lemma B.4, the inner conditional expectation of the first term on the right hand side of (B.33) can be written as

$$
\begin{align*}
\mathbf{X}^T \mathbb{E}_{0,\mathbf{I}_k} \left[ \frac{u'(\|\mathbf{z}\|)}{2\|\mathbf{z}\|} \mathbf{z}^T \mathbf{\Sigma}^{-1/2} \frac{\partial \mathbf{V}_P}{\partial \theta_j} \mathbf{\Sigma}^{-1/2} \mathbf{z} \mathbf{\Sigma}^{-1/2} \mathbf{z} \right] \\
&= \mathbf{X}^T \mathbb{E}_{0,\mathbf{I}_k} \left[ \frac{u'(\|\mathbf{z}\|) \|\mathbf{z}\|^2}{2} \right] \mathbb{E}_{0,\mathbf{I}_k} \left[ \mathbf{u}^T \mathbf{\Sigma}^{-1/2} \frac{\partial \mathbf{V}_P}{\partial \theta_j} \mathbf{\Sigma}^{-1/2} \mathbf{u} \mathbf{\Sigma}^{-1/2} \mathbf{u} \right],
\end{align*}
$$

Since the second term on the right hand side is the expectation with respect to a spherical density of an odd function of $\mathbf{u}$, this expectation is equal to zero due to Lemma B.4. Similarly, the second term on the right hand side of (B.33) has inner conditional expectation

$$ \mathbf{X}^T \mathbb{E}_{0,\mathbf{I}_k} \left[ u(\|\mathbf{z}\|) \mathbf{\Sigma}^{-1} \frac{\partial \mathbf{V}_P}{\partial \theta_j} \mathbf{\Sigma}^{-1/2} \mathbf{z} \right] = \mathbf{X}^T \mathbb{E}_{0,\mathbf{I}_k} \left[ u(\|\mathbf{z}\|) \|\mathbf{z}\| \right] \mathbf{\Sigma}^{-1} \frac{\partial \mathbf{V}_P}{\partial \theta_j} \mathbf{\Sigma}^{-1/2} \mathbb{E}_{0,\mathbf{I}_k} [\mathbf{u}] = \mathbf{0}, $$

due to Lemma B.4. It follows that

$$ \frac{\partial \Lambda_{\beta}(\xi_P)}{\partial \theta} = \mathbf{0}. $$
Next, using Lemma B.3 and (B.25), for all $j = 1, \dots, l$, consider

$$ 
\begin{aligned}
\frac{\partial \Lambda_{\theta,j}}{\partial \beta^T} &= \int \frac{\partial \Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi}_P)}{\partial \beta^T} dP(\mathbf{s}) \\
&= -\mathbb{E} \left[ \mathbb{E} \left[ \frac{u'(d)}{d} \mathbf{X}^T \mathbf{V}_P^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}_P) \cdot (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}_P)^T \mathbf{V}_P^{-1} \mathbf{H}_j \mathbf{V}_P^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}_P) \middle| \mathbf{X} \right] \right] \\
&\quad - \mathbb{E} \left[ \mathbb{E} \left[ u(d) \cdot 2 \mathbf{X}^T \mathbf{V}_P^{-1} \mathbf{H}_j \mathbf{V}_P^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}_P) \middle| \mathbf{X} \right] \right] \\
&\quad - \mathrm{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) \mathbb{E} \left[ \mathbb{E} \left[ \frac{\rho'(d)}{d} \mathbf{X}^T \mathbf{V}_P^{-1} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}_P) \middle| \mathbf{X} \right] \right].
\end{aligned} 
\quad (B.34) $$

According to Lemma B.4, the first term on the right hand side of (B.34) has inner conditional expectation

$$ 
\begin{aligned}
\mathbf{X}^T \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} \left[ \frac{u'(\|\mathbf{z}\|)}{\|\mathbf{z}\|} \mathbf{\Sigma}^{-1/2} \mathbf{z} \mathbf{z}^T \mathbf{\Sigma}^{-1/2} \mathbf{H}_j \mathbf{\Sigma}^{-1/2} \mathbf{z} \right] \\
= \mathbf{X}^T \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} \left[ u'(\|\mathbf{z}\|) \|\mathbf{z}\|^2 \right] \mathbf{\Sigma}^{-1/2} \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} \left[ \mathbf{u} \mathbf{u}^T \mathbf{\Sigma}^{-1/2} \mathbf{H}_j \mathbf{\Sigma}^{-1/2} \mathbf{u} \right].
\end{aligned} 
$$

Again, the second term on the right hand side is the expectation with respect to a spherical density of an odd function of $\mathbf{u}$, and is therefore equal to zero. Similarly, the inner expectation of the second term on the right hand side of (B.34) is equal to

$$ 2 \mathbf{X}^T \mathbf{\Sigma}^{-1} \mathbf{H}_j \mathbf{\Sigma}^{-1/2} \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} [u(\|\mathbf{z}\|) \mathbf{z}] = 2 \mathbf{X}^T \mathbf{\Sigma}^{-1} \mathbf{H}_j \mathbf{\Sigma}^{-1/2} \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} [u(\|\mathbf{z}\|) \|\mathbf{z}\|] \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} [\mathbf{u}] = \mathbf{0}, $$

and the inner expectation of the third term on the right hand side of (B.34) is equal to

$$ \mathbf{X}^T \mathbf{\Sigma}^{-1/2} \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} \left[ \frac{\rho'(\|\mathbf{z}\|)}{\|\mathbf{z}\|} \mathbf{z} \right] = \mathbf{X}^T \mathbf{\Sigma}^{-1/2} \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} \left[ \rho'(\|\mathbf{z}\|) \right] \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} [\mathbf{u}] = \mathbf{0}. $$

It follows that

$$ \frac{\partial \Lambda_{\theta}(\boldsymbol{\xi}_P)}{\partial \beta} = \mathbf{0}. $$

Finally, to determine $\partial \Lambda_{\theta,j}(\boldsymbol{\xi}_P)/\partial \theta_s$, note that when $\mathbf{V}$ is linear, we can write

$$ \Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi}) = -\mathrm{vec} \left( \mathbf{V}^{-1} \mathbf{L}_j \mathbf{V}^{-1} \right)^T \mathrm{vec} \left( \Psi_{\mathbf{V}}(\mathbf{s}, \boldsymbol{\xi}) \right) $$

where $\Psi_{\mathbf{V}}$ is defined in (7.8) and has the property that

$$ \int \Psi_{\mathbf{V}}(\mathbf{s}, \boldsymbol{\xi}_P) dP(\mathbf{s}) = \mathbf{0}, $$

when $\mathbf{y} | \mathbf{X}$ has an elliptically contoured density $f_{\mu, \mathbf{\Sigma}}$ with parameters $\mu = \mathbf{X}\boldsymbol{\beta}_P$ and $\mathbf{\Sigma} = \mathbf{V}_P$. This means that for each $j, s = 1, \dots, l$, we have

$$ \frac{\partial \Lambda_{\theta,j}(\boldsymbol{\xi}_P)}{\partial \theta_s} = \int \frac{\partial \Psi_{\theta,j}(\mathbf{s}, \boldsymbol{\xi}_P)}{\partial \theta_s} dP(\mathbf{s}) = -\mathrm{vec}(\mathbf{\Sigma}^{-1} \mathbf{L}_j \mathbf{\Sigma}^{-1})^T \mathrm{vec} \left( \int \frac{\partial \Psi_{\mathbf{V}}(\mathbf{s}, \boldsymbol{\xi}_P)}{\partial \theta_s} dP(\mathbf{s}) \right) $$

where $\Psi_{\mathbf{V}}$ is defined in (7.8). As before, we find

$$ 
\begin{aligned}
\int \frac{\partial \Psi_{\mathbf{V}}(\mathbf{s}, \boldsymbol{\xi}_P)}{\partial \theta_s} dP(\mathbf{s}) = & -\mathbb{E} \left[ \mathbb{E} \left[ \frac{k u'(d)}{2d} \mathbf{z}^T \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \mathbf{z} \cdot \mathbf{\Sigma}^{1/2} \mathbf{z} \mathbf{z}^T \mathbf{\Sigma}^{1/2} \middle| \mathbf{X} \right] \right] \\
& + \mathbb{E} \left[ \mathbb{E} \left[ \frac{v'(d)}{2d} \mathbf{z}^T \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \mathbf{z} \cdot \mathbf{\Sigma} \middle| \mathbf{X} \right] \right] \\
& - \mathbb{E} \left[ \mathbb{E} \left[ v(d) \mathbf{L}_s \middle| \mathbf{X} \right] \right].
\end{aligned} 
\quad (B.35) $$
The first term on the right hand side of (B.35) is equal to

$$ 
\begin{aligned}
& \mathbb{E}_{0, \mathbf{I}_k} \left[ \frac{ku'(\|\mathbf{z}\|)}{2\|\mathbf{z}\|} \mathbf{z}^T \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \mathbf{z} \mathbf{\Sigma}^{1/2} \mathbf{z} \mathbf{z}^T \mathbf{\Sigma}^{1/2} \right] \\
& = \mathbb{E}_{0, \mathbf{I}_k} \left[ \frac{ku'(\|\mathbf{z}\|)\|\mathbf{z}\|^3}{2} \right] \mathbb{E}_{0, \mathbf{I}_k} \left[ \mathbf{u}^T \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \mathbf{u} \mathbf{\Sigma}^{1/2} \mathbf{u} \mathbf{u}^T \mathbf{\Sigma}^{1/2} \right].
\end{aligned}
$$

Furthermore, we can write

$$ 
\begin{aligned}
& \mathrm{vec}(\mathbf{\Sigma}^{-1} \mathbf{L}_j \mathbf{\Sigma}^{-1})^T \mathrm{vec} \left( \mathbb{E}_{0, \mathbf{I}_k} \left[ \mathbf{u}^T \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \mathbf{u} \mathbf{\Sigma}^{1/2} \mathbf{u} \mathbf{u}^T \mathbf{\Sigma}^{1/2} \right] \right) \\
& = \mathrm{vec}(\mathbf{L}_j)^T (\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1}) \mathbb{E}_{0, \mathbf{I}_k} \left[ \mathrm{vec} \left( \mathbf{\Sigma}^{1/2} \mathbf{u} \mathbf{u}^T \mathbf{\Sigma}^{1/2} \right) \mathbf{u}^T \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \mathbf{u} \right] \\
& = \mathrm{vec}(\mathbf{L}_j)^T (\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1}) (\mathbf{\Sigma}^{1/2} \otimes \mathbf{\Sigma}^{1/2}) \mathbb{E}_{0, \mathbf{I}_k} \left[ \mathrm{vec} \left( \mathbf{u} \mathbf{u}^T \right) \mathrm{vec} \left( \mathbf{u} \mathbf{u}^T \right)^T \mathrm{vec} \left( \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \right) \right] \\
& = \mathrm{vec} \left( \mathbf{\Sigma}^{-1/2} \mathbf{L}_j \mathbf{\Sigma}^{-1/2} \right)^T \frac{1}{k(k+2)} \left( \mathbf{I}_{k^2} + \mathbf{K}_{k,k} + \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \right) \mathrm{vec} \left( \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \right),
\end{aligned}
$$

using Lemma B.4. Application of property (B.32) and the fact that $\mathrm{vec}(\mathbf{A})^T \mathrm{vec}(\mathbf{B}) = \mathrm{tr}(\mathbf{A}\mathbf{B})$, yields

$$ 
\begin{aligned}
& \mathrm{vec}(\mathbf{\Sigma}^{-1} \mathbf{L}_j \mathbf{\Sigma}^{-1})^T \mathrm{vec} \left( \mathbb{E}_{0, \mathbf{I}_k} \left[ \mathbf{u}^T \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \mathbf{u} \mathbf{\Sigma}^{1/2} \mathbf{u} \mathbf{u}^T \mathbf{\Sigma}^{1/2} \right] \right) \\
& = \frac{1}{k(k+2)} \left( 2\mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_j \mathbf{\Sigma}^{-1} \mathbf{L}_s) + \mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_j)\mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_s) \right).
\end{aligned}
$$

It follows that the first term on the right hand side of (B.35) leads to a first term in $\partial\Lambda_{\theta,j}(\xi_P)/\partial\theta_s$, which is equal to

$$ \frac{\mathbb{E}_{0, \mathbf{I}_k} [u'(\|\mathbf{z}\|)\|\mathbf{z}\|^3]}{2(k+2)} \left( 2\mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_j \mathbf{\Sigma}^{-1} \mathbf{L}_s) + \mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_j)\mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_s) \right). \quad (B.36) $$

The second term on the right hand side of (B.35) is equal to

$$ 
\begin{aligned}
& \mathbb{E}_{0, \mathbf{I}_k} \left[ \frac{v'(\|\mathbf{z}\|)}{2\|\mathbf{z}\|} \mathbf{z}^T \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \mathbf{z} \mathbf{\Sigma} \right] \\
& = \mathbb{E}_{0, \mathbf{I}_k} \left[ \frac{v'(\|\mathbf{z}\|)\|\mathbf{z}\|}{2} \right] \mathbb{E}_{0, \mathbf{I}_k} \left[ \mathbf{u}^T \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \mathbf{u} \right] \mathbf{\Sigma} \\
& = \mathbb{E}_{0, \mathbf{I}_k} \left[ \frac{v'(\|\mathbf{z}\|)\|\mathbf{z}\|}{2} \right] \mathrm{vec} \left( \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \right)^T \mathrm{vec} \left( \mathbb{E}_{0, \mathbf{I}_k} [\mathbf{u} \mathbf{u}^T] \right) \mathbf{\Sigma} \\
& = \mathbb{E}_{0, \mathbf{I}_k} \left[ \frac{v'(\|\mathbf{z}\|)\|\mathbf{z}\|}{2} \right] \mathrm{vec} \left( \mathbf{\Sigma}^{-1/2} \mathbf{L}_s \mathbf{\Sigma}^{-1/2} \right)^T \mathrm{vec} \left( \frac{1}{k} \mathbf{I}_k \right) \mathbf{\Sigma} \\
& = \frac{\mathbb{E}_{0, \mathbf{I}_k} [v'(\|\mathbf{z}\|)\|\mathbf{z}\|]}{2k} \mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_s) \mathbf{\Sigma},
\end{aligned}
$$

using Lemma B.4. This leads to a second term in $\partial\Lambda_{\theta,j}(\xi_P)/\partial\theta_s$, which is equal to

$$ 
\begin{aligned}
& - \frac{\mathbb{E}_{0, \mathbf{I}_k} [v'(\|\mathbf{z}\|)\|\mathbf{z}\|]}{2k} \mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_s) \mathrm{vec}(\mathbf{\Sigma}^{-1} \mathbf{L}_j \mathbf{\Sigma}^{-1})^T \mathrm{vec}(\mathbf{\Sigma}) \\
& = - \frac{\mathbb{E}_{0, \mathbf{I}_k} [v'(\|\mathbf{z}\|)\|\mathbf{z}\|]}{2k} \mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_s) \mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_j).
\end{aligned}
\quad (B.37) $$

The third term on the right hand side of (B.35) leads to a third term in $\partial\Lambda_{\theta,j}(\xi_P)/\partial\theta_s$, which is equal to

$$ \mathbb{E}_{0, \mathbf{I}_k} [v(\|\mathbf{z}\|)] \mathrm{vec}(\mathbf{\Sigma}^{-1} \mathbf{L}_j \mathbf{\Sigma}^{-1})^T \mathrm{vec}(\mathbf{L}_s) = \mathbb{E}_{0, \mathbf{I}_k} [v(\|\mathbf{z}\|)] \mathrm{tr}(\mathbf{\Sigma}^{-1} \mathbf{L}_j \mathbf{\Sigma}^{-1} \mathbf{L}_s). \quad (B.38) $$
We conclude that $\partial\Lambda_{\theta,j}(\xi_P)/\partial\theta_s$ consists of three terms given in (B.36), (B.37) and (B.38). This means that $\partial\Lambda_{\theta,j}(\xi_P)/\partial\theta_s$ has a term $\text{tr}(\Sigma^{-1}\mathbf{L}_j\Sigma^{-1}\mathbf{L}_s)$ with coefficient

$$ \frac{\mathbb{E}_{\mathbf{0},\mathbf{I}_k} [u'(\|\mathbf{z}\|) \|\mathbf{z}\|^3]}{(k+2)} + \mathbb{E}_{\mathbf{0},\mathbf{I}_k} [v(\|\mathbf{z}\|)] = \gamma_1, $$

and a term $\text{tr}(\Sigma^{-1}\mathbf{L}_s)\text{tr}(\Sigma^{-1}\mathbf{L}_j)$ with coefficient

$$ \frac{\mathbb{E}_{\mathbf{0},\mathbf{I}_k} [u'(\|\mathbf{z}\|) \|\mathbf{z}\|^3]}{2(k+2)} - \frac{\mathbb{E}_{\mathbf{0},\mathbf{I}_k} [v'(\|\mathbf{z}\|) \|\mathbf{z}\|]}{2k} = -\gamma_2, $$

where $\gamma_1$ and $\gamma_2$ are defined in (8.5), and where we use that $u'(s)s^3 = \rho''(s)s^2 - \rho'(s)s$ and $v(s) = \rho'(s)s - \rho(s) + b_0$. Finally, from the definition of $\mathbf{L}$ in (7.9) it follows that the $l \times l$ matrix with entries

$$ \begin{aligned} & \gamma_1 \text{tr}(\Sigma^{-1}\mathbf{L}_j\Sigma^{-1}\mathbf{L}_s) - \gamma_2 \text{tr}(\Sigma^{-1}\mathbf{L}_s)\text{tr}(\Sigma^{-1}\mathbf{L}_j) \\ &= \gamma_1 \text{vec}(\mathbf{L}_j)^T (\Sigma^{-1/2} \otimes \Sigma^{-1/2}) (\Sigma^{-1/2} \otimes \Sigma^{-1/2}) \text{vec}(\mathbf{L}_s) \\ & \quad - \gamma_2 \text{vec}(\mathbf{L}_j)^T (\Sigma^{-1/2} \otimes \Sigma^{-1/2}) \text{vec}(\mathbf{I}_k) \text{vec}(\mathbf{I}_k)^T (\Sigma^{-1/2} \otimes \Sigma^{-1/2}) \text{vec}(\mathbf{L}_s) \\ &= \gamma_1 \mathbf{L}^T (\Sigma^{-1} \otimes \Sigma^{-1}) \mathbf{L} - \gamma_2 \text{vec}(\mathbf{L}_j)^T \text{vec}(\Sigma^{-1}) \text{vec}(\Sigma^{-1})^T \text{vec}(\mathbf{L}_s), \end{aligned} $$

is the matrix

$$ \gamma_1 \mathbf{L}^T (\Sigma^{-1} \otimes \Sigma^{-1}) \mathbf{L} - \gamma_2 \mathbf{L}^T \text{vec}(\Sigma^{-1}) \text{vec}(\Sigma^{-1})^T \mathbf{L}. $$

This proves the lemma. $\square$

**Lemma B.5.** Suppose that $\rho$ satisfies (R3)-(R4). Let $\gamma_1$ and $\gamma_2$ defined in (8.5) and suppose that $\gamma_1 > 0$. Then the inverse of $\partial\Lambda_{\theta}(\xi(P))/\partial\theta^T$ exists and is given by

$$ a \left( \mathbf{L}^T (\Sigma^{-1} \otimes \Sigma^{-1}) \mathbf{L} \right)^{-1} + b \left( \mathbf{L}^T (\Sigma^{-1} \otimes \Sigma^{-1}) \mathbf{L} \right)^{-1} \mathbf{L}^T \text{vec}(\Sigma^{-1}) \text{vec}(\Sigma^{-1})^T \mathbf{L} \left( \mathbf{L}^T (\Sigma^{-1} \otimes \Sigma^{-1}) \mathbf{L} \right)^{-1} $$

where $a = 1/\gamma_1$ and $b = \gamma_2/(\gamma_1(\gamma_1 - k\gamma_2))$.

*Proof.* Together with Lemma 2, first write

$$ \begin{aligned} \frac{\partial \Lambda_{\theta}(\xi(P))}{\partial \theta^T} &= \gamma_1 \mathbf{L}^T (\Sigma^{-1} \otimes \Sigma^{-1}) \mathbf{L} - \gamma_2 \mathbf{L}^T \text{vec}(\Sigma^{-1}) \text{vec}(\Sigma^{-1})^T \mathbf{L} \\ &= \gamma_1 \mathbf{E}^T \mathbf{E} - \gamma_2 \mathbf{E}^T \text{vec}(\mathbf{I}_k) \text{vec}(\mathbf{I}_k)^T \mathbf{E}, \end{aligned} $$

where $\mathbf{E} = (\Sigma^{-1/2} \otimes \Sigma^{-1/2}) \mathbf{L}$. Since $\mathbf{V}(\theta(P)) = \Sigma$, by definition of $\mathbf{L}$, it follows that for $\theta(P) = (\theta_1, \dots, \theta_l)^T$,

$$ \mathbf{L} \theta(P) = \sum_{j=1}^{l} \theta_j \text{vec}(\mathbf{L}_j) = \text{vec} \left( \sum_{j=1}^{l} \theta_j \mathbf{L}_j \right) = \text{vec}(\Sigma). $$

This means that

$$ \mathbf{E} \theta(P) = (\Sigma^{-1/2} \otimes \Sigma^{-1/2}) \mathbf{L} \theta(P) = (\Sigma^{-1/2} \otimes \Sigma^{-1/2}) \text{vec}(\Sigma) = \text{vec}(\mathbf{I}_k). \quad (B.39) $$

Since $\mathbf{L}$ has full rank, also $\mathbf{E}$ has full rank. This means that $(\mathbf{E}^T\mathbf{E})^{-1}$ exists, and satisfies

$$ (\mathbf{E}^T\mathbf{E})^{-1}\mathbf{E}^T\text{vec}(\mathbf{I}_k) = \theta(P). \quad (B.40) $$
Now, write

$$
\begin{aligned}
& \gamma_1 \mathbf{E}^T \mathbf{E} - \gamma_2 \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} \\
& = \left( \gamma_1 \mathbf{I}_k - \gamma_2 \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \right) (\mathbf{E}^T \mathbf{E}).
\end{aligned}
$$

When we multiply the first matrix from the left with a matrix of the same type,

$$
\left( a \mathbf{I}_k + b \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \right) \left( \gamma_1 \mathbf{I}_k - \gamma_2 \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \right), \quad (B.41)
$$

then we find four terms. A term $\mathbf{I}_k$ with coefficient $a\gamma_1$, two terms

$$
\mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1},
$$

with coefficient $-a\gamma_2 + b\gamma_1$, and the term

$$
\mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1}
$$

with coefficient $-b\gamma_2$. Consider the scalar valued inner product in the middle

$$
\mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) = \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} \boldsymbol{\theta}(P) = \mathrm{vec}(\mathbf{I}_k)^T \mathrm{vec}(\mathbf{I}_k) = k, \quad (B.42)
$$

by application of (B.40) and then (B.39). It follows that the term with coefficient $-b\gamma_2$ reduces to

$k \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1}.$

Hence the matrix product in (B.41) is equal to

$$
a\gamma_1 \mathbf{I}_k + (-a\gamma_2 + b\gamma_1 - kb\gamma_2) \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \quad (B.43)
$$

When we multiply the same matrix from the right,

$$
\left( \gamma_1 \mathbf{I}_k - \gamma_2 \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \right) \left( a \mathbf{I}_k + b \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \right),
$$

we find the same result (B.43). This matrix is equal to $\mathbf{I}_k$ if and only if $a\gamma_1 = 1$ and $-a\gamma_2 + b\gamma_1 - kb\gamma_2 = 0$, or equivalently

$$
\begin{align}
a &= 1/\gamma_1 \nonumber \\
b &= \frac{a\gamma_2}{\gamma_1 - k\gamma_2} = \frac{\gamma_2}{\gamma_1(\gamma_1 - k\gamma_2)}, \tag{B.44}
\end{align}
$$

where we use that $\gamma_1 > 0$ and $\gamma_1 - k\gamma_2 = \mathbb{E}_{0,\mathbf{I}_k}[\rho'(\|\mathbf{z}\|)]/2 > 0$, due to (R3)-(R4). We conclude
that the inverse of $\gamma_1\mathbf{I}_k - \gamma_2\mathbf{E}^T\mathrm{vec}(\mathbf{I}_k)\mathrm{vec}(\mathbf{I}_k)^T\mathbf{E}(\mathbf{E}^T\mathbf{E})^{-1}$ exists and is equal to

$$
a \mathbf{I}_k + b \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1}
$$

with $a$ and $b$ given in (B.44). Hence, the inverse of the matrix $\gamma_1 \mathbf{E}^T \mathbf{E} - \gamma_2 \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E}$ is equal to

$$
\begin{align*}
& (\mathbf{E}^T \mathbf{E})^{-1} \left( \gamma_1 \mathbf{I}_k - \gamma_2 \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \right)^{-1} \\
&= (\mathbf{E}^T \mathbf{E})^{-1} \left( a \mathbf{I}_k + b \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1} \right) \\
&= a (\mathbf{E}^T \mathbf{E})^{-1} + b (\mathbf{E}^T \mathbf{E})^{-1} \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E} (\mathbf{E}^T \mathbf{E})^{-1}.
\end{align*}
$$

After inserting $\mathbf{E} = (\mathbf{\Sigma}^{-1/2} \otimes \mathbf{\Sigma}^{-1/2}) \mathbf{L}$, this finishes the proof. $\square$

<figure><img src="image_4.png" /></figure>
## Proof of Corollary 5

*Proof.* Since $\rho$ is strictly increasing on $[0, c_0]$, the function $u(s) = \rho'(s)/s > 0$, for $0 < s \le c_0$ and zero for $s > c_0$. This means that

$$ \alpha = \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} \left[ \left( 1 - \frac{1}{k} \right) \frac{\rho'(\|\mathbf{z}\|)}{\|\mathbf{z}\|} + \frac{1}{k} \rho''(\|\mathbf{z}\|) \right] \ge \frac{1}{k} \mathbb{E}_{\mathbf{0}, \mathbf{I}_k} [\rho''(\|\mathbf{z}\|)] > 0. $$

Furthermore, since $\mathbf{X}$ has full rank, the inverse of $\mathbb{E}[\mathbf{X}^T \mathbf{\Sigma}^{-1} \mathbf{X}]$ exists. It follows that the matrix $\partial\Lambda_{\beta}(\xi(P))/\partial\beta$ in (8.3) is non-singular. According to Lemma B.5, also $\partial\Lambda_{\theta}(\xi(P))/\partial\theta$ is non-singular. Together, with Lemma B.3 and Lemma 2, we conclude that $\partial\Lambda/\partial\xi$ is continuously differentiable with a non-singular derivative at $\xi(P)$, so that Theorem 5 applies. Together with Lemma 2, this implies that

$$ \begin{aligned} \mathrm{IF}(\mathbf{s}_0, \beta, P) &= - \left( \frac{\partial \Lambda_{\beta}(\xi(P))}{\partial \beta} \right)^{-1} \Psi_{\beta}(\mathbf{s}_0, \xi(P)) \\ &= \frac{u(d_0)}{\alpha} \left( \mathbb{E}[\mathbf{X}^T \mathbf{\Sigma}^{-1} \mathbf{X}] \right)^{-1} \mathbf{X}_0^T \mathbf{\Sigma}^{-1} (\mathbf{y}_0 - \mathbf{X}_0 \beta) \end{aligned} $$

where $d_0^2 = (\mathbf{y}_0 - \mathbf{X}_0\beta)^T \mathbf{\Sigma}^{-1}(\mathbf{y}_0 - \mathbf{X}_0\beta)$. From Theorem 5, together with Lemma 2, it also follows that

$$ \begin{aligned} \mathrm{IF}(\mathbf{s}_0, \theta, P) &= - \left( \frac{\partial \Lambda_{\theta}(\xi(P))}{\partial \theta} \right)^{-1} \Psi_{\theta}(\mathbf{s}_0, \xi(P)) \\ &= \left( \frac{\partial \Lambda_{\theta}(\xi(P))}{\partial \theta} \right)^{-1} \mathbf{L}^T (\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1}) \times \\ & \quad \times \mathrm{vec} \left( k u(d_0) (\mathbf{y}_0 - \mathbf{X}_0 \beta)(\mathbf{y}_0 - \mathbf{X}_0 \beta)^T - v(d_0) \mathbf{\Sigma} \right) \quad (\text{B.45}) \\ &= k u(d_0) \left( \frac{\partial \Lambda_{\theta}(\xi(P))}{\partial \theta} \right)^{-1} \mathbf{L}^T (\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1}) \mathrm{vec} \left( (\mathbf{y}_0 - \mathbf{X}_0 \beta)(\mathbf{y}_0 - \mathbf{X}_0 \beta)^T \right) \\ & \quad - v(d_0) \left( \frac{\partial \Lambda_{\theta}(\xi(P))}{\partial \theta} \right)^{-1} \mathbf{L}^T (\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1}) \mathrm{vec}(\mathbf{\Sigma}). \end{aligned} $$

Consider the first term on the right hand side of (B.45). We have that

$$ \begin{aligned} & (\mathbf{\Sigma}^{-1/2} \otimes \mathbf{\Sigma}^{-1/2}) \mathrm{vec} \left( (\mathbf{y}_0 - \mathbf{X}_0\beta)(\mathbf{y}_0 - \mathbf{X}_0\beta)^T \right) \\ &= \mathrm{vec} \left( \mathbf{\Sigma}^{-1/2} (\mathbf{y}_0 - \mathbf{X}_0\beta)(\mathbf{y}_0 - \mathbf{X}_0\beta)^T \mathbf{\Sigma}^{-1/2} \right) \end{aligned} $$

and from Lemma B.5,

$$ \begin{aligned} & \left( \frac{\partial \Lambda_{\theta}(\xi(P))}{\partial \theta} \right)^{-1} \mathbf{L}^T (\mathbf{\Sigma}^{-1/2} \otimes \mathbf{\Sigma}^{-1/2}) \\ &= a(\mathbf{E}^T \mathbf{E})^{-1} \mathbf{E}^T + b(\mathbf{E}^T \mathbf{E})^{-1} \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E}(\mathbf{E}^T \mathbf{E})^{-1} \mathbf{E}^T, \end{aligned} \quad (\text{B.46}) $$

where $\mathbf{E} = (\mathbf{\Sigma}^{-1/2} \otimes \mathbf{\Sigma}^{-1/2})\mathbf{L}$. This implies

$$ \begin{aligned} & \left( \frac{\partial \Lambda_{\theta}(\xi(P))}{\partial \theta} \right)^{-1} \mathbf{L}^T (\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1}) \mathrm{vec} \left( (\mathbf{y}_0 - \mathbf{X}_0\beta)(\mathbf{y}_0 - \mathbf{X}_0\beta)^T \right) \\ &= a(\mathbf{E}^T \mathbf{E})^{-1} \mathbf{E}^T \mathrm{vec} \left( \mathbf{\Sigma}^{-1/2}(\mathbf{y}_0 - \mathbf{X}_0\beta)(\mathbf{y}_0 - \mathbf{X}_0\beta)^T \mathbf{\Sigma}^{-1/2} \right) \\ & \quad + b(\mathbf{E}^T \mathbf{E})^{-1} \mathbf{E}^T \mathrm{vec}(\mathbf{I}_k) \mathrm{vec}(\mathbf{I}_k)^T \mathbf{E}(\mathbf{E}^T \mathbf{E})^{-1} \mathbf{E}^T \times \\ & \quad \times \mathrm{vec} \left( \mathbf{\Sigma}^{-1/2}(\mathbf{y}_0 - \mathbf{X}_0\beta)(\mathbf{y}_0 - \mathbf{X}_0\beta)^T \mathbf{\Sigma}^{-1/2} \right). \end{aligned} $$
The first term on the right hand side is equal to

$$a(\mathbf{L}^T(\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1})\mathbf{L})^{-1} \mathbf{L}^T(\mathbf{\Sigma}^{-1/2} \otimes \mathbf{\Sigma}^{-1/2}) \times \\ \times \mathrm{vec}(\mathbf{\Sigma}^{-1/2}(\mathbf{y}_0 - \mathbf{X}_0\beta)(\mathbf{y}_0 - \mathbf{X}_0\beta)^T \mathbf{\Sigma}^{-1/2})$$

and, with (B.39) and (B.40), the second term on the right hand side is equal to

$$\begin{align*}
b\boldsymbol{\theta}(P)\boldsymbol{\theta}(P)^T \mathbf{E}^T \mathrm{vec}(\mathbf{\Sigma}^{-1/2}(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})^T \mathbf{\Sigma}^{-1/2}) \\
&= b\boldsymbol{\theta}(P)\mathrm{vec}(\mathbf{I}_k)^T \mathrm{vec}(\mathbf{\Sigma}^{-1/2}(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})^T \mathbf{\Sigma}^{-1/2}) \\
&= b\boldsymbol{\theta}(P)\mathrm{tr}(\mathbf{\Sigma}^{-1/2}(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})^T \mathbf{\Sigma}^{-1/2}) \\
&= bd_0^2\boldsymbol{\theta}(P).
\end{align*}$$

It follows that the first term on the right hand side of (B.45) is equal to

$$aku(d_0)(\mathbf{L}^T(\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1})\mathbf{L})^{-1} \mathbf{L}^T(\mathbf{\Sigma}^{-1/2} \otimes \mathbf{\Sigma}^{-1/2}) \times \\ \times \mathrm{vec}(\mathbf{\Sigma}^{-1/2}(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})^T \mathbf{\Sigma}^{-1/2}) \\ + bku(d_0)d_0^2\boldsymbol{\theta}(P).$$

Next consider the second term on the right hand side of (B.45). We have that

$$(\mathbf{\Sigma}^{-1/2} \otimes \mathbf{\Sigma}^{-1/2})\mathrm{vec}(\mathbf{\Sigma}) = \mathrm{vec}(\mathbf{I}_k),$$

and with (B.46), together with (B.39) and (B.40),

$$\begin{align*}
\left(\frac{\partial \Lambda_{\boldsymbol{\theta}}(\boldsymbol{\xi}(P))}{\partial \boldsymbol{\theta}}\right)^{-1} \mathbf{L}^T(\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1})\mathrm{vec}(\mathbf{\Sigma}) \\
&= a(\mathbf{E}^T\mathbf{E})^{-1}\mathbf{E}^T\mathrm{vec}(\mathbf{I}_k) + b(\mathbf{E}^T\mathbf{E})^{-1}\mathbf{E}^T\mathrm{vec}(\mathbf{I}_k)\mathrm{vec}(\mathbf{I}_k)^T\mathbf{E}(\mathbf{E}^T\mathbf{E})^{-1}\mathbf{E}^T\mathrm{vec}(\mathbf{I}_k) \\
&= a\boldsymbol{\theta}(P) + b\boldsymbol{\theta}(P)\mathrm{vec}(\mathbf{I}_k)^T\mathbf{E}(\mathbf{E}^T\mathbf{E})^{-1}\mathbf{E}^T\mathrm{vec}(\mathbf{I}_k) \\
&= a\boldsymbol{\theta}(P) + b\boldsymbol{\theta}(P)\mathrm{vec}(\mathbf{I}_k)^T\mathbf{E}\boldsymbol{\theta}(P) \\
&= a\boldsymbol{\theta}(P) + b\boldsymbol{\theta}(P)\mathrm{vec}(\mathbf{I}_k)^T\mathrm{vec}(\mathbf{I}_k) \\
&= (a + bk)\boldsymbol{\theta}(P).
\end{align*}$$

It follows that the second term on the right hand side of (B.45) is equal to

$$-v(d_0)(a + bk)\boldsymbol{\theta}(P).$$

Putting things together, we find that $\mathrm{IF}(\mathbf{s}_0, \boldsymbol{\theta}, P)$ is equal to

$$aku(d_0)(\mathbf{L}^T(\mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1})\mathbf{L})^{-1} \mathbf{L}^T(\mathbf{\Sigma}^{-1/2} \otimes \mathbf{\Sigma}^{-1/2}) \times \\ \times \mathrm{vec}(\mathbf{\Sigma}^{-1/2}(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})^T \mathbf{\Sigma}^{-1/2}) \\ + (bku(d_0)d_0^2 - av(d_0) - bkv(d_0))\boldsymbol{\theta}(P).$$

Since $v(d_0) = u(d_0)d_0^2 - \rho(d_0) + b_0$, we have that

$$\begin{align*}
bku(d_0)d_0^2 - av(d_0) - bkv(d_0) &= -\left(\frac{v(d_0)}{\gamma_1} - \frac{k\gamma_2}{\gamma_1(\gamma_1 - k\gamma_2)}(\rho(d_0) - b_0)\right) \\
&= -\frac{u(d_0)d_0^2}{\gamma_1} + \left(\frac{1}{\gamma_1} + \frac{k\gamma_2}{\gamma_1(\gamma_1 - k\gamma_2)}\right)(\rho(d_0) - b_0) \\
&= -\frac{u(d_0)d_0^2}{\gamma_1} + \frac{\rho(d_0) - b_0}{\gamma_1 - k\gamma_2}
\end{align*}$$
We conclude that $\mathrm{IF}(\mathbf{s}_0, \boldsymbol{\theta}, P)$ is given by

$$ \frac{ku(d_0)}{\gamma_1} (\mathbf{L}^T(\boldsymbol{\Sigma}^{-1} \otimes \boldsymbol{\Sigma}^{-1})\mathbf{L})^{-1} \mathbf{L}^T \mathrm{vec}(\boldsymbol{\Sigma}^{-1}(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})(\mathbf{y}_0 - \mathbf{X}_0\boldsymbol{\beta})^T \boldsymbol{\Sigma}^{-1}) \\ + \left(-\frac{u(d_0)d_0^2}{\gamma_1} + \frac{\rho(d_0) - b_0}{\gamma_1 - k\gamma_2}\right) \boldsymbol{\theta}(P). $$

This proves the corollary.

□

## B.5 Proofs of Section 9

As preparation we first establish the following lemma, which is similar to Lemma 22 in [21].

**Lemma B.6.** Let $\rho(\cdot)$ be a real-valued function of bounded variation on $\mathbb{R}^+$. The class of all functions on $\mathbb{R}^p$ of the form

$$ \mathbf{s} = (\mathbf{y}, \mathbf{X}) \mapsto \rho(\|\mathbf{A}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\|) $$

with $\mathbf{A}$ ranging over all $k \times k$ matrices and $\boldsymbol{\beta}$ ranging over $\mathbb{R}^q$, has polynomial discrimination.

*Proof.* Consider the class of functions

$$ g_{\mathbf{A},\boldsymbol{\beta}}(\mathbf{s}, t) = \|\mathbf{A}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\|^2 - \rho^{-1}(t)^2. $$

According to Lemma 18 in [21], it suffices to show that the functions $g_{\mathbf{A},\boldsymbol{\beta}}(\cdot, \cdot)$ span a finite-dimensional vector space. In order to do so, write

$$ \|\mathbf{A}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\|^2 = \sum_{i=1}^{k} \sum_{j=1}^{k} (\mathbf{A}^T \mathbf{A})_{ij} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})_i (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})_j \\ = \sum_{i=1}^{k} \sum_{j=1}^{k} \sum_{s=1}^{k} a_{is} a_{sj} \left( y_i - \sum_{r=1}^{q} x_{ir} \beta_r \right) \left( y_j - \sum_{w=1}^{q} x_{jw} \beta_w \right) $$

This is a polynomial in $\mathbf{s} = (y_1, \dots, y_k, x_{11}, \dots, x_{kq})$, with coefficients in $\mathbb{R}$. This means that the class of functions $g_{\mathbf{A},\boldsymbol{\beta}}(\mathbf{s}, t) = \|\mathbf{A}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\|^2 - \rho^{-1}(t)^2$ forms a finite dimensional vector space. □

A useful first step is the following lemma.

**Lemma B.7.** Let $u: \mathbb{R} \to \mathbb{R}$ be a function of bounded variation. Let $\mathbf{s} = (\mathbf{y}, \mathbf{X}) = (s_1, \dots, s_p) \in \mathbb{R}^p$, and define

$$ g(\mathbf{s}, \boldsymbol{\beta}, \mathbf{V}) = u(\|\mathbf{V}^{-1/2}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\|) \quad \text{for } \boldsymbol{\beta} \in \mathbb{R}^q, \mathbf{V} \in \mathrm{PDS}(k). $$

Consider the classes of functions

$$ \mathcal{F} = \{g(\mathbf{s}, \boldsymbol{\beta}, \mathbf{V}) : \boldsymbol{\beta} \in \mathbb{R}^q, \mathbf{V} \in \mathrm{PDS}(k)\}, $$

$$ \mathcal{F}_a = \{g(\mathbf{s}, \boldsymbol{\beta}, \mathbf{V}) s_a : g \in \mathcal{F}\}, $$

$$ \mathcal{F}_{ab} = \{g(\mathbf{s}, \boldsymbol{\beta}, \mathbf{V}) s_a s_b : g \in \mathcal{F}\}, $$

for $a, b = 1, \dots, p$. Denote by $\mathcal{G}$, $\mathcal{G}_a$, and $\mathcal{G}_{ab}$, the corresponding classes of graphs of the functions in $\mathcal{F}$, $\mathcal{F}_a$, and $\mathcal{F}_{ab}$, respectively. Then $\mathcal{G}$, $\mathcal{G}_a$, and $\mathcal{G}_{ab}$, all have polynomial discrimination for $a, b = 1, \dots, p$.

*Proof.* Because the function $u(\cdot)$ is of bounded variation, it follows from Lemma B.6 that the class $\mathcal{G}$ has polynomial discrimination. To show the same for the class $\mathcal{G}_a$, suppose that $\mathcal{G}_a$ is not of polynomial discrimination. This means that for every integer $N$ there exists a set $V \subset \mathbb{R}^{p+1}$ of $N$ points, such that all subsets of $V$ can be written as $D_a \cap V$ for some $D_a \in \mathcal{G}_a$. Let $V = \{\mathbf{v}_1, \dots, \mathbf{v}_N\}$, where $\mathbf{v}_i = (\mathbf{s}_i, t_i)$, for $i = 1, \dots, N$. Then for every $(\mathbf{s}, t) \in V$ with $\mathbf{s} = (s_1, \dots, s_a, \dots, s_p)$, it must
hold that $s_a \neq 0$, otherwise this point cannot be separated from the other points by an element of the class $G_a$. Define the set $V_a = \{(s, t/s_a) : (s, t) \in V\}$. Note that for $D_a \in G_a$, we have

$$ (s, t) \in D_a \quad \Leftrightarrow \quad \begin{aligned} & 0 \le t \le g(s, \beta, \mathbf{V})s_a \text{ or } g(s, \beta, \mathbf{V})s_a \le t \le 0 \\ & 0 \le \frac{t}{s_a} \le g(s, \beta, \mathbf{V}) \text{ or } g(s, \beta, \mathbf{V}) \le \frac{t}{s_a} \le 0 \\ & (s, t/s_a) \in D, \end{aligned} $$

where $D \in G$. This implies that every subset of $V_a$ can be written as $D \cap V_a$, for some $D \in G$. However, this is in contradiction with the fact that $G$ has polynomial discrimination. We conclude that also $G_a$ has polynomial discrimination. A similar argument yields that $G_{ab}$ has polynomial discrimination. □

Lemma B.7 is comparable to Lemma 3 in [18], for similar classes of functions built from functions $g(\mathbf{y}, \mathbf{t}, \mathbf{C}) = u(\|\mathbf{C}^{-1/2}(\mathbf{y} - \mathbf{t})\|)$, where $\mathbf{t} \in \mathbb{R}^k$ and $\mathbf{C} \in \text{PDS}(k)$. With Lemma B.7 we can establish suitable bounds on the third term of (9.2). This is provided by the following key lemma. Once having established Lemma B.8, asymptotic normality can be derived easily from (9.2).

**Lemma B.8.** Let $\Psi = (\Psi_{\beta}, \Psi_{\theta})$ be defined in (7.5) and let $\xi_n = \xi(\mathbb{P}_n)$ and $\xi_P = \xi(P)$ be the solutions to minimization problems (3.1) and (3.5). Suppose that $\rho$ satisfies (R1)-(R4), such that $u(s)$ is of bounded variation, and suppose that $\mathbf{V}$ satisfies (V4). Suppose that $\xi_n \to \xi_P$, in probability, and that $\mathbb{E}\|s\|^2 < \infty$. Then

$$ \int (\Psi(\mathbf{s}, \xi_n) - \Psi(\mathbf{s}, \xi_P)) d(\mathbb{P}_n - P)(\mathbf{s}) = o_P(1/\sqrt{n}). \quad (B.47) $$

*Proof.* First write $\Psi_{\theta,j}(\mathbf{s}, \xi) = \Psi_{2,j}(\mathbf{s}, \xi) - \Psi_{3,j}(\mathbf{s}, \xi)$, for $j = 1, \dots, l$, where

$$ \begin{aligned} \Psi_{2,j}(\mathbf{s}, \xi) &= u(d)(\mathbf{y} - \mathbf{X}\beta)^T \mathbf{V}^{-1} \mathbf{H}_j \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X}\beta) \\ \Psi_{3,j}(\mathbf{s}, \xi) &= \text{tr} \left( \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \theta_j} \right) (\rho(d) - b_0), \end{aligned} \quad (B.48) $$

where $\mathbf{H}_j$ and $d = d(\mathbf{s}, \xi)$ are defined in (7.3) and (5.1). It suffices to show that

$$ \int (\Psi_{\beta}(\mathbf{s}, \xi_n) - \Psi_{\beta}(\mathbf{s}, \xi_P)) d(\mathbb{P}_n - P)(\mathbf{s}) = o_P(1/\sqrt{n}), \quad (B.49) $$

$$ \int (\Psi_{2,j}(\mathbf{s}, \xi_n) - \Psi_{2,j}(\mathbf{s}, \xi_P)) d(\mathbb{P}_n - P)(\mathbf{s}) = o_P(1/\sqrt{n}), \quad (B.50) $$

$$ \int (\Psi_{3,j}(\mathbf{s}, \xi_n) - \Psi_{3,j}(\mathbf{s}, \xi_P)) d(\mathbb{P}_n - P)(\mathbf{s}) = o_P(1/\sqrt{n}), \quad (B.51) $$

for $j = 1, \dots, l$.

To obtain (B.51), first write $\mathbf{M}_n = \mathbf{V}(\theta_n)^{-1}$ and $\mathbf{M}_P = \mathbf{V}(\theta_P)^{-1}$, so that $\mathbf{M}_n \to \mathbf{M}_P$, in probability, according to condition (V1). Decompose as follows

$$ \begin{aligned} \Psi_{3,j}(\mathbf{s}, \xi_n) - \Psi_{3,j}(\mathbf{s}, \xi_P) &= \text{tr} \left( \mathbf{M}_n^{-1} \frac{\partial \mathbf{V}(\theta_n)}{\partial \theta_j} \right) (\rho(d(\mathbf{s}, \xi_n)) - \rho(d(\mathbf{s}, \xi_P))) \\ &\quad + \left\{ \text{tr} \left( \mathbf{M}_n^{-1} \frac{\partial \mathbf{V}(\theta_n)}{\partial \theta_j} \right) - \text{tr} \left( \mathbf{M}_P^{-1} \frac{\partial \mathbf{V}(\theta_P)}{\partial \theta_j} \right) \right\} (\rho(d(\mathbf{s}, \xi_P)) - b_0). \end{aligned} \quad (B.52) $$

After integration with respect to $\mathbb{P}_n - P$, the first term on the right hand side of (B.52) becomes

$$ \text{tr} \left( \mathbf{M}_n^{-1} \frac{\partial \mathbf{V}(\theta_n)}{\partial \theta_j} \right) \int (\rho(d(\mathbf{s}, \xi_n)) - \rho(d(\mathbf{s}, \xi_P))) d(\mathbb{P}_n - P)(\mathbf{s}), \quad (B.53) $$
where with (V4),

$$ \mathrm{tr} \left( \mathbf{M}_n^{-1} \frac{\partial \mathbf{V}(\boldsymbol{\theta}_n)}{\partial \theta_j} \right) \rightarrow \mathrm{tr} \left( \mathbf{M}_P^{-1} \frac{\partial \mathbf{V}(\boldsymbol{\theta}_P)}{\partial \theta_j} \right) \qquad (\text{B.54}) $$

in probability. Furthermore, note that all functions $\rho(d(\cdot, \boldsymbol{\xi}))$, for $\boldsymbol{\xi} \in \mathbb{R}^q \times \mathbb{R}^l$ are members of the class

$$ \mathcal{F} = \left\{ \rho(\|\mathbf{V}^{-1/2}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\|) : \boldsymbol{\beta} \in \mathbb{R}^q, \mathbf{V} \in \mathrm{PDS}(k) \right\}. $$

From (R1)-(R2) it follows that the function $\rho$ is of bounded variation (being the sum of two monotone functions). According to Lemma B.7, the class $\mathcal{G}$, consisting of subgraphs of functions in the class $\mathcal{F}$, has polynomial discrimination. Moreover, the class $\mathcal{F}$ has a constant envelope. Then as a result of the empirical process theory developed in Pollard [23] (e.g., see Theorem 1 in [18]), it follows that for every $\delta > 0$,

$$ \sup_{f_1, f_2 \in [\delta]} \sqrt{n} \left| \int (f_1(\mathbf{s}) - f_2(\mathbf{s})) \, d(\mathbb{P}_n - P)(\mathbf{s}) \right| \rightarrow 0 $$

in probability, where

$$ [\delta] = \left\{ (f_1, f_2) : f_1, f_2 \in \mathcal{F} \text{ and } \int (f_1 - f_2)^2 dP \le \delta^2 \right\}. $$

Because $\boldsymbol{\xi}_n \rightarrow \boldsymbol{\xi}_P$ in probability one, the pair of functions $f_1(\mathbf{s}) = \rho(d(\mathbf{s}, \boldsymbol{\xi}_n))$ and $f_2(\mathbf{s}) = \rho(d(\mathbf{s}, \boldsymbol{\xi}_P)$ are in the set $[\delta]$, for sufficiently large $n$, with probability tending to one. It follows that

$$ \begin{aligned} & \sqrt{n} \left| \int \left( \rho(d(\mathbf{s}, \boldsymbol{\xi}_n) - \rho(d(\mathbf{s}, \boldsymbol{\xi}_P)) \right) d(\mathbb{P}_n - P)(\mathbf{s}) \right| \\ & \le \sup_{f_1, f_2 \in [\delta]} \sqrt{n} \left| \int (f_1(\mathbf{s}) - f_2(\mathbf{s})) d(\mathbb{P}_n - P)(\mathbf{s}) \right| \rightarrow 0, \end{aligned} \qquad (\text{B.55}) $$

in probability. Together with the fact that tr($\mathbf{M}_P^{-1} \partial \mathbf{V}(\boldsymbol{\theta}_P)/\partial \theta_j$) is bounded, this proves that (B.53) is of the order $o_P(1/\sqrt{n})$. For the second term on the right hand side of (B.52), we have that according to the central limit theorem

$$ \int (\rho(d(\mathbf{s}, \boldsymbol{\xi}_P)) - b_0) d(\mathbb{P}_n - P)(\mathbf{s}) = O_P(1/\sqrt{n}). $$

Together with (B.54) this implies that, after integration with respect to $\mathbb{P}_n - P$, the second term on the right hand side of (B.52) is of the order $o_P(1/\sqrt{n})$. This proves (B.51).

To obtain (B.49), decompose as follows

$$ \begin{aligned} & \Psi_{\boldsymbol{\beta}}(\mathbf{s}, \boldsymbol{\xi}_n) - \Psi_{\boldsymbol{\beta}}(\mathbf{s}, \boldsymbol{\xi}_P) \\ &= \left\{ u(d(\mathbf{s}, \boldsymbol{\xi}_n)) \mathbf{X}^T \mathbf{M}_n \mathbf{y} - u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \mathbf{X}^T \mathbf{M}_P \mathbf{y} \right\} \\ &\quad - \left\{ u(d(\mathbf{s}, \boldsymbol{\xi}_n)) \mathbf{X}^T \mathbf{M}_n \mathbf{X} \boldsymbol{\beta}_n - u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \mathbf{X}^T \mathbf{M}_P \mathbf{X} \boldsymbol{\beta}_P \right\}. \end{aligned} \qquad (\text{B.56}) $$

We will treat both terms on the right hand side separately. For the first term on the right hand side of (B.56), we write

$$ \left\{ u(d(\mathbf{s}, \boldsymbol{\xi}_n)) - u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \right\} \mathbf{X}^T \mathbf{M}_n \mathbf{y} + u(d(\mathbf{s}, \boldsymbol{\xi}_P) \mathbf{X}^T (\mathbf{M}_n - \mathbf{M}_P) \mathbf{y}. \qquad (\text{B.57}) $$

First consider the first term in (B.57). We consider each single element of the vector $\mathbf{X}^T \mathbf{M}_n \mathbf{y} \in \mathbb{R}^q$ separately. For $i = 1, \dots, q$ fixed, write

$$ (\mathbf{X}^T \mathbf{M}_n \mathbf{y})_i = \sum_{j=1}^{k} x_{ji} (\mathbf{M}_n \mathbf{y})_j = \sum_{j=1}^{k} x_{ji} \sum_{s=1}^{k} m_{n,js} y_s = \sum_{j=1}^{k} \sum_{s=1}^{k} x_{ji} y_s m_{n,js}. $$
Hence, after integration with respect to $\mathbb{P}_n - P$, the i-th component of the first term in (B.57) can be written as a finite sum with summands

$$ \int \left( u(d(\mathbf{s}, \boldsymbol{\xi}_n)) - u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \right) x_{ji} y_s d(\mathbb{P}_n - P)(\mathbf{s}) m_{n,js}, $$

for $i = 1, \dots, q$ and $j, s \in \{1, \dots, k\}$ fixed, where $m_{n,js} \to m_{P,js}$, Because $u$ is of bounded variation and since $\mathbf{s} = (\mathbf{y}, \mathbf{X}) = (y_1, \dots, y_k, x_{11}, \dots, x_{kq})$, all functions $u(d(\mathbf{s}, \boldsymbol{\xi}))x_{ji}y_s$, for $\boldsymbol{\xi} \in \mathbb{R}^q \times \mathbb{R}^l$, are members of the class

$$ \mathcal{F}_{ab} = \{ u(\|\mathbf{V}^{-1/2}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\|)s_a s_b : \boldsymbol{\beta} \in \mathbb{R}^q, \mathbf{V} \in \text{PDS}(k) \}. $$

According to Lemma B.7, the corresponding class of subgraphs has polynomial discrimination, so similar to (B.55) it follows that for each $i = 1, \dots, q$ and $j, s \in \{1, \dots, k\}$ fixed,

$$ \int \left( u(d(\mathbf{s}, \boldsymbol{\xi}_n)) - u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \right) x_{ji} y_s d(\mathbb{P}_n - P)(\mathbf{s}) = o_P(n^{-1/2}), $$

which means that

$$ \int \left\{ u(d(\mathbf{s}, \boldsymbol{\xi}_n)) - u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \right\} \mathbf{X}^T \mathbf{M}_P \mathbf{y} d(\mathbb{P}_n - P)(\mathbf{s}) = o_P(n^{-1/2}). \quad (B.58) $$

Next, consider the second term on the right hand side of (B.57). Similar to the first term, after integration with respect to $\mathbb{P}_n - P$, the i-th component of the second term on the right hand side of (B.57) can be written as a finite sum of summands

$$ \int u(d(\mathbf{s}, \boldsymbol{\xi}_P)) x_{ji} y_s d(\mathbb{P}_n - P)(\mathbf{s}) (m_{n,js} - m_{P,js}), $$

for $i = 1, \dots, q$, and $j, s \in \{1, \dots, k\}$ fixed. According to the central limit theorem, the integral is of the order $O_P(n^{-1/2})$, and since $m_{n,js} \to m_{P,js}$, in probability, it follows that the product is of the order $o_P(n^{-1/2})$. We conclude that

$$ \int u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \mathbf{X}^T (\mathbf{M}_n - \mathbf{M}_P) \mathbf{y} d(\mathbb{P}_n - P)(\mathbf{s}) = o_P(n^{-1/2}). \quad (B.59) $$

Putting together (B.58) and (B.59), it follows for the first term on the right hand side of (B.56) that

$$ \int \left\{ u(d(\mathbf{s}, \boldsymbol{\xi}_n)) \mathbf{X}^T \mathbf{M}_n \mathbf{y} - u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \mathbf{X}^T \mathbf{M}_P \mathbf{y} \right\} d(\mathbb{P}_n - P)(\mathbf{s}) = o_P(n^{-1/2}). \quad (B.60) $$

For the second term on the right hand side of (B.56), we write

$$ \begin{aligned} & u(d(\mathbf{s}, \boldsymbol{\xi}_n)) \mathbf{X}^T \mathbf{M}_n \mathbf{X} \boldsymbol{\beta}_n - u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \mathbf{X}^T \mathbf{M}_P \mathbf{X} \boldsymbol{\beta}_P \\ & = \left( u(d(\mathbf{s}, \boldsymbol{\xi}_n)) - u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \right) \mathbf{X}^T \mathbf{M}_n \mathbf{X} \boldsymbol{\beta}_n \\ & \quad + u(d(\mathbf{s}, \boldsymbol{\xi}_P)) \left( \mathbf{X}^T \mathbf{M}_n \mathbf{X} \boldsymbol{\beta}_n - \mathbf{X}^T \mathbf{M}_P \mathbf{X} \boldsymbol{\beta}_P \right). \end{aligned} \quad (B.61) $$

Consider the first term on the right hand side of (B.61). For $i = 1, \dots, q$ fixed, write

$$ \begin{aligned} (\mathbf{X}^T \mathbf{M}_n \mathbf{X} \boldsymbol{\beta}_n)_i &= \sum_{j=1}^k x_{ji} (\mathbf{M}_n \mathbf{X} \boldsymbol{\beta}_n)_j = \sum_{j=1}^k x_{ji} \sum_{s=1}^k m_{n,js} (\mathbf{X} \boldsymbol{\beta}_n)_s \\ &= \sum_{j=1}^k x_{ji} \sum_{s=1}^k m_{n,js} \sum_{t=1}^k x_{st} \boldsymbol{\beta}_{n,t} = \sum_{j=1}^k \sum_{s=1}^k \sum_{t=1}^k x_{ji} x_{st} m_{n,js} \boldsymbol{\beta}_{n,t}. \end{aligned} $$