Multivariate outlier detection based on a robust Mahalanobis distance with shrinkage estimators

TL;DR

This paper introduces a robust multivariate outlier detection method using shrinkage estimators for Mahalanobis distance, demonstrating improved accuracy and efficiency over existing techniques, especially with non-normal data.

Contribution

The paper proposes a novel robust Mahalanobis distance approach based on shrinkage estimators, with properties like affine equivariance and high computational efficiency.

Findings

High correct detection rates in simulations

Low false detection rates across datasets

Reduced computation time compared to existing methods

Abstract

A collection of robust Mahalanobis distances for multivariate outlier detection is proposed, based on the notion of shrinkage. Robust intensity and scaling factors are optimally estimated to define the shrinkage. Some properties are investigated, such as affine equivariance and breakdown value. The performance of the proposal is illustrated through the comparison to other techniques from the literature, in a simulation study and with a real dataset. The behavior when the underlying distribution is heavy-tailed or skewed, shows the appropriateness of the method when we deviate from the common assumption of normality. The resulting high correct detection rates and low false detection rates in the vast majority of cases, as well as the significantly smaller computation time shows the advantages of our proposal.

Tables

Table 1. Table 1 : Combinations of location and dispersion

Name	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
${\hat{\mu }}_R$	${\hat{\mu }}_{CCM}$	${\hat{\mu }}_{Sh\left(CCM\right)}$	${\hat{\mu }}_{Sh\left(CCM\right)}$	${\hat{\mu }}_{MM}$	${\hat{\mu }}_{Sh\left(MM\right)}$	${\hat{\mu }}_{Sh\left(MM\right)}$
${\hat{\Sigma }}_R$	${\hat{\Sigma }}_{Sh\left(CCM\right)}$	${\hat{\Sigma }}_{Sh\left(CCM\right)}$	${\hat{\Sigma }}_{Sh\left(Sh\left(CCM\right)\right)}$	${\hat{\Sigma }}_{Sh\left(MM\right)}$	${\hat{\Sigma }}_{Sh\left(MM\right)}$	${\hat{\Sigma }}_{Sh\left(Sh\left(MM\right)\right)}$

Table 2. Table 2 : Correct detection rates, with Normal distribution.

$\delta =5$		$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	1	1	0,9000	1	1	1	1	1	1	1	1
	0.2	0,8700	0,8700	0,5100	0,9500	0,9941	1	1	1	1	1	1
	0.3	0,0600	0,0600	0,9800	0,1500	0,5719	0,8766	0,8782	0,8782	0,9146	0,9090	0,9130
10	0.1	0,9900	0,9900	0,8600	1	1	1	1	1	1	1	1
	0.2	0,2800	0,2800	0,4600	0,9416	1	1	1	1	1	1	1
	0.3	0	0	0,9900	0,1612	0,7205	0,8774	0,8747	0,8750	0,9711	0,9672	0,9711
30	0.1	0,1900	0,1900	1	1	1	1	1	1	1	1	1
	0.2	0	0	0,1	1	1	1	1	1	1	1	1
	0.3	0	0	0,6100	0,0100	0,9407	0,5308	0,5275	0,5286	0,9990	0,9988	0,9991

Table 3. Table 3 : Correct detection rates, with Normal distribution.

$\delta =10$	$\alpha$	$\lambda =1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,8480	0,8465	0,9900	1	1	1	1	1	1	1	1
	0.3	0,2190	0,1976	0,9307	0,9591	0,9991	1	1	1	1	1	1
10	0.1	1	1	0,9800	1	1	1	1	1	1	1	1
	0.2	0,8623	0,8548	0,6558	1	1	1	1	1	1	1	1
	0.3	0,2280	0,2046	0,4618	0,9911	1	1	1	1	1	1	1
30	0.1	1	1	0,8919	1	1	1	1	1	1	1	1
	0.2	0,4879	0,4654	0,0125	1	1	1	1	1	1	1	1
	0.3	0,0810	0,0509	0,1087	1	1	1	1	1	1	1	1

Table 4. Table 4 : Simulation results for correlated data.

	MCD		Adj MCD		Kurtosis		OGK		COM		RMD-Shrinkage
$\alpha$	c	f	c	f	c	f	c	f	c	f	c	f
0.1	1	0,0397	1	0,0226	1	0,0371	1	0,0736	1	0,0025	1	0,0128
0.2	0,8659	0,0127	0,8565	0,0062	0,8771	0,0453	0,9792	0,0533	1	0,0011	1	0,0013
0.3	0,1504	0,0762	0,1238	0,0614	0,8186	0,0443	0,4780	0,0460	0,8302	0,0001	0,9274	0

Table 5. Table 5 : Correct and false detection rates of RMD-Shrinkage for transformed data.

			$\delta =5$		$\delta =10$
	$p$	$\alpha$	c	f	c	f
$\lambda =0.1$	5	0.1	1	0,0454	1	0,0455
		0.2	1	0,0155	1	0,0165
		0.3	0,9709	0,0034	1	0,0023
	10	0.1	1	0,0328	1	0,0252
		0.2	1	0,0088	1	0,0062
		0.3	0,9844	0,0023	1	0,0009
	30	0.1	1	0,0089	1	0,0074
		0.2	1	0,0006	1	0,0003
		0.3	1	0	1	0
$\lambda =1$	5	0.1	1	0,0451	1	0,0400
		0.2	1	0,0189	1	0,0120
		0.3	0,9344	0,0039	1	0,0046
	10	0.1	1	0,0282	1	0,0279
		0.2	1	0,0113	1	0,0071
		0.3	0,9872	0,0020	1	0,0010
	30	0.1	1	0,0093	1	0,0072
		0.2	1	0,0006	1	0,0004
		0.3	1	0	1	0

Table 6. Table 6 : Simulation results for breakdown value.

$n=1000$		Symmetric		Asymmetric
$p$	$\alpha$	c	f	c	f
10	0.1	1	0,0055	1	0,0047
	0.2	1	0,0001	1	0,0002
	0.3	1	0	1	0
	0.4	1	0	1	0
	0.45	1	0	1	0
30	0.1	1	0,0002	1	0,0002
	0.2	1	0	1	0
	0.3	1	0	1	0
	0.4	1	0	1	0
	0.45	1	0	1	0
50	0.1	1	0	1	0
	0.2	1	0	1	0
	0.3	1	0	1	0
	0.4	1	0	1	0
	0.45	1	0	1	0
80	0.1	1	0	1	0
	0.2	1	0	1	0
	0.3	1	0	1	0
	0.4	1	0	1	0
	0.45	1	0	1	0
100	0.1	1	0	1	0
	0.2	1	0	1	0
	0.3	1	0	1	0
	0.4	1	0	1	0
	0.45	1	0	1	0

Table 7. Table 7 : Computational times with Normal data.

$\delta =5$	$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMD-Shrinkage
5	0.1	1,0951	0,7670	0,1880	0,1087	0,0228	0,0096
	0.2	0,7619	0,7910	0,0499	0,0203	0,0088	0,0085
	0.3	0,7605	0,8304	0,0266	0,0196	0,0089	0,0074
	Mean	0,8725	0,7961	0,0882	0,0495	0,0135	0,0085
10	0.1	1,3184	0,9970	0,2191	0,1626	0,0247	0,0200
	0.2	1,0329	1,0477	0,1358	0,0793	0,0120	0,0118
	0.3	0,9685	1,0641	0,0482	0,0865	0,0128	0,0108
	Mean	1,1066	1,0363	0,1344	0,1095	0,0165	0,0142
30	0.1	6,2387	6,0934	0,7154	0,8969	0,2000	0,2206
	0.2	5,8676	6,3999	1,4635	0,8158	0,1687	0,1804
	0.3	5,9453	7,0405	1,6572	0,8407	0,1669	0,1674
	Mean	6,0172	6,5113	1,2787	0,8511	0,1785	0,1895

Table 8. Table 8 : Detected outliers.

Method	Out_in	Out_out	Out_total
MCD	72	4	76
Adj MCD	64	4	68
Kurtosis	158	4	162
OGK	144	4	148
COM	59	4	63
RMD-Shrinkage	25	3	28

Table 9. (a) Inside and outside the fences.

Method	Out_in	Out_out	Out_total
MCD	72	4	76
Adj MCD	64	4	68
Kurtosis	158	4	162
OGK	144	4	148
COM	59	4	63
RMD-Shrinkage	25	3	28

Table 10. (b) Inside the 50 % percent 50 50\% of the most central data.

Method	Out_in	Out_total
MCD	29	76
Adj MCD	27	68
Kurtosis	65	162
OGK	58	148
COM	20	63
RMD-Shrinkage	7	28

Table 11. Table 9 : Computational times for each methods with the WDBC dataset.

Method	MCD	Adj MCD	Kurtosis	OGK	COM	RMD-Shrinkage
Times in sec.	12,155	12,3378	6,3077	3,5325	0,3534	0,3299

Table 12. Table 10: False detection rates with Normal distribution α = 0 𝛼 0 \alpha=0 .

$p$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0,06440	0,04640	0,02630	0,08120	0,00480	0,03140	0,02820	0,02770	0,03100	0,02900	0,00316
10	0,11760	0,09830	0,07110	0,09580	0,00217	0,00245	0,00222	0,00215	0,00172	0,00161	0,00160
30	0,06276	0,03922	0,00804	0,09084	0,00008	0,00003	0,00003	0,00003	0,00003	0,00002	0,00001

Table 13. Table 11: Correct detection rates with Normal distribution.

$\delta =5$		$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	1	1	0,9000	1	1	1	1	1	1	1	1
	0.2	0,8700	0,8700	0,5100	0,9500	0,9941	1	1	1	1	1	1
	0.3	0,0600	0,0600	0,9800	0,1500	0,5719	0,8766	0,8782	0,8782	0,9146	0,9090	0,9130
10	0.1	0,9900	0,9900	0,8600	1	1	1	1	1	1	1	1
	0.2	0,2800	0,2800	0,4600	0,9416	1	1	1	1	1	1	1
	0.3	0	0	0,9900	0,1612	0,7205	0,8774	0,8747	0,8750	0,9711	0,9672	0,9711
30	0.1	0,1900	0,1900	1	1	1	1	1	1	1	1	1
	0.2	0	0	0,1	1	1	1	1	1	1	1	1
	0.3	0	0	0,6100	0,0100	0,9407	0,5308	0,5275	0,5286	0,9990	0,9988	0,9991
$\delta =5$	$\alpha$	$\lambda =1$
5	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,8578	0,8486	0,9602	0,9654	0,9975	1	1	1	1	1	1
	0.3	0,1955	0,1664	0,9336	0,5792	0,8735	0,8935	0,8947	0,8938	0,8740	0,8698	0,8755
