R/predict.mvout.R

In mvout: Robust Multivariate Outlier Detection

predict.mvout <-
  function(object, 
           x, 
           type = c("distance", "outlier", "scores"), 
           thresh = 0.01, 
           ...){
    # predictions from 'mvout' objects
    # Nathaniel E. Helwig
    # updated: 2025-05-21
    
    
    #########***#########   INITIAL CHECKS   #########***#########
    
    # check object
    if(!inherits(object, "mvout")){
      stop("Input 'object' should be of class 'mvout'")
    }
    
    # check n and p
    n <- length(object$mcd$mcd.wt)
    p <- nrow(object$mcd$cov)
    
    # check x
    nox <- FALSE
    if(missing(x)){
      nox <- TRUE
      x <- object$args$x
      if(is.null(x)) stop("Input 'x' must be provided when object$args$x is NULL")
    }
    x <- as.matrix(x)
    newn <- nrow(x)
    if(ncol(x) != p) stop("Input 'x' has the wrong number of columns")
    
    # check 'type'
    type <- as.character(type[1])
    type <- pmatch(type, c("distance", "outlier", "scores"))
    if(is.na(type)) stop("Invalid 'type' input. See help files.")
    type <- c("distance", "outlier", "scores")[type]
    
    # check 'type' w.r.t object$args$method
    if(type == "scores" && object$args$method == "none"){
      stop("Input 'type' cannot be 'scores' if object$args$method is 'none'")
    }
    
    # check 'thresh'
    thresh <- as.numeric(thresh[1])
    if(thresh <= 0 | thresh >= 1) stop("Input 'thresh' must satisfy: 0 < thresh < 1")
    
    
    #########***#########   ROBUST Z-SCORES   #########***#########
    
    # center (and scale?)
    if(object$args$standardize){
      x <- scale(x, center = object$mcd$center, scale = sqrt(diag(object$mcd$cov)))
      orig.cov <- object$mcd$cov
      object$mcd$cov <- object$mcd$cor
    } else {
      x <- scale(x, center = object$mcd$center, scale = FALSE)
    }
    
    
    #########***#########   MAHALANOBIS DISTANCE   #########***#########
    
    # directional distances?
    directional <- FALSE
    if(any(object$args$direction != "two.sided")){
      directional <- TRUE
      xorig <- x
      indx <- which(object$args$direction != "two.sided")
      for(j in indx){
        sgn <- ifelse(object$args$direction[j] == "less", -1, 1)
        x[,j] <- sgn * pmax(sgn * x[,j], 0)
      }
    }
    
    # robust distance calculation
    if(object$args$method == "none"){
      
      rdis <- mahalanobis(x, center = FALSE, cov = object$mcd$cov)
      
    } else if(object$args$method[1] == "princomp"){
      
      object$loadings <- tcrossprod(object$loadings, object$invrot)
      eigvals <- colSums(object$loadings^2)
      object$loadings <- object$loadings %*% diag(1/sqrt(eigvals), nrow = object$args$factors)
      scores <- x %*% object$loadings
      rdis <- mahalanobis(scores, center = FALSE, cov = diag(eigvals, nrow = object$args$factors))
      if(type == "scores"){
        if(directional) x <- xorig
        scores <- x %*% object$loadings %*% diag(1/sqrt(eigvals), nrow = object$args$factors)
      }
      p <- object$args$factors
      
    } else {
      
      object$loadings <- tcrossprod(object$loadings, object$invrot)
      covmat <- tcrossprod(object$loadings) + diag(object$uniquenesses)
      rdis <- mahalanobis(x, center = FALSE, cov = covmat)
      if(type == "scores"){
        if(directional) x <- xorig
        PsiInvLo <- diag(1 / object$uniquenesses, nrow = p, ncol = p) %*% object$loadings
        if(object$args$scores == "regression"){
          scores <- x %*% PsiInvLo %*% solve(diag(object$args$factors) + crossprod(object$loadings, PsiInvLo))
        } else {
          scores <- x %*% PsiInvLo %*% solve(crossprod(object$loadings, PsiInvLo))
        }
      }
      
    } # end if(method == "none")
    
    
    #########***#########   RETURN RESULT   #########***#########
    
    if(type == "distance"){
      return(rdis)
    } else if(type == "outlier"){
      cval.in <- qchisq(1 - thresh, df = p)
      cval.out <- qf(1 - thresh, df1 = p, df2 = n - p) * ((n - 1) * p / (n - p))
      if(nox){
        cval <- ifelse(object$mcd$mcd.wt, cval.in, cval.out)
      } else {
        cval <- cval.out
      }
      return(ifelse(rdis > cval, TRUE, FALSE))
    } else {
      scores <- scores %*% object$invrot
      rownames(scores) <- rownames(x)
      colnames(scores) <- paste0("Factor", 1:ncol(scores))
      return(scores)
    }
    
  } # predict.mvout