Simulate multivariate outliers that are hidden in 2D scatterplots

How could I simulate multivariate outliers that are "hidden" in all pairwise 2D scatterplots between the variables? By "hidden" I mean that they can't be seen (as obvious outliers) or detected visually in all pairwise 2D scatterplots (e.g. a scatterplot matrix). The outliers should be suspicious, however, using multivariate outlier detection techniques such as the Mahalanobis distance.

For simplicity, I would simulate the non-outliers from a multivariate normal distribution preferably with an arbitrary covariance matrix.

The purpose of this is in teaching: I want to illustrate that pairwise scatterplots are not always a good tool to detect multivariate outliers. For this purpose, it makes little sense to simulate more than, say, $$p=5$$ dimensions.

Here's one way to get an outlier or outliers, even with uncorrelated variables (it can sometimes be 'easier' to hide outliers with correlated ones though).

Pick a dimension, $$p$$. Say $$p=9$$ (You may want $$p> 9$$ to be sure to get obvious outliers in $$p$$-D and clearly-not-looking-like-outliers in 2D.)

Start with a population that's iid $$N(0,I)$$; let's then take a sample of 10000 points from that say. Now consider the limit of what would be an outlier in any of the bivariate margins where $$X_i=X_j=u$$ for some $$u>0$$ (e.g. visually, anything inside $$u=k/\sqrt{2}$$ might look okay for some $$k \sim 3$$ or so, YMMV); so let's say $$u=2.12$$. The idea is that $$k$$ sd's from the mean would be your cutoff for a univariate outlier.

Set something like that as an upper bound.

Now identify what would be an outlier in all dimensions combined ($$X_1=X_2=...=X_p = l$$, for some $$l>0$$); say $$k/\sqrt{p}$$. So if $$k\sim 3$$, and $$p=9$$, that would give $$l\sim 1$$-ish. Let's bump that up and say $$l=1.6$$ (to get some nice clear outliers; or push it higher if you like).

Generate some values $$z_i$$, such that $$z_i$$ is between $$l$$ and $$u$$. Let $$X_1=X_2=...=X_p=z_i$$. You can add some jitter and if you wish, flip the signs on some of the $$X$$'s.

What we're doing is generating some way in the direction of the "corners" of a box in p-dimensions such that we're inside a circle in every pair of dimensions but outside a ball in $$\mathbb{R}^p$$.

You'll want to play with the numbers until it suits you but this is the basic idea.

Here's an example. Here I used $$p=25$$ and generated 10000 iid standard normal vectors of length $$p$$, then added a single outlier at $$(1.625,1.625,1.625,...,1.625)'$$. I found the worst bivariate case in the sense of identifying the pairwise dimension where the Mahalanobis' distance of the outlier had the highest quantile. Here's a display showing that and the Mahalanobis distance for all p-dimensions:

The added outlier is marked with a red line segment under the histogram (around 2.3 on the left plot and around 8.13 on the right).

We see that the 2D worst case isn't so extreme we'd be likely to call it an outlier, but in the full 25 dimensions it sticks out the end more than all the most outlying cases of the original data, even with 10000 points.

• Thanks, that's great. Could this be easily extended to the case where the population is multivariate normal with an arbitrary covariance matrix instead of the identity matrix? – COOLSerdash Oct 16 '19 at 8:40
• If you're measuring outlyingness by say a Mahalanobis distance type measure, take the above setup and transform it to have an arbitrary covariance matrix. (It may be easier to visually hide bivariate outliers in some directions than others however.) – Glen_b -Reinstate Monica Oct 16 '19 at 22:18
• This is a good idea (+1). I developed code for it when the question appeared and discovered it doesn't work quite as well as one might hope. The problem is that many of these "corner" areas are also corners in the scatterplot matrix and therefore show up as obvious outliers anyway. The best solution I could think of was to generate a lot of such corner points and reject those that do show up as outliers in the SPM. That might be horribly inefficient in large dimensions. – whuber Oct 17 '19 at 10:47
• In large dimension you would be able to push down the bivariate upper bound. e.g. if d=100, something that's 1sd above the mean on every variable is 10 sd's from the center in d-space but is only 1.4 sd's from the center in every pairwise plot. With some care to choose $l$ and $u$ you shouldn't have to throw out many at all. – Glen_b -Reinstate Monica Oct 17 '19 at 12:15
• That's right--but it can take considerable experimentation. On the other hand, it's extremely frustrating to review SPMs for 100 variables! That suggests the question concerns much smaller dimensions than 100 (as confirmed by a recent edit). – whuber Oct 17 '19 at 13:48