7
$\begingroup$

Currently I’m trying to understand and use the Stahel-Donoho outlyngness measure.

But unfortunately I’ve got a problem in the part where one is taking the maximum over the set of directions. I found papers and tutorials about the Stahel-Donoho estimator, but their explanation is not good enough for me. In one paper, it has been suggested that one way to find this set of directions is by taking $p$ directions through 2 randomly chosen observations. In another paper [0], the authors use many directions through 2 randomly chosen observations. The original definition itself uses as direction $d$ the direction orthogonal to the hyperplane through $p$ randomly chosen data points. Then, the full set of directions is obtained by repeating this procedure for many random $p$ subsets (see page 3, 2nd paragraph from the bottom of this paper).

Unfortunately I cannot get any idea how to deal with this problem. Should I use random projection? The problem I'm having is that I don’t know how to obtain the directions $d$ and projections. I really stuck at the moment. Could someone kindly help me to understand or provide me a good reference for better understanding?

  • [0] Hubert, M., Rousseeuw P. J. and Vanden Branden, K. (2005). ROBPCA: A New Approach to Robust Principal Component Analysis. Technometrics Volume 47, Issue 1.
$\endgroup$
9
  • $\begingroup$ " In some papers, it has been mentioned that there is two way to obtain direction which is random and fixed direction." "In the literature it has been suggested that another possibility is taking p times a direction through 2 randomly chosen observations. " Can you cite your sources for those statements? $\endgroup$
    – user603
    Commented Nov 9, 2014 at 9:25
  • 1
    $\begingroup$ AN OUTLIER MAP FOR SUPPORT VECTOR MACHINE CLASSIFICATION[link]( arxiv.org/pdf/1009.5818.pdf) ROBPCA: A New Approach to Robust Principal Component Analysislink Combining Random and Specific Directions for Outlier Detection and Robust Estimation in High-Dimensional Multivariate Datalink $\endgroup$ Commented Nov 9, 2014 at 9:57
  • $\begingroup$ I have edited your question to make it more readable and include the links you provided. Can you confirm that I didn t alter its meaning? $\endgroup$
    – user603
    Commented Nov 9, 2014 at 10:55
  • 1
    $\begingroup$ okay, I know the answers to all these questions;) would you prefer a math answer or one using R? $\endgroup$
    – user603
    Commented Nov 9, 2014 at 15:29
  • 1
    $\begingroup$ Oh thank you very much. Would it possible for me to ask for both :D $\endgroup$ Commented Nov 9, 2014 at 15:47

1 Answer 1

10
$\begingroup$

Summary: the only 'correct' definition is the original one [0][1], the other ones were designed to solve a problem with it that happens in a specific context and would probably have been best called 'pseudo-Stahel-Donoho' distances because they don't, in general, have the same interpretation.

Now, I will show you how to obtain the standard SD distances in R, explain why they are computed the way they are, then explain what these pseudo-SD 'distances' are and why the need for them arose.

Given an $n$ by $p$ data matrix $X=\{x_1,\ldots,x_n\}$ whose entries lie in general (linear) position in $\mathbb{R}^p$, the original definition is to obtain a $p$-vector $d$ as the vector of coefficients of the hyperplane through $p$ data points chosen randomly out of your sample of $n$ observations. In R this is done so:

moar_directions<-function(x){
    n<-nrow(x)
    p<-ncol(x)
    P<-x[sample(n,p),,drop=FALSE]
    G<-rep(NA,p)
    E<-matrix(1,p,1)
    if((qrP<-qr(P))$rank==p) G<-solve(qrP,E)
    return(G)
}
n<-100
p<-5
x<-matrix(rnorm(n*p),nc=p)
moar_directions(x)

(see for example line 13 of robustbase::adjOutlyingness here, which you can check by tipping:

library(robustbase)
body(robustbase::adjOutlyingness) 

and [2]). Then, projecting the data unto the direction normal to $d$ is done so:

x%*%moar_directions(x)

The SD outlyingness of $x_i$ with respect to $X=\{x_1,\ldots,x_n\}$ on a single projection unto $d$ is ([0][1]):

