Package 'rrcov'

Description

The data on annual maximum streamflow at 104 gaging stations in the central Appalachia region of the United States contains the sample L-moments ratios (L-CV, L-skewness and L-kurtosis) as used by Hosking and Wallis (1997) to illustrate regional freqency analysis (RFA).

Usage

data(Appalachia)

Format

A data frame with 104 observations on the following 3 variables:

L-CV

L-coefficient of variation

L-skewness

L-coefficient of skewness

L-kurtosis

L-coefficient of kurtosis

Details

The sample L-moment ratios (L-CV, L-skewness and L-kurtosis) of a site are regarded as a point in three dimensional space.

Source

Hosking, J. R. M. and J. R. Wallis (1997), Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, p.175–185

References

Neykov, N.M., Neytchev, P.N., Van Gelder, P.H.A.J.M. and Todorov V. (2007), Robust detection of discordant sites in regional frequency analysis, Water Resources Research, 43, W06417, doi:10.1029/2006WR005322

Examples

data(Appalachia)

    # plot a matrix of scatterplots
    pairs(Appalachia,
          main="Appalachia data set",
          pch=21,
          bg=c("red", "green3", "blue"))

    mcd<-CovMcd(Appalachia)
    mcd
    plot(mcd, which="dist", class=TRUE)
    plot(mcd, which="dd", class=TRUE)

    ##  identify the discordant sites using robust distances and compare 
    ##  to the classical ones
    mcd <- CovMcd(Appalachia)
    rd <- sqrt(getDistance(mcd))
    ccov <- CovClassic(Appalachia)
    cd <- sqrt(getDistance(ccov))
    r.out <- which(rd > sqrt(qchisq(0.975,3)))
    c.out <- which(cd > sqrt(qchisq(0.975,3)))
    cat("Robust: ", length(r.out), " outliers: ", r.out,"\n")
    cat("Classical: ", length(c.out), " outliers: ", c.out,"\n")

Description

Produces a biplot from an object (derived from) Pca-class.

Usage

## S4 method for signature 'Pca'
biplot(x, choices=1L:2L, scale=1, ...)

Arguments

Side Effects

a plot is produced on the current graphics device.

Methods

biplot

signature(x = Pca): Plot a biplot, i.e. represent both the observations and variables of a matrix of multivariate data on the same plot. See also biplot.princomp.

References

Gabriel, K. R. (1971). The biplot graphical display of matrices with applications to principal component analysis. Biometrika, 58, 453–467.

See Also

Pca-class, PcaClassic, PcaRobust-class.

Examples

require(graphics)
biplot(PcaClassic(USArrests, k=2))

Description

The data set bus (Hettich and Bay, 1999) corresponds to a study in automatic vehicle recognition (see Maronna et al. 2006, page 213, Example 6.3)). This data set from the Turing Institute, Glasgow, Scotland, contains measures of shape features extracted from vehicle silhouettes. The images were acquired by a camera looking downward at the model vehicle from a fixed angle of elevation. Each of the 218 rows corresponds to a view of a bus silhouette, and contains 18 attributes of the image.

Usage

data(bus)

Format

A data frame with 218 observations on the following 18 variables:

V1

compactness

V2

circularity

V3

distance circularity

V4

radius ratio

V5

principal axis aspect ratio

V6

maximum length aspect ratio

V7

scatter ratio

V8

elongatedness

V9

principal axis rectangularity

V10

maximum length rectangularity

V11

scaled variance along major axis

V12

scaled variance along minor axis

V13

scaled radius of gyration

V14

skewness about major axis

V15

skewness about minor axis

V16

kurtosis about minor axis

V17

kurtosis about major axis

V18

hollows ratio

Source

Hettich, S. and Bay, S.D. (1999), The UCI KDD Archive, Irvine, CA:University of California, Department of Information and Computer Science, 'http://kdd.ics.uci.edu'

References

Maronna, R., Martin, D. and Yohai, V., (2006). Robust Statistics: Theory and Methods. Wiley, New York.

Examples

