R/plot.MIPCA.R

In missMDA: Handling Missing Values with Multivariate Data Analysis

plot.MIPCA <- function(x,choice="all",axes=c(1,2),new.plot=TRUE,main=NULL,level.conf=0.95, graph.type=c("ggplot","classic"), ...){

####
 procrustes <- function(amat, target, orthogonal = FALSE, translate = FALSE,
        magnify = FALSE) {
        for (i in nrow(amat):1) {
            if (any(is.na(amat)[i, ]) | any(is.na(target)[i,
                ])) {
                amat <- amat[-i, ]
                target <- target[-i, ]
            }
        }
        dA <- dim(amat)
        dX <- dim(target)
        if (length(dA) != 2 || length(dX) != 2)
            stop("arguments amat and target must be matrices")
        if (any(dA != dX))
            stop("dimensions of amat and target must match")
        if (length(attr(amat, "tmat")))
            stop("oblique loadings matrix not allowed for amat")       
if (orthogonal) {
            if (translate) {
                p <- dX[1]
                target.m <- (rep(1/p, p) %*% target)[, ]
                amat.m <- (rep(1/p, p) %*% amat)[, ]
                target.c <- scale(target, center = target.m,
                  scale = FALSE)
                amat.c <- scale(amat, center = amat.m, scale = FALSE)
                j <- svd(crossprod(target.c, amat.c))
            }
            else {
                amat.c <- amat
                j <- svd(crossprod(target, amat))
            }
       

            rot <- j$v %*% t(j$u)
            if (magnify)
                beta <- sum(j$d)/sum(amat.c^2)
            else beta <- 1

            B <- beta * amat.c %*% rot
            if (translate)
                B <- B + rep(as.vector(target.m), rep.int(p,
                  dX[2]))
   
       value <- list(rmat = B, tmat = rot, magnify = beta)
            if (translate)
                value$translate <- target.m - (rot %*% amat.m)[,
                  ]
    
  }

        else {
            b <- solve(amat, target)
            gamma <- sqrt(diag(solve(crossprod(b))))
            rot <- b * rep(gamma, rep.int(dim(b)[1], length(gamma)))
            B <- amat %*% rot
            fcor <- solve(crossprod(rot))
            value <- list(rmat = B, tmat = rot, correlation = fcor)
        }

        return(value)
    }
####
  res <- x
  if (!inherits(res, "MIPCA")) stop("non convenient data")
  ncp <- max(axes)
  reference <- FactoMineR::PCA(res$res.imputePCA,scale.unit=res$call$scale,graph=FALSE,ncp=ncp)
  rec.pca <- res$res.imputePCA
#  rec <- reconst(reference,ncp)
#  rec.pca <- as.matrix(res$call$X)
#  rec.pca[res$call$missing] <- rec[res$call$missing]
  
  res.var <- res.supp <- rec.pca
  res.procrustes <- reference$ind$coord[,1:ncp]
  res.dim <- as.matrix(res$res.imputePCA)

##for (i in 1:dim(res$res.MI)[3]){
## rec.pca <- res$res.MI[,,i]
for (i in 1:length(res$res.MI)){
 rec.pca <- res$res.MI[[i]]
 acpfin <- FactoMineR::PCA(rec.pca, scale.unit=res$call$scale,graph=FALSE,ncp=ncp)

 tourne <- procrustes(acpfin$ind$coord[,1:ncp], reference$ind$coord[,1:ncp],orthogonal = TRUE, translate = TRUE, magnify = TRUE)$rmat

 colnames(tourne) <- colnames(res.procrustes)
 colnames(rec.pca) <- colnames(res.supp)
 res.procrustes <- rbind.data.frame(res.procrustes,tourne)
 res.supp <- rbind.data.frame(res.supp,rec.pca)
 res.var <- cbind.data.frame(res.var,rec.pca)
 res.dim <- cbind.data.frame(res.dim,acpfin$ind$coord[,1:ncp])
}

####
if (!is.null(main)) title <- main
graph.type <- match.arg(graph.type[1],c("ggplot","classic"))
if (graph.type=="ggplot") graph <- list()
if ((choice=="all")|(choice=="ind.proc")){
  if ((new.plot)&!nzchar(Sys.getenv("RSTUDIO_USER_IDENTITY"))) dev.new()
  oo=FactoMineR::PCA(res.procrustes,ind.sup=c((nrow(res$call$X)+1):nrow(res.procrustes)),scale.unit=FALSE,graph=FALSE)
  oo$eig=reference$eig
  if (is.null(rownames(res$call$X))) rownames(res$call$X) <- 1:nrow(res$call$X)
  el=coord.ellipse(cbind.data.frame(as.factor(rep(rownames(res$call$X),res$call$nboot)),oo$ind.sup$coord[,axes]),level.conf=level.conf) 
  if (is.null(main)) title="Multiple imputation using Procrustes" 
  PlotIndProc <- plot(oo,axes=axes,col.ind.sup=rep(1:nrow(res$call$X),res$call$nboot),label="ind",ellipse=el,col.quali=1, title=title,invisible="ind.sup",new.plot=FALSE,graph.type=graph.type)
  if (graph.type=="ggplot"){
    print(PlotIndProc)
    graph$PlotIndProc <- PlotIndProc
  }
#  if (!is.null(add.tab)){
#    vrai = PCA(add.tab,graph=FALSE,scale=res$call$scale)
#    tourne <- procrustes(vrai$ind$coord[,axes], reference$ind$coord[,axes],orthogonal = TRUE, translate = TRUE, magnify = TRUE)$rmat
#    points(tourne[,axes],cex=0.9,col=2)
#  }
}

if ((choice=="all")|(choice=="dim")){
  if ((new.plot)&!nzchar(Sys.getenv("RSTUDIO_USER_IDENTITY"))) dev.new()
  colnames(res.dim)=paste("V",1:ncol(res.dim))
  ooo=FactoMineR::PCA(res.dim,quanti.sup=(ncol(res$call$X)+1):ncol(res.dim),scale.unit=res$call$scale,graph=FALSE)
  ooo$eig=reference$eig
  if (is.null(main)) title="Projection of the Principal Components"  
  PlotDim <- plot(ooo,choi="var",axes=axes,title=title,label="none",new.plot=FALSE,invisible="var",graph.type=graph.type)
  if (graph.type=="ggplot"){
    print(PlotDim)
    graph$PlotDim <- PlotDim
  }
}

if ((choice=="all")|(choice=="ind.supp")){
  if ((new.plot)&!nzchar(Sys.getenv("RSTUDIO_USER_IDENTITY"))) dev.new()
  oo=FactoMineR::PCA(res.supp,ind.sup=c((nrow(res$call$X)+1):nrow(res.supp)),scale.unit=res$call$scale,graph=FALSE,ncp=ncp)
  el=coord.ellipse(cbind.data.frame(as.factor(rep(rownames(res$call$X),res$call$nboot)),oo$ind.sup$coord),level.conf = level.conf,axes = axes)
  if (is.null(main)) title="Supplementary projection"    
  PlotIndSupp <- plot(oo,axes=axes,col.ind.sup=rep(1:nrow(res$call$X),res$call$nboot),label="ind",ellipse=el,col.quali=1,
    title=title,invisible="ind.sup",new.plot=FALSE, graph.type=graph.type)
  if (graph.type=="ggplot"){
    print(PlotIndSupp)
    graph$PlotIndSupp <- PlotIndSupp
  }
#  if (!is.null(add.tab)){
#    dele = PCA(rbind.data.frame(rec.pca,add.tab),ind.sup=c((nrow(res$call$X)+1):(2*nrow(res$call$X))),scale=scale,graph=FALSE)
#    points(dele$ind.sup$coord[,axes],col=2)
#  }
}

if ((choice=="all")|(choice=="var")){
  if ((new.plot)&!nzchar(Sys.getenv("RSTUDIO_USER_IDENTITY"))) dev.new()
  color = c("black", "red", "green3", "blue", "cyan", "magenta", 
            "darkgray", "darkgoldenrod", "darkgreen", "violet", 
            "turquoise", "orange", "lightpink", "lavender", "yellow", 
            "lightgreen", "lightgrey", "lightblue", "darkkhaki", 
            "darkmagenta", "darkolivegreen", "lightcyan", "darkorange", 
            "darkorchid", "darkred", "darksalmon", "darkseagreen", 
            "darkslateblue", "darkslategray", "darkslategrey", 
            "darkturquoise", "darkviolet", "lightgray", "lightsalmon", 
            "lightyellow", "maroon")
  colnames(res.var)=paste("V",1:ncol(res.var))
  colnames(res.var)[1:ncol(res$call$X)]=colnames(res$call$X)
  oo=FactoMineR::PCA(res.var,quanti.sup=c((ncol(res$call$X)+1):ncol(res.var)),scale.unit=res$call$scale,graph=FALSE)
  if (is.null(main)) title="Variable representation"    
  PlotVar <- plot(oo, axes=axes, choix = "var", title=title,invisible = "quanti.sup", col.hab = color[1:ncol(res$call$X)],new.plot=FALSE,graph.type=graph.type)
  if (graph.type=="classic") {
    for (k in 1:res$call$nboot) points(oo$quanti.sup$coord[((k-1)*ncol(res$call$X)+1):(k*ncol(res$call$X)),axes[1]], oo$quanti.sup$coord[((k-1)*ncol(res$call$X)+1):(k*ncol(res$call$X)),axes[2]], col = color[1:ncol(res$call$X)], pch = 15, cex = 0.3)
  } else {
    y <- variables <- NULL ## to avoid no visible binding for global variable
    dta <- cbind.data.frame(x=oo$quanti.sup$coord[,axes[1]],y=oo$quanti.sup$coord[,axes[2]],variables=rep(colnames(res$call$X),res$call$nboot))
    PlotVar <- PlotVar + geom_point(data=dta,aes(x=x,y=y,color=variables),alpha=0.5,shape=20,size=2)+guides(color = guide_legend(title=NULL,override.aes = list(alpha=1, size=4)))
    print(PlotVar)
    graph$PlotVar <- PlotVar
  }
}
  if (graph.type=="ggplot") return(graph)
}