## ----echo = FALSE-------------------------------------------------------------
knitr::opts_chunk$set(
 fig.width = 8 ,
 fig.height = 12,
 fig.align ='center'
)

## -----------------------------------------------------------------------------
library("mvdpd")
library("cellWise")
library("robustbase")
library("dplyr")
library("ggplot2")
library("knitr")

## -----------------------------------------------------------------------------
loo.analysis <- function(X, method) {
  p <- ncol(X)
  var_names <- colnames(X)
  eigen_res <- list()
  mu_res <- list()
  md_res <- list()
  ## Leave-One-Out
  for (i in 1:nrow(X)) {
    Xb <- X[-i,]
    res <- method(Xb)
    mu_res[[i]] <- data.frame(
      LeaveOut = i,
      Method = res$method,
      Parameter = var_names,
      Estimate = res$mu
    )
    rownames(mu_res[[i]]) <- var_names
  
    eigen_res[[i]] <- data.frame(
      LeaveOut = i,
      Method = res$method,
      Parameter = paste0("e", 1:p),
      Eigenvalues = eigen(res$Sigma)$values
    )
    rownames(eigen_res[[i]]) <- paste0("e", 1:p)
  
    md_res[[i]] <- data.frame(
      LeaveOut = i,
      Method = res$method,
      Observation = paste0("Obs", seq_len(nrow(X))),
      Distance = sqrt(mahalanobis(X,
        center = res$mu,
        cov = res$Sigma))
    )
  }  
  res <- list(mu=mu_res, eigen=eigen_res, mahalanobis=md_res)
  return(res)
}

## -----------------------------------------------------------------------------
ML <- function(X, ...) {
  list(mu=colMeans(X), Sigma=cov(X), method="ML/MCL")
}

## -----------------------------------------------------------------------------
MCD <- function(X, ...) { 
  mcd <- covMcd(X, ...)
  list(mu=mcd$center, Sigma=mcd$cov, method="MCD")
}

## -----------------------------------------------------------------------------
CELLMCD <- function(X, ...) { 
  cellmcd <- cellMCD(X, checkPars=list(silent=TRUE), ...)
  list(mu=cellmcd$mu, Sigma=cellmcd$S, method="CellMCD")
}

## -----------------------------------------------------------------------------
MDPD <- function(X, beta, ...) { 
  mdpd <- mvnormDPD(X, beta=beta, method="multivariate", ...)
  list(mu=mdpd$mu, Sigma=mdpd$Sigma, method=paste0("MDPD(", beta, ")"))
}

## -----------------------------------------------------------------------------
CDPD <- function(X, beta, ...) { 
  cdpd <- mvnormDPD(X, beta=beta, method="composite", ...)
  list(mu=cdpd$mu, Sigma=cdpd$Sigma, method=paste0("CDPD(", beta, ")"))
}

## -----------------------------------------------------------------------------
perform.analysis <- function(X) {  
  resML <- loo.analysis(X, method=ML)
  resMCD <- loo.analysis(X, method=MCD)
  resCELLMCD <- loo.analysis(X, method=CELLMCD)
  resMDPD1 <- loo.analysis(X, method=function(X) MDPD(X, beta=0.1))
  resMDPD3 <- loo.analysis(X, method=function(X) MDPD(X, beta=0.3))
  resMDPD5 <- loo.analysis(X, method=function(X) MDPD(X, beta=0.5))
  resCDPD1 <- loo.analysis(X, method=function(X) CDPD(X, beta=0.1))
  resCDPD3 <- loo.analysis(X, method=function(X) CDPD(X, beta=0.3))
  resCDPD5 <- loo.analysis(X, method=function(X) CDPD(X, beta=0.5))
  mu_df <- bind_rows(resML$mu, resMCD$mu, resCELLMCD$mu, 
    resMDPD1$mu, resMDPD3$mu, resMDPD5$mu,
    resCDPD1$mu, resCDPD3$mu, resCDPD5$mu)
  eigen_df <- bind_rows(resML$eigen, resMCD$eigen, resCELLMCD$eigen, 
    resMDPD1$eigen, resMDPD3$eigen, resMDPD5$eigen,
    resCDPD1$eigen, resCDPD3$eigen, resCDPD5$eigen)

  md_df <- bind_rows(resML$mahalanobis, 
    resMCD$mahalanobis, resCELLMCD$mahalanobis, 
    resMDPD1$mahalanobis, resMDPD3$mahalanobis, resMDPD5$mahalanobis,
    resCDPD1$mahalanobis, resCDPD3$mahalanobis, resCDPD5$mahalanobis)
  
  mu_var <- mu_df %>%
    group_by(Method, Parameter) %>%
    summarise(
      Variance = var(Estimate),
      .groups = "drop"
    )

  eigen_var <- eigen_df %>%
    group_by(Method, Parameter) %>%
    summarise(
      Variance = var(Eigenvalues),
      .groups = "drop"
    )

  md_var <- md_df %>%
    group_by(Method, Observation) %>%
    summarise(
      Variance = var(Distance),
      .groups = "drop"
    )
  res <- list(mu_df=mu_df, mu_var=mu_var,
    eigen_df=eigen_df, eigen_var=eigen_var,
    md_df=md_df,  md_var=md_var)
  return(res)
}