## Reproduce Table 6.3 from Maronna et al. (2006), page 213
    data(bus)
    bus <- as.matrix(bus)
    
    ## calculate MADN for each variable
    xmad <- apply(bus, 2, mad)      
    cat("\nMin, Max of MADN: ", min(xmad), max(xmad), "\n")


    ## MADN vary between 0 (for variable 9) and 34. Therefore exclude 
    ##  variable 9 and divide the remaining variables by their MADNs.
    bus1 <- bus[, -9]
    madbus <- apply(bus1, 2, mad)
    bus2 <- sweep(bus1, 2, madbus, "/", check.margin = FALSE)

    ## Compute classical and robust PCA (Spherical/Locantore, Hubert, MCD and OGK)    
    pca  <- PcaClassic(bus2)
    rpca <- PcaLocantore(bus2)
    pcaHubert <- PcaHubert(bus2, k=17, kmax=17, mcd=FALSE)
    pcamcd <- PcaCov(bus2, cov.control=CovControlMcd())
    pcaogk <- PcaCov(bus2, cov.control=CovControlOgk())

    ev    <- getEigenvalues(pca)
    evrob <- getEigenvalues(rpca)
    evhub <- getEigenvalues(pcaHubert)
    evmcd <- getEigenvalues(pcamcd)
    evogk <- getEigenvalues(pcaogk)

    uvar    <- matrix(nrow=6, ncol=6)
    svar    <- sum(ev)
    svarrob <- sum(evrob)
    svarhub <- sum(evhub)
    svarmcd <- sum(evmcd)
    svarogk <- sum(evogk)
    for(i in 1:6){
        uvar[i,1] <- i
        uvar[i,2] <- round((svar - sum(ev[1:i]))/svar, 3)
        uvar[i,3] <- round((svarrob - sum(evrob[1:i]))/svarrob, 3)
        uvar[i,4] <- round((svarhub - sum(evhub[1:i]))/svarhub, 3)
        uvar[i,5] <- round((svarmcd - sum(evmcd[1:i]))/svarmcd, 3)
        uvar[i,6] <- round((svarogk - sum(evogk[1:i]))/svarogk, 3)
    }
    uvar <- as.data.frame(uvar)
    names(uvar) <- c("q", "Classical","Spherical", "Hubert", "MCD", "OGK")
    cat("\nBus data: proportion of unexplained variability for q components\n")
    print(uvar)
 
    ## Reproduce Table 6.4 from Maronna et al. (2006), page 214
    ##
    ## Compute classical and robust PCA extracting only the first 3 components
    ## and take the squared orthogonal distances to the 3-dimensional hyperplane
    ##
    pca3 <- PcaClassic(bus2, k=3)               # classical
    rpca3 <- PcaLocantore(bus2, k=3)            # spherical (Locantore, 1999)
    hpca3 <- PcaHubert(bus2, k=3)               # Hubert
    dist <- pca3@od^2
    rdist <- rpca3@od^2
    hdist <- hpca3@od^2

    ## calculate the quantiles of the distances to the 3-dimensional hyperplane
    qclass  <- round(quantile(dist, probs = seq(0, 1, 0.1)[-c(1,11)]), 1)
    qspc <- round(quantile(rdist, probs = seq(0, 1, 0.1)[-c(1,11)]), 1)
    qhubert <- round(quantile(hdist, probs = seq(0, 1, 0.1)[-c(1,11)]), 1)
    qq <- cbind(rbind(qclass, qspc, qhubert), round(c(max(dist), max(rdist), max(hdist)), 0))
    colnames(qq)[10] <- "Max"
    rownames(qq) <- c("Classical", "Spherical", "Hubert")
    cat("\nBus data: quantiles of distances to hiperplane\n")
    print(qq)

    ## 
    ## Reproduce Fig 6.1 from Maronna et al. (2006), page 214
    ## 
    cat("\nBus data: Q-Q plot of logs of distances to hyperplane (k=3) 
    \nfrom classical and robust estimates. The line is the identity diagonal\n")
    plot(sort(log(dist)), sort(log(rdist)), xlab="classical", ylab="robust")
    lines(sort(log(dist)), sort(log(dist)))

Description

This data set is based on the bushfire data set which was used by Campbell (1984) to locate bushfire scars - see bushfire in package robustbase. The original dataset contains satelite measurements on five frequency bands, corresponding to each of 38 pixels.

Usage

data(bushmiss)

Format

A data frame with 190 observations on 6 variables.

The original data set consists of 38 observations in 5 variables. Based on it four new data sets are created in which some of the data items are replaced by missing values with a simple "missing completely at random " mechanism. For this purpose independent Bernoulli trials are realized for each data item with a probability of success 0.1, 0.2, 0.3, 0.4, where success means that the corresponding item is set to missing. The obtained five data sets, including the original one (each with probability of a data item to be missing equal to 0, 0.1, 0.2, 0.3 and 0.4 which is reflected in the new variable MPROB) are merged. (See also Beguin and Hulliger (2004).)

Source

Maronna, R.A. and Yohai, V.J. (1995) The Behavoiur of the Stahel-Donoho Robust Multivariate Estimator. Journal of the American Statistical Association 90, 330–341.

Beguin, C. and Hulliger, B. (2004) Multivariate outlier detection in incomplete survey data: the epidemic algorithm and transformed rank correlations. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 127, 2, 275–294.

Examples

## The following code will result in exactly the same output
##  as the one obtained from the original data set
data(bushmiss)
bf <- bushmiss[bushmiss$MPROB==0,1:5]
plot(bf)
covMcd(bf)


## Not run: 
##  This is the code with which the missing data were created:
##
##  Creates a data set with missing values (for testing purposes)
##  from a complete data set 'x'. The probability of
##  each item being missing is 'pr' (Bernoulli trials).
##
getmiss <- function(x, pr=0.1)
{
    n <- nrow(x)
    p <- ncol(x)
    done <- FALSE
    iter <- 0
    while(iter <= 50){
        bt <- rbinom(n*p, 1, pr)
        btmat <- matrix(bt, nrow=n)
        btmiss <- ifelse(btmat==1, NA, 0)
        y <- x+btmiss
        if(length(which(rowSums(nanmap(y)) == p)) == 0)
            return (y)
        iter <- iter + 1
    }
    y
}

## End(Not run)

Description

A data frame containing 11 variables with different dimensions of 111 cars

Usage

data(Cars)

Format

A data frame with 111 observations on the following 11 variables.

length

a numeric vector

wheelbase

a numeric vector

width

a numeric vector

height

a numeric vector

front.hd

a numeric vector

rear.hd

a numeric vector

front.leg

a numeric vector

rear.seating

a numeric vector

front.shoulder

a numeric vector

rear.shoulder

a numeric vector

luggage

a numeric vector

Source

Consumer reports. (April 1990). http://backissues.com/issue/Consumer-Reports-April-1990, pp. 235–288.

References

Chambers, J. M. and Hastie, T. J. (1992). Statistical models in S. Cole, Pacific Grove, CA: Wadsworth and Brooks, pp. 46–47.

M. Hubert, P. J. Rousseeuw, K. Vanden Branden (2005), ROBPCA: A new approach to robust principal components analysis, Technometrics, 47, 64–79.

Examples

data(Cars)

## Plot a pairwise scaterplot matrix
    pairs(Cars[,1:6])

    mcd <- CovMcd(Cars[,1:6])    
    plot(mcd, which="pairs")
    
## Start with robust PCA
    pca <- PcaHubert(Cars, k=ncol(Cars), kmax=ncol(Cars))
    pca

## Compare with the classical PCA
    prcomp(Cars)

## or  
    PcaClassic(Cars, k=ncol(Cars), kmax=ncol(Cars))
    
## If you want to print the scores too, use
    print(pca, print.x=TRUE)

## Using the formula interface
    PcaHubert(~., data=Cars, k=ncol(Cars), kmax=ncol(Cars))

## To plot the results:

    plot(pca)                    # distance plot
    pca2 <- PcaHubert(Cars, k=4)  
    plot(pca2)                   # PCA diagnostic plot (or outlier map)
    
## Use the standard plots available for prcomp and princomp
    screeplot(pca)    # it is interesting with all variables    
    biplot(pca)       # for biplot we need more than one PCs
    
## Restore the covraiance matrix     
    py <- PcaHubert(Cars, k=ncol(Cars), kmax=ncol(Cars))
    cov.1 <- py@loadings %*% diag(py@eigenvalues) %*% t(py@loadings)
    cov.1

Description

The data on annual precipitation totals for the North Cascades region contains the sample L-moments ratios (L-CV, L-skewness and L-kurtosis) for 19 sites as used by Hosking and Wallis (1997), page 53, Table 3.4, to illustrate screening tools for regional freqency analysis (RFA).

Usage

data(Cascades)

Format

A data frame with 19 observations on the following 3 variables.

L-CV

L-coefficient of variation

L-skewness

L-coefficient of skewness

L-kurtosis

L-coefficient of kurtosis

Details

The sample L-moment ratios (L-CV, L-skewness and L-kurtosis) of a site are regarded as a point in three dimensional space.

Source

Hosking, J. R. M. and J. R. Wallis (1997), Regional Frequency Analysis: An Approach Based on L-moments. Cambridge University Press, p. 52–53

References

Neykov, N.M., Neytchev, P.N., Van Gelder, P.H.A.J.M. and Todorov V. (2007), Robust detection of discordant sites in regional frequency analysis, Water Resources Research, 43, W06417, doi:10.1029/2006WR005322

Examples

data(Cascades)

    # plot a matrix of scatterplots
    pairs(Cascades,
          main="Cascades data set",
          pch=21,
          bg=c("red", "green3", "blue"))

    mcd<-CovMcd(Cascades)
    mcd
    plot(mcd, which="dist", class=TRUE)
    plot(mcd, which="dd", class=TRUE)

    ##  identify the discordant sites using robust distances and compare 
    ##  to the classical ones
    rd <- sqrt(getDistance(mcd))
    ccov <- CovClassic(Cascades)
    cd <- sqrt(getDistance(ccov))
    r.out <- which(rd > sqrt(qchisq(0.975,3)))
    c.out <- which(cd > sqrt(qchisq(0.975,3)))
    cat("Robust: ", length(r.out), " outliers: ", r.out,"\n")
    cat("Classical: ", length(c.out), " outliers: ", c.out,"\n")

Description

The class Cov represents an estimate of the multivariate location and scatter of a data set. The objects of class Cov contain the classical estimates and serve as base for deriving other estimates, i.e. different types of robust estimates.

Objects from the Class

Objects can be created by calls of the form new("Cov", ...), but the usual way of creating Cov objects is a call to the function Cov which serves as a constructor.

Slots

call:

Object of class "language"

cov:

covariance matrix

center:

location

n.obs:

number of observations used for the computation of the estimates

mah:

mahalanobis distances

det:

determinant

flag:

flags (FALSE if suspected an outlier)

method:

a character string describing the method used to compute the estimate: "Classic"

singularity:

a list with singularity information for the covariance matrix (or NULL of not singular)

X:

data

Methods

getCenter

signature(obj = "Cov"): location vector

getCov

signature(obj = "Cov"): covariance matrix

getCorr

signature(obj = "Cov"): correlation matrix

getData

signature(obj = "Cov"): data frame

getDistance

signature(obj = "Cov"): distances

getEvals

signature(obj = "Cov"): Computes and returns the eigenvalues of the covariance matrix

getDet

signature(obj = "Cov"): Computes and returns the determinant of the covariance matrix (or 0 if the covariance matrix is singular)

getShape

signature(obj = "Cov"): Computes and returns the shape matrix corresponding to the covariance matrix (i.e. the covariance matrix scaled to have determinant =1)

getFlag

signature(obj = "Cov"): Flags observations as outliers if the corresponding mahalanobis distance is larger then qchisq(prob, p) where prob defaults to 0.975.

isClassic

signature(obj = "Cov"): returns TRUE by default. If necessary, the robust classes will override

plot

signature(x = "Cov"): plot the object

show

signature(object = "Cov"): display the object

summary

signature(object = "Cov"): calculate summary information

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

showClass("Cov")

Description

Computes the classical estimates of multivariate location and scatter. Returns an S4 class CovClassic with the estimated center, cov, Mahalanobis distances and weights based on these distances.

Usage

CovClassic(x, unbiased=TRUE)
    Cov(x, unbiased=TRUE)

Arguments

Value

An object of class "CovClassic".

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

Cov-class, CovClassic-class

Examples

data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
cv <- CovClassic(hbk.x)
cv
summary(cv)
plot(cv)

Description

The class CovClassic represents an estimate of the multivariate location and scatter of a data set. The objects of class CovClassic contain the classical estimates.

Objects from the Class

Objects can be created by calls of the form new("CovClassic", ...), but the usual way of creating CovClassic objects is a call to the function CovClassic which serves as a constructor.

Slots

call:

Object of class "language"

cov:

covariance matrix

center:

location

n.obs:

number of observations used for the computation of the estimates

mah:

mahalanobis distances

method:

a character string describing the method used to compute the estimate: "Classic"

singularity:

a list with singularity information for the ocvariance matrix (or NULL of not singular)

X:

data

Methods

getCenter

signature(obj = "CovClassic"): location vector

getCov

signature(obj = "CovClassic"): covariance matrix

getCorr

signature(obj = "CovClassic"): correlation matrix

getData

signature(obj = "CovClassic"): data frame

getDistance

signature(obj = "CovClassic"): distances

getEvals

signature(obj = "CovClassic"): Computes and returns the eigenvalues of the covariance matrix

plot

signature(x = "CovClassic"): plot the object

show

signature(object = "CovClassic"): display the object

summary

signature(object = "CovClassic"): calculate summary information

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
cv <- CovClassic(hbk.x)
cv
summary(cv)
plot(cv)

Description

The class "CovControl" is a VIRTUAL base control class for the derived classes representing the control parameters for the different robust methods

Arguments

Objects from the Class

A virtual Class: No objects may be created from it.

Methods

No methods defined with class "CovControl" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Description

This function will create a control object CovControlMcd containing the control parameters for CovMcd

Usage

CovControlMcd(alpha = 0.5, nsamp = 500, scalefn=NULL, maxcsteps=200, 
seed = NULL, trace= FALSE, use.correction = TRUE)

Arguments

Value

A CovControlMcd object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMcd", alpha=0.75)
    ctrl2 <- CovControlMcd(alpha=0.75)

    data(hbk)
    CovMcd(hbk, control=ctrl1)

Description

This class extends the CovControl class and contains the control parameters for "CovMcd"

Objects from the Class

Objects can be created by calls of the form new("CovControlMcd", ...) or by calling the constructor-function CovControlMcd.

Slots

alpha:

numeric parameter controlling the size of the subsets over which the determinant is minimized, i.e., alpha*n observations are used for computing the determinant. Allowed values are between 0.5 and 1 and the default is 0.5.

nsamp

number of subsets used for initial estimates or "best", "exact" or "deterministic". Default is nsamp = 500. For nsamp="best" exhaustive enumeration is done, as long as the number of trials does not exceed 5000. For "exact", exhaustive enumeration will be attempted however many samples are needed. In this case a warning message will be displayed saying that the computation can take a very long time.

For "deterministic", the deterministic MCD is computed; as proposed by Hubert et al. (2012) it starts from the $h$ most central observations of six (deterministic) estimators.

scalefn

function to compute a robust scale estimate or character string specifying a rule determining such a function.

maxcsteps

maximal number of concentration steps in the deterministic MCD; should not be reached.

seed:

starting value for random generator. Default is seed = NULL

use.correction:

whether to use finite sample correction factors. Default is use.correction=TRUE.

trace, tolSolve:

from the "CovControl" class.

Extends

Class "CovControl", directly.

Methods

restimate

signature(obj = "CovControlMcd"): the generic function restimate allows the different methods for robust estimation to be used polymorphically - this function will call CovMcd passing it the control object and will return the obtained CovRobust object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMcd", alpha=0.75)
    ctrl2 <- CovControlMcd(alpha=0.75)

    data(hbk)
    CovMcd(hbk, control=ctrl1)

Description

This function will create a control object CovControlMest containing the control parameters for CovMest

Usage

CovControlMest(r = 0.45, arp = 0.05, eps = 0.001, maxiter = 120)

Arguments

Value

A CovControlMest object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMest", r=0.4)
    ctrl2 <- CovControlMest(r=0.4)

    data(hbk)
    CovMest(hbk, control=ctrl1)

Description

This class extends the CovControl class and contains the control parameters for CovMest

Objects from the Class

Objects can be created by calls of the form new("CovControlMest", ...) or by calling the constructor-function CovControlMest.

Slots

r:

a numeric value specifying the required breakdown point. Allowed values are between (n - p)/(2 * n) and 1 and the default is 0.45

arp:

a numeric value specifying the asympthotic rejection point, i.e. the fraction of points receiving zero weight (see Rocke (1996)). Default is 0.05

eps:

a numeric value specifying the relative precision of the solution of the M-estimate. Defaults to 1e-3

maxiter:

maximum number of iterations allowed in the computation of the M-estimate. Defaults to 120

trace, tolSolve:

from the "CovControl" class.

Extends

Class "CovControl", directly.

Methods

restimate

signature(obj = "CovControlMest"): the generic function restimate allowes the different methods for robust estimation to be used polymorphically - this function will call CovMest passing it the control object and will return the obtained CovRobust object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMest", r=0.4)
    ctrl2 <- CovControlMest(r=0.4)

    data(hbk)
    CovMest(hbk, control=ctrl1)

Description

This function will create a control object CovControlMMest containing the control parameters for CovMMest

Usage

CovControlMMest(bdp = 0.5, eff=0.95, maxiter = 50, sest=CovControlSest(),
        trace = FALSE, tolSolve = 1e-7)

Arguments

Value

A CovControlSest object.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMMest", bdp=0.25)
    ctrl2 <- CovControlMMest(bdp=0.25)

    data(hbk)
    CovMMest(hbk, control=ctrl1)

Description

This class extends the CovControl class and contains the control parameters for CovMMest

Objects from the Class

Objects can be created by calls of the form new("CovControlMMest", ...) or by calling the constructor-function CovControlMMest.

Slots

bdp

a numeric value specifying the required breakdown point. Allowed values are between 0.5 and 1 and the default is bdp=0.5.

eff

a numeric value specifying the required efficiency for the MM estimates. Default is eff=0.95.

sest

an CovControlSest object containing control parameters for the initial S-estimate.

maxiter

maximum number of iterations allowed in the computation of the MM-estimate. Default is maxiter=50.

trace, tolSolve:

from the "CovControl" class. tolSolve is used as a convergence tolerance for the MM-iteration.

Extends

Class "CovControl", directly.

Methods

restimate

signature(obj = "CovControlMMest"): the generic function restimate allowes the different methods for robust estimation to be used polymorphically - this function will call CovMMest passing it the control object and will return the obtained CovRobust object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMMest", bdp=0.25)
    ctrl2 <- CovControlMMest(bdp=0.25)

    data(hbk)
    CovMMest(hbk, control=ctrl1)

Description

This function will create a control object CovControlMrcd containing the control parameters for CovMrcd

Usage

CovControlMrcd(alpha = 0.5, h=NULL, maxcsteps=200, rho=NULL, 
    target=c("identity", "equicorrelation"), maxcond=50,
    trace= FALSE)

Arguments

Value

A CovControlMrcd object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMrcd", alpha=0.75)
    ctrl2 <- CovControlMrcd(alpha=0.75)

    data(hbk)
    CovMrcd(hbk, control=ctrl1)

Description

This class extends the CovControl class and contains the control parameters for "CovMrcd"

Objects from the Class

Objects can be created by calls of the form new("CovControlMrcd", ...) or by calling the constructor-function CovControlMrcd.

Slots

alpha:

numeric parameter controlling the size of the subsets over which the determinant is minimized, i.e., alpha*n observations are used for computing the determinant. Allowed values are between 0.5 and 1 and the default is 0.5.

h

the size of the subset (can be between ceiling(n/2) and n). Normally NULL and then it h will be calculated as h=ceiling(alpha*n). If h is provided, alpha will be calculated as alpha=h/n.

maxcsteps

maximal number of concentration steps in the deterministic MCD; should not be reached.

rho

regularization parameter. Normally NULL and will be estimated from the data.

target

structure of the robust positive definite target matrix: a) "identity": target matrix is diagonal matrix with robustly estimated univariate scales on the diagonal or b) "equicorrelation": non-diagonal target matrix that incorporates an equicorrelation structure (see (17) in paper).