10	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,9016	0,8935	0,7212	0,9993	1	1	1	1	1	1	1
	0.3	0,2375	0,2080	0,5935	0,6108	0,9505	0,8846	0,8838	0,8838	0,9608	0,9581	0,9621
30	0.1	1	1	0,8816	1	1	1	1	1	1	1	1
	0.2	0,4461	0,4232	0,0154	1	1	1	1	1	1	1	1
	0.3	0,0823	0,0532	0,1483	0,9772	0,9990	0,7142	0,7035	0,7059	1	1	1

Table 14. Table 12: Correct detection rates with Normal distribution.

$\delta =10$		$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	1	1	0,9200	1	1	1	1	1	1	1	1
	0.2	0,8200	0,8200	0,6500	1	1	1	1	1	1	1	1
	0.3	0,1000	0,1000	1	0,7400	1	1	1	1	1	1	1
10	0.1	1	1	0,9200	1	1	1	1	1	1	1	1
	0.2	0,7200	0,7200	0,4400	1	1	1	1	1	1	1	1
	0.3	0,0500	0,0500	0,9700	0,7500	1	1	1	1	1	1	1
30	0.1	0,8800	0,8800	1	1	1	1	1	1	1	1	1
	0.2	0	0	0,1200	1	1	1	1	1	1	1	1
	0.3	0	0	0,5400	0,9300	1	1	1	1	1	1	1
$\delta =10$	$\alpha$	$\lambda =1$
5	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,8480	0,8465	0,9900	1	1	1	1	1	1	1	1
	0.3	0,2190	0,1976	0,9307	0,9591	0,9991	1	1	1	1	1	1
10	0.1	1	1	0,9800	1	1	1	1	1	1	1	1
	0.2	0,8623	0,8548	0,6558	1	1	1	1	1	1	1	1
	0.3	0,2280	0,2046	0,4618	0,9911	1	1	1	1	1	1	1
30	0.1	1	1	0,8919	1	1	1	1	1	1	1	1
	0.2	0,4879	0,4654	0,0125	1	1	1	1	1	1	1	1
	0.3	0,0810	0,0509	0,1087	1	1	1	1	1	1	1	1

Table 15. Table 13: False detection rates with Normal distribution δ = 5 𝛿 5 \delta=5 .

		$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	0,0327	0,0184	0,0336	0,0640	0,0040	0,0080	0,0070	0,0073	0,0069	0,0067	0,0076
	0.2	0,0265	0,0171	0,0512	0,0600	0,0028	0,0015	0,0015	0,0012	0,0013	0,0013	0,0013
	0.3	0,1504	0,1247	0,0373	0,1641	0,0015	0	0	0	0	0	0
10	0.1	0,0735	0,0529	0,1167	0,0823	0,0027	0,0055	0,0057	0,0048	0,0025	0,0024	0,0027
	0.2	0,1566	0,1330	0,2803	0,0667	0,0023	0,0002	0,0002	0,0001	0,0001	0,0001	0,0001
	0.3	0,2485	0,2206	0,0767	0,2502	0,0010	0	0	0	0	0	0
30	0.1	0,0767	0,0507	0,0078	0,0699	5E-05	0,0007	0,0006	0,0006	0,0001	4,4E-05	0,0001
	0.2	0,1079	0,0784	0,0565	0,0552	3E-05	0	0	0	0	0	0
	0.3	0,1491	0,1150	0,0547	0,5030	3E-05	0	0	0	0	0	0
	$\alpha$	$\lambda =1$
5	0.1	0,0339	0,0198	0,0341	0,0636	0,0033	0,0084	0,0073	0,0069	0,0075	0,0071	0,0077
	0.2	0,0110	0,0055	0,0347	0,0507	0,0017	0,0014	0,0012	0,0012	0,0010	0,0009	0,0012
	0.3	0,0374	0,0260	0,0379	0,0383	0,0001	0	0	0	0	0	0
10	0.1	0,0704	0,0505	0,0917	0,0842	0,0023	0,0047	0,0042	0,0041	0,0027	0,0022	0,0029
	0.2	0,0270	0,0171	0,0860	0,0610	0,0010	0,0002	0,0002	0,0002	0,0001	0,0001	0,0001
	0.3	0,0851	0,0706	0,0782	0,0457	0,0003	0	0	0	0	0	0
30	0.1	0,0347	0,0115	0,0081	0,0713	0,0001	0,0008	0,0007	0,0007	0,0001	0,0001	0,0001
	0.2	0,0357	0,0198	0,0066	0,0544	0	0	0	0	0	0	0
	0.3	0,0552	0,0343	0,0077	0,0366	0	0	0	0	0	0	0

Table 16. Table 14: False detection rates with Normal distribution δ = 10 𝛿 10 \delta=10 .

		$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	0,0309	0,0154	0,0301	0,0681	0,0022	0,0090	0,0076	0,0066	0,0066	0,0058	0,0069
	0.2	0,0337	0,0248	0,0414	0,0485	0,0036	0,0005	0,0004	0,0004	0,0011	0,0008	0,0007
	0.3	0,1434	0,1183	0,0309	0,0603	0,0031	0	0	0	0	0	0
10	0.1	0,0665	0,0464	0,1147	0,0772	0,0017	0,0041	0,0047	0,0042	0,0034	0,0030	0,0035
	0.2	0,0809	0,0647	0,2827	0,0606	0,0022	0,0002	0,0001	0,0001	0	0	0
	0.3	0,2353	0,2088	0,0934	0,0895	0,0008	0	0	0	0	0	0
30	0.1	0,0415	0,0174	0,0081	0,0715	4E-05	0,0011	0,0010	0,0011	0,0001	0,0001	0,0001
	0.2	0,1072	0,0776	0,0481	0,0561	3E-05	0	0	0	0	0	0
	0.3	0,1500	0,1163	0,0578	0,0623	0,0001	0	0	0	0	0	0
	$\alpha$	$\lambda =1$
5	0.1	0,0332	0,0174	0,0297	0,0682	0,0021	0,0076	0,0068	0,0069	0,0058	0,0052	0,0062
	0.2	0,0136	0,0058	0,0352	0,0535	0,0023	0,0011	0,0012	0,0011	0,0008	0,0008	0,0007
	0.3	0,0385	0,0289	0,0372	0,0397	0,0007	0	0	0	0,0001	0	0
10	0.1	0,0668	0,0469	0,0907	0,0771	0,0020	0,0036	0,0036	0,0031	0,0025	0,0021	0,0026
	0.2	0,0278	0,0170	0,0930	0,0629	0,0009	0,0002	0,0001	0,0002	0	0	0
	0.3	0,0856	0,0694	0,0685	0,0440	0,0004	0	0	0	0	0	0
30	0.1	0,0351	0,0122	0,0077	0,0735	4E-05	0,0009	0,0008	0,0009	0,0001	2,21E-05	0,0001
	0.2	0,0334	0,0181	0,0053	0,0535	0	0	0	0	0	0	0
	0.3	0,0577	0,0368	0,0085	0,0388	0	0	0	0	0	0	0

Table 17. Table 15: Computational times with Normal data δ = 5 𝛿 5 \delta=5 and λ = 1 𝜆 1 \lambda=1 .

$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMD-Shrinkage
5	0.1	0,8547	0,8078	0,0959	0,0181	0,0070	0,0029
	0.2	1,2146	0,7129	0,0763	0,0186	0,0061	0,0034
	0.3	1,0064	0,7544	0,0612	0,0176	0,0063	0,0025
	Mean	1,0252	0,7584	0,0778	0,0181	0,0065	0,0030
10	0.1	1,0090	1,1250	0,1592	0,0793	0,0113	0,0047
	0.2	1,0135	1,0448	0,1679	0,0623	0,0100	0,0025
	0.3	1,0335	1,0595	0,1515	0,0612	0,0091	0,0009
	Mean	1,0187	1,0765	0,1595	0,0676	0,0101	0,0027
30	0.1	6,5263	6,2629	0,4788	0,9623	0,2530	0,2752
	0.2	6,1268	6,2031	0,5737	0,8317	0,1700	0,2139
	0.3	5,9068	6,1767	0,5034	0,9472	0,1791	0,2101
	Mean	6,1866	6,2142	0,5186	0,9137	0,2007	0,2331

Table 18. Table 16: Computational times with Normal data δ = 10 𝛿 10 \delta=10 and λ = 0.1 𝜆 0.1 \lambda=0.1 .

$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMD-Shrinkage
5	0.1	0,7704	0,7425	0,1090	0,0195	0,0107	0,0026
	0.2	0,6878	0,6534	0,0573	0,0195	0,0079	0,0033
	0.3	0,6990	0,7434	0,0294	0,0268	0,0072	0,0054
	Mean	0,7191	0,7131	0,0652	0,0219	0,0086	0,0038
10	0.1	1,0631	1,1115	0,1297	0,0824	0,0101	0,0066
	0.2	1,1940	0,9625	0,1179	0,0719	0,0080	0,0041
	0.3	1,0673	0,9902	0,0756	0,0663	0,0097	0,0047
	Mean	1,1081	1,0214	0,1077	0,0735	0,0093	0,0051
30	0.1	6,0350	6,0913	0,7347	0,8205	0,1701	0,1697
	0.2	6,4506	6,1627	0,7385	0,7491	0,1837	0,1572
	0.3	6,2114	6,1146	1,1310	0,7342	0,1714	0,1281
	Mean	6,2324	6,1229	0,8681	0,7679	0,1750	0,1517

Table 19. Table 17: Computational times with Normal data δ = 10 𝛿 10 \delta=10 and λ = 1 𝜆 1 \lambda=1 .

$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMD-Shrinkage
5	0.1	0,7956	0,8301	0,0814	0,0176	0,0078	0,0037
	0.2	0,7939	0,7684	0,0836	0,0189	0,0124	0,0019
	0.3	0,8614	0,7207	0,0662	0,0183	0,0049	0,0014
	Mean	0,8170	0,7731	0,0770	0,0183	0,0084	0,0023
10	0.1	0,9990	1,0609	0,1350	0,0634	0,0117	0,0047
	0.2	1,0917	1,1028	0,1613	0,0682	0,0093	0,0049
	0.3	1,0111	1,1860	0,1610	0,0766	0,0089	0,0025
	Mean	1,0340	1,1166	0,1524	0,0694	0,0100	0,0040
30	0.1	5,7161	5,6622	0,5731	0,7563	0,1693	0,1220
	0.2	5,6394	5,6993	0,5821	0,7191	0,1654	0,1083
	0.3	5,7104	5,7975	0,9465	0,7193	0,1561	0,1138
	Mean	5,6886	5,7197	0,7006	0,7316	0,1636	0,1147

