R/MVRr.r

In jedazard/MVR: Mean-Variance Regularization

##########################################################################################################################################
# MVR
##########################################################################################################################################

##########################################################################################################################################
# 1. END-USER REGULARIZATION & VARIANCE STABILIZATION FUNCTIONS
##########################################################################################################################################

##########################################################################################################################################
################
# Usage         :
################
#                   mvr(data,
#                       block=rep(1,nrow(data)),
#                       tolog=FALSE,
#                       nc.min=1,
#                       nc.max=30,
#                       probs=seq(0, 1, 0.01),
#                       B=100,
#                       parallel=FALSE,
#                       conf=NULL,
#                       verbose=TRUE,
#                       seed=NULL)
#
################
# Description   :
################
#
################
# Arguments     :
################
#
################
# Values        :
################
#
################
#
##########################################################################################################################################

mvr <- function(data,
                block=rep(1,nrow(data)),
                tolog=FALSE,
                nc.min=1,
                nc.max=30,
                probs=seq(0, 1, 0.01),
                B=100,
                parallel=FALSE,
                conf=NULL,
                verbose=TRUE,
                seed=NULL) {

   if (any(table(block) == 1))
      stop("All group sample sizes must be greater than 1!")

   p <- ncol(data)
   if (is.null(colnames(data))) colnames(data) <- paste("v", 1:p, sep="")

   block <- factor(x=block, levels=unique(as.character(block)))
   lev <- levels(block)
   ng <- nlevels(block)
   def <- vector(mode="list", length=ng)
   tab <- numeric(ng)
   for (g in 1:ng) {
      def[[g]] <- which(block == lev[g])
      tab[g] <- length(def[[g]])
   }

   X.raw <- as.matrix(data)
   if (tolog) {
      if (!is.empty(which(data <= 0))) {
         data[data <= 0] <- 1
         warning("Negative or null values have been changed to 1 before taking log2-transformation \n")
      }
      data <- log2(data)
   }

   if (parallel) {     
      if (conf$type == "SOCKET") {
         cl <- parallel::makeCluster(spec=conf$spec,
                                     type="PSOCK",
                                     homogeneous=conf$homo,
                                     outfile=conf$outfile,
                                     verbose=conf$verbose)
      } else if (conf$type == "MPI") {
        cl <- parallel::makeCluster(spec=conf$spec,
                                     type="MPI",
                                     homogeneous=conf$homo,
                                     outfile=conf$outfile,
                                     verbose=conf$verbose)
      } else {
        stop("Unrecognized cluster type\n")
      }
   } else {
        cl <- NULL
   }

   MVR <- vector(mode="list", length=ng)
   M <- matrix(data=NA, nrow=p, ncol=ng, dimnames=list(colnames(data), lev))
   for (g in 1:ng) {
      if (verbose) cat("Computing optimal cluster configuration of group ", g, " ... \n", sep="")
      MVR[[g]] <- MeanVarReg(data=data[def[[g]], , drop=FALSE],
                             nc.min=nc.min,
                             nc.max=nc.max,
                             probs=probs,
                             B=B,
                             parallel=parallel,
                             cl=cl,
                             verbose=verbose,
                             seed=seed)
      M[,g] <- MVR[[g]]$membership
      if (verbose) cat("group : ", g, " done \n")
   }

   if (parallel) {
      parallel::stopCluster(cl)  
   }

   clus <- merging.cluster(M)
   vec.av <- numeric(p)
   vec.sd <- numeric(p)
   for (k in unique(clus)) {
      pma <- pmatch(x=colnames(data), table=names(clus))
      ind <- (clus == k)[pma]
      vec.av[ind] <- mean(pooled.mean(x=data[, ind, drop=FALSE], block=block))
      vec.sd[ind] <- mean(pooled.sd(x=data[, ind, drop=FALSE], block=block))
   }
   X.mvr <- scale(x=data, center=vec.av, scale=vec.sd)
   X.mvr <- as.matrix(X.mvr)
   
   cat("Finished!\n")

   return(structure(list(Xraw = X.raw,
                         Xmvr = X.mvr,
                         centering = vec.av,
                         scaling = vec.sd,
                         MVR = MVR,
                         block = block,
                         tolog = tolog,
                         nc.min = nc.min,
                         nc.max = nc.max,
                         probs = probs,
                         seed = seed),
                    class = "MVR"))

}
##########################################################################################################################################