$$O(x_i,d,X)=\frac{|x_i'd-\mbox{med}_i(x_i'd)|}{\mbox{mad}_i(x_i'd)}\;\;\;(0)$$

which in R is obtained as:

w<-x%*%moar_directions(x)
abs(w-median(w))/mad(w)

And computing the SD outlyingness as the maximum outlyingness over all the members of $B_p^n$, where $B_p^n$ is the set of all ${n\choose p}$ such directions $d$ in a $n$ by $p$ data matrix whose entries lie in general (linear) position in $\mathbb{R}^p$. Often, $|B_p^n|={n\choose p}$ will be too large, and it will not be possible to consider all its members in which case one can sample randomly $K$ directions from it to compute an approximation to the SD outlyigness of $x_i$ w.r.t. $X$:

$$O(x_i,X)=\underset{\{d_k\}_{k=1}^K\in B_p^n}{\max}\frac{|x_i'd_k-\mbox{med}_i(x_i'd_k)|}{\mbox{mad}_i(x_i'd_k)}$$

which is R is obtained as:

library(matrixStats)
K<-100
AllThoseProj<-matrix(NA,n,K)
for(i in 1:K){
    w<-moar_directions(x)
    AllThoseProj[,i]<-abs(w-median(w))/mad(w)
}
rowMaxs(AllThoseProj)

(I will use for a couple of lines K=100 as my definition of 'many' then latter explain what is meant by many in this context) To motivate the original definition, notice[3] that:

$$d(x_i,X)=\underset{d\in B_p}{\max}\frac{|x_i'd_k-\mbox{mean}_i(x_i'd_k)|}{\mbox{sd}_i(x_i'd_k)}\;\;(1)$$

where

$$d^2(x_i,X)=(x_i-\mbox{mean}_i(x_i))'\mbox{Cov}_i(x_i)^{-1}(x_i-\mbox{mean}_i(x_i))$$

is the vector of squared Mahalanobis distances. So the original SD outlyingness index was designed as a way to compute a consistent estimator of $d^2(x_i)$ with 50% breakdown point (see [1]) and equality (1) only holds when the directions $d$ are defined as above.

Now, I will delve a bit more on what is meant by 'many' in ' Often, $|B_p^n|$ will be too large, and ... one can sample randomly many directions'

The notion of many that should be used here depends on the objective being pursued. If you are using the SD to approximate the vector of $d(x_i,X)$'s, then $K$ in the high hundreds should already give a good approximation. If you are using the SD as a robust alternative to the vector of $d(x_i,X)$'s, then a much higher value of $K$ will be necessary (see page 13 of [0] for a discussion of this).

Now, I will delve a bit more on these pseudo-SD distances.

The problem with the way $d$ is defined above is that it is only uniquely defined in settings where $n>p$. In the last decade a lot of research started to be done on high dimensional robustness. The idea built up incrementally to search for outliers through projection pursuit (as was done with the SD distances) but using a type of projection that would also 'work' when $p>n$. There are many such notions, but, gradually, a consensus emerged around the the idea of using directions through two points (see [4] for a late implementation of this idea). Directions through two points are defined as (see [4], and line 52 of body(rrcov:::extradir) in package rrcov for example):

moar_directions_2points<-function(x){
        n<-nrow(x)
        p<-ncol(x)
        P<-x[sample(n,2),,drop=FALSE]
        G<-rep(NA,p)
        E<-G
        E<-P[1,]-P[2,]
        N<-sqrt(crossprod(E))
        if(N>1e-8)  G<-E/N        
        return(G)
}

Now, if we substitute in equation (0) the original definition of the directions $d$ by the ones above, the new outlyingness index (let's call it pseudo-SD) is no longer (except when $p=2$) a consistent estimator of the vector of $d(x_i,X)$'s, so it is a bit hard to interpret what it is. On the other hand, it can still be computed when $p>n$ in which case the vector of $d(x_i)$'s is not even defined anyway. The fact that we lost consistency by using directions through 2 data points also means we no longer have a target (a quantity we are trying to estimate, like we did with the the vector of $d(x_i,X)$'s in the case of the SD) so the question of what is meant by 'many' here is not really addressable. rrcov for example uses 'many' (in the case of directions through two points) to mean 250.

References:

  • [0] Stahel W. (1981). Breakdown of Covariance Estimators. Research Report 31, Fachgrupp fur Statistik, E.T.H. Zurich.
  • [1] Donoho. D.L. (1982). Breakdown properties of multivariate location estimators Ph.D. Qualifying Paper Harvard University.
  • [2] Hubert, M. and Van der Veeken, S. (2007). Outlier detection for skewed data. Journal of chemometrics vol:22 issue:3-4 pages:235-246.
  • [3] Rousseeuw, P.J. and Leroy, A.M. (1987). Robust Regression and Outlier Detection. Wiley, New York.
  • [4] Hubert, M., Rousseeuw P. J. and Vanden Branden, K. (2005). ROBPCA: A New Approach to Robust Principal Component Analysis. Technometrics Volume 47, Issue 1.
$\endgroup$
7
  • 1
    $\begingroup$ Thank you so much for your detailed answer.Your explanation is very helpful. $\endgroup$ Commented Nov 10, 2014 at 20:27
  • $\begingroup$ According to you for high dimensional data ($p>n$), I should use the second function (moar_directions_2points). In ( stat.osu.edu/~statgen/joul_win2010/Robpca05.pdf) it has been suggested that if ${n\choose p}> 250$ we should take randomly 250 directions form $B$. So I think in second step (sample randomly $k$ directions) I should chose $k=250$ instead of $k=100$. Am I right? $\endgroup$ Commented Nov 15, 2014 at 9:51
  • $\begingroup$ Moreover would you please suggest me some references which is simply explain the direction. I’ve tried to understand it from PCA concept but unfortunately much success has not been achieved so far! Thank you $\endgroup$ Commented Nov 15, 2014 at 9:56
  • $\begingroup$ Is it possible to comment on lines 7:10 of the moar_directions_2points funcion? $\endgroup$ Commented Nov 15, 2014 at 11:42
  • 1
    $\begingroup$ first question: read the last 2 sentences in my answer, they address your first question. second question: yes, see [4], third question: as I wrote, the function does the same thing as line 52 of rrcov:::extradir but feel free to use the code in that function if mine is not clear. $\endgroup$
    – user603
    Commented Nov 16, 2014 at 20:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.