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.
 A: 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.