##########################################################################################################################################
# 2. END-USER REGULARIZED TESTS-STATISTICS FUNCTIONS
##########################################################################################################################################

##########################################################################################################################################
################
# Usage         :
################
#                   mvrt.test(data,
#                             obj=NULL,
#                             block,
#                             tolog=FALSE,
#                             nc.min=1,
#                             nc.max=30,
#                             pval=FALSE,
#                             replace=FALSE,
#                             n.resamp=100,
#                             parallel=FALSE,
#                             conf=NULL,
#                             verbose=TRUE,
#                             seed=NULL)
#
################
# Description   :
################
#
################
# Arguments     :
################
#
################
# Values        :
################
#
################
#
##########################################################################################################################################

mvrt.test <- function(data,
                      obj=NULL,
                      block,
                      tolog=FALSE,
                      nc.min=1,
                      nc.max=30,
                      pval=FALSE,
                      replace=FALSE,
                      n.resamp=100,
                      parallel=FALSE,
                      conf=NULL,
                      verbose=TRUE,
                      seed=NULL) {

   if (is.null(obj)) {
      block <- factor(x=block, levels=unique(as.character(block)))
   } else {
      block <- factor(x=obj$block, levels=unique(as.character(obj$block)))
      data <- obj$Xraw
      nc.min <- obj$nc.min
      nc.max <- obj$nc.max
      tolog <- obj$tolog
   }

   if (any(table(block) == 1))
      stop("All group sample sizes must be greater than 1!")

   if (nlevels(block) <= 1)
      stop("Number of sample groups less than 2. You must have at least two sample groups!")

   p <- ncol(data)

   if (is.null(colnames(data))) colnames(data) <- paste("v", 1:p, sep="")

   if (tolog) {
      if (!is.empty(which(data <= 0))) {
         data[data <= 0] <- 1
         warning("Negative or null values have been changed to 1 before taking log2-transformation \n")
      }
      data <- log2(data)
   }

   lev <- levels(block)
   ng <- nlevels(block)
   def <- vector(mode="list", length=ng)
   tab <- numeric(ng)
   for (g in 1:ng) {
      def[[g]] <- which(block == lev[g])
      tab[g] <- length(def[[g]])
   }

   if (parallel) {
      if (conf$type == "SOCKET") {
         cl <- parallel::makeCluster(spec=conf$spec,
                                     type="PSOCK",
                                     homogeneous=conf$homo,
                                     outfile=conf$outfile,
                                     verbose=conf$verbose)
      } else if (conf$type == "MPI") {
        cl <- parallel::makeCluster(spec=conf$spec,
                                     type="MPI",
                                     homogeneous=conf$homo,
                                     outfile=conf$outfile,
                                     verbose=conf$verbose)
      } else {
        stop("Unrecognized cluster type\n")
      }
   } else {
        cl <- NULL
   }

   t.reg <- mvrt(x=data,
                 obj=obj,
                 lev=lev,
                 tab=tab,
                 ng=ng,
                 def=def,
                 nc.min=nc.min,
                 nc.max=nc.max,
                 parallel=parallel,
                 cl=cl,
                 verbose=verbose,
                 seed=seed)

   if (pval == FALSE) {
      p.value <- NULL
   } else {
      if (verbose == TRUE) {
         cat("Computing t-statistics null distributions ...\n")
      }
      resamp <- 1:n.resamp
      if (!parallel) {
         if (is.null(seed)) {
            digits <- getOption("digits")
            seed <- runif(n=n.resamp, min=1, max=2) * 10^(digits-2)
         } else {
            seed <- (0:(n.resamp-1)) + seed
         }
         stat.bo <- statresamp(X=data,
                               block=block,
                               nc.min=nc.min,
                               nc.max=nc.max,
                               resamp=resamp,
                               replace=replace,
                               parallel=parallel,
                               verbose=verbose,
                               seed=seed)
      } else {
         clusterSetRNGStream(cl=cl, iseed=seed)
         stat.cl <- clusterApply(cl=cl, x=resamp, fun=statresamp,
                                 X=data,
                                 block=block,
                                 nc.min=nc.min,
                                 nc.max=nc.max,
                                 replace=replace,
                                 parallel=parallel,
                                 verbose=verbose,
                                 seed=NULL)
         stat.bo <- matrix(data=NA, nrow=0, ncol=p)
         for (b in 1:n.resamp) {
            stat.bo <- rbind(stat.bo, stat.cl[[b]])
         }
      }
      if (verbose == TRUE) {
         cat("Computation of p-values ...\n")
      }
      if (!is.null(seed)) {
         set.seed(seed)
      }
      p.value <- numeric(p)
      for (j in 1:p) {
         x <- sum(abs(stat.bo[,j]) > abs(t.reg[j]), na.rm=TRUE)
         if (replace == FALSE) {
            p.value[j] <- statmod::permp(x=x, n1=tab[1], n2=tab[2], nperm=n.resamp)
         } else {
            p.value[j] <- x/n.resamp
         }
      }
   }

   if (parallel) {
      parallel::stopCluster(cl)
   }

   cat("Finished!\n")

   return(structure(list(statistic = t.reg,
                         p.value = p.value,
                         seed = seed),
                    class = "MVRT"))

}
##########################################################################################################################################





