Description Usage Arguments Details Value Note Author(s) References See Also Examples
View source: R/adjoutlyingness.R
For an n * p data matrix (or data frame) x
,
compute the “outlyingness” of all n observations.
Outlyingness here is a generalization of the DonohoStahel
outlyingness measure, where skewness is taken into account via the
medcouple, mc()
.
1 2 3 4 5  adjOutlyingness(x, ndir = 250, p.samp = p, clower = 4, cupper = 3,
alpha.cutoff = 0.75, coef = 1.5,
qr.tol = 1e12, keep.tol = 1e12,
only.outlyingness = FALSE, maxit.mult = max(100, p),
trace.lev = 0)

x 
a numeric 
ndir 
positive integer specifying the number of directions that should be searched. 
p.samp 
the sample size to use for finding good random
directions, must be at least 
clower, cupper 
the constant to be used for the lower and upper
tails, in order to transform the data towards symmetry. You can set

alpha.cutoff 
number in (0,1) specifying the quantiles
(α, 1α) which determine the “outlier”
cutoff. The default, using quartiles, corresponds to the definition
of the medcouple ( 
coef 
positive number specifying the factor with which the
interquartile range ( 
qr.tol 
positive tolerance to be used for 
keep.tol 
positive tolerance to determine which of the sample
direction should be kept, namely only those for which
x * B is larger than 
only.outlyingness 
logical indicating if the final outlier determination should be skipped. In that case, a vector is returned, see ‘Value:’ below. 
maxit.mult 
integer factor; 
trace.lev 
an integer, if positive allows to monitor the direction search. 
FIXME: Details in the comment of the Matlab code; also in the reference(s).
The method as described can be useful as preprocessing in FASTICA (http://www.cis.hut.fi/projects/ica/fastica/; see also the R package fastICA.
If only.outlyingness
is true, a vector adjout
,
otherwise, as by default, a list with components
adjout 
numeric of 
cutoff 
cutoff for “outlier” with respect to the adjusted
outlyingnesses, and depending on 
nonOut 
logical of 
The result is random as it depends on the sample of
ndir
directions chosen; hence set.seed()
yourself
for reproducibility!
Till Aug/Oct. 2014, the default values for clower
and cupper
were
accidentally reversed, and the signs inside exp(.)
where swapped
in the (now corrected) two expressions
1 2 
already in the code from Antwerpen (‘mcrsoft/adjoutlingness.R’), contrary to the published reference.
Further, the original algorithm had not been scaleequivariant in the direction construction, which has been amended in 201410 as well.
The results, including diagnosed outliers, therefore have changed, typically slightly, since robustbase version 0.920.
Guy Brys; help page and improvements by Martin Maechler
Brys, G., Hubert, M., and Rousseeuw, P.J. (2005) A Robustification of Independent Component Analysis; Journal of Chemometrics, 19, 1–12.
Hubert, M., Van der Veeken, S. (2008) Outlier detection for skewed data; Journal of Chemometrics 22, 235–246.
For the uptodate reference, please consult http://wis.kuleuven.be/stat/robust
the adjusted boxplot, adjbox
and the medcouple,
mc
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51  ## An Example with bad condition number and "border case" outliers
dim(longley)
set.seed(1) ## result is random!
ao1 < adjOutlyingness(longley)
## which are outlying ?
which(!ao1$nonOut) ## one: "1948"  for this seed! (often: none)
stopifnot(all(ao1$nonOut[2]))
## An Example with outliers :
dim(hbk)
set.seed(1)
ao.hbk < adjOutlyingness(hbk)
str(ao.hbk)
hist(ao.hbk $adjout)## really two groups
table(ao.hbk$nonOut)## 14 outliers, 61 nonoutliers:
## outliers are :
which(! ao.hbk$nonOut) # 1 .. 14  but not for all random seeds!
## here, they are the same as found by (much faster) MCD:
cc < covMcd(hbk)
stopifnot(all(cc$mcd.wt == ao.hbk$nonOut))
## This is revealing: About 12 cases, where outliers are *not* == 1:14
## but needs almost 1 [sec] per call:
if(interactive()) {
for(i in 1:30) {
print(system.time(ao.hbk < adjOutlyingness(hbk)))
if(!identical(iout < which(!ao.hbk$nonOut), 1:14)) {
cat("Outliers:\n"); print(iout)
}
}
}
## "Central" outlyingness: *not* calling mc() anymore, since 20141211:
trace(mc)
out < capture.output(
oo < adjOutlyingness(hbk, clower=0, cupper=0)
)
untrace(mc)
stopifnot(length(out) == 0)
## A rankdeficient case
T < tcrossprod(data.matrix(toxicity))
try(adjOutlyingness(T, maxit. = 20, trace.lev = 2)) # fails and recommends:
T. < fullRank(T)
aT < adjOutlyingness(T.)
plot(sort(aT$adjout, decreasing=TRUE), log="y")
plot(T.[,9:10], col = (1:2)[1 + (aT$adjout > 10000)])
## .. (not conclusive; directions are random, more 'ndir' makes a difference!)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.