Table 20. Table 18: False detection rates with Student-t distribution with 3 d.f, α = 0 𝛼 0 \alpha=0 .

$p$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0,12860	0,11060	0,17670	0,21320	0,07730	0,20290	0,20070	0,20030	0,17950	0,17660	0,18370
10	0,17560	0,15640	0,31610	0,24330	0,08410	0,19040	0,19550	0,19570	0,11780	0,11220	0,12760
30	0,15820	0,13530	0,20380	0,28760	0,08270	0,16660	0,16430	0,16490	0,11220	0,11130	0,11400

Table 21. Table 19: Correct detection rates with Student-t with 3 d.f. and δ = 5 𝛿 5 \delta=5 .

	$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	0,9300	0,9300	0,4900	1	1	1	1	1	1	1	1
	0.2	0,1800	0,1800	0,2902	0,5300	0,6872	1	1	1	0,9927	0,9915	0,9915
	0.3	0,0007	0,0007	0,6890	0,0007	0,0411	0,6573	0,6446	0,6441	0,6878	0,6758	0,7072
10	0.1	0,5107	0,5107	0,3407	0,9800	1	1	1	1	1	1	1
	0.2	0,0009	0,0009	0,1809	0,3709	0,6576	1,0000	1,0000	1,0000	0,9997	0,9997	1
	0.3	0,0100	0,0100	0,6200	0,0203	0,0400	0,6500	0,6525	0,6522	0,8021	0,7947	0,8137
30	0.1	0,0003	0,0003	0,2000	1	1	1	1	1	1	1	1
	0.2	0,0004	0,0004	0,1703	0,1905	0,6055	1	1	1	1	1	1
	0.3	0	0	1	0,0002	0	0,1689	0,1638	0,1649	0,7437	0,7417	0,7628
$p$	$\alpha$	$\lambda =1$
5	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,7445	0,7263	0,8974	0,9177	0,9155	0,9970	0,9952	0,9952	0,9922	0,9919	0,9924
	0.3	0,1844	0,1591	0,8421	0,4498	0,4732	0,8442	0,8485	0,8457	0,8411	0,8360	0,8377
10	0.1	0,9826	0,9819	0,9492	1	1	1	1	1	1	1	1
	0.2	0,5344	0,5237	0,5320	0,9341	0,9629	1	1	1	1	0,9996	0,9996
	0.3	0,1864	0,1637	0,6358	0,3756	0,4901	0,8602	0,8571	0,8575	0,9126	0,9063	0,9126
30	0.1	0,9923	0,9917	0,4932	1	1	1	1	1	1	1	1
	0.2	0,1823	0,1581	0,2148	1	1	1	1	1	1	1	1
	0.3	0,1658	0,1403	0,5935	0,4170	0,8556	0,6581	0,6485	0,6501	0,9670	0,9664	0,9689

Table 22. Table 20: Correct detection rates with Student-t with 3 d.f. and δ = 10 𝛿 10 \delta=10 .

	$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	0,9900	0,9900	0,6600	1	1	1	1	1	1	1	1
	0.2	0,5400	0,5400	0,4700	0,9600	1	1	1	1	1	1	1
	0.3	0,0400	0,0400	0,8985	0,3803	0,7211	0,9900	0,9900	0,9900	1	1	0,9900
10	0.1	0,9300	0,9300	0,4900	1	1	1	1	1	1	1	1
	0.2	0,1000	0,1000	0,2200	0,9300	1	1	1	1	1	1	1
	0.3	0	0	0,7804	0,2200	0,7709	1	1	1	1	1	1
30	0.1	0,0200	0,0200	0,2594	1	1	1	1	1	1	1	1
	0.2	0	0	0,3098	1	1	1	1	1	1	1	1
	0.3	0,0001	0,0001	1	0,0003	0,7862	1	1	1	1	1	1
$p$	$\alpha$	$\lambda =1$
5	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,8890	0,8859	0,9580	1	1	1	1	1	1	1	1
	0.3	0,2324	0,2127	0,9470	0,7598	0,9791	0,9998	0,9998	0,9998	1	1	1
10	0.1	1	1	0,9830	1	1	1	1	1	1	1	1
	0.2	0,7050	0,6964	0,5019	1	1	1	1	1	1	1	1
	0.3	0,2563	0,2314	0,6768	0,8453	0,9914	1	1	1	1	1	1
30	0.1	1	1	0,5614	1	1	1	1	1	1	1	1
	0.2	0,2407	0,2166	0,2203	1	1	1	1	1	1	1	1
	0.3	0,1640	0,1394	0,7130	0,9624	1	1	1	1	1	1	1

Table 23. Table 21: False detection rates with Student-t with 3 d.f. and δ = 5 𝛿 5 \delta=5 .

		$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	0,0901	0,0741	0,2020	0,1840	0,0641	0,1305	0,1293	0,1289	0,1155	0,1153	0,1208
	0.2	0,1591	0,1381	0,3055	0,2014	0,0494	0,0741	0,0743	0,0748	0,0668	0,0656	0,0694
	0.3	0,2237	0,1972	0,2046	0,4273	0,0654	0,0299	0,0300	0,0300	0,0298	0,0298	0,0310
10	0.1	0,1724	0,1528	0,3689	0,2149	0,0703	0,1876	0,1851	0,1874	0,1289	0,1259	0,1377
	0.2	0,2460	0,2224	0,4709	0,2743	0,0534	0,0896	0,0882	0,0890	0,0643	0,0636	0,0687
	0.3	0,2978	0,2723	0,3106	0,5624	0,0833	0,0411	0,0411	0,0410	0,0375	0,0364	0,0387
30	0.1	0,1775	0,1523	0,3154	0,2626	0,0717	0,3472	0,3450	0,3452	0,1563	0,1557	0,1581
	0.2	0,2116	0,1831	0,4937	0,4232	0,0551	0,1459	0,1454	0,1452	0,0779	0,0777	0,0785
	0.3	0,2562	0,2235	0,1669	0,8151	0,0991	0,0418	0,0418	0,0417	0,0401	0,0400	0,0446
	$\alpha$	$\lambda =1$
5	0.1	0,0809	0,0662	0,2027	0,1903	0,0664	0,1361	0,1333	0,1346	0,1194	0,1183	0,1232
	0.2	0,0619	0,0500	0,1846	0,1613	0,0513	0,0800	0,0799	0,0796	0,0711	0,0702	0,0727
	0.3	0,1098	0,0951	0,1537	0,1340	0,0362	0,0375	0,0370	0,0372	0,0384	0,0380	0,0394
10	0.1	0,1209	0,1032	0,2985	0,2166	0,0769	0,1922	0,1910	0,1927	0,1414	0,1383	0,1497
	0.2	0,1235	0,1073	0,2940	0,1805	0,0538	0,0986	0,0960	0,0962	0,0710	0,0701	0,0754
	0.3	0,1705	0,1523	0,2096	0,1712	0,0360	0,0392	0,0394	0,0397	0,0361	0,0355	0,0375
30	0.1	0,1039	0,0819	0,1931	0,2587	0,0694	0,3466	0,3454	0,3460	0,1547	0,1540	0,1570
	0.2	0,1514	0,1291	0,1888	0,2311	0,0557	0,1524	0,1519	0,1520	0,0781	0,0778	0,0788
	0.3	0,1525	0,1304	0,1681	0,2248	0,0392	0,0445	0,0447	0,0444	0,0364	0,0363	0,0368

Table 24. Table 22: False detection rates with Student-t with 3 d.f. and δ = 10 𝛿 10 \delta=10 .

		$\lambda =0.1$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	0,0831	0,0669	0,2094	0,1955	0,0636	0,1378	0,1352	0,1368	0,1234	0,1222	0,1273
	0.2	0,1090	0,0935	0,2544	0,1602	0,0442	0,0693	0,0690	0,0687	0,0629	0,0609	0,0647
	0.3	0,2158	0,1904	0,1578	0,2988	0,0365	0,0330	0,0327	0,0331	0,0294	0,0295	0,0314
10	0.1	0,1272	0,1091	0,3392	0,2118	0,0737	0,1874	0,1851	0,1839	0,1316	0,1285	0,1387
	0.2	0,2362	0,2148	0,4324	0,2050	0,0603	0,1033	0,1023	0,1013	0,0776	0,0761	0,0818
	0.3	0,3045	0,2798	0,2388	0,4548	0,0346	0,0321	0,0320	0,0330	0,0253	0,0250	0,0277
30	0.1	0,1785	0,1533	0,3069	0,2618	0,0702	0,3474	0,3449	0,3457	0,1545	0,1530	0,1572
	0.2	0,2129	0,1841	0,4603	0,2327	0,0529	0,1384	0,1382	0,1379	0,0647	0,0641	0,0659
	0.3	0,2611	0,2282	0,1588	0,7518	0,0391	0,0387	0,0386	0,0387	0,0258	0,0258	0,0264
	$\alpha$	$\lambda =1$
5	0.1	0,0843	0,0660	0,1778	0,1854	0,0636	0,1368	0,1353	0,1361	0,1221	0,1213	0,1264
	0.2	0,0454	0,0343	0,1855	0,1658	0,0506	0,0774	0,0779	0,0788	0,0670	0,0666	0,0684
	0.3	0,0964	0,0818	0,1508	0,1287	0,0264	0,0287	0,0276	0,0274	0,0269	0,0269	0,0272
10	0.1	0,1151	0,0973	0,2771	0,2101	0,0716	0,1870	0,1859	0,1872	0,1325	0,1295	0,1403
	0.2	0,0892	0,0768	0,2694	0,1759	0,0531	0,0896	0,0872	0,0887	0,0641	0,0635	0,0682
	0.3	0,1494	0,1337	0,2122	0,1437	0,0328	0,0371	0,0360	0,0360	0,0305	0,0300	0,0318
30	0.1	0,1027	0,0803	0,1994	0,2635	0,0710	0,3423	0,3417	0,3420	0,1523	0,1509	0,1552
	0.2	0,1451	0,1233	0,1937	0,2297	0,0551	0,1453	0,1453	0,1455	0,0685	0,0679	0,0699
	0.3	0,1501	0,1279	0,1633	0,1997	0,0409	0,0390	0,0388	0,0389	0,0244	0,0243	0,0250

Table 25. Table 23: False detection rates with Exponential distribution α = 0 𝛼 0 \alpha=0 .