##########################################################################################################################################
# 3. END-USER DIAGNOSTIC PLOTS FOR QUALITY CONTROL
##########################################################################################################################################

##########################################################################################################################################
################
# Usage         :
################
#                   cluster.diagnostic(obj,
#                                      span=0.75,
#                                      degree=2,
#                                      family="gaussian",
#                                      title="Cluster Diagnostic Plots",
#                                      device=NULL,
#                                      file="Cluster Diagnostic Plots",
#                                      path=getwd()
#                                      horizontal=FALSE,
#                                      width=8.5,
#                                      height=11, ...)
#
################
# Description   :
################
#
################
# Arguments     :
################
#
################
# Values        :
################
#
################
#
##########################################################################################################################################

cluster.diagnostic <- function(obj,
                               span=0.75,
                               degree=2,
                               family="gaussian",
                               title="Cluster Diagnostic Plots",
                               device=NULL,
                               file="Cluster Diagnostic Plots",
                               path=getwd(),
                               horizontal=FALSE,
                               width=8.5,
                               height=11, ...) {

   block <- factor(x=obj$block, levels=unique(as.character(obj$block)))
   lev <- levels(block)
   ng <- nlevels(block)

   clusterplot <- function(obj, ng, title, span, degree, family, ...) {
      par(mfcol=c(3, ng), oma=c(0, 0, 3, 0), mar=c(4, 3, 3, 1), mgp=c(2, 0.5, 0), xpd=FALSE)

      nc.min <- obj$nc.min
      nc.max <- obj$nc.max
      probs <- obj$probs
      n <- nrow(obj$Xmvr)
      p <- ncol(obj$Xmvr)

      X <- matrix(data=rnorm(n=n*p, mean=0, sd=1), ncol=p)
      mu.null       <- pooled.mean(x=X, block=block)
      mu.null.quant <- quantile(x=mu.null, probs=probs)
      sd.null       <- (pooled.sd(x=X, block=block) /sqrt(n-ng)) * sqrt(rchisq(n=p, df=n-ng))
      sd.null.quant <- quantile(x=sd.null, probs=probs)

      for (g in 1:ng) {
         gap <- obj$MVR[[g]]$gap
         sde <- obj$MVR[[g]]$sde
         nc <- obj$MVR[[g]]$nc
         low <- loess(gap ~ as.numeric(nc.min:nc.max), span=span, degree=degree, family=family)$fitted
         plot(nc.min:nc.max, gap, type="n", axes=FALSE, main=paste("Group ", g, sep=""), cex.main=1,
              xlim=range(nc.min, nc.max), ylim=range(0, max(gap + sde), min(gap - sde)),
              xlab="Number of clusters", ylab="Sim Statistic")
         lines(nc.min:nc.max, gap, type="b", pch=3, col=1, lty=1)
         lines(low, col=2, lty=5)
         arrows(nc.min:nc.max, gap + sde, nc.min:nc.max, gap - sde, angle=90, code=3, length=0.1, col=1)
         arrows(nc, min(gap) + max(gap)/5, nc, min(gap) + max(gap)/10, length=0.05, angle=20, code=2, col=2, lwd=2)
         text(x=nc, y=min(gap) + max(gap)/5, labels=nc, pos=3, col=2)
         axis(side=1, at=seq(from=nc.min, to=nc.max, by=1), pos=0)
         axis(side=2, at=NULL, pos=nc.min)
         legend("topright", inset=0.05, legend=c("Sim Stat", "Lowess"), col=c(1,2), lty=c(1,5), cex=0.7)

         mu.quant <- obj$MVR[[g]]$mu.quant
         sd.quant <- obj$MVR[[g]]$sd.quant
         w <- which((nc.min:nc.max) <= obj$MVR[[g]]$nc)

         plot(probs, mu.null.quant, type="l", lty=1, col="green", lwd=4, cex.main=1,
              main=paste("Group ", g, sep=""), xlab="Percentiles", ylab="Transformed Pooled Means")
         matplot(probs, t(mu.quant), type="l", lty=4, col="black", lwd=1, add=TRUE)
         matplot(probs, t(mu.quant)[,w], type="l", lty=4, col="red", lwd=1, add=TRUE)
         segments(x0=0.5, x1=0.5, y0=min(mu.null.quant, mu.quant), y1=0, col="black", lty=2, lwd=0.5)
         segments(x0=0.5, x1=0, y0=0, y1=0, col="black", lty=2, lwd=0.5)

         plot(probs, sd.null.quant, type="l", lty=1, col="green", lwd=4, cex.main=1,
              main=paste("Group ", g, sep=""), xlab="Percentiles", ylab="Transformed Pooled Standard Deviations")
         matplot(probs, t(sd.quant), type="l", lty=4, col="black", lwd=1, add=TRUE)
         matplot(probs, t(sd.quant)[,w], type="l", lty=4, col="red", lwd=1, add=TRUE)
         segments(x0=0.5, x1=0.5, y0=min(sd.null.quant, sd.quant), y1=1, col="black", lty=2, lwd=0.5)
         segments(x0=0.5, x1=0, y0=1, y1=1, col="black", lty=2, lwd=0.5)
      }
      mtext(text=title, cex=1, side=3, outer=TRUE)
   }

   if (is.null(device)) {
      dev.new(width=width, height=height, title="Cluster Diagnostic Plots", noRStudioGD = TRUE)
      clusterplot(obj=obj, ng=ng, title=title, span=span, degree=degree, family=family)
   } else if (device == "PS") {
      path <- normalizePath(path=paste(path, "/", sep=""), winslash="\\", mustWork=FALSE)
      file <- paste(file, ".ps", sep="")
      cat("\nOUTPUT: \n")
      cat("Filename : ", file, "\n")
      cat("Directory: ", path, "\n")
      postscript(file=paste(path, file, sep=""), width=width, height=height, onefile=TRUE, horizontal=horizontal)
      clusterplot(obj=obj, ng=ng, title=title, span=span, degree=degree, family=family)
      dev.off()
   } else if (device == "PDF") {
      path <- normalizePath(path=paste(path, "/", sep=""), winslash="\\", mustWork=FALSE)
      file <- paste(file, ".pdf", sep="")
      cat("\nOUTPUT: \n")
      cat("Filename : ", file, "\n")
      cat("Directory: ", path, "\n")
      pdf(file=paste(path, file, sep=""), width=width, height=height, onefile=TRUE, paper=ifelse(test=horizontal, yes="USr", no="US"))
      clusterplot(obj=obj, ng=ng, title=title, span=span, degree=degree, family=family)
      dev.off()
   } else {
      stop("Currently allowed display device are \"PS\" (Postscript) or \"PDF\" (Portable Document Format) \n")
   }

}
##########################################################################################################################################