maxcond

maximum condition number allowed (see step 3.4 in algorithm 1).

trace, tolSolve:

from the "CovControl" class.

Extends

Class "CovControl", directly.

Methods

restimate

signature(obj = "CovControlMrcd"): the generic function restimate allows the different methods for robust estimation to be used polymorphically - this function will call CovMrcd passing it the control object and will return the obtained CovRobust object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

"CovControlMcd"

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMrcd", alpha=0.75)
    ctrl2 <- CovControlMrcd(alpha=0.75)

    data(hbk)
    CovMrcd(hbk, control=ctrl1)

Description

This function will create a control object CovControlMve containing the control parameters for CovMve

Usage

CovControlMve(alpha = 0.5, nsamp = 500, seed = NULL, trace= FALSE)

Arguments

Value

A CovControlMve object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMve", alpha=0.75)
    ctrl2 <- CovControlMve(alpha=0.75)

    data(hbk)
    CovMve(hbk, control=ctrl1)

Description

This class extends the CovControl class and contains the control parameters for "CovMve"

Objects from the Class

Objects can be created by calls of the form new("CovControlMve", ...) or by calling the constructor-function CovControlMve.

Slots

alpha:

numeric parameter controlling the size of the subsets over which the determinant is minimized, i.e., alpha*n observations are used for computing the determinant. Allowed values are between 0.5 and 1 and the default is 0.5.

nsamp:

number of subsets used for initial estimates or "best" or "exact". Default is nsamp = 500. For nsamp="best" exhaustive enumeration is done, as long as the number of trials does not exceed 5000. For "exact", exhaustive enumeration will be attempted however many samples are needed. In this case a warning message will be displayed saying that the computation can take a very long time.

seed:

starting value for random generator. Default is seed = NULL

trace, tolSolve:

from the "CovControl" class.

Extends

Class "CovControl", directly.

Methods

restimate

signature(obj = "CovControlMve"): the generic function restimate allowes the different methods for robust estimation to be used polymorphically - this function will call CovMve passing it the control object and will return the obtained CovRobust object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlMve", alpha=0.75)
    ctrl2 <- CovControlMve(alpha=0.75)

    data(hbk)
    CovMve(hbk, control=ctrl1)

Description

This function will create a control object CovControlOgk containing the control parameters for CovOgk

Usage

CovControlOgk(niter = 2, beta = 0.9, mrob = NULL, 
vrob = .vrobGK, smrob = "scaleTau2", svrob = "gk")

Arguments

Details

If the user does not specify a scale and covariance function to be used in the computations or specifies one by using the arguments smrob and svrob (i.e. the names of the functions as strings), a native code written in C will be called which is by far faster than the R version.

If the arguments mrob and vrob are not NULL, the specified functions will be used via the pure R implementation of the algorithm. This could be quite slow.

Value

A CovControlOgk object

Author(s)

Valentin Todorov [email protected]

References

Maronna, R.A. and Zamar, R.H. (2002) Robust estimates of location and dispersion of high-dimensional datasets; Technometrics 44(4), 307–317.

Yohai, R.A. and Zamar, R.H. (1998) High breakdown point estimates of regression by means of the minimization of efficient scale JASA 86, 403–413.

Gnanadesikan, R. and John R. Kettenring (1972) Robust estimates, residuals, and outlier detection with multiresponse data. Biometrics 28, 81–124.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlOgk", beta=0.95)
    ctrl2 <- CovControlOgk(beta=0.95)

    data(hbk)
    CovOgk(hbk, control=ctrl1)

Description

This class extends the CovControl class and contains the control parameters for "CovOgk"

Objects from the Class

Objects can be created by calls of the form new("CovControlOgk", ...) or by calling the constructor-function CovControlOgk.

Slots

niter

number of iterations, usually 1 or 2 since iterations beyond the second do not lead to improvement.

beta

coverage parameter for the final reweighted estimate

mrob

function for computing the robust univariate location and dispersion - defaults to the tau scale defined in Yohai and Zamar (1998)

vrob

function for computing robust estimate of covariance between two random vectors - defaults the one proposed by Gnanadesikan and Kettenring (1972)

smrob

A string indicating the name of the function for computing the robust univariate location and dispersion - defaults to scaleTau2 - the scale 'tau' function defined in Yohai and Zamar (1998)

svrob

A string indicating the name of the function for computing robust estimate of covariance between two random vectors - defaults to gk, the one proposed by Gnanadesikan and Kettenring (1972).

trace, tolSolve:

from the "CovControl" class.

Extends

Class "CovControl", directly.

Methods

restimate

signature(obj = "CovControlOgk"): the generic function restimate allowes the different methods for robust estimation to be used polymorphically - this function will call CovOgk passing it the control object and will return the obtained CovRobust object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlOgk", beta=0.95)
    ctrl2 <- CovControlOgk(beta=0.95)

    data(hbk)
    CovOgk(hbk, control=ctrl1)

Description

This function will create a control object CovControlSde containing the control parameters for CovSde

Usage

CovControlSde(nsamp = 0, maxres = 0, tune = 0.95, eps = 0.5, prob = 0.99,
    seed = NULL, trace = FALSE, tolSolve = 1e-14)

Arguments

Value

A CovControlSde object.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlSde", nsamp=2000)
    ctrl2 <- CovControlSde(nsamp=2000)

    data(hbk)
    CovSde(hbk, control=ctrl1)

Description

This class extends the CovControl class and contains the control parameters for CovSde

Objects from the Class

Objects can be created by calls of the form new("CovControlSde", ...) or by calling the constructor-function CovControlSde.

Slots

nsamp

a positive integer giving the number of resamples required

maxres

a positive integer specifying the maximum number of resamples to be performed including those that are discarded due to linearly dependent subsamples.

tune

a numeric value between 0 and 1 giving the fraction of the data to receive non-zero weight. Default is tune = 0.95.

prob

a numeric value between 0 and 1 specifying the probability of high breakdown point; used to compute nsamp when nsamp is omitted. Default is prob = 0.99.

eps

a numeric value between 0 and 0.5 specifying the breakdown point; used to compute nsamp when nresamp is omitted. Default is eps = 0.5.

seed

starting value for random generator. Default is seed = NULL.

trace, tolSolve:

from the "CovControl" class.

Extends

Class "CovControl", directly.

Methods

restimate

signature(obj = "CovControlSde"): ...

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlSde", nsamp=2000)
    ctrl2 <- CovControlSde(nsamp=2000)

    data(hbk)
    CovSde(hbk, control=ctrl1)

Description

This function will create a control object CovControlSest containing the control parameters for CovSest

Usage

CovControlSest(bdp = 0.5, arp = 0.1, eps = 1e-5, maxiter = 120,
        nsamp = 500, seed = NULL, trace = FALSE, tolSolve = 1e-14, method= "sfast")

Arguments

Value

A CovControlSest object.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlSest", bdp=0.4)
    ctrl2 <- CovControlSest(bdp=0.4)

    data(hbk)
    CovSest(hbk, control=ctrl1)

Description

This class extends the CovControl class and contains the control parameters for CovSest

Objects from the Class

Objects can be created by calls of the form new("CovControlSest", ...) or by calling the constructor-function CovControlSest.

Slots

bdp

a numeric value specifying the required breakdown point. Allowed values are between (n - p)/(2 * n) and 1 and the default is bdp=0.45.

arp

a numeric value specifying the asympthotic rejection point (for the Rocke type S estimates), i.e. the fraction of points receiving zero weight (see Rocke (1996)). Default is arp=0.1.

eps

a numeric value specifying the relative precision of the solution of the S-estimate (bisquare and Rocke type). Default is to eps=1e-5.

maxiter

maximum number of iterations allowed in the computation of the S-estimate (bisquare and Rocke type). Default is maxiter=120.

nsamp

the number of random subsets considered. Default is nsamp = 500.

seed

starting value for random generator. Default is seed = NULL.

method

Which algorithm to use: 'sfast'=FAST-S, 'surreal'=Ruppert's SURREAL algorithm, 'bisquare'=Bisquare S-estimation with HBDP start or 'rocke' for Rocke type S-estimates

trace, tolSolve:

from the "CovControl" class.

Extends

Class "CovControl", directly.

Methods

restimate

signature(obj = "CovControlSest"): the generic function restimate allowes the different methods for robust estimation to be used polymorphically - this function will call CovSest passing it the control object and will return the obtained CovRobust object

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