$p$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0,15070	0,13310	0,29900	0,26540	0,21340	0,16060	0,15850	0,16090	0,11820	0,11580	0,11280
10	0,18710	0,16890	0,36000	0,26580	0,38930	0,13638	0,15970	0,15910	0,13454	0,13920	0,13609
30	0,14292	0,11952	0,13558	0,27938	0,44708	0,13324	0,13266	0,13246	0,18204	0,18092	0,18464

Table 26. Table 24: Correct detection rates with Exponential distribution.

		$\delta =5$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,8578	0,8486	0,9602	0,9654	0,9975	1	1	1	0,9993	0,9993	0,9993
	0.3	0,1955	0,1664	0,8336	0,5792	0,8735	0,8935	0,8947	0,8938	0,8740	0,8698	0,8755
10	0.1	0,9983	0,9976	0,9974	0,9983	0,9975	0,9993	0,9993	0,9993	0,9995	0,9995	0,9999
	0.2	0,9671	0,9649	0,9893	0,9984	0,9923	0,9992	0,9992	0,9992	0,9984	0,9984	0,9986
	0.3	0,7821	0,7792	0,9847	0,9887	0,9832	0,9957	0,9958	0,9957	0,9951	0,9951	0,9958
30	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,9985	0,9985	0,9987	1	1	1	1	1	1	1	1
	0.3	0,7829	0,7550	1	1	1	1	1	1	1	1	1
$p$	$\alpha$	$\delta =10$
5	0.1	0,9991	0,9970	1	1	0,9991	0,9991	0,9991	0,9991	1	1	1
	0.2	0,9770	0,9733	0,9991	0,9991	0,9927	0,9961	0,9961	0,9961	0,9952	0,9952	0,9997
	0.3	0,8185	0,7952	0,9900	0,9966	0,9853	0,9916	0,9919	0,9919	0,9906	0,9903	0,9911
10	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,9756	0,9756	1	1	1	1	1	1	1	1	1
	0.3	0,7905	0,7902	1	1	1	1	1	1	1	1	1
30	0.1	1	1	1	1	1	1	1	1	1	1	1
	0.2	0,9987	0,9987	1	1	1	1	1	1	1	1	1
	0.3	0,7719	0,7601	1	1	1	1	1	1	1	1	1

Table 27. Table 25: False detection rates with Exponential distribution.

		$\delta =5$
$p$	$\alpha$	MCD	Adj MCD	Kurtosis	OGK	COM	RMDv1	RMDv2	RMDv3	RMDv4	RMDv5	RMDv6
5	0.1	0,0339	0,0198	0,0341	0,0636	0,0033	0,0084	0,0073	0,0069	0,0075	0,0071	0,0077
	0.2	0,0110	0,0055	0,0347	0,0507	0,0017	0,0014	0,0012	0,0012	0,0010	0,0009	0,0012
	0.3	0,0374	0,0260	0,0379	0,0383	0,0001	0,0000	0,0000	0,0000	0,0000	0,0000	0,0000
10	0.1	0,1194	0,1018	0,3285	0,2299	0,0683	0,2380	0,2362	0,2372	0,2302	0,2258	0,2442
	0.2	0,0315	0,0235	0,2634	0,1838	0,0466	0,1176	0,1183	0,1182	0,1338	0,1302	0,1479
	0.3	0,0021	0,0016	0,1890	0,1445	0,0258	0,0575	0,0576	0,0576	0,0722	0,0681	0,0803
30	0.1	0,0887	0,0639	0,1276	0,2371	0,0335	0,4021	0,4020	0,4015	0,2726	0,2697	0,2781
	0.2	0,0239	0,0099	0,1390	0,2022	0,0202	0,1751	0,1750	0,1754	0,1531	0,1504	0,1583
	0.3	0,0001	0,0000	0,1196	0,1734	0,0111	0,0525	0,0528	0,0524	0,0551	0,0541	0,0580
$p$	$\alpha$	$\delta =10$
5	0.1	0,0842	0,0679	0,2846	0,2198	0,0799	0,1567	0,1572	0,1557	0,1659	0,1631	0,1708
	0.2	0,0271	0,0176	0,2408	0,1824	0,0551	0,0743	0,0752	0,0754	0,1015	0,0994	0,1039
	0.3	0,0011	0,0006	0,1969	0,1451	0,0288	0,0302	0,0299	0,0289	0,0460	0,0443	0,0470
10	0.1	0,1148	0,0961	0,3201	0,2284	0,0620	0,2103	0,2121	0,2131	0,2015	0,1963	0,2161
	0.2	0,0357	0,0261	0,2647	0,1841	0,0404	0,0895	0,0906	0,0890	0,0942	0,0912	0,1032
	0.3	0,0028	0,0022	0,1873	0,1474	0,0240	0,0298	0,0293	0,0309	0,0411	0,0388	0,0463
30	0.1	0,0863	0,0621	0,1343	0,2416	0,0334	0,3536	0,3530	0,3536	0,2394	0,2371	0,2450
	0.2	0,0253	0,0116	0,1179	0,2043	0,0213	0,1053	0,1056	0,1058	0,0913	0,0899	0,0944
	0.3	0,0000	0,0000	0,1240	0,1662	0,0108	0,0127	0,0128	0,0130	0,0146	0,0144	0,0156

Equations

$$M{D}_i=(({\mathbf{x}}_i-\hat{\mathbf{\mu }}){\hat{\Sigma }}^{-1}({\mathbf{x}}_i-\hat{\mathbf{\mu }}{)}^T{)}^{1/2},$$

$$RM{D}_i:=(({\mathbf{x}}_i-{\hat{\mathbf{\mu }}}_R){\hat{\Sigma }}_R^{-1}({\mathbf{x}}_i-{\hat{\mathbf{\mu }}}_R{)}^T{)}^{1/2},$$

$${\hat{E}}_{Sh}=(1-\eta )\hat{E}+\eta \hat{T}$$

$${\hat{\mathbf{\mu }}}_{CCM}=(median({\mathbf{x}}_{\cdot 1}),...,median({\mathbf{x}}_{\cdot p})),$$

$${\hat{\mathbf{\mu }}}_{MM}=argmi{n}_{{\mathbf{x}}_m,m\in \{1,...,n\}}\frac{1}{n}\sum _{i=1}^n||{\mathbf{x}}_m-{\mathbf{x}}_i|{|}_1$$

$${\hat{\mathbf{\mu }}}_{Sh(CCM)}=(1-\eta ){\hat{\mathbf{\mu }}}_{CCM}+\eta {\nu }_{\mathbf{\mu }}\mathbf{e}$$

$$\begin{array}{cc}{\text{min}}_{{\nu }_{\mathbf{\mu }},\eta }& E\left[{\Vert {\hat{\mathbf{\mu }}}_{Sh(CCM)}-\mathbf{\mu }\Vert }_2^2\right]\\ \text{s.t.}& {\hat{\mathbf{\mu }}}_{Sh(CCM)}=(1-\eta ){\hat{\mathbf{\mu }}}_{CCM}+\eta {\nu }_{\mathbf{\mu }}\mathbf{e},\end{array}$$

$${\nu }_{\mathbf{\mu }}$$

$${\sigma }_{{{\hat{\mathbf{\mu }}}_{CCM}}_j}^2=Var({{\hat{\mathbf{\mu }}}_{CCM}}_j)=\frac{\pi }{2n}{\sigma }_{{\mathbf{x}}_{\cdot \mathbf{j}}}^2$$

$$E\left[{\Vert {\hat{\mathbf{\mu }}}_{CCM}-\mathbf{\mu }\Vert }_2^2\right]$$

$${\hat{\mathbf{\mu }}}_{Sh(MM)}=(1-\eta ){\hat{\mathbf{\mu }}}_{MM}+\eta {\nu }_{\mathbf{\mu }}\mathbf{e}$$

$$\begin{array}{cc}{\text{min}}_{{\nu }_{\mathbf{\mu }},\eta }& E\left[{\Vert {\hat{\mathbf{\mu }}}_{Sh(MM)}-\mathbf{\mu }\Vert }_2^2\right]\\ \text{s.t.}& {\hat{\mathbf{\mu }}}_{Sh(MM)}=(1-\eta ){\hat{\mathbf{\mu }}}_{MM}+\eta {\nu }_{\mathbf{\mu }}\mathbf{e},\end{array}$$

$${\nu }_{\mathbf{\mu }}$$

$$E\left[{\Vert {\hat{\mathbf{\mu }}}_{MM}-\mathbf{\mu }\Vert }_2^2\right]=E\left[\sum _{j=1}^p({{\hat{\mathbf{\mu }}}_{MM}}_j-{\mathbf{\mu }}_j{)}^2\right]=\sum _{j=1}^p{\sigma }_{{{\hat{\mathbf{\mu }}}_{MM}}_j}^2$$

$${\hat{\mathbf{\mu }}}_{MM}\sim N_p\left(\mathbf{\mu },\frac{1}{n}{\hat{A}}^{-1}\hat{B}{\hat{A}}^{-1}\right),$$

$$E\left[{\Vert {\hat{\mathbf{\mu }}}_{MM}-\mathbf{\mu }\Vert }_2^2\right]=trace\left(\frac{1}{n}{\hat{A}}^{-1}\hat{B}{\hat{A}}^{-1}\right)$$

$$MAD(X)=median(|X-median(X)|),$$

$$MAD(X)={\sigma }_X{\Phi }^{-1}(3/4),$$

$${\sigma }_X^2=2.198\cdot MA{D}^2(X)$$

$$COM(X,Y)=med((X-med(X))(Y-med(Y)))$$

$$COM(\mathbf{x})=(COM({\mathbf{x}}_{\cdot j},{\mathbf{x}}_{\cdot t}))j,t=1,...,p$$

$${\hat{S}}_{CCM}=2.198\cdot COM(\mathbf{x})$$

$${\hat{\Sigma }}_{Sh}=(1-\eta )\hat{E}+\eta \hat{T},$$

$$trace({\hat{S}}_{CCM})$$

$${\hat{\Sigma }}_{Sh(CCM)}=(1-\eta ){\hat{S}}_{CCM}+\eta {\nu }_{\Sigma }I$$

$$\begin{array}{cc}{\text{min}}_{{\nu }_{\Sigma },\eta }& E\left[{\Vert {\hat{\Sigma }}_{Sh}-\Sigma \Vert }^2\right]\\ \text{s.t.}& {\hat{\Sigma }}_{Sh}=(1-\eta ){\hat{S}}_{CCM}+\eta {\nu }_{\Sigma }I,\end{array}$$

$${\nu }_{\Sigma }$$

$$m_A=m(X_A)=Am(X),S_A=S(X_A)=S(X)A^T$$

$$P=\left[\begin{array}{cc}{P}_1& 0\\ 0& P_2\end{array}\right]$$