##########################################################################################################################################
################
# Usage         :
################
#                   target.diagnostic(obj,
#                                     title="Target Moments Diagnostic Plots",
#                                     device=NULL,
#                                     file="Target Moments Diagnostic Plots",
#                                     path=getwd(),
#                                     horizontal=FALSE,
#                                     width=8.5,
#                                     height=6.5, ...)
#
################
# Description   :
################
#
################
# Arguments     :
################
#
################
# Values        :
################
#
################
#
##########################################################################################################################################

target.diagnostic <- function(obj,
                              title="Target Moments Diagnostic Plots",
                              device=NULL,
                              file="Target Moments Diagnostic Plots",
                              path=getwd(),
                              horizontal=FALSE,
                              width=8.5,
                              height=6.5, ...) {

   targetplot <- function(obj, title, ...) {
      par(mfrow=c(2, 3), oma=c(0, 0, 3, 0), mar=c(4, 3, 3, 1), mgp=c(2, 0.5, 0), xpd=FALSE)

      block <- factor(x=obj$block, levels=unique(as.character(obj$block)))
      lev <- levels(block)
      ng <- nlevels(block)
      n <- nrow(obj$Xmvr)
      p <- ncol(obj$Xmvr)
      data.raw <- obj$Xraw
      data.mvr <- obj$Xmvr

      omu <- pooled.mean(x=data.mvr, block=block)
      osd <- pooled.sd(x=data.mvr, block=block)
      X <- matrix(data=rnorm(n=n*p, mean=0, sd=1), ncol=p)
      emu <- pooled.mean(x=X, block=block)
      esd <- (pooled.sd(x=X, block=block) / sqrt(n-ng)) * sqrt(rchisq(n=p, df=n-ng))

      if (p > 50) {
         i <- sample(x=1:p, size=50)
      } else {
         i <- 1:p
      }
      p.tt.mu <- t.test(x=omu[i], y=emu[i], alternative="two.sided")$p.value
      p.tt.sd <- t.test(x=osd[i], y=esd[i], alternative="two.sided")$p.value
      p.ks.mu <- ks.test(x=omu, y=emu, alternative="two.sided")$p.value
      p.ks.sd <- ks.test(x=osd, y=esd, alternative="two.sided")$p.value

      dens <- density(pooled.mean(data.raw, block=block), na.rm=TRUE)
      plot(dens, type="l", col=1, lty=1, cex.main=1,
           xlim=range(0, dens$x), ylim=range(0, dens$y),
           xlab="Means", ylab="Density", main="Untransformed")
      abline(v=0, col=2, lty=2)
      abline(v=mean(pooled.mean(data.raw, block=block), na.rm=TRUE), col=1, lty=4)
      legend("topright", inset=0.05, legend="observed", lty=4, col=1, cex=0.7)

      dens.std <- density(pooled.mean(data.mvr, block=block), na.rm=TRUE)
      plot(dens.std, type="l", col=1, lty=1, cex.main=1,
           xlim=range(0, dens.std$x), ylim=range(0, dens.std$y),
           xlab="Means", ylab="Density", main="Mean-Variance Regularization")
      abline(v=0, col=2, lty=2)
      abline(v=mean(pooled.mean(data.mvr, block=block), na.rm=TRUE), col=1, lty=4)
      legend(x="topright", inset=0.05, legend=c("expected", "observed", paste("p=", round(p.tt.mu, 4), sep="")),
             lty=c(2,4,NA), col=c(2,1,NA), cex=0.7)

      qqplot(x=emu, y=omu, col=1, pch=".", cex=3, cex.main=1,
             main="QQ plot",
             xlab="Theoretical (Normal) Pooled Mean Quantiles",
             ylab="Observed Transformed Pooled Mean Quantiles")
      segments(x0=0, x1=0, y0=min(omu), y1=0, col="black", lty=2, lwd=0.5)
      segments(x0=min(emu), x1=0, y0=0, y1=0, col="black", lty=2, lwd=0.5)
      abline(a=0, b=1, lty=1, lwd=0.5, col="red")
      legend(x="topleft", inset=0.05, legend=paste("p=", round(p.ks.mu, 4), sep=""), cex=0.7)

      dens <- density(pooled.sd(data.raw, block=block), na.rm=TRUE)
      plot(dens, type="l", col=1, lty=1, cex.main=1,
           xlim=range(1, dens$x), ylim=range(0, dens$y),
           xlab="Standard Deviations", ylab="Density", main="Untransformed")
      abline(v=1, col=2, lty=2)
      abline(v=mean(pooled.sd(data.raw, block=block), na.rm=TRUE), col=1, lty=4)
      legend(x="topright", inset=0.05, legend="observed", lty=4, col=1, cex=0.7)

      dens.std <- density(pooled.sd(data.mvr, block=block), na.rm=TRUE)
      plot(dens.std, type="l", col=1, lty=1, cex.main=1,
           xlim=range(1, dens.std$x), ylim=range(0, dens.std$y),
           xlab="Standard Deviations", ylab="Density", main="Mean-Variance Regularization")
      abline(v=1, col=2, lty=2)
      abline(v=mean(pooled.sd(data.mvr, block=block), na.rm=TRUE), col=1, lty=4)
      legend(x="topright", inset=0.05, legend=c("expected", "observed", paste("p=", round(p.tt.sd, 4), sep="")),
             lty=c(2,4,NA), col=c(2,1,NA), cex=0.7)

      qqplot(x=esd, y=osd, col=1, pch=".", cex=3, cex.main=1,
             main="QQ plot",
             xlab="Theoretical (Chi) Pooled Std. Dev. Quantiles",
             ylab="Observed Transformed Pooled Std. Dev. Quantiles")
      segments(x0=1, x1=1, y0=min(osd), y1=1, col="black", lty=2, lwd=0.5)
      segments(x0=min(esd), x1=1, y0=1, y1=1, col="black", lty=2, lwd=0.5)
      abline(a=0, b=1, lty=1, lwd=0.5, col="red")
      legend(x="topleft", inset=0.05, legend=paste("p=", round(p.ks.sd, 4), sep=""), cex=0.7)

      mtext(text=title, cex=1, side=3, outer=TRUE)
   }

   if (is.null(device)) {
      dev.new(width=width, height=height, title="Target Moments Diagnostic Plots", noRStudioGD = TRUE)
      targetplot(obj=obj, title=title)
   } else if (device == "PS") {
      path <- normalizePath(path=paste(path, "/", sep=""), winslash="\\", mustWork=FALSE)
      file <- paste(file, ".ps", sep="")
      cat("\nOUTPUT: \n")
      cat("Filename : ", file, "\n")
      cat("Directory: ", path, "\n")
      postscript(file=paste(path, file, sep=""), width=width, height=height, onefile=TRUE, horizontal=horizontal)
      targetplot(obj=obj, title=title)
      dev.off()
   } else if (device == "PDF") {
      path <- normalizePath(path=paste(path, "/", sep=""), winslash="\\", mustWork=FALSE)
      file <- paste(file, ".pdf", sep="")
      cat("\nOUTPUT: \n")
      cat("Filename : ", file, "\n")
      cat("Directory: ", path, "\n")
      pdf(file=paste(path, file, sep=""), width=width, height=height, onefile=TRUE, paper=ifelse(test=horizontal, yes="USr", no="US"))
      targetplot(obj=obj, title=title)
      dev.off()
   } else {
      stop("Currently allowed display device are \"PS\" (Postscript) or \"PDF\" (Portable Document Format) \n")
   }

}
##########################################################################################################################################




