## The result will be (almost) identical to the raw MCD
## (since we do not do reweighting of MRCD)
##
data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
c0 <- CovMcd(hbk.x, alpha=0.75, use.correction=FALSE)
cc <- CovMrcd(hbk.x, alpha=0.75)
cc$rho
all.equal(c0$best, cc$best)
all.equal(c0$raw.center, cc$center)
all.equal(c0$raw.cov/c0$raw.cnp2[1], cc$cov/cc$cnp2)
summary(cc)
## the following three statements are equivalent
c1 <- CovMrcd(hbk.x, alpha = 0.75)
c2 <- CovMrcd(hbk.x, control = CovControlMrcd(alpha = 0.75))
## direct specification overrides control one:
c3 <- CovMrcd(hbk.x, alpha = 0.75,
control = CovControlMrcd(alpha=0.95))
c1
if (FALSE) {
## This is the first example from Boudt et al. (2020). The first variable is
## the dependent one, which we remove and remain with p=226 NIR absorbance spectra
data(octane)
octane <- octane[, -1] # remove the dependent variable y
n <- nrow(octane)
p <- ncol(octane)
## Compute MRCD with h=33, which gives approximately 15 percent breakdown point.
## This value of h was found by Boudt et al. (2020) using a data driven approach,
## similar to the Forward Search of Atkinson et al. (2004).
## The default value of h would be 20 (i.e. alpha=0.5)
out <- CovMrcd(octane, h=33)
out$rho
## Please note that in the paper is indicated that the obtained rho=0.1149, however,
## this value of rho is obtained if the parameter maxcond is set equal to 999 (this was
## the default in an earlier version of the software, now the default is maxcond=50).
## To reproduce the result from the paper, change the call to CovMrcd() as follows
## (this will not influence the results shown further):
## out <- CovMrcd(octane, h=33, maxcond=999)
## out$rho
robpca = PcaHubert(octane, k=2, alpha=0.75, mcd=FALSE)
(outl.robpca = which(robpca@flag==FALSE))
# Observations flagged as outliers by ROBPCA:
# 25, 26, 36, 37, 38, 39
# Plot the orthogonal distances versus the score distances:
pch = rep(20,n); pch[robpca@flag==FALSE] = 17
col = rep('black',n); col[robpca@flag==FALSE] = 'red'
plot(robpca, pch=pch, col=col, id.n.sd=6, id.n.od=6)
## Plot now the MRCD mahalanobis distances
pch = rep(20,n); pch[!getFlag(out)] = 17
col = rep('black',n); col[!getFlag(out)] = 'red'
plot(out, pch=pch, col=col, id.n=6)
}
Run the code above in your browser using DataLab