Weighted likelihood estimation of multivariate location and scatter

TL;DR

This paper introduces a new weighted likelihood estimation method for multivariate location and scatter that uses Mahalanobis distances for weighting, improving robustness and avoiding high-dimensional density estimation issues.

Contribution

It proposes a Mahalanobis distance-based weighting scheme for weighted likelihood estimation, enhancing robustness and computational efficiency in multivariate analysis.

Findings

Effective outlier detection rules developed

Robust dimensionality reduction techniques introduced

Method demonstrated through numerical and real data studies

Abstract

A novel approach to obtain weighted likelihood estimates of multivariate location and scatter is discussed. A weighting scheme is proposed that is based on the distribution of the Mahalanobis distances rather than the distribution of the data at the assumed model. This strategy allows to avoid the curse of dimensionality affecting non-parametric density estimation, that is involved in the construction of the weights through the Pearson residuals Markatou et al (1998). Then, weighted likelihood based outlier detection rules and robust dimensionality reduction techniques are developed. The effectiveness of the methodology is illustrated through some numerical studies and real data examples.

Tables

Table 1. Table 1: Diabetes data: misclassification rates for different rules based on all the data and on leave-one-out cross validation.

		ALL	CV
	MLE	0.131	0.131
	WLEA	0.124	0.131
LDA	WLEB	0.076	0.103
	MCDA	0.124	0.131
	MCDB	0.083	0.117
	MLE	0.076	0.110
QDA	WLE	0.076	0.103
	MCD	0.083	0.103

Equations

$$d(y;\mu ,\Sigma )=\sqrt{(y-\mu {)}^{{}^{\top }}{\Sigma }^{-1}(y-\mu )}.$$

$$\sum _{i=1}^n{w}_{i}s(y_i;\theta )=0,$$

$$w=w(y;\theta ,{\hat{M}}_n)=\frac{{\left[A(\delta (y;\theta ,{\hat{M}}_n))+1\right]}^{+}}{\delta (y;\theta ,{\hat{M}}_n)+1},$$

$$\delta (y;\theta ,{\hat{M}}_n)=\frac{{\hat{m}}_n(y)}{m^{\ast }(y;\theta )}-1$$

$${\hat{m}}_n(y)=\frac{1}{n}\sum _{i=1}^{n}k(y;y_i,h)$$

$$m^{\ast }(y;\theta )={\int }_{\mathcal{Y}}k(y;t,h)m(t;\theta )dt$$