##########################################################################################################################################
################
# Usage         :
################
#                   stabilization.diagnostic(obj,
#                                            span=0.5,
#                                            degree=2,
#                                            family="gaussian",
#                                            title="Stabilization Diagnostic Plots",
#                                            device=NULL,
#                                            file="Stabilization Diagnostic Plots",
#                                            path=getwd(),
#                                            horizontal=FALSE,
#                                            width=7,
#                                            height=5, ...)
#
################
# Description   :
################
#
################
# Arguments     :
################
#
################
# Values        :
################
#
################
#
##########################################################################################################################################

stabilization.diagnostic <- function(obj,
                                     span=0.5,
                                     degree=2,
                                     family="gaussian",
                                     title="Stabilization Diagnostic Plots",
                                     device=NULL,
                                     file="Stabilization Diagnostic Plots",
                                     path=getwd(),
                                     horizontal=FALSE,
                                     width=7,
                                     height=5, ...) {

   stabplot <- function(obj, title, span, degree, family, ...) {
      par(mfrow=c(1, 2), oma=c(0, 0, 3, 0), mar=c(4, 3, 3, 1), mgp=c(2, 0.5, 0), xpd=FALSE)

      block <- factor(x=obj$block, levels=unique(as.character(obj$block)))
      data.raw <- obj$Xraw
      data.mvr <- obj$Xmvr

      ord <- order(pooled.mean(x=data.raw, block=block))
      psd <- pooled.sd(x=data.raw, block=block)[ord]
      low <- loess(psd ~ as.numeric(1:ncol(data.raw)), span=span, degree=degree, family=family)$fitted
      plot(x=1:ncol(data.raw), y=psd, type="p",
           main="Variance-Mean Plot of Untransformed Data",
           xlab="rank(pooled mean)", ylab="pooled sd", col=2, pch=".", cex=3, cex.main=0.7)
      lines(low, col=1, lty=5, lwd=2)

      ord <- order(pooled.mean(x=data.mvr, block=block))
      psd <- pooled.sd(x=data.mvr, block=block)[ord]
      low <- loess(psd ~ as.numeric(1:ncol(data.mvr)), span=span, degree=degree, family=family)$fitted
      plot(x=1:ncol(data.mvr), y=psd, type="p",
           main="Variance-Mean Plot of MVR-transformed Data",
           xlab="rank(pooled mean)", ylab="pooled sd", col=2, pch=".", cex=3, cex.main=0.7)
      lines(low, col=1, lty=5, lwd=2)

      mtext(text=title, cex=1, side=3, outer=TRUE)
   }

   if (is.null(device)) {
      dev.new(width=width, height=height, title="Stabilization Diagnostic Plots", noRStudioGD = TRUE)
      stabplot(obj=obj, title=title, span=span, degree=degree, family=family)
   } else if (device == "PS") {
      path <- normalizePath(path=paste(path, "/", sep=""), winslash="\\", mustWork=FALSE)
      file <- paste(file, ".ps", sep="")
      cat("\nOUTPUT: \n")
      cat("Filename : ", file, "\n")
      cat("Directory: ", path, "\n")
      postscript(file=paste(path, file, sep=""), width=width, height=height, onefile=TRUE, horizontal=horizontal)
      stabplot(obj=obj, title=title, span=span, degree=degree, family=family)
      dev.off()
   } else if (device == "PDF") {
      path <- normalizePath(path=paste(path, "/", sep=""), winslash="\\", mustWork=FALSE)
      file <- paste(file, ".pdf", sep="")
      cat("\nOUTPUT: \n")
      cat("Filename : ", file, "\n")
      cat("Directory: ", path, "\n")
      pdf(file=paste(path, file, sep=""), width=width, height=height, onefile=TRUE, paper=ifelse(test=horizontal, yes="USr", no="US"))
      stabplot(obj=obj, title=title, span=span, degree=degree, family=family)
      dev.off()
   } else {
      stop("Currently allowed display device are \"PS\" (Postscript) or \"PDF\" (Portable Document Format) \n")
   }

}
##########################################################################################################################################