## -----------------------------------------------------------------------------
plot.results <- function(object) {
  # Boxplot of the location estimates
  mu_df_gg <- ggplot(object$mu_df, aes(x=Method, y=Estimate, fill=Parameter)) +
    geom_boxplot(position = position_dodge(0.8), width = 0.7) +
    scale_x_discrete(name="Method", limits=c("CellMCD", "ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")) +
    scale_y_continuous(name=expression(paste("Leave-one-out ",hat(mu)))) + 
    theme(
      axis.text = element_text(size = 12, face="bold"),
      axis.title.x = element_text(size = 14, face="bold"),
      axis.title.y = element_text(size = 16, face="bold"),
      legend.text = element_text(size = 12, face="bold"),
      legend.title = element_text(size = 14, face="bold")
    )

  # Boxplot of variances of the location estimates
  mu_var_gg <- ggplot(object$mu_var, aes(x = Method, y = Variance)) +
    geom_boxplot(fill = "grey85", width = 0.6) +
    geom_point(aes(color = Parameter), size = 3, 
      position = position_jitter(width = 0.08)) +
    scale_x_discrete(name="Method",
      limits=c("CellMCD", "ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")) +
    scale_y_continuous(name=expression(paste("Variances of leave-one-out ",
      hat(mu)))) + 
    theme(
      axis.text = element_text(size = 12, face="bold"),
      axis.title.x = element_text(size = 14, face="bold"),
      axis.title.y = element_text(size = 16, face="bold"),
      legend.text = element_text(size = 12, face="bold"),
      legend.title = element_text(size = 14, face="bold")
    )

  # Boxplot of the eigenvalues of Scatter estimates
  eigen_df_gg <- ggplot(object$eigen_df, aes(x=Method,
    y=Eigenvalues, fill=Parameter)) +
    geom_boxplot(position = position_dodge(0.8), width = 0.7) +
    scale_x_discrete(name="Method", limits=c("CellMCD", "ML/MCL",
      "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")) +
    scale_y_continuous(name=expression(paste("Leave-one-out eigenvalues of ",
      hat(Sigma)))) + 
    theme(
      axis.text = element_text(size = 12, face="bold"),
      axis.title.x = element_text(size = 14, face="bold"),
      axis.title.y = element_text(size = 16, face="bold"),
      legend.text = element_text(size = 12, face="bold"),
      legend.title = element_text(size = 14, face="bold")
    )

  # Boxplot of variances of the eigenvalues of Scatter estimates
  eigen_var_gg <- ggplot(object$eigen_var, aes(x = Method, y = Variance)) +
    geom_boxplot(fill = "grey85", width = 0.6) +
    geom_point(aes(color = Parameter), size = 3, 
      position = position_jitter(width = 0.08)) +
    scale_x_discrete(name="Method", limits=c("CellMCD",
      "ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")) +
    scale_y_continuous(name=expression(paste("Variances of leave-one-out 
      eigenvalues of ",hat(Sigma)))) + 
    theme(
      axis.text = element_text(size = 12, face="bold"),
      axis.title.x = element_text(size = 14, face="bold"),
      axis.title.y = element_text(size = 16, face="bold"),
      legend.text = element_text(size = 12, face="bold"),
      legend.title = element_text(size = 14, face="bold")
    )
  res <- list(mu_df_gg=mu_df_gg, eigen_df_gg=eigen_df_gg,
    mu_var_gg=mu_var_gg, eigen_var_gg=eigen_var_gg)
  return(res)
}

## -----------------------------------------------------------------------------
data(alcohol)
X <- as.matrix(alcohol)
X <- transfo(X)$Y
resAlcohol <- perform.analysis(X)
plotAlcohol <- plot.results(resAlcohol)  

## ----fig=TRUE-----------------------------------------------------------------
plotAlcohol$mu_df_gg

## ----fig=TRUE-----------------------------------------------------------------
plotAlcohol$mu_var_gg

## ----fig=TRUE-----------------------------------------------------------------
plotAlcohol$eigen_df_gg

## ----fig=TRUE-----------------------------------------------------------------
plotAlcohol$eigen_var_gg

## -----------------------------------------------------------------------------
data(milk)
X <- as.matrix(milk)
X <- transfo(X)$Y
resMilk <- perform.analysis(X)
plotMilk <- plot.results(resMilk)  

## ----fig=TRUE-----------------------------------------------------------------
plotMilk$mu_df_gg

## ----fig=TRUE-----------------------------------------------------------------
plotMilk$mu_var_gg

## ----fig=TRUE-----------------------------------------------------------------
plotMilk$eigen_df_gg

## ----fig=TRUE-----------------------------------------------------------------
plotMilk$eigen_var_gg

## -----------------------------------------------------------------------------
data(bushfire)
X <- as.matrix(bushfire)
X <- transfo(X)$Y
resBushfire <- perform.analysis(X)
plotBushfire <- plot.results(resBushfire)  

## ----fig=TRUE-----------------------------------------------------------------
plotBushfire$mu_df_gg

## ----fig=TRUE-----------------------------------------------------------------
plotBushfire$mu_var_gg

## ----fig=TRUE-----------------------------------------------------------------
plotBushfire$eigen_df_gg

## ----fig=TRUE-----------------------------------------------------------------
plotBushfire$eigen_var_gg

## -----------------------------------------------------------------------------
data(toxicity)
X <- as.matrix(toxicity)
X <- transfo(X)$Y
betas <- c(0, 0.1, 0.3, 0.5)
d <- ncol(X)
pairs <- which(upper.tri(matrix(0, d, d)), arr.ind=TRUE)
result <-matrix(NA, nrow = d*(d+3)/2, ncol=length(betas))
for (i in seq_along(betas)){
  res <- mvnormDPD(X, betas[i], method = "composite")
  cor_matrix <- res$rho
  result[,i] <- c(res$mu, res$sigma,
    cor_matrix[upper.tri(cor_matrix, diag = FALSE)])
}
colnames(result) <- c("ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")
rownames(result) <- c(paste0("mu_", 1:10), paste0("sigma2_", 1:10), paste0("rho_", pairs[,1], pairs[,2]))

## -----------------------------------------------------------------------------
kable(result)

## -----------------------------------------------------------------------------
data(toxicity)
X <- as.matrix(toxicity)
betas <- c(0, 0.1, 0.3, 0.5)
d <- ncol(X)
pairs <- which(upper.tri(matrix(0, d, d)), arr.ind=TRUE)
result <-matrix(NA, nrow = d*(d+3)/2, ncol=length(betas))
for (i in seq_along(betas)){
  res <- mvnormDPD(X, betas[i], method = "composite")
  cor_matrix <- res$rho
  result[,i] <- c(res$mu, res$sigma,
    cor_matrix[upper.tri(cor_matrix, diag = FALSE)])
}
colnames(result) <- c("ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")
rownames(result) <- c(paste0("mu_", 1:10), paste0("sigma2_", 1:10), paste0("rho_", pairs[,1], pairs[,2]))  

## -----------------------------------------------------------------------------
kable(result)