R/T2_qualitative.R

In T2Qv: Control Qualitative Variables

#' @import ca
#' @import highcharter
#' @import dplyr
#'
globalVariables(c("var","qchisq","sd","abline"))

#' @title Multivariate control chart for qualitative variables
#'
#' @description Multivariate control chart T2 Hotelling applicable for qualitative variables.
#' @param base Data set
#' @param IndK Character with the name of the column that specifies the partition of the data set in k tables.
#' @param dim Dimension taken for reduction. Initial dimension - 1 is recommended.
#' @param interactive If it is TRUE, the graph will be shown interactively. If FALSE, the graph is displayed flat. FALSE is the default.
#' @param alpha Type I error, it is recommended to reach this value by using the ARL.
#' @return A control chart made with the T2 hotelling statistic, applied to detect anomalies in any of the K tables obtained with the specification of \code{IndK}. The control limit of the graph is obtained from the number of dimensions \code{dim} and the type I error \code{alpha}.
#' @examples
#' data(Datak10Contaminated)
#' T2_qualitative(Datak10Contaminated,"GroupLetter",9, FALSE,0.0027)
#' @export
T2_qualitative <- function(base, IndK, dim, interactive=FALSE, alpha=0.0027){

  Table <- list()
  Ind <- base%>% pull(IndK)
  groupFactor=as.factor(Ind)

  for (i in 1:length(levels(groupFactor))){
    Table[[i]] <- subset(base, Ind==as.character(levels(groupFactor)[i]))
  }



  for (i in 1:length(levels(groupFactor))){
    Table[[i]] <- Table[[i]][,!names(Table[[i]]) %in% IndK, drop = F]}

  mjcatable <- function(base){
    MJCA <- mjca(base,nd = 2)
    BURT <- MJCA$Burt
    P <- BURT/sum(as.matrix(BURT))
    r <- apply(P, 1, sum)
    c <- apply(P, 2, sum)
    D_r_inv_squ <- diag(1/sqrt(r))
    D_c_inv_squ <- diag(1/sqrt(c))
    S <- D_r_inv_squ%*%(P-r%*%t(c))%*%D_c_inv_squ
    SVD <- svd(S)
    U <- SVD$u
    V <- SVD$v
    F <- D_r_inv_squ%*%U
    C <- D_c_inv_squ%*%V
    C <- C[,1:dim]
    return(C)
  }
  # mjcatable <- function(base, dime){

  #    MJCA <- mjca(base,nd = dime)
  #    Col <- data.frame(MJCA$colcoord[,1:dime])
  #    rownames(Col) <- MJCA$levelnames
  #    return(Col)
  # }
  colcoor <- list()



  for (i in 1:length(levels(groupFactor))){
    colcoor[[i]] <- mjcatable(Table[[i]])
  }

  NormalizacionAFM <- function(mjca){
    SVD <- svd(mjca)
    ACP1_VALUE <- (SVD$d[1])
    return(mjca/ACP1_VALUE)
  }
  coornorm <- list()
  for (i in 1:length(levels(groupFactor))){
    coornorm[[i]] <- NormalizacionAFM(abs(colcoor[[i]]))
    colnames(coornorm[[i]]) <- paste0("V",1:ncol(coornorm[[i]]))
  }


  colnum <- c()
  for (i in 1:length(coornorm)){
    colnum[i] <- ncol(coornorm[[i]])
  }

  General <- data.frame(coornorm[[1]][,1:min(colnum)])

  for (i in 2:length(levels(groupFactor))){
    General <- rbind(General,coornorm[[i]][,1:min(colnum)])
  }

  mu00 <- apply(General,2, 'median')
  muii <- list()
  n <- list()
  for (i in 1:length(levels(groupFactor))){
    muii[[i]] <- apply(coornorm[[i]][,1:min(colnum)],2,'median')
    n[[i]] <- nrow(coornorm[[i]])
  }

  t2 <- list()
  sigma <- var(General)

  for (i in 1:length(levels(groupFactor))){
    t2[[i]]=n[[i]]*(t(muii[[i]]-mu00)%*%solve(sigma)%*%(muii[[i]]-mu00))
  }

  T2 <- as.data.frame(as.matrix(t2))
  DtGraph <- data.frame(table=seq(1:nrow(T2)),hote=as.numeric(as.matrix(T2$V1)))
  p=ncol(base)
  m=length(colcoor)

  #alpha2 <- 1-(1-alpha)^dim
  alpha2 <- alpha

  LC <- qchisq(p=alpha2,df=dim, lower.tail = FALSE)

  # LC <-( (p*(m+1)*(m-1))/(m*(m-p)))*qf(alpha,p,(m-p))
  if (max(DtGraph$hote)>LC){
    YLIM=max(DtGraph$hote)+sd(DtGraph$hote)
  } else {
    YLIM=LC+sd(DtGraph$hote)
  }

  Categories <- levels(groupFactor)
  if (interactive==FALSE){
    plot(DtGraph, type='l', main="Multivariate Control Chart", sub=paste0("UCL = ",round(LC,2), ", alpha = ",alpha,", ARL = ",round(1/alpha)),xlab="k Table", ylab="T2 Hotelling",ylim=c(0,YLIM))+
      abline(h=LC, col="blue")
    message('If you want the interactive chart specify interactive = TRUE')

  }else{
    DtGraph$hote=round(DtGraph$hote,3)
    highchart()%>%
      hc_add_series(DtGraph, type='line',hcaes( y='hote'), name="T2 Hotelling", color="#1C3F63")%>%
      hc_title(text="Multivariate Control Chart")%>%
      hc_subtitle(text=paste0("UCL = ",round(LC,2), ", alpha = ",alpha,", ARL = ",round(1/alpha)))%>%
      hc_xAxis(categories=Categories)%>%
      hc_yAxis(max=YLIM,
               plotLines = list(list(
                 value = LC,
                 color = '#821D1D',
                 width = 3,
                 zIndex = 4,
                 label = list(text = "",
                              style = list( color = '#821D1D', fontWeight = 'bold' )))))
    # %>%
    #   hc_annotations(
    #     list(labelOptions = list(y = 35, x = 0, backgroundColor = '#E6EEFF', borderColor = "#1D4B5E"),
    #          labels = list(
    #            list(style = list(color = '#1D4B5E', fontSize = 8),
    #                 useHTML = TRUE,
    #                 point = list(x = 1, y = LC+2*sd(DtGraph$hote), xAxis = 0, yAxis = 0),text = paste0("UL = ",round(LC,2), "<br/> alpha = ",alpha,"<br/> ARL = ",round(1/alpha)))
    #          )
    #     )
    #   )
  }
}