##########################################################################################################################################
################
# Usage         :
################
#                   normalization.diagnostic(obj,
#                                            pal,
#                                            title="Normalization Diagnostic Plots",
#                                            device=NULL,
#                                            file="Normalization Diagnostic Plots",
#                                            path=getwd(),
#                                            horizontal=FALSE,
#                                            width=7,
#                                            height=8, ...)
#
################
# Description   :
################
#
################
# Arguments     :
################
#
################
# Values        :
################
#
################
#
##########################################################################################################################################

normalization.diagnostic <- function(obj,
                                     pal,
                                     title="Normalization Diagnostic Plots",
                                     device=NULL,
                                     file="Normalization Diagnostic Plots",
                                     path=getwd(),
                                     horizontal=FALSE,
                                     width=7,
                                     height=8, ...) {

   normplot <- function(obj, title, ...) {
      par(mfrow=c(2, 2), oma=c(0, 0, 3, 0), mar=c(4, 3, 3, 1), mgp=c(2, 0.5, 0), xpd=FALSE)

      block <- factor(x=obj$block, levels=unique(as.character(obj$block)))
      lev <- levels(block)
      ng <- nlevels(block)
      def <- vector(mode="list", length=ng)
      tab <- numeric(ng)
      for (g in 1:ng) {
         def[[g]] <- which(block == lev[g])
         tab[g] <- length(def[[g]])
      }
      n <- nrow(obj$Xmvr)
      p <- ncol(obj$Xmvr)
      data.raw <- obj$Xraw
      data.mvr <- obj$Xmvr
      x.tick <- seq(from=min(1:n), to=max(1:n), by=1)
      y.tick <- seq(from=min(1:p), to=max(1:p), by=min(c(p,100)))

      boxplot(as.data.frame(t(data.raw)), col=rep(2:(ng+1), tab),
              names=block, notch=FALSE, cex.axis=0.7, cex.main=1, las=2, log="",
              main="Boxplot of Untransformed Data", xlab="Samples", ylab="intensity")

      boxplot(as.data.frame(t(data.mvr)), col=rep(2:(ng+1), tab),
              names=block, notch=FALSE, cex.axis=0.7, cex.main=1, las=2, log="",
              main="Boxplot of MVR-transformed Data", xlab="Samples", ylab="intensity")

      image(x=1:n, y=1:p, z=data.raw[,rev(1:p)], col=pal, axes=F, asp=NA, xlab="", ylab="")
      title(main="Heatmap of Untransformed Data", cex.main=1, line=2)
      axis(side=1, at=(1:n)[x.tick], labels=colnames(data.raw)[x.tick], pos=-0.5, cex.axis=0.5, las=2)
      axis(side=2, at=(1:p)[p+1-y.tick], labels=as.character((1:p)[y.tick]), pos=(1:n)[1]-0.5, cex.axis=0.5, las=2)
      mtext(text="Variables", side=2, outer=F, cex=0.7, col=1, line=2)
      mtext(text="Samples", side=1, outer=F, cex=0.7, col=1, line=2)
      box(which="plot")

      image(x=1:n, y=1:p, z=data.mvr[,rev(1:p)], col=pal, axes=F, asp=NA, xlab="", ylab="")
      title(main="Heatmap of MVR-transformed Data", cex.main=1, line=2)
      axis(side=1, at=(1:n)[x.tick], labels=colnames(data.mvr)[x.tick], pos=-0.5, cex.axis=0.5, las=2)
      axis(side=2, at=(1:p)[p+1-y.tick], labels=as.character((1:p)[y.tick]), pos=(1:n)[1]-0.5, cex.axis=0.5, las=2)
      mtext(text="Variables", side=2, outer=F, cex=0.7, col=1, line=2)
      mtext(text="Samples", side=1, outer=F, cex=0.7, col=1, line=2)
      box(which="plot")

      mtext(text=title, cex=1, side=3, outer=TRUE)
   }

   if (is.null(device)) {
      dev.new(width=width, height=height, title="Normalization Diagnostic Plots", noRStudioGD = TRUE)
      normplot(obj=obj, title=title)
   } else if (device == "PS") {
      path <- normalizePath(path=paste(path, "/", sep=""), winslash="\\", mustWork=FALSE)
      file <- paste(file, ".ps", sep="")
      cat("\nOUTPUT: \n")
      cat("Filename : ", file, "\n")
      cat("Directory: ", path, "\n")
      postscript(file=paste(path, file, sep=""), width=width, height=height, onefile=TRUE, horizontal=horizontal)
      normplot(obj=obj, title=title)
      dev.off()
   } else if (device == "PDF") {
      path <- normalizePath(path=paste(path, "/", sep=""), winslash="\\", mustWork=FALSE)
      file <- paste(file, ".pdf", sep="")
      cat("\nOUTPUT: \n")
      cat("Filename : ", file, "\n")
      cat("Directory: ", path, "\n")
      pdf(file=paste(path, file, sep=""), width=width, height=height, onefile=TRUE, paper=ifelse(test=horizontal, yes="USr", no="US"))
      normplot(obj=obj, title=title)
      dev.off()
   } else {
      stop("Currently allowed display device are \"PS\" (Postscript) or \"PDF\" (Portable Document Format) \n")
   }

}
##########################################################################################################################################




##########################################################################################################################################
# 4. OTHER END-USER FUNCTIONS
##########################################################################################################################################

##########################################################################################################################################
#################
#Usage         :
################
#                   MVR.news(...)
#
################
# Description   :
################
#
################
# Arguments     :
################
#
################
# Values        :
################
#
##########################################################################################################################################

MVR.news <- function(...) {
   newsfile <- file.path(system.file(package="MVR"), "NEWS")
   file.show(newsfile)
}
##########################################################################################################################################