$$A_{\tau}(\delta)=\left\{\begin{array}[]{lc}\tau\left((\delta+1)^{1/\tau}-1\right)&\tau<\infty\\ \log(\delta+1)&\tau\rightarrow\infty\end{array}\right.$$

$$\underset{y}{\sup }|w(y,\hat{\theta },{\hat{M}}_n)-1|\stackrel{p}{\to }0$$

$$d^2(Y,\mu ,\Sigma )\sim {\chi }_p^2$$

$${\hat{\Sigma }}_u=\frac{{\sum }_{i=1}^n(y_i-\hat{\mu })(y_i-\hat{\mu }{)}^{{}^{\top }}{\hat{w}}_i}{\gamma {\sum }_{i=1}^n{\hat{w}}_i},\gamma =1-\frac{{\sum }_{i=1}^n{\hat{w}}_i^2}{({\sum }_{i=1}^n{\hat{w}}_i{)}^2},$$

$${\hat{d}}_i^2=d^2(y_i,\hat{\mu },{\hat{\Sigma }}_u)$$

$${\hat{m}}_n(d_i^2)$$

$$\delta ({\hat{d}}_i^2;\hat{\mu },{\hat{\Sigma }}_u,{\hat{M}}_n)=\frac{{\hat{m}}_n(d_i^2)}{m^{\ast }(d_i^2)}-1$$

$${\hat{w}}_i=\frac{{\left[A(\delta ({\hat{d}}_i^2,\hat{\mu },{\hat{\Sigma }}_u,{\hat{M}}_n))+1\right]}^{+}}{\delta (d_i^2,\hat{\mu },{\hat{\Sigma }}_u,{\hat{M}}_n)+1}$$

$$k(y;t,h)=\frac{1}{h}k\left(\frac{y-t}{h}\right)+\frac{1}{h}k\left(\frac{y+t}{h}\right)$$

$$k(y;t,h)=\Gamma (y;t/h+1,h)$$

$$p(x;p)=\frac{1}{2^{p/2}\Gamma (p/2)}\text{exp}\left[\frac{1}{2}(px-\text{exp}(x))\right],x\in \mathbb{R}.$$

$$d^2(Y_i,\hat{\mu },{\hat{\Sigma }}_u)\stackrel{d}{\to }\frac{(n-1{)}^2}{n}\text{Beta}\left(\frac{p}{2},\frac{n-p}{2}\right).$$

$${\text{Pr}}_{\hat{\theta }}\left[\delta (Y;\hat{\theta },{\hat{M}}_n)<-0.95\right]$$

$$\log {\hat{\pi }}_j-\frac{1}{2}{d}^2(y,{\hat{\mu }}_j,{\hat{\Sigma }}_p)$$

$$\log {\hat{\pi }}_j-\frac{1}{2}{d}^2(y,{\hat{\mu }}_j,{\hat{\Sigma }}_j)$$

$${\hat{\Sigma }}_p^A=\frac{{\sum }_{j=1}^k{\gamma }_j{\omega }_j{\hat{\Sigma }}_{uj}}{{\sum }_{j=1}^k{\gamma }_j{\omega }_j},{\omega }_j=\sum _{i=1}^{n_j}{\hat{w}}_{ij}$$

$$A-A_{(-1)}=\frac{\omega }{{\omega }_{(-i)}}(Y_i-\hat{\mu })(Y_i-\hat{\mu }{)}^{{}^{\top }}{w}_i.$$

$$A$$

$$A_{(-i)}$$

$$A-A_{(-i)}=Y_i{Y}_i^{{}^{\top }}{w}_i-\omega \hat{\mu }{\hat{\mu }}^{{}^{\top }}+{\omega }_{(-i)}{\hat{\mu }}_{(-i)}{\hat{\mu }}_{(-i)}^{{}^{\top }},$$

$$\omega \hat{\mu }={\omega }_{(-i)}{\hat{\mu }}_{(-i)}+Y_i{w}_i$$

$$R_{(-i)}=\frac{|A_{(-i)}|}{|A|}.$$

$$R_{(-i)}$$

$$(Y_i-\hat{\mu })\stackrel{d}{\to }{N}_p\left(0,\frac{n-1}{n}\Sigma \right)$$

$$\sqrt{\frac{\omega w_i}{{\omega }_{(-i)}^2{\gamma }_{(-i)}}}\stackrel{p}{\to }\sqrt{\frac{n}{(n-1)(n-2)}}.$$

Full text

**Keywords: Dimensionality Reduction; Discriminant Analysis; Mahalanobis distance; Multivariate Normal; Outlier detection; Pearson residuals; PCA; Robustness; Weighted Likelihood

Mathematics Subject Classification (2000): MSC 62F35; MSC 62G35; MSC 62H25; MSC 62H30**

1 Introduction

Several multivariate techniques are based on the the assumption of multivariate normality and the use of the sample mean vector and covariance matrix. It is well known that small departures from model assumptions may invalidate classical analysis completely [Maronna et al., 2006, Huber and Ronchetti, 2009]. Such departures result in data inadequacies that are typically observed in the form of several outliers. Outliers can be defined as observations that are highly unlikely to occur under the assumed model [Markatou et al., 1998]. In other words, outliers contaminate the data with respect to (w.r.t.) the postulated model. Contamination in the data may have dramatic effects on all those techniques based on multivariate estimation of location and scatter, such as Principal Component Analysis or Discriminant Analysis, for instance. On the contrary, by supplying robust estimates of multivariate location and covariance, one could rely on multivariate techniques that are resistant to contamination [Hubert et al., 2008]. Furthermore, the appropriate use of robust estimators may also lead to detect outliers, find unexpected structures in the data and explore the types of occurred departures. There is a growing literature on robust multivariate estimation. The reader is pointed to the book by Farcomeni and Greco [2016] for a recent account on multivariate settings.

Robust estimates of multivariate location and covariance are obtained by attaching a weight to each data point in order to bound the effect of possible outliers on the resulting fit. Weights are determined according to an outlyingness measure, that is a measure of the distance of the multivariate data point from the robust fit. In summary, we can consider three main classes of estimators.

Estimators based on hard trimming: weights are 0-1 and outliers are trimmed. The final estimate is based on a subset of the original data points, whose size is tuned by the user. The Minimum Covariance Determinant (MCD) is undoubtedly one of the most popular techniques [Rousseeuw, 1985, Croux and Haesbroeck, 1999]. 2. 2.

Estimators based on adaptive hard trimming: outliers are trimmed but the final sample is determined adaptively by the data. The main tool is represented by the Forward Search (FS, see Riani et al. [2009], Atkinson and Riani [2012] for a recent account). 3. 3.

Estimators based on soft trimming: outliers are down-weighted, with weights varying in $[0,1]$, and the final estimate consists of a weighted mean and weighted covariance matrix. This feature characterizes M-estimators and related methods such as S-estimators [Lopuhaa, 1989] and MM-estimators [Salibian-Barrera et al., 2006], but also the weighted likelihood estimator (WLE, Markatou et al. [1998]). In this class we also include those methods stemming from projection of multivariate data onto univariate directions, as the Stahel-Donoho estimator, for instance.

The weighting strategy that characterizes the MCD, the FS and M-estimation is based on the inspection of the Mahalanobis distances. Let $y=(y_{1},y_{2},\ldots,y_{p})^{{}^{\top}}$ denote a $p$-variate observation, $p>1$, sampled from a multivariate normal model, $Y\sim N_{p}(\mu,\Sigma)$ with mean vector $\mu=(\mu_{1},\mu_{2},\ldots,\mu_{p})^{{}^{\top}}$ and $p\times p$ covariance matrix $\Sigma$. The Mahalanobis distance is

[TABLE]

Let $(\hat{\mu},\hat{\Sigma})$ be a robust estimate of location and scatter, then data points are discarded or down-weighted according to their distance $d(y;\hat{\mu},\hat{\Sigma})$ from the robust fit: the larger the robust distance the closer to zero the weight and more likely the point will be treated as an outlier.

In a different fashion, the computation of the WLE is not based on such robust distances, but outlyingness is measured according to the agreement between the data, summarized by a non parametric density estimate, and the assumed multivariate normal model. Actually, the weighting scheme based on the computation of a multivariate density estimate becomes troublesome for large dimensions, because of the curse of dimensionality [Huber, 1985, Scott and Wand, 1991]. With growing dimensions the data are more sparse and kernel density estimation may become unfeasible. To the best of our knowledge, non parametric kernel estimation is implemented in a statistical software as R, up to three dimensions. The method by Duong [2007] is an exception since it allows to get a non parametric kernel estimate up to six dimensions. The reader is pointed to Deng and Wickham [2011] for a comparison of several density estimation methods available from R.

This feature represents a serious limitation of the weighted likelihood methodology in a multivariate framework. Such a restriction is much more annoying since all the other multivariate estimators that we have mentioned so far are well behaved in large dimensions, at least up to $p=100$. It is worth to stress here that we only consider the case where the sample size $n$ is larger than the dimension $p$.

In this paper, a novel approach to overcome this hindrance is presented. We introduce a weighting algorithm that is still based on non parametric density estimation, but now it is driven from (robust) distances rather than from the data. Hence, the new algorithm handles a univariate kernel density estimate. We obtain multivariate estimates of location and covariance that are consistent and fully efficient at the assumed multivariate normal model, robust w.r.t. the presence of outliers, with weights that depend on robust distances as well as for the MCD, the FS and M-type estimators.

The rest of the paper is structured as follows. Some background on the weighted likelihood is given in Section 2. The new weighting algorithm is introduced in Section 3, whereas an outlier detection rule is given in Section 4. Some numerical studies are given in Section 5 and real data examples concerning estimation, outlier detection, principal component analysis and discriminant analysi are discussed in Section 6.

2 Background

Let $y=(y_{1},\cdots,y_{n})$ be a random sample from a random variable $Y$ with unknown probability (density) function $m(y;\theta)$, $\theta\in\Theta\subset\mathbb{R}^{p}$, with $p\geq 1$ and let $\hat{M}_{n}$ be the empirical distribution function. A weighted likelihood estimate is defined as the root of the weighted likelihood estimating equation

[TABLE]

where $s(y_{i};\theta)$ denotes the $i$-th contribution to the score function and the weight $w_{i}$ is defined as

[TABLE]

where $[\cdot]^{+}$ denotes the positive part. The function $\delta(y;\theta,\hat{M}_{n})$ is the Pearson residual function [Markatou et al., 1998]

[TABLE]

and $A(\cdot)$ is the Residual Adjustment Function (RAF, Lindsay [1994], Park et al. [2002]). The Pearson residuals are evaluated by comparing a non parametric density estimate

[TABLE]

based on the kernel $k(\cdot;\cdot)$ with bandwidth $h$, and a smoothed version of the model density

[TABLE]

based on the same kernel function. Model smoothing leads to the desired asymptotic behavior of the weights, that will be described below, but, in finite samples, for large sample sizes relative to the dimensionality of the parameter space, its presence/absence does not have an important impact on the estimation process.

The RAF plays the role to bound the effect of large residuals on the fitting procedure, as well as the Huber and Tukey-bisquare function bound large distances in M-type estimation. Here, we consider the families of RAF based on the Power Divergence Measure

[TABLE]

Special cases are maximum likelihood ($\tau=1$, as the weights become all equal to one), Hellinger distance ($\tau=2$), Kullback–Leibler divergence ($\tau\rightarrow\infty$) and Neyman’s Chi–Square ($\tau=-1$). An alternative is represented by the families of RAF based on the Generalized Kullback-Leibler divergence (see Cressie and Read [1984, 1988], Park and Basu [2003] and references therein).

When the model is correctly specified, the Pearson residual function evaluated at the true parameter value converges almost surely to zero, whereas, otherwise, for each value of the parameters, large Pearson residuals detect regions where the observation is unlikely to occur under the assumed model. Hence, those observations lying in such regions are attached a weight that decreases with increasing Pearson residual. Large Pearson residuals and small weights will correspond to data points that are likely to be outliers.

Under classical regularity assumptions regarding the model, the kernel, the RAF and the weight function and a correctly specified model,

the WLE $\hat{\theta}$ is consistent and first order efficient [Markatou et al., 1998] 2. 2.

the robustness weights satisfy

[TABLE]

[Agostinelli, 2002, Agostinelli and Greco, 2013] 3. 3.

the weighted likelihood ratio test and its asymptotically equivalent versions share the same first order asymptotic properties of their genuine likelihood counterparts [Agostinelli and Markatou, 2001] .

3 Weighted likelihood estimation based on robust distances

At the multivariate normal model, the Mahalanobis distance given in (1) satisfies

[TABLE]

at the true parameter values. In order to define a set of weights whose computation does not need the evaluation of a multivariate kernel density estimate and does not suffer from any problem due to large dimensionality, we suggest to focus on the distribution of squared distances rather than observations. In other words, the weighting scheme will be based on Person residuals aiming at measuring the degree of agreement between a univariate kernel density estimate based on the vector of squared distances and their underlying $\chi^{2}_{p}$ distribution at the assumed multivariate normal model. Then, this strategy leads to down-weight those observations that exhibit a large distance from the robust fit.

The behavior of the Pearson residual function and the resulting weight function are exemplified in Figure 1. The true underlying model for the squared distances is assumed to be an $\epsilon$-contaminated model of the form $m(x)=(1-\epsilon)\chi^{2}_{p}(x)+\epsilon\chi^{2}_{p}(x,c)$, where the perturbing component is a non-central $\chi^{2}_{p}$ distribution with non centrality parameter $c$. The mixture model is shown in the left panel, with $p=2,c=5,\epsilon=0.05$. Large squared distances are likely to occur under the contaminating component and are expected to be down-weighted at the $\chi^{2}_{2}$ distribution. The middle panel displays the (asymptotic) Pearson residual function at the $\chi^{2}_{2}$ model: actually, it takes large values at large distances and, hence, detect a region where outlying distances are likely to occur. The weight function based on the Hellinger distance RAF is given in the right panel: it clearly decreases at large distances. The vertical dashed line in the third panel gives the $0.975$- level quantile of the $\chi^{2}_{2}$ distribution: this is the quantile commonly used to declare a large distance and detect outliers in robust multivariate estimation.

The WLE of multivariate location and scatter $(\hat{\mu},\hat{\Sigma})$ is a weighted mean and weighted covariance matrix with data dependent weights. It is a common practice to consider an unbiased weighted likelihood estimates of the covariance matrix, that can be defined as

[TABLE]

where $\hat{w}_{i}=w(\hat{d}_{i}^{2},\hat{\mu},\hat{\Sigma}_{u},\hat{M}(d_{i}^{2}))$, $\hat{d}_{i}=d(y_{i},\hat{\mu},\hat{\Sigma}_{u})$, , $i=1,2,\ldots,n$. Actually, this is the approach currently implemented in the R function cov.wt to get an unbiased estimate of scatter. When all the weights are equal to one, then $(n-1)$ appears in the denominator.

The computation of $(\hat{\mu},\hat{\Sigma}_{u})$ yields an iterative procedure, as illustrated in Algorithm 1. At each iteration, based on the current values $(\hat{\mu},\hat{\Sigma}_{u})$ robust distances are obtained. Then, their non parametric density estimate is fitted based on the chosen kernel and Pearson residuals and weights are updated. This algorithm shares the main features of the iterative procedure developed to obtain weighted likelihood estimates in linear regression [Agostinelli and Markatou, 1998] and generalized linear models [Alqallaf and Agostinelli, 2016]. Actually, in Algorithm 1 at each iteration squared distances and their non-parametric density estimate are updated, whereas the (smoothed) model is held fixed.

Some care is needed in the construction of the kernel density estimate $\hat{m}_{n}(d_{i}^{2})$ since it is not expected to allocate any weight before zero otherwise it will be biased at the boundary [Karunamuni and Alberts, 2005]. In the development of Algorithm 1 four methods are suggested to come through this issue. The first three are designed to obtain an unbiased kernel density estimate over $(0,\infty)$, whereas the fourth is based on the distribution of log-transformed squared distances, moving the problem over the whole real line.

The reflection technique [Silverman, 1986] is based on data augmentation by adding the reflections of all the points in the boundary. Then, it is possible to implement any method originally designed for the whole real line. A reflection kernel can be defines as follows

[TABLE]

where $k(\cdot)$ is a symmetric and differentiable probability density: the reflection of a normal density leads to a folded normal kernel. 2. 2.

A kernel density estimate over $(0,\infty)$ can be also obtained by first log-transforming the squared distances, fitting a non parametric density estimate over the whole real line, i.e. $\hat{m}_{n}(\log d^{2})$, and then back-transforming the fitted density to $(0,\infty)$, i.e. $\hat{m}_{n}(d^{2})=\frac{1}{d^{2}}\hat{m}_{n}(e^{\log d^{2}})$, [Bowman and Azzalini, 1997]. The corresponding smoothed model can be obtained according to the same scheme. In this paper, we make use of the code available from the R-package sm. 3. 3.

The Gamma kernel [Chen, 2000, Jones and Henderson, 2007]

[TABLE]

where $\Gamma(t;a,b)$ denotes the probability density function of a Gamma variate with shape parameter $a$ and scale parameter $b$. This is an appealing alternative that does not involve any data transformation. The role of $y$ and $t$ may also be switched. 4. 4.

Log-transformed squared distances are distributed according to a $\log\chi^{2}_{p}$ model whose probability density function is

[TABLE]

Then, Pearson residuals and weights can be evaluated on this new scale by comparing the fitted kernel density based on log-transformed squared distances with the $\log\chi^{2}_{p}$ distribution.

The smoothing parameter $h$ indexing the kernel function $k(\cdot;\cdot)$ involved in the construction of the Pearson residuals regulates the robustness/efficiency trade-off of the weighted likelihood methodology. An appealing feature is that $h$ may be set independently from the scale of the model. Its value may be determined in order to achieve a fixed expected downweighting level [Markatou et al., 1998], that is the expected number of observations that will be deleted under the specified model, a fixed asymptotic weight for a given point mass contamination $\epsilon$ at the outlying point $y$ [Agostinelli and Markatou, 2001], but also in an adaptive fashion by monitoring the empirical downweighting level $(1-\bar{\omega})$, with $\bar{\omega}=n^{-1}\sum_{i=1}^{n}\hat{w}_{i}$ [Greco, 2016].

Algorithm 1 may be initialized by drawing a large number of random subsets of fixed dimension $(p+1)$. The sample mean and covariance matrix are evaluated over each subsample and used as starting values [Markatou et al., 1998]. A deterministic solution to set initial values can be also implemented, stemming from that described in Hubert et al. [2012]. Strategies to select the best root are given in Markatou et al. [1998] and Agostinelli [2006].

4 The distribution of robust distances

The availability of robust estimates of location and scatter allows to activate some procedures designed to identify multivariate outliers. Actually, the use of robust estimates in place of the sample vector mean and covariance matrix avoids masking and swamping effects in outlier detection: there is masking whenever an outlier is not detected because of the presence of similar outliers, swamping when a genuine observation is flagged as an outlier.

The problem of outlier detection consists in testing the $n$ null hypotheses that each data point is a realization of a multivariate normal distribution $N_{p}(\mu,\Sigma)$. The detection rule will depend on the (asymptotic) distribution of the squared robust distances $d^{2}(y;\hat{\mu},\hat{\Sigma}_{u})$. A common approach to define cut-off values to flag outliers is based on the $\chi^{2}_{p}$ distribution to approximate the distribution of squared robust distances. A more accurate distributional result may be used after the computation of the MCD estimator [Cerioli, 2010], but not in the case of M-type estimation. A rule of thumb is based on the $0.975$-level quantile of the reference distribution. The outliers detection process could also be designed to take into account multiplicity arguments in the simultaneous testing of all the $n$ data points. For instance, cut-off values can be based on a $(1-\gamma)$-level quantile such that the simultaneous testing of all the data points corresponds to a global nominal level $\alpha$, that is $\gamma=1-(1-\alpha)^{1/n}$ [Cerioli, 2010] or by controlling the overall level of the simultaneous testing procedure by the False Discovery Rate [Cerioli and Farcomeni, 2011].

Here, we state a result concerning the distribution of robust distances at the postulated multivariate normal model based on the WLE, that resembles, asymptotically, the classical one on the Mahalanobis distance evaluated at the unbiased MLE [Gnanadesikan and Kettenring, 1972].

Proposition 1 Let $\hat{\theta}=(\hat{\mu},\hat{\Sigma}_{u})$ be the unbiased WLE of multivariate location and scatter. Assume that: (i) the model is correctly specified, that is $m(y)=m(y;\theta_{0})$ for some $\theta_{0}=(\mu_{0},\Sigma_{0})\in\Theta$; (ii) $\hat{\theta}$ is a consistent estimator of $\theta_{0}$; (iii) $\sup_{y}\left|w(y;\hat{\theta},\hat{M}_{n})-1\right|\stackrel{{\scriptstyle p}}{{\longrightarrow}}0$. Then

[TABLE]

The proof is given in the Appendix. It follows the guidelines of the classical proof based on the unbiased MLE that has been revised in Ververidis and Kotropoulos [2008]. It is worth to stress that the same result does not hold for M-type estimation since Huber and Tukey’s bisquare weights do not share the asymptotic behavior of the weights in (2) at the postulated model. A close result has been established in the case of the MCD [Cerioli, 2010].

5 Numerical studies

In this section we investigate the finite sample behavior of the newly proposed WLE of multivariate location and scatter through some numerical studies. The strategies outlined in Section 3 to compute Pearson residuals and the corresponding weights are all considered: folded normal kernel (WLEa), log and back transform (WLEb), log transform with $\log\chi^{2}_{p}$ (WLEc), gamma kernel (WLEd). The weights are based on the Hellinger distance RAF. The multivariate WLE has been also compared with the deterministic MCD and the S-estimator (with Rocke type weights, that has been designed to work properly for large dimensions), evaluated by using the functions from the R package rrcov. The WLE runs on a deterministic algorithm [Hubert et al., 2012], that is based on six initial solutions. For each starting point, Pearson residuals are evaluated and the solution with the lowest fitted probability

[TABLE]

is selected [Agostinelli, 2006]. Then, this initial proposal is iteratively updated according to Algorithm 1.

Several combinations of $(n,p)$ have been taken into account. Data have been generated according to a multivariate normal $N_{p}(0,I)$ with uncorrelated components and unit variance with point mass contamination, that is a percentage $\epsilon$ of outliers is driven by a multivariate normal model $N_{p}(ka,\delta I)$ with $k=1,2,3,\ldots,K$ and $\delta=0.01$ ($k=0$ gives the uncontaminated scenario). When $p\leq 10$, contamination is designed to affect all dimensions, $a=(1,1,\ldots,1)^{{}^{\top}}$, wheres when $p>10$, outliers only contaminate the first five dimensions $a=(1,1,1,1,1,0,\ldots,0)^{{}^{\top}}$. We show results corresponding to a contamination level $\epsilon=20\%$ and two data configurations: in the first $n=100$ and $p=10$, in the second $n=500$ and $p=50$. The most distant point mass contamination has been located at $K=5$ in the first case and at $K=10$ in the second one. All numerical studies are based on 1000 Monte Carlo trials.

The following performance measures were considered:

$||\hat{\mu}||$ 2. 2.

$\log\frac{\textrm{trace}(\hat{\Sigma}_{u})}{p}$ 3. 3.

$\log_{10}\textrm{cond}(\hat{\Sigma}_{u})$ 4. 4.

$||\hat{\Sigma}_{u_{jj}}-1||,\ j=1,2,\ldots,p$ 5. 5.

$||\hat{\Sigma}_{u_{jh}}||,\ j,h=1,2,\ldots,p,j\neq h$ 6. 6.

computational time (in seconds on a 3,4 GHz Intel Core i5).

All of them are expected to be as close as possible to zero.

Figure 2 displays the average performance measures for $n=100$, $p=10$. The WLE provides very accurate results and an appealing behavior compared to the S-estimator and the MCD, whatever the chosen weighting scheme, both for location and covariance estimation. It is worth noting that all the estimators provide less accurate results when outliers are not located at large distances. The WLE becomes less accurate at $k=2$ only and outperforms both the S-estimator and the MCD. The former suffers from contamination still at $k=3$, whereas the latter exhibits the desired performance at $k=2$, as well, but the WLE is to be preferred.

The results corresponding to the second data configuration, with $n=500$ and $p=50$ are showed in Figure 3. In the same fashion as above, contamination is such that all estimators may suffer from lack of robustness when outliers are not located at large distances. The WLE still exhibits a satisfactory behavior. All the estimators share the same features until $k=3$. After that, the MCD shows the desired redescending behavior, whereas the WLE is no more affected by contamination at $k=5$ and the S-estimator does not protect inference from contamination for all considered values of $k$. The fact that the WLE performs remarkably better than the S-estimator is a noticeable result, in that both of them fall in the general category of soft-trimming estimators, as stated in the Introduction.

The sixth panels in both Figures show the computational time. One needs to keep in mind that the comparison with the MCD and S- estimators is unfair, since the WLE is still based on an unoptimized R code, that will be soon available from package wle. Actually, computational time for the WLE remains in a feasible range, even when $p=50$. In the first scenario, the use of folded normal leads to save computational time w.r.t. the other kernels, especially for larger $k$. With growing dimensionality, the reflection kernel is still to be preferred but the procedure based on log and back transform is also competitive.

6 Real data examples

In this section we provide some real data examples concerning multivariate estimation of location and scatter, outlier detection, principal component analysis and discriminant analysis. The proposed weighted likelihood methodology is also compared with other popular robust multivariate tools.

6.1 Multivariate estimation

The StarsCYG data give the effective temperature at the surface and the light intensity, both on a log scale, of 47 stars in the star cluster CYG OB1. Five stars are clear outliers (the points 11, 20, 30, 34 correspond to giant stars that do not lie on the main sequence) and the point 7 also does not share the correlation structure of the remaining 43 stars. Figure 4 displays 0.975-level tolerance ellipses stemming from the proposed WLE, the WLE based on a bivariate kernel density estimate (WLEmulti), the MCD (with $50\%$ breakdown point), the MM-estimator (with $50\%$ breakdown point and $95\%$ shape efficiency) and the MLE. The weighted likelihood contours are based on the asymptotic result given in (3), whereas that stemming from the MCD is based on the distributional result given in Cerioli [2010] and that derived from the MM-estimator is based on the $\chi^{2}_{2}$ distribution. The classical tolerance ellipse is also based on the scaled Beta distribution. All the methods we outlined to compute the weights for the WLE gave very similar results and hence, only the one based on log and back transform of distances is shown. The fitted robust ellipses do not exhibit any significant difference and are all able to catch the correlation structure in the main sequence of stars: the WLE gives a correlation of 0.680, the MCD gives 0.655 and the MM gives 0.691. Moreover, the newly proposed WLE behaves not dissimilarly from the WLE based on the bivariate kernel. On the contrary, the MLE leads to inflated variability and negative correlation.

Figure 5 gives the final robustness weights corresponding to the WLE. All the outliers are given a weight that is very close to zero. The selected smoothing parameter leads to an empirical downweighting level that is about $11\%$ (that is 5 observations over 47).

6.2 Outliers detection

The Auto data give information on technical and insurance characteristics of $n=195$ cars collected in 1985 by the Insurance Institute for Highway Safety, for a total of $p=15$ variables. The car are of two types: running on a gasoline or diesel engine. The are only 20 cars running on diesel that may be identified as outliers w.r.t. the others but several outliers of different nature may arise corresponding, for instance, to cars with peculiar technical features or deserving specific insurance conditions. Figure 6 gives the robust distances corresponding to each car stemming from the WLE. The solid line gives the cut-off value at the 0.975 level, whereas the dashed line gives the cut-off value at the $(1-\gamma)\approx 0.9998702$ level quantile of the scaled Beta distribution in (3). The latter threshold has been computed to take into account multiplicity, in such a way that the simultaneous testing of all the data points correspond to a nominal level 0.025.

The group of cars running on diesel is clearly characterized by the largest distances and is well separated from the remaining cars. The inspection of Figure 6 also unveils some other outlying cars that may exhibit peculiar characteristics. The group of cars running on diesel and the other outliers are clearly spotted by the QQ-plot in Figure 7. The different nature of the several outliers that we have identified can be investigated further by exploring the distance-distance plot in Figure 8. The robust distances based on the WLE are compared with the classical distances based on the MLE. An important feature of such plot is that the cut-off values are determined according to the same scaled Beta distribution, hence being the same on both axes. It is worth noting that the group of car running on a diesel engine would not have been detected by looking at the classical distances based on the MLE. The results driven by the use of the MCD, S and MM estimators are very similar.

6.3 Principal Component Analysis

Principal Component Analysis (PCA) is undoubtedly the most popular technique for dimension reduction. The data are projected onto a lower dimensional sub-space so that they are as spread out as possible. This feature allows to express the covariance structure of the data by means of a small number of new variables (the principal components). These new variables are obtained as linear combination of the original set of variables and are orthogonal each other. The coefficients of the linear combinations are given by the eigenvectors of the covariance matrix and each component accounts for an amount of total variability proportional to the corresponding eigenvalue. PCA is clearly sensitive to the occurrence of outliers, that, in particular, may inflate the variability accounted for by the first components hence leading to wrongly rotated loadings. One approach to robust PCA is based on the eigen-decomposition of a robust estimate of covariance. Here, we employ the WLE to perform a robust PCA on the Auto data. The same example has been discussed in Farcomeni and Greco [2016]. Let us consider the first $k=3$ components. The percentage of explained variance from standard PCA is $76.5\%$ whereas the robust analysis gives a smaller value of $73.6\%$. In order to better explain the deleterious effect of outliers on standard PCA and the effectiveness of our weighted approach, Figure 9 displays the pairwise score-plots based on the first three components. The group of cars running on diesel is clearly spotted by the robust components in the left panels, whereas this does not happens in the right panel. A typical effect due to the presence of outliers can be seen in the last panel: the second and third component from standard PCA still show a linear trend and only the effect of the outlying cars leads to a null correlation.

Robust PCA is an effective tool in outlier detection when the dimensionality is not of a manageable size. The usual tool is an outlier map, displayed in Figure 10, that is obtained by plotting for each data point its score and orthogonal distance: the group of outlying cars is clearly separated from the rest but also other outlying points are visible. Guidelines to find the cut-off values are given in Hubert et al. [2005].

We only mention here, that the WLE of multivariate location and scatter could be used in the technique developed by Greco and Farcomeni [2016] to obtain sparse and robust PCA.

6.4 Discriminant Analysis

Discriminant analysis is concerned with the problem of assigning data to one of several groups. The observations within each group are assumed to arise from a multivariate normal distribution. In linear discriminant analysis (LDA) it is assumed homogeneity of the covariance matrices, whereas in quadratic discriminant analysis (QDA) the groups are allowed to have different scatters. Let $\pi_{j},j=1,2,\ldots,k$ denote the prior probabilities. The linear discriminant rule classifies observations by maximizing

[TABLE]

and the quadratic discriminant rule classifies observations by maximizing

[TABLE]

where $\hat{\pi}_{j}=\frac{n_{j}}{\sum_{j=1}^{k}n_{j}}$ is an estimate of prior probabilities to be used when prior information is not available, $\hat{\mu}_{j}$ is an estimate of the group vector mean, $\hat{\Sigma}_{p}$ is a pooled estimate of the common scatter and, $\hat{\Sigma}_{j}$ is an estimate of the group scatter matrix. Actually, an appealing approach to define a discriminant function that is not prone to contamination in the data is based on robust estimates of location and common covariance matrix [Hubert and Van Driessen, 2004, He and Fung, 2000]. Here, we apply weighted likelihood to perform robust LDA and QDA. In particular, we consider two different strategies to obtain a robust pooled estimate of the covariance matrix, in a fashion similar to what happens when using the MCD estimator [Todorov and Pires, 2007]. By paralleling the standard technique, the first estimate (WLEA) is obtained by averaging the unbiased estimates from each group as follows:

[TABLE]

The second estimate (WLEB) can be obtained after the following steps. First center the data from each group by a robust estimate of location $\hat{\mu}_{j0}$; one could use the L1 (spatial) median, for instance. Then, evaluate the WLE $(\hat{\mu}_{p},\hat{\Sigma}_{p}^{B})$ from all the centered data and update the group vector means as $\hat{\mu}_{j}^{B}=\hat{\mu}_{j0}+\hat{\mu}_{p}$. The latter approach needs only one robust estimate of covariance rather than one for each group as in the former one. Nevertheless, an alternative, even if slightly more demanding, still consists in running Algorithm 1 for each group and centering the data by using the WLE of location from each of them in the first step.

Let us apply weighted LDA and QDA to the Diabetes data. These data consists of three measurement of plasma, glucose, insuline and sspg, made on 145 non-obese adult patients classified into three groups: normal subjects, chemical diabetes and overt diabetes.

The data, along with the fitted groups according to LDA based on WLEB and QDA stemming from group-wise WLEs, are displayed in Figure 11. The fitted groups appears as 0.975-level tolerance ellipses. It is worth noting the differences among the two techniques concerning, in particular, the peculiar nature of the overt diabetes group. Actually, the nature of correlation between glucose and sspg and insuline and sspg in the third group is different from what happens in the other two groups. The entries in Table 1 give the estimates of the misclassification rate based on all the data (ALL) and on leave-one-out cross validation (CV) based on the WLE, MLE and MCD, for LDA and QDA. The use of robust estimates of multivariate location and scatter improves classification accuracy over the standard approach based on the MLE and the WLE performs satisfactory compared to the MCD. In particular, both LDA based on WLEB and QDA stemming from the group-wise WLEs lead to the same results.

Appendix

Lemma 1 Let $A=\omega\hat{\Sigma}=\omega\gamma\hat{\Sigma}_{u}$, $A_{(-i)}=\omega_{(-i)}\hat{\Sigma}_{(-i)}=\omega_{(-i)}\gamma_{(-i)}\hat{\Sigma}_{u(-i)}$, with $\omega=\sum_{i=1}^{n}w_{i}$ and the subscript $(-i)$ denote that the $i$-th observation has been removed. Then

[TABLE]

Proof. It is known that:

[TABLE]

Then

[TABLE]

and by replacing

[TABLE]

the equality stated in (4) is obtained.

Proof of Proposition 1. Let us consider the random quantity

[TABLE]

By using Lemma 1, we find that

[TABLE]

where $X_{i}=\sqrt{\frac{\omega w_{i}}{\omega^{2}_{(-i)}\gamma_{(-i)}}}(Y_{i}-\hat{\mu})$. At the assumed multivariate normal model, by applying the asymptotic results on the distribution of the WLE and the behavior of the weights, we have that

[TABLE]

and

[TABLE]

Then, from the Slutsky theorem, $X_{i}\stackrel{{\scriptstyle d}}{{\rightarrow}}\frac{1}{\sqrt{n-2}}N(0,\Sigma)$.

In a similar fashion, at the assumed multivariate normal model, we can state that $\hat{\Sigma}_{u(-1)}\stackrel{{\scriptstyle d}}{{\rightarrow}}W_{p}(\Sigma,n-1)$. By using the fact that, conditioning on the weights, $X_{i}$ and $\hat{\Sigma}_{u(-i)}$ are independent, since $\hat{\Sigma}_{u(-i)}$ does not involve $Y_{i}$, we have that, asymptotically, the distribution of $X_{i}^{{}^{\top}}\hat{\Sigma}_{u(-i)}^{-1}X_{i}$ is the Hotelling distribution. Then, by using standard results, it is possible to derive the asymptotic distribution of $R_{(-i)}$, that is

[TABLE]

According to Lemma 1, we are also allowed to write

[TABLE]

and

[TABLE]

Since $(1-R_{(-i)})\stackrel{{\scriptstyle d}}{{\rightarrow}}\textrm{Beta}\left(\frac{p}{2},\frac{n-p-1}{2}\right)$ and $\frac{\omega_{(-i)}\gamma}{w_{i}}\stackrel{{\scriptstyle p}}{{\rightarrow}}\frac{(n-1)^{2}}{n}$, then, again by Slutsky theorem, the result stated in (3) is obtained and the proof is concluded.