## the following two statements are equivalent
    ctrl1 <- new("CovControlSest", bdp=0.4)
    ctrl2 <- CovControlSest(bdp=0.4)

    data(hbk)
    CovSest(hbk, control=ctrl1)

Description

Computes a robust multivariate location and scatter estimate with a high breakdown point, using the ‘Fast MCD’ (Minimum Covariance Determinant) estimator.

Usage

CovMcd(x,
       raw.only=FALSE, alpha=control@alpha, nsamp=control@nsamp,
       scalefn=control@scalefn, maxcsteps=control@maxcsteps,
       initHsets=NULL, save.hsets=FALSE,
       seed=control@seed, trace=control@trace,
       use.correction=control@use.correction,
       control=CovControlMcd(), ...)

Arguments

Details

This function computes the minimum covariance determinant estimator of location and scatter and returns an S4 object of class CovMcd-class containing the estimates. The implementation of the function is similar to the existing R function covMcd() which returns an S3 object. The MCD method looks for the $h (> n/2)$ observations (out of $n$) whose classical covariance matrix has the lowest possible determinant. The raw MCD estimate of location is then the average of these $h$ points, whereas the raw MCD estimate of scatter is their covariance matrix, multiplied by a consistency factor and a finite sample correction factor (to make it consistent at the normal model and unbiased at small samples). Both rescaling factors are returned also in the vector raw.cnp2 of length 2. Based on these raw MCD estimates, a reweighting step is performed which increases the finite-sample efficiency considerably - see Pison et al. (2002). The rescaling factors for the reweighted estimates are returned in the vector cnp2 of length 2. Details for the computation of the finite sample correction factors can be found in Pison et al. (2002). The finite sample corrections can be suppressed by setting use.correction=FALSE. The implementation in rrcov uses the Fast MCD algorithm of Rousseeuw and Van Driessen (1999) to approximate the minimum covariance determinant estimator.

Value

An S4 object of class CovMcd-class which is a subclass of the virtual class CovRobust-class.

Author(s)

Valentin Todorov [email protected]

References

P. J. Rousseeuw and A. M. Leroy (1987) Robust Regression and Outlier Detection. Wiley.

P. J. Rousseeuw and K. van Driessen (1999) A fast algorithm for the minimum covariance determinant estimator. Technometrics 41, 212–223.

M. Hubert, P. Rousseeuw and T. Verdonck (2012) A deterministic algorithm for robust location and scatter. Journal of Computational and Graphical Statistics 21(3), 618–637.

Pison, G., Van Aelst, S., and Willems, G. (2002), Small Sample Corrections for LTS and MCD, Metrika, 55, 111-123.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

cov.rob from package MASS

Examples

data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
CovMcd(hbk.x)
cD <- CovMcd(hbk.x, nsamp = "deterministic")
summary(cD)

## the following three statements are equivalent
c1 <- CovMcd(hbk.x, alpha = 0.75)
c2 <- CovMcd(hbk.x, control = CovControlMcd(alpha = 0.75))
## direct specification overrides control one:
c3 <- CovMcd(hbk.x, alpha = 0.75,
             control = CovControlMcd(alpha=0.95))
c1

Description

This class, derived from the virtual class "CovRobust" accomodates MCD Estimates of multivariate location and scatter computed by the ‘Fast MCD’ algorithm.

Objects from the Class

Objects can be created by calls of the form new("CovMcd", ...), but the usual way of creating CovMcd objects is a call to the function CovMcd which serves as a constructor.

Slots

alpha:

Object of class "numeric" - the size of the subsets over which the determinant is minimized (the default is (n+p+1)/2)

quan:

Object of class "numeric" - the number of observations on which the MCD is based. If quan equals n.obs, the MCD is the classical covariance matrix.

best:

Object of class "Uvector" - the best subset found and used for computing the raw estimates. The size of best is equal to quan

raw.cov:

Object of class "matrix" the raw (not reweighted) estimate of location

raw.center:

Object of class "vector" - the raw (not reweighted) estimate of scatter

raw.mah:

Object of class "Uvector" - mahalanobis distances of the observations based on the raw estimate of the location and scatter

raw.wt:

Object of class "Uvector" - weights of the observations based on the raw estimate of the location and scatter

raw.cnp2:

Object of class "numeric" - a vector of length two containing the consistency correction factor and the finite sample correction factor of the raw estimate of the covariance matrix

cnp2:

Object of class "numeric" - a vector of length two containing the consistency correction factor and the finite sample correction factor of the final estimate of the covariance matrix.

iter, crit, wt:

from the "CovRobust" class.

call, cov, center, n.obs, mah, method, singularity, X:

from the "Cov" class.

Extends

Class "CovRobust", directly. Class "Cov", by class "CovRobust".

Methods

No methods defined with class "CovMcd" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovMcd, Cov-class, CovRobust-class

Examples

showClass("CovMcd")

Description

Computes constrained M-Estimates of multivariate location and scatter based on the translated biweight function (‘t-biweight’) using a High breakdown point initial estimate. The default initial estimate is the Minimum Volume Ellipsoid computed with CovMve. The raw (not reweighted) estimates are taken and the covariance matrix is standardized to determinant 1.

Usage

covMest(x, cor=FALSE, r = 0.45, arp = 0.05, eps=1e-3,
    maxiter=120, control, t0, S0)

Arguments

.

Details

Rocke (1996) has shown that the S-estimates of multivariate location and scatter in high dimensions can be sensitive to outliers even if the breakdown point is set to be near 0.5. To mitigate this problem he proposed to utilize the translated biweight (or t-biweight) method with a standardization step consisting of equating the median of rho(d) with the median under normality. This is then not an S-estimate, but is instead a constrained M-estimate. In order to make the smooth estimators to work, a reasonable starting point is necessary, which will lead reliably to a good solution of the estimator. In covMest the MVE computed by CovMve is used, but the user has the possibility to give her own initial estimates.

Value

An object of class "mest" which is basically a list with the following components. This class is "derived" from "mcd" so that the same generic functions - print, plot, summary - can be used. NOTE: this is going to change - in one of the next revisions covMest will return an S4 class "mest" which is derived (i.e. contains) form class "cov".

Note

The psi, rho and weight functions for the M estimation are encapsulated in a virtual S4 class PsiFun from which a PsiBwt class, implementing the translated biweight (t-biweight), is dervied. The base class PsiFun contains also the M-iteration itself. Although not documented and not accessibale directly by the user these classes will form the bases for adding other functions (biweight, LWS, etc.) as well as S-estimates.

Author(s)

Valentin Todorov [email protected],