$$P_1=\left[\begin{array}{ccc}1& 0.95& 0.3\\ 0.95& 1& 0.1\\ 0.3& 0.1& 1\end{array}\right],P_2=\left[\begin{array}{ccc}1& -0.499& -0.499\\ -0.499& 1& -0.499\\ -0.499& -0.499& 1\end{array}\right]$$

Full text

Keywords: multivariate distance, robust location and covariance matrix estimation, comedian matrix, multivariate $L_{1}$-median.

1 Introduction

The detection of outliers in multivariate data is an important task in Statistics, since that kind of data can distort any statistical procedure (Tarr et al. [2016]). The task of detecting multivariate outliers can be useful in various fields such as quality control, medicine, finance, image analysis, and chemistry (Vargas N [2003], Brettschneider et al. [2008], Hubert et al. [2008], Hubert and Debruyne [2010], Perrotta and Torti [2010] and Choi et al. [2016]). The concept of outlier is not well-defined in the literature since different authors tend to have varying definitions. Although they are generally defined as observations resulting from a secondary process, which differ from the background distribution. Thus, the outliers will differ from the main bulk of the data. They do not need to be especially high or low for all the variables in the data set, that is why the task of identifying multivariate outliers with the classical univariate methods commonly fail. In the multivariate sense, there must be considered both the distance of an observation from the centroid of the data, and the shape of the data. The Mahalanobis distance [Mahalanobis, 1936] is a well-known measure which takes it into account. For multivariate gaussian data, the distribution of the squared Mahalanobis distance, $MD^{2}$, is known [Gnanadesikan and Kettenring, 1972] to be chi-squared with $p$ (the dimension of the data, the number of variables) degrees of freedom, i.e. $\chi_{p}^{2}$. Then, the adopted rule for identifying the outliers is selecting the threshold as the 0.975 quantile of the $\chi_{p}^{2}$.

However, outliers need not necessarily have large MD values (Masking problem) and not all observations with large MD values are necessarily outliers (Swamping problem) [Hadi, 1992]. The problems of masking and swamping arise due to the influence of outliers on classical location and scatter estimates (sample mean and sample covariance matrix), which implies that the estimated distance will not be robust to outliers. The solution is to consider robust estimators of centrality and covariance matrix to obtain a robust Mahalanobis distance (RMD). Many robust estimators for location and covariance have been introduced in the literature [Maronna and Yohai, 1976]. Rousseeuw [1985] proposed the Minimum Covariance Determinant (MCD) estimator based on the computation of the ellipsoid with the smallest volume or with the smallest covariance determinant that would encompass at least half of the data points. The procedure required naive subsampling for minimizing the objective function of the MCD, but an improvement much more effective, the Fast-MCD, was introduced by Rousseeuw and Driessen [1999] and a code is available in MATLAB (Verboven and Hubert [2005]). Unfortunately Fast-MCD still requires substantial running times for large $p$, because the number of candidate solutions grows exponentially with the dimension $p$ of the sample and, as a consequence, the procedure becomes computationally expensive for even moderately sized problems.

On the other hand, the squared RMD distributional fit usually breaks down, i.e. it does not necessarily have to follow a chi-squared distribution when you deviate from the gaussian distribution. Thus, determining exact cutoff values for outlying distances continues to be a difficult problem and it has found much attention because no universally applicable method has been proposed. Despite this fact, the $\chi^{2}_{p;0.975}$ quantile is often considered as threshold for recognizing outliers in the robust distance case, but this approach may have some drawbacks. Evidence of this behavior is now well documented even in moderately large samples, especially when the number of variables increases (Becker and Gather [1999], Hardin and Rocke [2005], Cerioli et al. [2009] and Cerioli et al. [2008]). It is crucial to determine the threshold for the distances in order to decide whether an observation is an outlier. Filzmoser et al. [2005] proposed to use an adjusted quantile, instead of the classical choice of the $\chi^{2}_{p;0.975}$ quantile. The adjusted threshold is estimated adaptively from the data, but their proposal is defined for a specific robust Mahalanobis distance, the one based on the MCD estimator. Let us call this method Adj MCD. Peña and Prieto [2001] and Peña and Prieto [2007] proposed an algorithm called Kurtosis, based on the analysis of the projections of the sample points onto a certain set of directions obtained by maximizing and minimizing the kurtosis coefficient of the projections, and some random directions generated by a stratified sampling scheme. With the combination of random and specific directions, the authors proposed a powerful procedure for robust estimation and outlier detection. However, this procedure has some drawbacks when the dimension $p$ of the sample space grows, and in presence of correlation between the variables, the method looses power [Marcano and Fermín, 2013]. Maronna and Zamar [2002] proposed the Orthogonalized Gnanadesikan-Kettenring (OGK) estimator. It was the result of applying a general method to the pairwise robust scatter matrix from Gnanadesikan and Kettenring [1972], in order to obtain a positive-definite scatter matrix. On the other hand, a reweighing step can be used to identify outliers, where atypical observations get weight 0 and normal observations get weight 1. Sajesh and Srinivasan [2012] proposed the Comedian method (COM) to detect outliers from multivariate data based on the comedian matrix estimator from Falk [1997]. The method is found to be efficient under various simulation scenarios and suitable in high-dimensional data. Furthermore, there are several real scenarios where the number of variables is high in which outlier detection is very important. For example, medical imaging datasets often contain deviant observations due to pre-processing artifacts or large intrinsic inter-subject variability (Lazar [2008], Lindquist [2008], Monti [2011], Poline and Brett [2012]), in biological and financial studies (Chen et al. [2010] and Zeng et al. [2015]), and also in geochemical data, because of their complex nature (Reimann and Filzmoser [2000], Templ et al. [2008]).

In this article, a collection of RMD’s are proposed for outlier detection especially in high dimension. They are based on considering different combinations of robust estimators of location and covariance matrix. Two basic options are considered for the location parameter: a component-wise median and the $L_{1}$ multivariate median (Gower [1974], Brown [1983], Dodge [1987], Small [1990]). A notion called shrinkage estimator (Ledoit and Wolf [2003a], Ledoit and Wolf [2003b], Ledoit and Wolf [2004], DeMiguel et al. [2013], Gao [2016], citesun2018portfolio, Steland [2018]) is considered, which is aimed to reduce estimation error. The shrinkage is applied to both of the previous mentioned location estimators. As for the covariance matrix, the options basically consists on a shrinkage estimator over special cases of comedian matrices (Hall and Welsh [1985], Falk [1997]), which are based on a location parameter that will be estimated using a robust estimator of centrality in a way that a RMD can be obtained with meaningful combinations of both location and covariance matrix estimators. Simulation results demonstrates the satisfactory practical performance of our proposal, especially when the number of variables grows. The computational cost is studied by both simulations and a real dataset example.

The paper is organized as follows. Section 2 describes the shrinkage estimators both for the location and the covariance matrix, and the proposed combinations of these estimators in order to define a RMD. Section 3 shows a simulation study with contaminated multivariate gaussian data and when we deviate from the gaussian assumption, e.g. with skewed or heavy-tailed data, to compare the proposal with the other robust approaches: MCD, Adj MCD, Kurtosis, OGK and COM. To investigate the properties of affine equivariance and breakdown value, other simulation scenarios are proposed with correlated data, transformed data and large contaminated data. Section 4 shows the behavior with a real dataset example. Finally, Section 5 provides some conclusions.

2 A robust Mahalanobis distance based on shrinkage estimator

The classical Mahalanobis distance is defined for every $p-$dimensional observation $\mathbf{x}_{i}$ of the multivariate sample $\{\mathbf{x}_{1},...,\mathbf{x}_{n}\}$, as:

[TABLE]

where $\hat{\boldsymbol{\mu}}$ is the estimated multivariate location (sample mean) and $\hat{\Sigma}$ is the estimated covariance matrix (sample covariance matrix).

Since the problem with this definition is that the classical estimates of location and covariance matrix are often highly influenced by the presence of outliers [Rousseeuw and Van Zomeren, 1990], the solution is to consider robust estimates of centrality and covariance matrix, i.e. resistant against the influence of outlying observations, giving rise to a robust Mahalanobis distance, defined as:

[TABLE]

where $\hat{\boldsymbol{\mu}}_{R}$ and $\hat{\Sigma}_{R}$ are robust estimators of centrality and covariance matrix, respectively.

We propose to use a notion which is frequently used in finance and portfolio optimization, known as shrinkage (Equation 3). It is widely used in those fields because its good performance for “large $p$ small $n$” problems (see Couillet and McKay [2014], Chen et al. [2011] and Steland [2018]), although we focus on data with $n>p$. This estimator $\hat{E}_{Sh}$ relies on the fact that “shrinking” an estimator $\hat{E}$ of a parameter $\theta$ towards a target estimator $\hat{T}$, would help to reduce the estimation error, because although the shrinkage target is usually biased, it also contains less variance than the estimator $\hat{E}$. Therefore, under general conditions, there exists a shrinkage intensity $\eta$, so the resulting shrinkage estimator would contain less estimation error than $\hat{E}$ [James and Stein, 1961].

[TABLE]

The main advantage of using a shrinkage estimator is to obtain a trade-off between bias and variance. This approach can be applied to estimate both the location and dispersion parameters obtaining different meaningful combinations to define robust Mahalanobis distances. In the case of covariance matrices shrinkage has the additional advantage that it is always positive definite and well conditioned.

2.1 Location parameter

Let $\textbf{x}=\{\mathbf{x}_{\cdot 1},...,\mathbf{x}_{\cdot p}\}$ be the $n\times p$ data matrix with $n$ being the sample size and $p$ the number of variables. Based on the fact that the $median$ is a better choice in terms of robustness, we start by considering as a location estimator the component-wise median:

[TABLE]

where $median$ denotes the univariate median and $(\mathbf{x}_{\cdot j})=(x_{1j},...,x_{nj})^{T}$ for all $j=1,...,p$ is the $j$-th column of $\mathbf{x}$.

Another option is to consider a multivariate median $\hat{\boldsymbol{\mu}}_{MM}$ called $L_{1}-$median which is a robust and highly efficient estimator of central tendency (Lopuhaa and Rousseeuw [1991], Vardi and Zhang [2000], Oja [2010]). It is defined as:

[TABLE]

DeMiguel et al. [2013] proposed a shrinkage estimator over the sample mean, towards a scaled vector of ones as the target. In the same way we propose to study shrinkage estimators for both (4) and (5). Consider $\nu_{\boldsymbol{\mu}}\mathbf{e}$ as the target estimator $\hat{T}$ in (3), where $\mathbf{e}$ is the $p-$dimensional vector of ones, and consider $\hat{\boldsymbol{\mu}}_{CCM}$ as the sample estimator $\hat{E}$. Then, the shrinkage estimator over the component-wise median is:

[TABLE]

The scaling factor $\nu_{\boldsymbol{\mu}}$ and the intensity $\eta$ should minimize the expected quadratic loss, that is:

[TABLE]

where $\left\|\mathbf{x}\right\|_{2}^{2}=\sum_{j=1}^{p}x_{j}^{2}$.

Proposition 1

The solution of the problem in (7) is:

[TABLE]

See the proof in Section 1 from the Supplementary Material. Note that the denominator in the above expression (8) has no problem when estimating, but the numerator is not straightforward because $\boldsymbol{\mu}$ is unknown. Then, it is necessary to provide another expression for the numerator. Chu [1955] investigated the distribution for the sample median estimator and obtained the following result about the variance in presence of normality. Fix $j$, for $j\in\{1,...,p\}$:

[TABLE]

Therefore, the numerator in the expression (8) for determining the $\eta$ in Proposition 1 is:

[TABLE]

We need to estimate $\sigma_{\mathbf{x}_{\cdot j}}^{2}$ robustly and we will do so as explained in the next sub-section with property (23).

On the other hand, consider $\nu_{\boldsymbol{\mu}}\mathbf{e}$ again as the target estimator $\hat{T}$ and consider $\hat{\boldsymbol{\mu}}_{MM}$ as the sample estimator $\hat{E}$, in (3). Then, the shrinkage estimator over the multivariate $L_{1}-$median is:

[TABLE]

The scaling factor $\nu_{\boldsymbol{\mu}}$ and the intensity $\eta$ should minimize the expected quadratic loss:

[TABLE]

where $\left\|\mathbf{x}\right\|_{2}^{2}=\sum_{j=1}^{p}x_{j}^{2}$.

Proposition 2

The solution of the problem in (12) is:

[TABLE]

The proof is also in Section 1 from the Supplementary Material. As in the previous case, the denominator in the $\eta$ expression (13) can be described as:

[TABLE]

Bose and Chaudhuri [1993], Bose [1995] and Möttönen et al. [2010] investigated the asymptotic distribution for the $L_{1}-$median, and they obtained the following result in presence of normality:

[TABLE]

where $\hat{A}(\mathbf{x}_{i})=\frac{1}{||\mathbf{x}_{i}||_{2}}\left(I_{p}-\frac{\mathbf{x}_{i}\mathbf{x}_{i}^{T}}{||\mathbf{x}_{i}||^{2}_{2}}\right)$ and $\hat{B}(\mathbf{x}_{i})=\frac{\mathbf{x}_{i}\mathbf{x}_{i}^{T}}{||\mathbf{x}_{i}||_{2}^{2}}$, with $\mathbf{x}_{i}\in\mathbb{R}^{p}$, for each $i=1,...,n$.

The numerator in the expression (13) can be obtained with the above property:

[TABLE]

2.2 Dispersion parameter

Based on the median concept, which is a robust measure of location, one can define a robust measure of dispersion for a random variable $X$, which is the Median Absolute Deviation (MAD) from the data’s median:

[TABLE]

Falk [1997] showed the following relation, assuming normality, between the $MAD$ and the standard deviation $\sigma_{X}$:

[TABLE]

where $\Phi$ denotes the standard normal cdf. Taking the square in (17) we obtain a relation between the variance $\sigma_{X}^{2}$ and $MAD^{2}(X)$:

[TABLE]

Extending the idea of the $MAD$, a robust measure of dependence between two random variables $X$ and $Y$ is the comedian (Falk [1997]):

[TABLE]

The comedian generalizes the MAD, because $COM(X,X)=MAD^{2}(X)$, and also has the highest possible breakdown point (Falk [1997]). An important fact is that the comedian parallels the covariance, but the latter requires the existence of the first two moments of the two random variables, whereas the comedian always exists. Other known properties of the comedian are that it is symmetric, location invariant and scale equivariant. Furthermore, Hall and Welsh [1985] discussed about the strong consistency and asymptotic normality of the MAD, and Falk [1997] established similar results for the comedian.

Finally, a comedian matrix can be defined based on a multivariate version of (19). Let $\textbf{x}=\{\mathbf{x}_{\cdot 1},...,\mathbf{x}_{\cdot p}\}$ be the $n\times p$ data matrix with $n$ being the sample size and $p$ the number of variables. Then the comedian matrix is defined as:

[TABLE]

Note that from relation described in (18), one can consider the adjusted comedian:

[TABLE]

Note that $\hat{S}_{CCM}$ is a robust alternative for the covariance matrix, but in general it is not positive (semi-) definite (see Falk [1997]). Since we need this property for inverting the covariance matrix in a Mahalanobis distance, we propose a shrinkage over $\hat{S}_{CCM}$, because of its advantage of providing always a positive definite and well-conditioned matrix. Therefore, if a shrinkage estimator is considered in (3) for the dispersion parameter:

[TABLE]

we propose to use in (22), the estimator $\hat{E}=\hat{S}_{CCM}$.

Recall the previous sub-section 2.1 in which we needed to provide a robust estimator for $\sigma_{\mathbf{x}_{\cdot j}}^{2}$, for each $j=1,...,p$ (Equation 10) note that, because of the relation in (18):

[TABLE]

Thus, when considering a shrinkage estimator of the component-wise median, in order to estimate the variance of $\boldsymbol{\hat{\mu}}_{CCM}$ needed in the expression (8) for the shrinkage intensity $\eta$, and according to the relation (10), we propose to estimate $\sum_{j=1}^{p}\sigma^{2}_{\mathbf{x}_{\cdot j}}$ using the $trace(\hat{S}_{CCM})$.

About the shrinkage target $\hat{T}$, several choices have been proposed in the literature. For example, Ledoit and Wolf [2003b] proposed a weighted average of the sample covariance matrix and a single-index covariance matrix. Ledoit and Wolf [2003a] proposed selecting the shrinkage target as a “constant correlation matrix”, whose correlations are set equal to the average of all sample correlations. Finally, Ledoit and Wolf [2004] proposed to use a multiple of the identity matrix as the shrinkage target. The authors proved that the resulting shrinkage covariance matrix is well-conditioned, even if the sample covariance matrix is not. There is also another approach introduced by DeMiguel et al. [2013]. The authors proposed a shrinkage estimator both for the covariance matrix and its inverse. The estimators were constructed as a convex combination of the sample covariance matrix or its inverse, respectively, and a scaled shrinkage target, which they consider the scaled identity matrix as Ledoit and Wolf [2004]. Therefore, we propose to use as shrinkage target $\hat{T}=\nu_{\Sigma}I$. Thus (22) results in:

[TABLE]

Lastingly, the scaling parameter $\nu_{\Sigma}$ and the shrinkage intensity parameter $\eta$ in (24) need to be estimated. They both are chosen to minimize the expected quadratic loss as in Ledoit and Wolf [2004]:

[TABLE]

where $\left\|A\right\|^{2}=trace(AA^{T})/p$.

Proposition 3

The solution of the problem (25) is:

[TABLE]

The proof can be found in Section 1 from the Supplementary Material. In practice, we propose to estimate the numerator of the expression for $\eta$ as Ledoit and Wolf [2003a], Ledoit and Wolf [2003b] and Ledoit and Wolf [2004], but considering $\hat{S}_{CCM}$ instead of the sample covariance matrix, as the estimator of $\Sigma$.

Note that the comedian matrix depends on centered data considering the component-wise median $\boldsymbol{\hat{\mu}}_{CCM}$. A special case of comedian matrix can be defined if the data are centered using a different location estimator. We propose to center the data using the other location estimators described in Subsection 2.1, i.e. the multivariate $L_{1}-$median $\boldsymbol{\hat{\mu}}_{MM}$, and the shrinkage estimators $\boldsymbol{\hat{\mu}}_{Sh(CCM)}$ and $\boldsymbol{\hat{\mu}}_{Sh(MM)}$. We will consider shrinkages over those special comedian matrices.

$\hat{\Sigma}_{Sh(MM)}=(1-\eta)\hat{S}_{MM}+\eta\nu_{\Sigma}I$, with for $j,t=1,...,p$:

$\hat{S}_{MM}=2.198\cdot COM_{MM}(\mathbf{x})=2.198\cdot(med((\mathbf{x}_{\cdot j}-(\boldsymbol{\hat{\mu}}_{MM})_{j})(\mathbf{x}_{\cdot t}-(\boldsymbol{\hat{\mu}}_{MM})_{t}))$ 2. 2.

$\hat{\Sigma}_{Sh(Sh(CCM))}=(1-\eta)\hat{S}_{Sh(CCM)}+\eta\nu_{\Sigma}I$, with for $j,t=1,...,p$:

$\hat{S}_{Sh(CCM)}=2.198\cdot COM_{Sh(CCM)}(\mathbf{x})=2.198\cdot(med((\mathbf{x}_{\cdot j}-(\boldsymbol{\hat{\mu}}_{Sh(CCM)})_{j})(\mathbf{x}_{\cdot t}-(\boldsymbol{\hat{\mu}}_{Sh(CCM)})_{t}))$ 3. 3.

$\hat{\Sigma}_{Sh(Sh(MM))}=(1-\eta)\hat{S}_{Sh(MM)}+\eta\nu_{\Sigma}I$, with for $j,t=1,...,p$:

$\hat{S}_{Sh(MM)}=2.198\cdot COM_{Sh(MM)}(\mathbf{x})=2.198\cdot(med((\mathbf{x}_{\cdot j}-(\boldsymbol{\hat{\mu}}_{Sh(MM)})_{j})(\mathbf{x}_{\cdot t}-(\boldsymbol{\hat{\mu}}_{Sh(MM)})_{t}))$

The optimal expression for the parameters $\eta$ and $\nu_{\Sigma}$ in the above cases is analogous to the Proposition 3, but considering in each case the sample estimator as the corresponding special comedian matrix.

2.3 Proposed Robust Mahalanobis Distances

A robust Mahalanobis distance can be defined as in (2), for each of the following 6 possible combinations for the location and the dispersion estimators (see Table 1). Note that they are meaningful combinations because the shrinkage estimator of dispersion is made upon a special comedian matrix closely based on the location estimator jointly considered for defining the RMD.

For all our proposed combinations, the threshold considered to detect the outliers is the $\chi^{2}_{p;0.975}$ quantile.

3 Simulation results

3.1 Normal distribution

A simulation study is performed considering a $p-$dimensional random variable $X$ following a contaminated multivariate normal distribution given as a mixture of normals of the form $(1-\alpha)N(\mathbf{0},I)+\alpha N(\delta\mathbf{e},\lambda I)$, where $\mathbf{e}$ denotes the $p-$dimensional vector of ones. This model is analogous to the one used by Rousseeuw and Driessen [1999], Peña and Prieto [2001], Filzmoser et al. [2005], Peña and Prieto [2007], Maronna and Zamar [2002] and Sajesh and Srinivasan [2012]. This experiment has been conducted for different values of the sample-space dimension $p=5,10,30$, and the chosen sample size in relation to the dimension was $n=100,100,500$, respectively. The contamination levels were $\alpha=0,0.1,0.2,0.3$, the distance of the outliers $\delta=5,10$ and the concentration of the contamination $\lambda=0.1,1$. For each set of values, 100 random sample repetitions have been generated.

For the methods mentioned in previous sections some measures are studied: the correct detection rates (c) and the false detection rates (f). The method MCD refers to the RMD based on the MCD estimator and with the classical threshold, the method Adj MCD refers to the latter distance considering the adjusted quantile of Filzmoser et al. [2005], the method Kurtosis refers to the Peña and Prieto [2007] approach, the method OGK refers to the Orthogonalized Gnanadesikan-Kettenring method proposed by Maronna and Zamar [2002] and COM is the Comedian method proposed by Sajesh and Srinivasan [2012]. We have also presented the results for the collection RMDv1-RMDv6 proposed in Table 1. All simulations were performed in Matlab. Section 2 in the Supplementary Material shows the tables corresponding to all simulation scenarios. Here we show only the most significant and representative results.

Nevertheless, the tables show general outcomes. For example, Adj MCD, actually improves MCD with respect to the “f”, lowering it, and in most cases maintaining the same “c”. Although, in other cases it also slightly lowers the “c”. On the other hand, the false detection rates in case of no contamination are sufficiently low for all methods, but our proposed collection shows the lowest values especially in high dimension, actually here the best performance is observed for RMDv6. With certain percent of contamination, the worst behavior of our proposed methods is when dimension is low and the highest percentage of outliers are considered to be near the center of the data. This matter can be seen in Table 11 which corresponds to the correct detection rates. In this case Kurtosis has better performance. Although, this only happens in two cases, and in all other cases MCD, Adj MCD, Kurtosis and OGK are the ones with the worst behavior. Meanwhile, COM is a good competitor, but the best overall performance is made by RMDv6 especially in high dimension.

Another situation is when outliers are far from the center of the data, i.e. $\delta=10$. This scenario is shown in Table 12. It is clear that our proposed methods lead to the best performance, achieving $100\%$ of correct detection rate, for all dimension and percentage of contamination considered. Other tables about the “c” can be found in the Supplementary Material, as well as the false detection rate tables, which show that in the vast majority of cases our proposal have an “f” value equal to zero and when not, a value very close to zero, which is what is desirable.

3.2 $t_{3}$-distribution

In order to check the behavior of the methods when the distribution deviates from normality, a simulation study is performed considering a $p-$dimensional random variable $X$ following a contaminated multivariate $t$-distribution with $3$ degrees of freedom of the form $(1-\alpha)T_{3}(\mathbf{0},I)+\alpha T_{3}(\delta\mathbf{e},\lambda I)$. The first parameter of the notation of $T_{3}(\cdot,\cdot)$ refers to the mean and the second one to the covariance matrix. The parameters for the contamination are the same considered above and the same measures “c” and “f” are studied. All the results can be founded in the tables from Section 2 in the Supplementary Material. It should be noted the unsatisfactory behavior of our competitors with respect to the “c” especially in high dimension or with large contamination level, meanwhile in most cases we attain a $100\%$ correct detection rate. With respect to the “f” value, all methods show non-zero “f” values, and the best performance is showed by COM and our proposed methods.

3.3 Exponential distribution

We considered also a $p-$dimensional random variable $X$ following a contaminated multivariate exponential distribution given as a mixture $(1-\alpha)Exp(\mathbf{0})+\alpha Exp(\delta\mathbf{e})$. The parameter of the notation $Exp(\cdot)$ refers to the mean. This case is analogous to the previous ones, with the difference that only the schemes associated with the distance of the outliers are considered. The tables with the results can be found in Section 2 in the Supplementary Material, and it can be seen that our methods achieve $100\%$ of correct detection rate in the majority of cases, except in some cases when dimension is low. Actually in all cases the best performance is showed by RMDv6. On the other hand, all methods show more or less the same “f” value.

3.4 Summary and selection of one of our proposed distances

In the simulation study, for each contamination scheme we have also calculated a measure called F-score (Goutte and Gaussier [2005], Sokolova et al. [2006], Powers [2011]), often used in Engineering, which is a measure of a test’s accuracy. Its expression is F-score$=2PR/(P+R)$, where $P$ is called precision and $R$ is known as the recall. The precision $P$ is the number of correct detected outliers divided by the total number of detected outliers, and the recall $R$ is the number of correct detected outliers divided by the real total number of outliers. Thus, this measure provides a trade-off between the two desired outcomes: a high rate of correctly identified outliers and a low rate of observations mislabel as outliers. The results are not included in the paper for avoiding large extension, but the method with the overall classification between the top 3 best positions ranking with respect to the F-score, is method RMDv6.

It is clear the out-performance of our proposed methods with gaussian data, especially in high dimension and even when we deviate from the normality assumption, for example when considering skewed and heavy-tailed distributions like the multivariate $t_{3}$-distribution and the multivariate exponential distribution. From all of our six proposed robust distances, the one that shows the best results in the vast majority of cases is RMDv6. Thus, we decided to select it as the best one in the matter of performance, and from now on we will refer to it as RMD-Shrinkage. Now we will proceed to study some properties like the behavior under correlated data, the affine equivariance, the breakdown value, and the computational time.

3.5 Correlation and affine equivariance

Consider $X=\{\mathbf{x_{1}},...,\mathbf{x_{n}}\}$ and a pair of multivariate location and covariance estimators $(m,S)$. In general, these estimators are called affine equivariant if for any nonsingular matrix $A$ it holds that:

[TABLE]

The affine transformation of $X$ is $X_{A}=\{A\mathbf{x_{1}},...,A\mathbf{x_{n}}\}$. Affine equivariance implies that the estimator transforms well under any nonsingular reparametrization of the space of the $\mathbf{x}_{i}$. The data might for instance be rotated, translated or rescaled (for example through a change of the measurement units).

The method RMD-Shrinkage is ultimately based on not affine equivariant estimators which are the $L_{1}$-median (Lopuhaa and Rousseeuw [1991]) and the comedian matrix. However, the $L_{1}$-median is orthogonal equivariant, i.e. it satisfies Equation (26) with $A$ any orthogonal matrix $(A^{\prime}=A^{-1})$. This implies that the $L_{1}$-median transforms appropriately under all transformations that preserve Euclidean distances (such as translations, rotations and reflections). About the comedian matrix, which always exists, it is symmetric, location invariant and scale equivariant (Falk [1997]), i.e. $COM(X,\mathbf{a}Y+\mathbf{b})=\mathbf{a}COM(X,Y)=\mathbf{a}COM(Y,X)$. Since the proposed method is not affine equivariant, it is important to investigate the behavior under correlated data. Devlin et al. [1981] used a correlation matrix $P$ for generating Monte Carlo data from different distributions of moderate dimension $p=6$. The matrix has the form:

[TABLE]

[TABLE]

The reason for the selection of the matrix $P$ is because the dimension is large enough to study multivariate estimators and the range of correlation values is large. This way the differences in the abilities of the methods to detect outliers from highly correlated data can be observed. For the simulations, $n=100$ observations were generated from a mixture of Normals $(1-\alpha)N(\mathbf{0},P)+\alpha N(5\mathbf{e},P)$. The contamination level $\alpha=10\%,20\%,30\%$.

Table 4 shows that the correct and false detection rates of MCD, Adj MCD, Kurtosis and OGK are worse than that of our proposal . On the other hand, COM show more or less the same behavior in case of $10\%$ and $20\%$ of contamination, and slightly worse than our proposal when the contamination level increases to $30\%$. The methods MCD, Adj MCD and Kurtosis are affine equivariant, while OGK and COM are not. Hence, the proposed procedure RMD-Shrinkage is more efficient than other affine and not affine equivariant methods in case of correlated datasets. Also the false detection rate is very low even in this case of presence of correlation.

Affine equivariance of the estimators is equivalent to say that the robust Mahalanobis distance is affine invariant:

[TABLE]

Maronna and Zamar [2002] and Sajesh and Srinivasan [2012] proposed to investigate the lack of equivariance with transformed data, by simulations. We study the same for our proposal. They propose to generate random matrices as $A=TD$, where $T$ is a random orthogonal matrix and $D=diag(u_{1},...,u_{p})$, where the $u_{j}$’s are independent and uniformly distributed in $(0,1)$. Then, the proposed simulations consist on affinely transform each generated data matrix $X$ in each repetition, by applying the random matrix of transformation $A$ to $X$, in order to obtain $X_{A}$.

The contamination scheme consist in a mixture of normals $(1-\alpha)N(\mathbf{0},I)+\alpha N(\delta\mathbf{e},\lambda I)$. The dimension $p=5,10,30$, with sample size $n=100,100,500$ respectively, the contamination level $\alpha=10\%,20\%,30\%$, the distance of the outliers $\delta=5,10$, and the concentration of the contamination $\lambda=0.1,1$. Table 5 shows the obtained results about the correct and false detection rates.

As it can be observed, even under affine transformations, RMD-Shrinkage is able to detect all the outliers, except for a few cases (in bold type) that corresponds to large contamination level ($30\%$) in case of outliers close to the center of the distribution. However, it can be noted that these cases improve in performance when dimension increases.

3.6 Breakdown value

The maximum proportion of outliers that the estimator can safely tolerate is known as the breakdown value. For an outlier detection method, it can be defined as the maximum $\alpha^{*}$ outliers that the procedure can successfully detect, so that if $\alpha>\alpha^{*}$ the method will fail to identify most of the true outliers and it will falsely detect many inliers, reducing drastically the“c” and inflating the “f”. Thus, it is necessary to use the correct and false detection rates for studying the breakdown value of the outlier detection procedure.

Analogously as Sajesh and Srinivasan [2012], we consider two forms of contamination: $\alpha$ percent symmetric, for which the $i$th observation is multiplied by $100i$, and $\alpha$ percent asymmetric, for which the $i$th observation is replaced by $(100i)\mathbf{e}$, $i=1,...,n\alpha$, where $\mathbf{e}=(1,...,1)$. In the first case the outliers are symmetrically distributed, and asymmetrically in the second case.

The dimensions considered are $p=10,30,50,80,100$ and the sample size $n=1000$. The contamination level $\alpha=10\%,20\%,30\%,40\%,45\%$. Table 6 gives the resulting correct and false detection rates for both forms of contamination. For all the values of $p$ and $\alpha$ the method attains a $100\%$ correct detection rate. The maximum false detection rate for symmetric case is $0.55\%$ and for the asymmetric case is $0.47\%$. These results, even for high values of $\alpha$, indicates that RMD-Shrinkage can robustly detect large amount of contamination.

3.7 Computational times

Table 7 show the resulting computational times in seconds for the Normal case when outliers are close to the center of the data and they are concentrated. The other tables can be founded in the Supplementary Material.

On average, the fastest methods are COM and RMD-Shrinkage with very similar computational times. When compared to the MCD and its adjusted version Adj MCD, the latters are much more slower than our proposal. The Kurtosis and OGK are not that slower but they show worse times than ours. Thus, RMD-Shrinkage shows competitive computational times.

4 Real dataset

The proposed RMD’s are applied to a real dataset to evaluate their performance. The following dataset was taken from the UCI Knowledge Discovery in Databases Archive [Bay, 1999]. Specifically, we have chosen the Breast Cancer Wisconsin (Diagnostic) Data Set (WDBC). Features are computed from a digitized image of a fine needle aspirate of a breast mass. They describe 30 characteristics of the cell nuclei present in the image, for 569 samples, from which 357 are benign and 212 malign. We propose to study only the 357 benign data, as Maronna and Zamar [2002]. Therefore, this example has dimension $p=30$ and sample size $n=357$. We applied each method for detecting outliers and we retained the results, along with the computational times.

In order to interpret the outcome, we show the standardized data (after the detection) only for better visualization aim. We have also plotted the multivariate $L_{1}$ median and a kind of “multivariate boxplot”, which is based on the idea from Sun and Genton [2011] method, but for finite dimensional. What the “box” would be is constructed sorting the data according to their $L_{1}$ depth value. The corresponding $Q_{1}$ and $Q_{3}$ “quartiles” delimiting the “box” are in fact the minimum and maximum values for each coordinate taking only into account the $50\%$ of the most central data. Thus, the “fences” can be constructed with the same approach $F_{1}=Q_{1}-1.5RI$ and $F_{2}=Q_{3}+1.5RI$, where the “interquartile range” is $RI=Q_{3}-Q_{1}$. Then, we can look for each method’s result how many detected outliers are inside the “fences” for all their coordinates, and how many are outside the “fences”. Figure 4 shows the data in blue color plotted in parallel coordinates (Inselberg and Dimsdale [1990], Wegman [1990], Inselberg [2009]), the “box” delimiting the $50\%$ of most central data in yellow color, the “fences” in red and the multivariate $L_{1}-$median in cyan.

Table 8(a) shows the detected outliers by each method. Outside the “fences” there are 3 or 4 for all the methods. Also, the method Kurtosis detected 162 outliers out of the 357 data. More or less like OGK, which detected 148. Furthermore, our method RMD-Shrinkage is the one that label less amount of data as outliers. Table 8(b) shows how many outliers belong to the $50\%$ of the most central data, according to the $L_{1}$-median. We can investigate the shape of the detected outliers inside the “multivariate box”, in order to see if they are similar or near to the median, or if they have a distinct shape. Figure 11 shows the shape of some of our competitors detected outliers that belong to the $50\%$ of the most central data. In cyan color is the multivariate median, in yellow color the “box” and in blue color the detected outlier. For all of our competitors there seems to be some outliers having a shape very similar to the multivariate median or close to it for all the values of its components, leading us to think that maybe they are detecting too much. However, in the figure associated with RMD-Shrinkage, we can see that all outliers are quite different than the multivariate median, in fact, they might be shape outliers. For a final argument, we can say that all the outliers that belong to the “multivariate box” which are detected by the method RMD-Shrinkage, are actually detected by our competitors (the intersection matches), so this also makes us think that our method detects just enough.

Table 9 shows the computational times for each method in the task of detecting outliers with this example of real dataset. The results demonstrates that our competitors are much more slower than our proposal, except for COM which has a similar computational time.

5 Conclusions

Correct detection of outliers in the multivariate case is well-known to be a very important task for thorough data analysis. In order to reach that goal properly, it is necessary to consider the shape of the data and its structure in the multivariate space. That is the reason why the Mahalanobis distance approach is frequently used for the task of identifying the outliers. Different robust Mahalanobis distances can be defined according to the selected robust location and dispersion estimators. There are various robust estimators in the literature that have been considered in this paper. A collection of different combinations of robust location and covariance matrix estimators based on the notion of shrinkage is proposed, in order to define with each combination a robust Mahalanobis distance for the outlier detection problem. The performance of the proposed RMD’s and the others from the literature is shown through a simulation study. It can be concluded that our competitors do not have a very desirable over-all behavior, especially in high dimension. The proposed RMD’s have the ability to discover outliers with high precision in the vast majority of cases in the simulations, with gaussian data and with skewed or heavy-tailed distributions. RMD-Shrinkage is the most competitive version, as the simulation results showed. That is the reason why it is selected and some properties are investigated. The behavior under correlated and transformed data shows that RMD-Shrinkage is approximately affine equivariant. With highly contaminated data it is shown that the approach has high breakdown value even in high dimension. There is also evidence of its inexpensive computational time. A real dataset example is also studied, in which the results bear out with the latter conclusions.

The results presented in this article emphasize the advantages of using shrinkage estimators for the location and dispersion in the definition of a robust Mahalanobis distance. It remains to be examined whether the proposal could be improved by adapting the adjusted quantile to the proposed robust distances. It could also be an interesting matter to study, whether the use of the different definitions of “depth” in the literature (Tukey [1975], Liu et al. [1990], Serfling [2002], Chen et al. [2009], Agostinelli and Romanazzi [2011], Paindaveine and Van Bever [2013]), could improve the performance of the approach, as it is known that depth is a robust measure for location.

6 Supplementary Material

6.1 Proofs

6.1.1 Proof of Theorem 1.

The optimization problem is:

[TABLE]

where $\left\|\mathbf{x}\right\|_{2}^{2}=\sum_{j=1}^{p}x_{j}^{2}$ and the associated inner product is: $\langle x,y\rangle=\sum_{j=1}^{p}x_{j}y_{j}$.

The objective function is equivalent to:

[TABLE]

The latter element in the above expression is equal to zero because $E(\boldsymbol{\hat{\mu}}_{CCM})=\boldsymbol{\mu}$ (see Chu [1995]). Then, the optimization problem (29) reduces to minimize:

[TABLE]

In order to find the optimal $\nu_{\boldsymbol{\mu}}$, it is necessary to minimize only the right element of the above expression.

[TABLE]

Then, with respect to the scaling parameter, the first order optimality condition give:

[TABLE]

Thus:

[TABLE]

Estimating $\boldsymbol{\mu}$ with $\boldsymbol{\hat{\mu}}_{CCM}$, we obtain:

[TABLE]

In (30), with respect to the shrinkage intensity parameter $\eta$, the first order optimality condition give:

[TABLE]

Hence:

[TABLE]

6.1.2 Proof of Theorem 2.

The optimization problem is:

[TABLE]

where $\left\|x\right\|_{2}^{2}=\sum_{j=1}^{p}x_{j}^{2}$.

Similarly to the previous demonstration, we can consider the following expression for the objective function:

[TABLE]

The expectation of the inner product is equal to zero because Bose and Chaudhuri [1993] and Bose [1995] investigated the asymptotic distribution for the $L_{1}-$median, and they obtained the following result about the covariance matrix in presence of normality:

[TABLE]

Then, the optimization problem (31) reduces to minimize:

[TABLE]

Then, the optimal parameter $\nu_{\boldsymbol{\mu}}$ can be found minimizing only the right element of the above expression, which is the only one depending on that parameter.

[TABLE]

The associated first order optimality condition give:

[TABLE]

Therefore:

[TABLE]

In practice, we propose to estimate $\boldsymbol{\mu}$ with $\boldsymbol{\hat{\mu}}_{MM}$. Thus:

[TABLE]

With respect to the shrinkage intensity parameter $\eta$, the first order optimality condition associated to (32), give:

[TABLE]

Hence:

[TABLE]

6.1.3 Proof of Theorem 3.

The optimization problem is:

[TABLE]

where $\left\|A\right\|^{2}=trace(AA^{T})/p$, and the associated inner product is $\langle A_{1},A_{2}\rangle=trace(A_{1}A_{2}^{T})/p$.

Analogous to the previous Propositions, the objective function in the above minimization problem (33) can be seen as:

[TABLE]

In this case, note that the latter element in the above expression is equal to zero because $E(\hat{S}_{CCM})=\Sigma$. Hence, the optimization problem (33) reduces to minimize the following expression:

[TABLE]

The optimal $\nu_{\Sigma}$ can be obtained by minimizing only the right element of the above expression, because it is the only one depending on that parameter. Also, note that:

[TABLE]

Then, the first order optimality condition with respect to the scaling parameter, give:

[TABLE]

Therefore:

[TABLE]

In practice, we propose to estimate $\Sigma$ with $\hat{S}_{CCM}$, thus:

[TABLE]

In (34), with respect to the shrinkage intensity parameter $\eta$, the first order optimality condition give:

[TABLE]

6.2 Tables

6.2.1 Normal distribution

Table 10 shows the false detection rates when there is no contamination. The Tables 11-12 show the correct detection rates (c) and Tables 13-14 the false detection rates (f), for each method, corresponding to the simulations with multivariate Normal distribution, for contamination levels $\alpha=0.1,0.2,0.3$, dimension $p=5,10,30$, distance of the outliers $\delta=5,10$ and concentration of the contamination $\lambda=0.1,1$.

6.2.2 Computational times

6.3 Multivariate Student-t distribution with 3 degrees of freedom

Table 18 shows the false detection rates when there is no contamination. Tables 19-20 show the correct detection rates (c) and Tables 21-22 the false detection rates (f), for each method, corresponding to the simulations with multivariate Student-t distribution with 3 degrees of freedom, for contamination levels $\alpha=0.1,0.2,0.3$, dimension $p=5,10,30$, distance of the outliers $\delta=5,10$ and concentration of the contamination $\lambda=0.1,1$.

6.4 Multivariate Exponential distribution

Table 23 shows the false detection rates when there is no contamination. Table 24 shows the correct detection rates (c) and Table 25 the false detection rates (f), for each method, corresponding to the simulations with multivariate Exponential distribution, for contamination levels $\alpha=0.1,0.2,0.3$, dimension $p=5,10,30$ and distance of the outliers $\delta=5,10$.

6.5 Figures

The following are figures corresponding to the real dataset example.