(some code from C. Becker - http://www.sfb475.uni-dortmund.de/dienst/de/content/struk-d/bereicha-d/tpa1softw-d.html)

References

D.L.Woodruff and D.M.Rocke (1994) Computable robust estimation of multivariate location and shape on high dimension using compound estimators, Journal of the American Statistical Association, 89, 888–896.

D.M.Rocke (1996) Robustness properties of S-estimates of multivariate location and shape in high dimension, Annals of Statistics, 24, 1327-1345.

D.M.Rocke and D.L.Woodruff (1996) Identification of outliers in multivariate data Journal of the American Statistical Association, 91, 1047–1061.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

covMcd

Description

Computes constrained M-Estimates of multivariate location and scatter based on the translated biweight function (‘t-biweight’) using a High breakdown point initial estimate as defined by Rocke (1996). The default initial estimate is the Minimum Volume Ellipsoid computed with CovMve. The raw (not reweighted) estimates are taken and the covariance matrix is standardized to determinant 1.

Usage

CovMest(x, r = 0.45, arp = 0.05, eps=1e-3,
        maxiter=120, control, t0, S0, initcontrol)

Arguments

Details

Rocke (1996) has shown that the S-estimates of multivariate location and scatter in high dimensions can be sensitive to outliers even if the breakdown point is set to be near 0.5. To mitigate this problem he proposed to utilize the translated biweight (or t-biweight) method with a standardization step consisting of equating the median of rho(d) with the median under normality. This is then not an S-estimate, but is instead a constrained M-estimate. In order to make the smooth estimators to work, a reasonable starting point is necessary, which will lead reliably to a good solution of the estimator. In CovMest the MVE computed by CovMve is used, but the user has the possibility to give her own initial estimates.

Value

An object of class CovMest-class which is a subclass of the virtual class CovRobust-class.

Note

The psi, rho and weight functions for the M estimation are encapsulated in a virtual S4 class PsiFun from which a PsiBwt class, implementing the translated biweight (t-biweight), is dervied. The base class PsiFun contains also the M-iteration itself. Although not documented and not accessibale directly by the user these classes will form the bases for adding other functions (biweight, LWS, etc.) as well as S-estimates.

Author(s)

Valentin Todorov [email protected],

(some code from C. Becker - http://www.sfb475.uni-dortmund.de/dienst/de/content/struk-d/bereicha-d/tpa1softw-d.html)

References

D.L.Woodruff and D.M.Rocke (1994) Computable robust estimation of multivariate location and shape on high dimension using compound estimators, Journal of the American Statistical Association, 89, 888–896.

D.M.Rocke (1996) Robustness properties of S-estimates of multivariate location and shape in high dimension, Annals of Statistics, 24, 1327-1345.

D.M.Rocke and D.L.Woodruff (1996) Identification of outliers in multivariate data Journal of the American Statistical Association, 91, 1047–1061.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

covMcd, Cov-class, CovMve, CovRobust-class, CovMest-class

Examples

library(rrcov)
data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
CovMest(hbk.x)

## the following four statements are equivalent
c0 <- CovMest(hbk.x)
c1 <- CovMest(hbk.x, r = 0.45)
c2 <- CovMest(hbk.x, control = CovControlMest(r = 0.45))
c3 <- CovMest(hbk.x, control = new("CovControlMest", r = 0.45))

## direct specification overrides control one:
c4 <- CovMest(hbk.x, r = 0.40,
             control = CovControlMest(r = 0.25))
c1
summary(c1)
plot(c1)

Description

This class, derived from the virtual class "CovRobust" accomodates constrained M-Estimates of multivariate location and scatter based on the translated biweight function (‘t-biweight’) using a High breakdown point initial estimate (Minimum Covariance Determinant - ‘Fast MCD’)

Objects from the Class

Objects can be created by calls of the form new("CovMest", ...), but the usual way of creating CovMest objects is a call to the function CovMest which serves as a constructor.

Slots

vt:

Object of class "vector" - vector of weights (v)

iter, crit, wt:

from the "CovRobust" class.

call, cov, center, n.obs, mah, method, singularity, X:

from the "Cov" class.

Extends

Class "CovRobust", directly. Class "Cov", by class "CovRobust".

Methods

No methods defined with class "CovMest" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovMest, Cov-class, CovRobust-class

Examples

showClass("CovMest")

Description

Computes MM-Estimates of multivariate location and scatter starting from an initial S-estimate

Usage

CovMMest(x, bdp = 0.5, eff = 0.95, eff.shape=TRUE, maxiter = 50, 
        trace = FALSE, tolSolve = 1e-7, control)

Arguments

Details

Computes MM-estimates of multivariate location and scatter starting from an initial S-estimate.

Value

An S4 object of class CovMMest-class which is a subclass of the virtual class CovRobust-class.

Author(s)

Valentin Todorov [email protected]

References

Tatsuoka, K.S. and Tyler, D.E. (2000). The uniqueness of S and M-functionals under non-elliptical distributions. Annals of Statistics 28, 1219–1243

M. Salibian-Barrera, S. Van Aelstt and G. Willems (2006). Principal components analysis based on multivariate MM-estimators with fast and robust bootstrap. Journal of the American Statistical Association 101, 1198–1211.

R. A. Maronna, D. Martin and V. Yohai (2006). Robust Statistics: Theory and Methods. Wiley, New York.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

library(rrcov)
data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
CovMMest(hbk.x)

## the following four statements are equivalent
c0 <- CovMMest(hbk.x)
c1 <- CovMMest(hbk.x, bdp = 0.25)
c2 <- CovMMest(hbk.x, control = CovControlMMest(bdp = 0.25))
c3 <- CovMMest(hbk.x, control = new("CovControlMMest", bdp = 0.25))

## direct specification overrides control one:
c4 <- CovMMest(hbk.x, bdp = 0.40,
             control = CovControlMMest(bdp = 0.25))
c1
summary(c1)
plot(c1)

## Deterministic MM-estmates
CovMMest(hbk.x, control=CovControlMMest(sest=CovControlSest(method="sdet")))

Description

This class, derived from the virtual class "CovRobust" accomodates MM Estimates of multivariate location and scatter.

Objects from the Class

Objects can be created by calls of the form new("CovMMest", ...), but the usual way of creating CovSest objects is a call to the function CovMMest which serves as a constructor.

Slots

det, flag, iter, crit:

from the "CovRobust" class.

c1

tuning parameter of the loss function for MM-estimation (depend on control parameters eff and eff.shape). Can be computed by the internal function .csolve.bw.MM(p, eff, eff.shape=TRUE). For the tuning parameters of the underlying S-estimate see the slot sest and "CovSest".

sest

an CovSest object containing the initial S-estimate.

call, cov, center, n.obs, mah, method, singularity, X:

from the "Cov" class.

Extends

Class "CovRobust", directly. Class "Cov", by class "CovRobust".

Methods

No methods defined with class "CovMMest" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovMMest, Cov-class, CovRobust-class

Examples

showClass("CovMMest")

Description

Computes a robust multivariate location and scatter estimate with a high breakdown point, using the Minimum Regularized Covariance Determonant (MRCD) estimator.

Usage

CovMrcd(x,
       alpha=control@alpha, 
       h=control@h,
       maxcsteps=control@maxcsteps,
       initHsets=NULL, save.hsets=FALSE,
       rho=control@rho,
       target=control@target,
       maxcond=control@maxcond,
       trace=control@trace,
       control=CovControlMrcd())

Arguments

Details

This function computes the minimum regularized covariance determinant estimator (MRCD) of location and scatter and returns an S4 object of class CovMrcd-class containing the estimates. Similarly like the MCD method, MRCD looks for the $h (> n/2)$ observations (out of $n$) whose classical covariance matrix has the lowest possible determinant, but replaces the subset-based covariance by a regularized covariance estimate, defined as a weighted average of the sample covariance of the h-subset and a predetermined positive definite target matrix. The Minimum Regularized Covariance Determinant (MRCD) estimator is then the regularized covariance based on the h-subset which makes the overall determinant the smallest. A data-driven procedure sets the weight of the target matrix (rho), so that the regularization is only used when needed.

Value

An S4 object of class CovMrcd-class which is a subclass of the virtual class CovRobust-class.

Author(s)

Kris Boudt, Peter Rousseeuw, Steven Vanduffel and Tim Verdonk. Improved by Joachim Schreurs and Iwein Vranckx. Adapted for rrcov by Valentin Todorov [email protected]

References

Kris Boudt, Peter Rousseeuw, Steven Vanduffel and Tim Verdonck (2020) The Minimum Regularized Covariance Determinant estimator, Statistics and Computing, 30, pp 113–128 doi:10.1007/s11222-019-09869-x.

Mia Hubert, Peter Rousseeuw and Tim Verdonck (2012) A deterministic algorithm for robust location and scatter. Journal of Computational and Graphical Statistics 21(3), 618–637.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovMcd

Examples

## 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

## Not run: 

##  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)

## End(Not run)

Description

This class, derived from the virtual class "CovRobust" accomodates MRCD Estimates of multivariate location and scatter computed by a variant of the ‘Fast MCD’ algorithm.

Objects from the Class

Objects can be created by calls of the form new("CovMrcd", ...), but the usual way of creating CovMrcd objects is a call to the function CovMrcd which serves as a constructor.

Slots

alpha:

Object of class "numeric" - the size of the subsets over which the determinant is minimized (the default is (n+p+1)/2)

quan:

Object of class "numeric" - the number of observations on which the MCD is based. If quan equals n.obs, the MCD is the classical covariance matrix.

best:

Object of class "Uvector" - the best subset found and used for computing the raw estimates. The size of best is equal to quan

cnp2:

Object of class "numeric" - containing the consistency correction factor of the estimate of the covariance matrix.

icov:

The inverse of the covariance matrix.

rho:

The estimated regularization parameter.

target:

The estimated target matrix.

crit:

from the "CovRobust" class.

call, cov, center, n.obs, mah, method, X:

from the "Cov" class.

Extends

Class "CovRobust", directly. Class "Cov", by class "CovRobust".

Methods

No methods defined with class "CovMrcd" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovMrcd, Cov-class, CovRobust-class, CovMcd-class

Examples

showClass("CovMrcd")

Description

Computes a robust multivariate location and scatter estimate with a high breakdown point, using the ‘MVE’ (Minimum Volume Ellipsoid) estimator.

Usage

CovMve(x, alpha = 1/2, nsamp = 500, seed = NULL, trace = FALSE, control)

Arguments

Details

This function computes the minimum volume ellipsoid estimator of location and scatter and returns an S4 object of class CovMve-class containing the estimates.

The approximate estimate is based on a subset of size alpha*n with an enclosing ellipsoid of smallest volume. The mean of the best found subset provides the raw estimate of the location, and the rescaled covariance matrix is the raw estimate of scatter. The rescaling of the raw covariance matrix is by median(dist)/qchisq(0.5, p) and this scale factor is returned in the slot raw.cnp2. Currently no finite sample corrction factor is applied. The Mahalanobis distances of all observations from the location estimate for the raw covariance matrix are calculated, and those points within the 97.5 under Gaussian assumptions are declared to be good. The final (reweightd) estimates are the mean and rescaled covariance of the good points. The reweighted covariance matrix is rescaled by 1/pgamma(qchisq(alpha, p)/2, p/2 + 1)/alpha (see Croux and Haesbroeck, 1999) and this scale factor is returned in the slot cnp2.

The search for the approximate solution is made over ellipsoids determined by the covariance matrix of p+1 of the data points and applying a simple but effective improvement of the subsampling procedure as described in Maronna et al. (2006), p. 198. Although there exists no formal proof of this improvement (as for MCD and LTS), simulations show that it can be recommended as an approximation of the MVE.

Value

An S4 object of class CovMve-class which is a subclass of the virtual class CovRobust-class.

Note

Main reason for implementing the MVE estimate was that it is the recommended initial estimate for S estimation (see Maronna et al. (2006), p. 199) and will be used by default in CovMest (after removing the correction factors from the covariance matrix and rescaling to determinant 1).

Author(s)

Valentin Todorov [email protected] and Matias Salibian-Barrera [email protected]

References

P. J. Rousseeuw and A. M. Leroy (1987) Robust Regression and Outlier Detection. Wiley.

C. Croux and G. Haesbroeck (1999). Influence function and efficiency of the minimum covariance determinant scatter matrix estimator. Journal of Multivariate Analysis, 71, 161–190.

R. A. Maronna, D. Martin and V. Yohai (2006). Robust Statistics: Theory and Methods. Wiley, New York.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

cov.rob from package MASS

Examples

data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
CovMve(hbk.x)

## the following three statements are equivalent
c1 <- CovMve(hbk.x, alpha = 0.75)
c2 <- CovMve(hbk.x, control = CovControlMve(alpha = 0.75))
## direct specification overrides control one:
c3 <- CovMve(hbk.x, alpha = 0.75,
             control = CovControlMve(alpha=0.95))
c1

Description

This class, derived from the virtual class "CovRobust" accomodates MVE Estimates of multivariate location and scatter computed by the ‘Fast MVE’ algorithm.

Objects from the Class

Objects can be created by calls of the form new("CovMve", ...), but the usual way of creating CovMve objects is a call to the function CovMve which serves as a constructor.

Slots

alpha:

Object of class "numeric" - the size of the subsets over which the volume of the ellipsoid is minimized (the default is (n+p+1)/2)

quan:

Object of class "numeric" - the number of observations on which the MVE is based. If quan equals n.obs, the MVE is the classical covariance matrix.

best:

Object of class "Uvector" - the best subset found and used for computing the raw estimates. The size of best is equal to quan

raw.cov:

Object of class "matrix" the raw (not reweighted) estimate of location

raw.center:

Object of class "vector" - the raw (not reweighted) estimate of scatter

raw.mah:

Object of class "Uvector" - mahalanobis distances of the observations based on the raw estimate of the location and scatter

raw.wt:

Object of class "Uvector" - weights of the observations based on the raw estimate of the location and scatter

raw.cnp2:

Object of class "numeric" - a vector of length two containing the consistency correction factor and the finite sample correction factor of the raw estimate of the covariance matrix

cnp2:

Object of class "numeric" - a vector of length two containing the consistency correction factor and the finite sample correction factor of the final estimate of the covariance matrix.

iter, crit, wt:

from the "CovRobust" class.

call, cov, center, n.obs, mah, method, singularity, X:

from the "Cov" class.

Extends

Class "CovRobust", directly. Class "Cov", by class "CovRobust".

Methods

No methods defined with class "CovMve" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovMve, Cov-class, CovRobust-class

Examples

showClass("CovMve")

Description

Computes a robust multivariate location and scatter estimate with a high breakdown point, using the pairwise algorithm proposed by Marona and Zamar (2002) which in turn is based on the pairwise robust estimator proposed by Gnanadesikan-Kettenring (1972).

Usage

CovOgk(x, niter = 2, beta = 0.9, control)

Arguments

Details

The method proposed by Marona and Zamar (2002) allowes to obtain positive-definite and almost affine equivariant robust scatter matrices starting from any pairwise robust scatter matrix. The default robust estimate of covariance between two random vectors used is the one proposed by Gnanadesikan and Kettenring (1972) but the user can choose any other method by redefining the function in slot vrob of the control object CovControlOgk. Similarly, the function for computing the robust univariate location and dispersion used is the tau scale defined in Yohai and Zamar (1998) but it can be redefined in the control object.

The estimates obtained by the OGK method, similarly as in CovMcd are returned as 'raw' estimates. To improve the estimates a reweighting step is performed using the coverage parameter beta and these reweighted estimates are returned as 'final' estimates.

Value

An S4 object of class CovOgk-class which is a subclass of the virtual class CovRobust-class.

Note

If the user does not specify a scale and covariance function to be used in the computations or specifies one by using the arguments smrob and svrob (i.e. the names of the functions as strings), a native code written in C will be called which is by far faster than the R version.

If the arguments mrob and vrob are not NULL, the specified functions will be used via the pure R implementation of the algorithm. This could be quite slow.

See CovControlOgk for details.

Author(s)

Valentin Todorov [email protected] and Kjell Konis [email protected]

References

Maronna, R.A. and Zamar, R.H. (2002) Robust estimates of location and dispersion of high-dimensional datasets; Technometrics 44(4), 307–317.

Yohai, R.A. and Zamar, R.H. (1998) High breakdown point estimates of regression by means of the minimization of efficient scale JASA 86, 403–413.

Gnanadesikan, R. and John R. Kettenring (1972) Robust estimates, residuals, and outlier detection with multiresponse data. Biometrics 28, 81–124.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovMcd, CovMest

Examples

data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
CovOgk(hbk.x)

## the following three statements are equivalent
c1 <- CovOgk(hbk.x, niter=1)
c2 <- CovOgk(hbk.x, control = CovControlOgk(niter=1))

## direct specification overrides control one:
c3 <- CovOgk(hbk.x, beta=0.95,
             control = CovControlOgk(beta=0.99))
c1

x<-matrix(c(1,2,3,7,1,2,3,7), ncol=2)
##  CovOgk(x)   - this would fail because the two columns of x are exactly collinear.
##              In order to fix it, redefine the default 'vrob' function for example
##              in the following way and pass it as a parameter in the control
##              object.
cc <- CovOgk(x, control=new("CovControlOgk",
                            vrob=function(x1, x2, ...)
                            {
                                r <- .vrobGK(x1, x2, ...)
                                if(is.na(r))
                                    r <- 0
                                r
                            })
)
cc

Description

This class, derived from the virtual class "CovRobust" accomodates OGK Estimates of multivariate location and scatter computed by the algorithm proposed by Marona and Zamar (2002).

Objects from the Class

Objects can be created by calls of the form new("CovOgk", ...), but the usual way of creating CovOgk objects is a call to the function CovOgk which serves as a constructor.

Slots

raw.cov:

Object of class "matrix" the raw (not reweighted) estimate of covariance matrix

raw.center:

Object of class "vector" - the raw (not reweighted) estimate of the location vector

raw.mah:

Object of class "Uvector" - mahalanobis distances of the observations based on the raw estimate of the location and scatter

raw.wt:

Object of class "Uvector" - weights of the observations based on the raw estimate of the location and scatter

iter, crit, wt:

from the "CovRobust" class.

call, cov, center, n.obs, mah, method, singularity, X:

from the "Cov" class.

Extends

Class "CovRobust", directly. Class "Cov", by class "CovRobust".

Methods

No methods defined with class "CovOgk" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovMcd-class, CovMest-class

Examples

showClass("CovOgk")

Description

Computes a robust multivariate location and scatter estimate with a high breakdown point, using one of the available estimators.

Usage

CovRobust(x, control, na.action = na.fail)

Arguments

Details

This function simply calls the restimate method of the control object control. If a character string naming an estimator is specified, a new control object will be created and used (with default estimation options). If this argument is missing or a character string 'auto' is specified, the function will select the robust estimator according to the size of the dataset. If there are less than 1000 observations and less than 10 variables or less than 5000 observations and less than 5 variables, Stahel-Donoho estimator will be used. Otherwise, if there are less than 50000 observations either bisquare S-estimates (for less than 10 variables) or Rocke type S-estimates (for 10 to 20 variables) will be used. In both cases the S iteration starts at the initial MVE estimate. And finally, if there are more than 50000 observations and/or more than 20 variables the Orthogonalized Quadrant Correlation estimator (CovOgk with the corresponding parameters) is used.

Value

An object derived from a CovRobust object, depending on the selected estimator.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
CovRobust(hbk.x)
CovRobust(hbk.x, CovControlSest(method="bisquare"))

Description

CovRobust is a virtual base class used for deriving the concrete classes representing different robust estimates of multivariate location and scatter. Here are implemeted the standard methods common for all robust estimates like show, summary and plot. The derived classes can override these methods and can define new ones.

Objects from the Class

A virtual Class: No objects may be created from it.

Slots

iter:

number of iterations used to compute the estimates

crit:

value of the criterion function

wt:

weights

call, cov, center, n.obs, mah, method, singularity, X:

from the "Cov" class.

Extends

Class "Cov", directly.

Methods

isClassic

signature(obj = "CovRobust"): Will return FALSE, since this is a 'Robust' object

getMeth

signature(obj = "CovRobust"): Return the name of the particular robust method used (as a character string)

show

signature(object = "CovRobust"): display the object

plot

signature(x = "CovRobust"): plot the object

getRaw

signature(obj = "CovRobust"): Return the object with the reweighted estimates replaced by the raw ones (only relevant for CovMcd, CovMve and CovOgk)

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

Cov-class, CovMcd-class, CovMest-class, CovOgk-class

Examples

data(hbk)
     hbk.x <- data.matrix(hbk[, 1:3])
     cv <- CovMest(hbk.x)               # it is not possible to create an object of
                                        # class CovRobust, since it is a VIRTUAL class
     cv
     summary(cv)                        # summary method for class CovRobust
     plot(cv)                           # plot method for class CovRobust

Description

Compute a robust estimate of location and scale using the Stahel-Donoho projection based estimator

Usage

CovSde(x, nsamp, maxres, tune = 0.95, eps = 0.5, prob = 0.99, 
seed = NULL, trace = FALSE, control)

Arguments

Details

The projection based Stahel-Donoho estimator posses very good statistical properties, but it can be very slow if the number of variables is too large. It is recommended to use this estimator if n <= 1000 and p<=10 or n <= 5000 and p<=5. The number of subsamples required is calculated to provide a breakdown point of eps with probability prob and can reach values larger than the larger integer value - in such case it is limited to .Machine$integer.max. Of course you could provide nsamp in the call, i.e. nsamp=1000 but this will not guarantee the required breakdown point of th eestimator. For larger data sets it is better to use CovMcd or CovOgk. If you use CovRobust, the estimator will be selected automatically according on the size of the data set.

Value

An S4 object of class CovSde-class which is a subclass of the virtual class CovRobust-class.

Note

The Fortran code for the Stahel-Donoho method was taken almost with no changes from package robust which in turn has it from the Insightful Robust Library (thanks to by Kjell Konis).

Author(s)

Valentin Todorov [email protected] and Kjell Konis [email protected]

References

R. A. Maronna and V.J. Yohai (1995) The Behavior of the Stahel-Donoho Robust Multivariate Estimator. Journal of the American Statistical Association 90 (429), 330–341.

R. A. Maronna, D. Martin and V. Yohai (2006). Robust Statistics: Theory and Methods. Wiley, New York.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
CovSde(hbk.x)

## the following four statements are equivalent
c0 <- CovSde(hbk.x)
c1 <- CovSde(hbk.x, nsamp=2000)
c2 <- CovSde(hbk.x, control = CovControlSde(nsamp=2000))
c3 <- CovSde(hbk.x, control = new("CovControlSde", nsamp=2000))

## direct specification overrides control one:
c4 <- CovSde(hbk.x, nsamp=100,
             control = CovControlSde(nsamp=2000))
c1
summary(c1)
plot(c1)

## Use the function CovRobust() - if no estimation method is
##  specified, for small data sets CovSde() will be called
cr <- CovRobust(hbk.x)
cr

Description

This class, derived from the virtual class "CovRobust" accomodates Stahel-Donoho estimates of multivariate location and scatter.

Objects from the Class

Objects can be created by calls of the form new("CovSde", ...), but the usual way of creating CovSde objects is a call to the function CovSde which serves as a constructor.

Slots

iter, crit, wt:

from the "CovRobust" class.

call, cov, center, n.obs, mah, method, singularity, X:

from the "Cov" class.

Extends

Class "CovRobust", directly. Class "Cov", by class "CovRobust".

Methods

No methods defined with class "CovSde" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovSde, Cov-class, CovRobust-class

Examples

showClass("CovSde")

Description

Computes S-Estimates of multivariate location and scatter based on Tukey's biweight function using a fast algorithm similar to the one proposed by Salibian-Barrera and Yohai (2006) for the case of regression. Alternativley, the Ruppert's SURREAL algorithm, bisquare or Rocke type estimation can be used.

Usage

CovSest(x, bdp = 0.5, arp = 0.1, eps = 1e-5, maxiter = 120,
        nsamp = 500, seed = NULL, trace = FALSE, tolSolve = 1e-14,
        scalefn, maxisteps=200, 
        initHsets = NULL, save.hsets = missing(initHsets) || is.null(initHsets),
        method = c("sfast", "surreal", "bisquare", "rocke", "suser", "sdet"), 
        control, t0, S0, initcontrol)

Arguments

Details

Computes multivariate S-estimator of location and scatter. The computation will be performed by one of the following algorithms:

FAST-S

An algorithm similar to the one proposed by Salibian-Barrera and Yohai (2006) for the case of regression

SURREAL

Ruppert's SURREAL algorithm when method is set to 'surreal'

BISQUARE

Bisquare S-Estimate with method set to 'bisquare'

ROCKE

Rocke type S-Estimate with method set to 'rocke'

Except for the last algorithm, ROCKE, all other use Tukey biweight loss function. The tuning parameters used in the loss function (as determined by bdp) are returned in the slots cc and kp of the result object. They can be computed by the internal function .csolve.bw.S(bdp, p).

Value

An S4 object of class CovSest-class which is a subclass of the virtual class CovRobust-class.

Author(s)

Valentin Todorov [email protected], Matias Salibian-Barrera [email protected] and Victor Yohai [email protected]. See also the code from Kristel Joossens, K.U. Leuven, Belgium and Ella Roelant, Ghent University, Belgium.

References

M. Hubert, P. Rousseeuw and T. Verdonck (2012) A deterministic algorithm for robust location and scatter. Journal of Computational and Graphical Statistics 21(3), 618–637.

M. Hubert, P. Rousseeuw, D. Vanpaemel and T. Verdonck (2015) The DetS and DetMM estimators for multivariate location and scatter. Computational Statistics and Data Analysis 81, 64–75.

H.P. Lopuhaä (1989) On the Relation between S-estimators and M-estimators of Multivariate Location and Covariance. Annals of Statistics 17 1662–1683.

D. Ruppert (1992) Computing S Estimators for Regression and Multivariate Location/Dispersion. Journal of Computational and Graphical Statistics 1 253–270.

M. Salibian-Barrera and V. Yohai (2006) A fast algorithm for S-regression estimates, Journal of Computational and Graphical Statistics, 15, 414–427.

R. A. Maronna, D. Martin and V. Yohai (2006). Robust Statistics: Theory and Methods. Wiley, New York.

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

Examples

library(rrcov)
data(hbk)
hbk.x <- data.matrix(hbk[, 1:3])
cc <- CovSest(hbk.x)
cc

## summry and different types of plots
summary(cc)                         
plot(cc)                            
plot(cc, which="dd")
plot(cc, which="pairs")
plot(cc, which="xydist")

## the following four statements are equivalent
c0 <- CovSest(hbk.x)
c1 <- CovSest(hbk.x, bdp = 0.25)
c2 <- CovSest(hbk.x, control = CovControlSest(bdp = 0.25))
c3 <- CovSest(hbk.x, control = new("CovControlSest", bdp = 0.25))

## direct specification overrides control one:
c4 <- CovSest(hbk.x, bdp = 0.40,
             control = CovControlSest(bdp = 0.25))
c1
summary(c1)
plot(c1)

## Use the SURREAL algorithm of Ruppert
cr <- CovSest(hbk.x, method="surreal")
cr

## Use Bisquare estimation
cr <- CovSest(hbk.x, method="bisquare")
cr

## Use Rocke type estimation
cr <- CovSest(hbk.x, method="rocke")
cr

## Use Deterministic estimation
cr <- CovSest(hbk.x, method="sdet")
cr

Description

This class, derived from the virtual class "CovRobust" accomodates S Estimates of multivariate location and scatter computed by the ‘Fast S’ or ‘SURREAL’ algorithm.

Objects from the Class

Objects can be created by calls of the form new("CovSest", ...), but the usual way of creating CovSest objects is a call to the function CovSest which serves as a constructor.

Slots

iter, crit, wt:

from the "CovRobust" class.

iBest, nsteps, initHsets:

parameters for deterministic S-estimator (the best initial subset, number of concentration steps to convergence for each of the initial subsets, and the computed initial subsets, respectively).

cc, kp

tuning parameters used in Tukey biweight loss function, as determined by bdp. Can be computed by the internal function .csolve.bw.S(bdp, p).

call, cov, center, n.obs, mah, method, singularity, X:

from the "Cov" class.

Extends

Class "CovRobust", directly. Class "Cov", by class "CovRobust".

Methods

No methods defined with class "CovSest" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

CovSest, Cov-class, CovRobust-class

Examples

showClass("CovSest")

Description

The data set contains five measurements made on 145 non-obese adult patients classified into three groups.

The three primary variables are glucose intolerance (area under the straight line connecting glucose levels), insulin response to oral glucose (area under the straight line connecting insulin levels) and insulin resistance (measured by the steady state plasma glucose (SSPG) determined after chemical suppression of endogenous insulin secretion). Two additional variables, the relative weight and fasting plasma glucose, are also included.

Reaven and Miller, following Friedman and Rubin (1967), applied cluster analysis to the three primary variables and identified three clusters: "normal", "chemical diabetic", and "overt diabetic" subjects. The column group contains the classifications of the subjects into these three groups, obtained by current medical criteria.

Usage

data(diabetes)

Format

A data frame with the following variables:

rw

relative weight, expressed as the ratio of actual weight to expected weight, given the person's height.

fpg

fasting plasma glucose level.

glucose

area under plasma glucose curve after a three hour oral glucose tolerance test (OGTT).

insulin

area under plasma insulin curve after a three hour oral glucose tolerance test (OGTT).

sspg

Steady state plasma glucose, a measure of insulin resistance.

group

the type of diabetes: a factor with levels normal, chemical and overt.

Source

Reaven, G. M. and Miller, R. G. (1979). An attempt to define the nature of chemical diabetes using a multidimensional analysis. Diabetologia 16, 17–24. Andrews, D. F. and Herzberg, A. M. (1985). Data: A Collection of Problems from Many Fields for the Student and Research Worker, Springer-Verlag, Ch. 36.

References

Reaven, G. M. and Miller, R. G. (1979). An attempt to define the nature of chemical diabetes using a multidimensional analysis. Diabetologia 16, 17–24.

Friedman, H. P. and Rubin, J. (1967). On some invariant criteria for grouping data. Journal of the American Statistical Association 62, 1159–1178.

Hawkins, D. M. and McLachlan, G. J., 1997. High-breakdown linear discriminant analysis. Journal of the American Statistical Association 92 (437), 136–143.

Examples

data(diabetes)
(cc <- Linda(group~insulin+glucose+sspg, data=diabetes))
(pr <- predict(cc))

Description

The Fish Catch data set contains measurements on 159 fish caught in the lake Laengelmavesi, Finland.

Usage

data(fish)

Format

A data frame with 159 observations on the following 7 variables.

Weight

Weight of the fish (in grams)

Length1

Length from the nose to the beginning of the tail (in cm)

Length2

Length from the nose to the notch of the tail (in cm)

Length3

Length from the nose to the end of the tail (in cm)

Height

Maximal height as % of Length3

Width

Maximal width as % of Length3

Species

Species

Details

The Fish Catch data set contains measurements on 159 fish caught in the lake Laengelmavesi, Finland. For the 159 fishes of 7 species the weight, length, height, and width were measured. Three different length measurements are recorded: from the nose of the fish to the beginning of its tail, from the nose to the notch of its tail and from the nose to the end of its tail. The height and width are calculated as percentages of the third length variable. This results in 6 observed variables, Weight, Length1, Length2, Length3, Height, Width. Observation 14 has a missing value in variable Weight, therefore this observation is usually excluded from the analysis. The last variable, Species, represents the grouping structure: the 7 species are 1=Bream, 2=Whitewish, 3=Roach, 4=Parkki, 5=Smelt, 6=Pike, 7=Perch. This data set was also analyzed in the context of robust Linear Discriminant Analysis by Todorov (2007), Todorov and Pires (2007).

Source

Journal of Statistical Education, Fish Catch Data Set, [https://jse.amstat.org/datasets/fishcatch.dat.txt] accessed November, 2023.

References

Todorov, V. (2007 Robust selection of variables in linear discriminant analysis, Statistical Methods and Applications, 15, 395–407, doi:10.1007/s10260-006-0032-6.

Todorov, V. and Pires, A.M. (2007) Comparative performance of several robust linear discriminant analysis methods, REVSTAT Statistical Journal, 5, 63–83.

Examples

data(fish)

    # remove observation #14 containing missing value
    fish <- fish[-14,]

    # The height and width are calculated as percentages 
    #   of the third length variable
    fish[,5] <- fish[,5]*fish[,4]/100
    fish[,6] <- fish[,6]*fish[,4]/100
 
    # plot a matrix of scatterplots
    pairs(fish[1:6],
          main="Fish Catch Data",
          pch=21,
          bg=c("red", "green3", "blue", "yellow", "magenta", "violet", 
          "turquoise")[unclass(fish$Species)])

Description

A data set that contains the spectra of six different cultivars of the same fruit (cantaloupe - Cucumis melo L. Cantaloupensis group) obtained from Colin Greensill (Faculty of Engineering and Physical Systems, Central Queensland University, Rockhampton, Australia). The total data set contained 2818 spectra measured in 256 wavelengths. For illustrative purposes are considered only three cultivars out of it, named D, M and HA with sizes 490, 106 and 500, respectively. Thus the data set thus contains 1096 observations. For more details about this data set see the references below.

Usage

data(fruit)

Format

A data frame with 1096 rows and 257 variables (one grouping variable – cultivar – and 256 measurement variables).

Source

Colin Greensill (Faculty of Engineering and Physical Systems, Central Queensland University, Rockhampton, Australia).

References

Hubert, M. and Van Driessen, K., (2004). Fast and robust discriminant analysis. Computational Statistics and Data Analysis, 45(2):301–320. doi:10.1016/S0167-9473(02)00299-2.

Vanden Branden, K and Hubert, M, (2005). Robust classification in high dimensions based on the SIMCA Method. Chemometrics and Intelligent Laboratory Systems, 79(1-2), pp. 10–21. doi:10.1016/j.chemolab.2005.03.002.

Hubert, M, Rousseeuw, PJ and Verdonck, T, (2012). A Deterministic Algorithm for Robust Location and Scatter. Journal of Computational and Graphical Statistics, 21(3), pp 618–637. doi:10.1080/10618600.2012.672100.

Examples

data(fruit)
 table(fruit$cultivar)

Description

Accessor methods to the slots of objects of classCov and its subclasses

Usage

getCenter(obj)
getCov(obj)
getCorr(obj)
getData(obj)
getDistance(obj)
getEvals(obj)
getDet(obj)
getShape(obj)
getFlag(obj, prob=0.975)
getMeth(obj)
isClassic(obj)
getRaw(obj)

Arguments

Methods

obj = "Cov"

generic functions - see getCenter, getCov, getCorr, getData, getDistance, getEvals, getDet, getShape, getFlag, isClassic

obj = "CovRobust"

generic functions - see getCenter, getCov, getCorr, getData, getDistance, getEvals, getDet, getShape, getFlag, getMeth, isClassic

Description

A simple function to calculate the points of a confidence ellipsoid, by default dist=qchisq(0.975, 2)

Usage

getEllipse(loc = c(0, 0), cov = matrix(c(1, 0, 0, 1), ncol = 2), crit = 0.975)

Arguments

Value

A matrix with two columns containing the calculated points.

Author(s)

Valentin Todorov, [email protected]

Examples

data(hbk)
cc <- cov.wt(hbk)
e1 <- getEllipse(loc=cc$center[1:2], cov=cc$cov[1:2,1:2])
e2 <- getEllipse(loc=cc$center[1:2], cov=cc$cov[1:2,1:2], crit=0.99)
plot(X2~X1, data=hbk,
    xlim=c(min(X1, e1[,1], e2[,1]), max(X1,e1[,1], e2[,1])),
    ylim=c(min(X2, e1[,2], e2[,2]), max(X2,e1[,2], e2[,2])))
lines(e1, type="l", lty=1, col="red")
lines(e2, type="l", lty=2, col="blue")
legend("topleft", legend=c(0.975, 0.99), lty=1:2, col=c("red", "blue"))

Description

Accessor methods to the slots of objects of class Pca and its subclasses

Arguments

Methods

obj = "Pca"

Accessors for object of class Pca

obj = "PcaRobust"

Accessors for object of class PcaRobust

obj = "PcaClassic"

Accessors for object of class PcaClassic

Description

The hemophilia data set contains two measured variables on 75 women, belonging to two groups: n1=30 of them are non-carriers (normal group) and n2=45 are known hemophilia A carriers (obligatory carriers).

Usage

data(hemophilia)

Format

A data frame with 75 observations on the following 3 variables.

AHFactivity

AHF activity

AHFantigen

AHF antigen

gr

group - normal or obligatory carrier

Details

Originally analized in the context of discriminant analysis by Habemma and Hermans (1974). The objective is to find a procedure for detecting potential hemophilia A carriers on the basis of two measured variables: X1=log10(AHV activity) and X2=log10(AHV-like antigen). The first group of n1=30 women consists of known non-carriers (normal group) and the second group of n2=45 women is selected from known hemophilia A carriers (obligatory carriers). This data set was also analyzed by Johnson and Wichern (1998) as well as, in the context of robust Linear Discriminant Analysis by Hawkins and McLachlan (1997) and Hubert and Van Driessen (2004).

Source

Habemma, J.D.F, Hermans, J. and van den Broek, K. (1974) Stepwise Discriminant Analysis Program Using Density Estimation in Proceedings in Computational statistics, COMPSTAT'1974 (Physica Verlag, Heidelberg, 1974, pp 101–110).

References

Johnson, R.A. and Wichern, D. W. Applied Multivariate Statistical Analysis (Prentice Hall, International Editions, 2002, fifth edition)

Hawkins, D. M. and McLachlan, G.J. (1997) High-Breakdown Linear Discriminant Analysis J. Amer. Statist. Assoc. 92 136–143.

Hubert, M., Van Driessen, K. (2004) Fast and robust discriminant analysis, Computational Statistics and Data Analysis, 45 301–320.

Examples

data(hemophilia)
plot(AHFantigen~AHFactivity, data=hemophilia, col=as.numeric(as.factor(gr))+1)
##
## Compute robust location and covariance matrix and 
## plot the tolerance ellipses
(mcd <- CovMcd(hemophilia[,1:2]))
col <- ifelse(hemophilia$gr == "carrier", 2, 3) ## define clours for the groups
plot(mcd, which="tolEllipsePlot", class=TRUE, col=col)

Description

”This radar data was collected by a system in Goose Bay, Labrador. This system consists of a phased array of 16 high-frequency antennas with a total transmitted power on the order of 6.4 kilowatts. The targets were free electrons in the ionosphere. "good" radar returns are those showing evidence of some type of structure in the ionosphere. "bad" returns are those that do not; their signals pass through the ionosphere. Received signals were processed using an autocorrelation function whose arguments are the time of a pulse and the pulse number. There were 17 described by 2 attributes per pulse number, corresponding to the complex values returned by the function resulting from the complex electromagnetic signal.” [UCI archive]

Usage

data(ionosphere)

Format

A data frame with 351 rows and 33 variables: 32 measurements and one (the last, Class) grouping variable: 225 'good' and 126 'bad'.

The original dataset at UCI contains 351 rows and 35 columns. The first 34 columns are features, the last column contains the classification label of 'g' and 'b'. The first feature is binary and the second one is only 0s, one grouping variable - factor with labels 'good' and 'bad'.

Source

Source: Space Physics Group; Applied Physics Laboratory; Johns Hopkins University; Johns Hopkins Road; Laurel; MD 20723

Donor: Vince Sigillito ([email protected])

The data have been taken from the UCI Repository Of Machine Learning Databases at https://archive.ics.uci.edu/ml/datasets/ionosphere

This data set, with the original 34 features is available in the package mlbench and a different data set (refering to the same UCI repository) is available in the package dprep (archived on CRAN).

References

Sigillito, V. G., Wing, S. P., Hutton, L. V., and Baker, K. B. (1989). Classification of radar returns from the ionosphere using neural networks. Johns Hopkins APL Technical Digest, 10, 262-266.

Examples

data(ionosphere)
 ionosphere[, 1:6] |> pairs()

Description

Returns TRUE if the covariance matrix contained in a Cov-class object (or derived from) is singular.

Usage

## S4 method for signature 'Cov'
isSingular(obj)

Arguments

Methods

isSingular

signature(x = Cov): Check if a covariance matrix (object of class Cov-class) is singular.

See Also

Cov-class, CovClassic, CovRobust-class.

Examples

data(hbk)
cc <- CovClassic(hbk)
isSingular(cc)

Description

The class Lda serves as a base class for deriving all other classes representing the results of classical and robust Linear Discriminant Analisys methods

Objects from the Class

A virtual Class: No objects may be created from it.

Slots

call:

the (matched) function call.

prior:

prior probabilities used, default to group proportions

counts:

number of observations in each class

center:

the group means

cov:

the common covariance matrix

ldf:

a matrix containing the linear discriminant functions

ldfconst:

a vector containing the constants of each linear discriminant function

method:

a character string giving the estimation method used

X:

the training data set (same as the input parameter x of the constructor function)

grp:

grouping variable: a factor specifying the class for each observation.

covobj:

object of class "Cov" containing the estimate of the common covariance matrix of the centered data. It is not NULL only in case of method "B".

control:

object of class "CovControl" specifying which estimate and with what estimation options to use for the group means and common covariance (or NULL for classical linear discriminant analysis)

Methods

predict

signature(object = "Lda"): calculates prediction using the results in object. An optional data frame or matrix in which to look for variables with which to predict. If omitted, the training data set is used. If the original fit used a formula or a data frame or a matrix with column names, newdata must contain columns with the same names. Otherwise it must contain the same number of columns, to be used in the same order.

show

signature(object = "Lda"): prints the results

summary

signature(object = "Lda"): prints summary information

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

LdaClassic, LdaClassic-class, LdaRobust-class

Examples

showClass("Lda")

Description

Performs a linear discriminant analysis and returns the results as an object of class LdaClassic (aka constructor).

Usage

LdaClassic(x, ...)


## Default S3 method:
LdaClassic(x, grouping, prior = proportions, tol = 1.0e-4, ...)

Arguments

Value

Returns an S4 object of class LdaClassic

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

Lda-class, LdaClassic-class,

Examples

## Example anorexia
library(MASS)
data(anorexia)

## rrcov: LdaClassic()
lda <- LdaClassic(Treat~., data=anorexia)
predict(lda)@classification

## MASS: lda()
lda.MASS <- lda(Treat~., data=anorexia)
predict(lda.MASS)$class

## Compare the prediction results of MASS:::lda() and LdaClassic()
all.equal(predict(lda)@classification, predict(lda.MASS)$class)

Description

Contains the results of a classical Linear Discriminant Analysis

Objects from the Class

Objects can be created by calls of the form new("LdaClassic", ...) but the usual way of creating LdaClassic objects is a call to the function LdaClassic which serves as a constructor.

Slots

call:

The (matched) function call.

prior:

Prior probabilities used, default to group proportions

counts:

number of observations in each class

center:

the group means

cov:

the common covariance matrix

ldf:

a matrix containing the linear discriminant functions

ldfconst:

a vector containing the constants of each linear discriminant function

method:

a character string giving the estimation method used

X:

the training data set (same as the input parameter x of the constructor function)

grp:

grouping variable: a factor specifying the class for each observation.

Extends

Class "Lda", directly.

Methods

No methods defined with class "LdaClassic" in the signature.

Author(s)

Valentin Todorov [email protected]

References

Todorov V & Filzmoser P (2009), An Object Oriented Framework for Robust Multivariate Analysis. Journal of Statistical Software, 32(3), 1–47. doi:10.18637/jss.v032.i03.

See Also

LdaRobust-class, Lda-class, LdaClassic

Examples

showClass("LdaClassic")

Description

Performs robust linear discriminant analysis by the projection-pursuit approach - proposed by Pires and Branco (2010) - and returns the results as an object of class LdaPP (aka constructor).

Usage

LdaPP(x, ...)
## S3 method for class 'formula'
LdaPP(formula, data, subset, na.action, ...)
## Default S3 method:
LdaPP(x, grouping, prior = proportions, tol = 1.0e-4,
                 method = c("huber", "mad", "sest", "class"),
                 optim = FALSE,
                 trace=FALSE, ...)

Arguments

Details

Currently the algorithm is implemented only for binary classification and in the following will be assumed that only two groups are present.

The PP algorithm searches for low-dimensional projections of higher-dimensional data where a projection index is maximized. Similar to the original Fisher's proposal the squared standardized distance between the observations in the two groups is maximized. Instead of the sample univariate mean and standard deviation (T,S) robust alternatives are used. These are selected through the argument method and can be one of

huber

the pair (T,S) are the robust M-estimates of location and scale

mad

(T,S) are the Median and the Median Absolute Deviation

sest

the pair (T,S) are the robust S-estimates of location and scale

class

(T,S) are the mean and the standard deviation.

The first approximation A1 to the solution is obtained by investigating a finite number of candidate directions, the unit vectors defined by all pairs of points such that one belongs to the first group and the other to the second group. The found solution is stored in the slots raw.ldf and raw.ldfconst.

The second approximation A2 (optional) is performed by a numerical optimization algorithm using A1 as initial solution. The Nelder and Mead method implemented in the function optim is applied. Whether this refinement will be used is controlled by the argument optim. If optim=TRUE the result of the optimization is stored into the slots ldf and ldfconst. Otherwise these slots are set equal to raw.ldf and raw.ldfconst.

Value

Returns an S4 object of class LdaPP-class

Warning

Still an experimental version! Only binary classification is supported.

Author(s)

Valentin Todorov [email protected] and Ana Pires [email protected]

References

Pires, A. M. and A. Branco, J. (2010) Projection-pursuit approach to robust linear discriminant analysis Journal Multivariate Analysis, Academic Press, Inc., 101, 2464–2485.

See Also

Linda, LdaClassic

Examples

##
## Function to plot a LDA separation line
##
lda.line <- function(lda, ...)
{
    ab <- lda@ldf[1,] - lda@ldf[2,]
    cc <- lda@ldfconst[1] - lda@ldfconst[2]
    abline(a=-cc/ab[2], b=-ab[1]/ab[2],...)
}

data(pottery)
x <- pottery[,c("MG", "CA")]
grp <- pottery$origin
col <- c(3,4)
gcol <- ifelse(grp == "Attic", col[1], col[2])
gpch <- ifelse(grp == "Attic", 16, 1)

##
## Reproduce Fig. 2. from Pires and branco (2010)
##
plot(CA~MG, data=pottery, col=gcol, pch=gpch)

## Not run: 

ppc <- LdaPP(x, grp, method="class", optim=TRUE)
lda.line(ppc, col=1, lwd=2, lty=1)

pph <- LdaPP(x, grp, method="huber",optim=TRUE)
lda.line(pph, col=3, lty=3)

pps <- LdaPP(x, grp, method="sest", optim=TRUE)
lda.line(pps, col=4, lty=4)

ppm <- LdaPP(x, grp, method="mad", optim=TRUE)
lda.line(ppm, col=5, lty=5)

rlda <- Linda(x, grp, method="mcd")
lda.line(rlda, col=6, lty=1)

fsa <- Linda(x, grp, method="fsa")
lda.line(fsa, col=8, lty=6)

## Use the formula interface:
##
LdaPP(origin~MG+CA, data=pottery)       ## use the same two predictors
LdaPP(origin~., data=pottery)           ## use all predictor variables

##
## Predict method
data(pottery)
fit <- LdaPP(origin~., data = pottery)
predict(fit)

## End(Not run)