R/plot_imprecise.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

#' @aliases 
#' plot, plot.imprecise, plot.summary.imprecise
#' 
#' @title 
#' Object Visualization in Imprecise Inferential Framework
#'
#' @description 
#' \code{plot} use xy-coordinate information produced from \code{\link{iprior}}
#' for an imprecise prior and the information produced from \code{\link{update}}
#' for an imprecise posterior.
#'
#' @param x 
#' an object of class of \code{imprecise} or \code{imprecise.summary}.
#'
#' @param ... other passing arguments to \code{print}
#' 
#' @details Visualization of Imprecise Prior and Posterior
#' \itemize{
#' 
#'   \item{circle2d}{if xreg=TRUE, it is suggested to draw a polygon
#' or scatter plot in 2 dimensional space.
#' If xreg=FALSE, there is no available visiualization method; then,
#' use a perspective plot (or level plots) with contours.}
#' 
#'   \item{circle3d}{if xreg=TRUE, it is suggested to draw a
#' 3-dimensional scatter plot.
#' If xreg=FALSE, no methods since a number of hyperparameters is
#' limited to 2 not 3.}
#' 
#'   \item{eqns2d}{if xreg=TRUE, it would be good to draw a
#' two-dimensional polygon or scatter plot.
#' If xreg=FALSE, there is no available visualization method; then,
#' use a perspective plot (or level plots) with contours. }
#' 
#'   \item{eqns3d}{if xreg=TRUE, it would be good to generate a
#' three-dimensional scatter plot.
#' If xreg=FALSE, no methods since a number of hyperparameters is
#' limited to 2 not 3.}
#'
#'   \item{others}{if a number of dimensions is more than 4, there is
#' no particular visualization metho since one more additional
#' dimension is required to present a magnitude of imprecision.
#' In this case, it would be good to generate a pairwise plot with
#' contours.}
#' }
#'
#' @note 
#' The function \code{inout} in the package \code{mgcv} is used. 
#'
#' @seealso 
#' \code{\link{in.out}}, \code{\link{chull}}, \code{\link{panel.3dwire}}
#'
#' @keywords 
#' convex, polygon, polyhedren, circle, sphere
#' 
#' @author Chel Hee Lee <\email{gnustats@@gmail.com}>
#'
#' @method plot imprecise
#' @S3method plot imprecise
plot.imprecise <- function(x, ...){

  # sanity check
  stage <- x$stage
  # stopifnot(stage != "model")
  
  xtms <- do.call(rbind, x$xtms)
  m0shape <- x$m0shape
  
  ## put the messages 
  ## 'pairwise' visualization is required if dimension is greater than 3.
  
  fn.3d <- function(...){
    # Scatterplot in a single 3D panel
    stopifnot(m0shape %in% c("eqns3d", "sphere3d"))
    xtms <- as.data.frame(xtms)
    names(xtms) <- c("x", "y", "z")
    tobj <- cloud(z~x*y, data=xtms, pretty=TRUE,
      scales=list(cex=1.5, arrows=FALSE, col="black"),
      par.settings = list(axis.line = list(col = "transparent")),
      par.box=list(col='black', lwd=0.25), 
      ...)
    print(tobj)
    
    # pairwise plot
    
  }
  
  fn.2d <- function(...){
    stopifnot(m0shape %in% c("eqns2d", "circle2d"))
    plot(x=xtms, ...)
    polygon(xtms[chull(xtms),], col="azure2", border="darkblue", lty="dashed", lwd=3, ...)
    points(xtms, col="red", pch=19, cex=0.5)
    text(x=xtms, labels=rownames(xtms), pos=4)
  }

  switch(m0shape, 
    "eqns2d"=fn.2d(...),
    "circle2d"=fn.2d(...),
    "eqns3d"=fn.3d(...),
    "sphere3d"=fn.3d(...))
  invisible(x)
}
NULL

#' @rdname plot.imprecise
#' @method plot summary.imprecise
#' @S3method plot summary.imprecise
plot.summary.imprecise <- function(x, ...){
  
  xreg <- x$xreg
  m0shape <- x$m0shape
  xnames <- if(xreg) colnames(x$X) else NULL

  fn.sphere3dxreg <- function(...){
    # TODO: pairwise scatterplot
    stopifnot(xreg, length(xnames)==3)
    est <- as.data.frame(x$est)
    tobj <- splom(est, type="p", pch=19, cex=0.5, col="darkblue", ...)
    print(tobj)
    
    # TODO: Scatterplot in a single panel
#    stopifnot(m0shape %in% c("eqns3d", "sphere3d"))
#    xtms <- as.data.frame(xtms)
#    names(xtms) <- c("x", "y", "z")
#    tobj <- cloud(z~x*y, data=xtms, pretty=TRUE,
#      scales=list(cex=1.5, arrows=FALSE, col="black"),
#      par.settings = list(axis.line = list(col = "transparent")),
#      par.box=list(col='black', lwd=0.25), 
#      ...)
#    print(tobj)
  }
  
  fn.eqns3dxreg <- function(...){
    stopifnot(xreg, length(xnames)==3)
    est <- na.omit(as.data.frame(x$est))
    tobj <- splom(est, type="p", col="darkblue", ...)
    print(tobj)
  }
  
  fn.2dxreg <- function(M0.rm=FALSE, ref.lines=FALSE, ...){ 
    ## this function covers both 'circle2d' and 'eqns2d'.
  
    stopifnot(xreg, length(xnames)==2)
    xi <- if(is.null(x$xi)) NULL else x$xi
    
    # keep only the names of extreme points (M0) 
    m0xtms <- as.data.frame(do.call(rbind, x$xtms))
    c0idx <- chull(m0xtms)
    m0names <- character(nrow(m0xtms))
    m0names[c0idx] <- paste("x", c0idx, sep="")

    m1xtms <- as.data.frame(x$est)
    if(nrow(m0xtms) != nrow(m1xtms)) stop("Different lengths of 'm0xtms' and 'm1xtms'")
    m1names <- m0names

    infs <- do.call(c, lapply(m1xtms, min))
    sups <- do.call(c, lapply(m1xtms, max))
    idx.inf <- do.call(c, lapply(m1xtms, which.min))
    idx.sup <- do.call(c, lapply(m1xtms, which.max))
    infnames <- supnames <- character(nrow(m1xtms))
    infnames[idx.inf] <- paste("x", idx.inf, sep=".")
    supnames[idx.sup] <- paste("x", idx.sup, sep=".")
    
    if(!M0.rm) plot(x=m0xtms, type="n", ...) else plot(x=m1xtms, ...)
    if(!M0.rm){
      polygon(m0xtms[c0idx,], col=alpha("azure2", 0.5), border="darkblue", lty="dashed", lwd=2)
      points(m0xtms, col="gray", cex=0.75)
      text(m0xtms, labels=m0names, pos=3, cex=0.75)
    }
    polygon(m1xtms[chull(m1xtms),], col="yellow1", border="darkgreen", lty="solid", lwd=2)
    points(m1xtms, col="red1", cex=0.75, pch=19)
    text(m1xtms, infnames, pos=2, cex=0.75)
    text(m1xtms, supnames, pos=4, cex=0.75)
    
    mle <- x$init$cfs
    points(x=mle[1], y=mle[2], pch="*", cex=2)
    
    if (ref.lines) abline(v=infs[1], lty="dashed", col="darkred")
    if (ref.lines) abline(h=infs[2], lty="dashed", col="darkred")
    if (ref.lines) abline(v=sups[1], lty="dashed", col="darkblue")
    if (ref.lines) abline(h=sups[2], lty="dashed", col="darkblue")
    
    if(!is.null(xi)){ # for prediction
      mu <- x$mui

      mu.inf <- which.min(mu)
      mu.sup <- which.max(mu)

      mu.inf.beta <- m1xtms[mu.inf,]
      mu.sup.beta <- m1xtms[mu.sup,]

      abline(a=log(mu[mu.sup])/xi[2], b=-xi[1]/xi[2], lwd=2)
      abline(a=log(mu[mu.inf])/xi[2], b=-xi[1]/xi[2], lwd=2)

      points(x=mu.inf.beta[1], y=mu.inf.beta[2], col="black", cex=2)
      text(x=mu.inf.beta[1], y=mu.inf.beta[2], pos=2, labels=mu.inf)

      points(x=mu.sup.beta[1], y=mu.sup.beta[2], col="black", cex=2)
      text(x=mu.sup.beta[1], y=mu.sup.beta[2], pos=4, labels=mu.sup)
    }
  
  } # end of fn.2dxreg
  
  fn.eqns2d <- function(bnd, clevels=10, M0.rm=FALSE, ...){
    ## without covarites
    stopifnot(!xreg)
    if(missing(bnd)) bnd <- NULL
    # print(bnd)
    # message("'bnd' is printed")
    
    if(is.null(bnd)){
    
      xtms <- do.call(rbind, x$xtms)
      theta <- x$est
      xyz <- as.data.frame(cbind(xtms, theta))
      names(xyz) <- c("x", "y", "z")
      ds <- xyz
    
    } 
    
    if(!is.null(bnd)){
    
#			message("working with given bounds")
    
      grid <- do.call(rbind, x$xtms)
      # print(grid)
      bnd <- do.call(rbind, bnd$xtms)
      # print(bnd)
      xtms <- bnd[chull(bnd),]
      # xtms <- do.call(rbind, bnd$xtms)
      # print(xtms)
      
      idx <- in.out(bnd=xtms, x=grid)
      inner <- grid[idx,]
      theta <- x$est
      theta <- theta[idx]
      innerx <- as.data.frame(cbind(inner, theta))
      ds <- innerx
      names(ds) <- c("x", "y", "z")

    }
    
    # clevels <- seq(from=min(theta), to=max(theta), length.out=10)

    panel.3d.wireframe <- function(x, y, z, xlim, ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled, rot.mat, distance, drape=drape, ...){

      clines <- contourLines(x, y, matrix(z, nrow = length(x), byrow = TRUE), nlevels = clevels) 

      for (ll in clines){
        m <- ltransform3dto3d(rbind(ll$x, ll$y, zlim.scaled[1]), rot.mat, distance)
        panel.lines(m[1,], m[2,])
        panel.text(x=m[1,1], y=m[2,1], lab=ll$level, cex=0.5)
      } 

      xtms <- xtms[chull(xtms),]
      xtms <- rbind(xtms, xtms[1,])
      refBox <- data.frame(x=xtms[,1], y=xtms[,2], z=zlim[1])
      refBox$minz <- min(ds$z)
      refBox$maxz <- max(ds$z)

      panel.3dwire(x=x, y=y, z=z, xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled, zlim.scaled=zlim.scaled, xlim=xlim, ylim=ylim, zlim=zlim, drape=TRUE, rot.mat, distance, ...)
      
      # TODO: 
      # panel.3dscatter(x=ds[,1], y=ds[,2], z=ds[,3], xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled, zlim.scaled=zlim.scaled, xlim=xlim, ylim=ylim, zlim=zlim, .scale=TRUE, rot.mat, distance, type="l", ...)

      # Reference box for contour plot 
      xcoords <- xlim.scaled[1] + diff(xlim.scaled)*(refBox$x - xlim[1])/diff(xlim)
      ycoords <- ylim.scaled[1] + diff(ylim.scaled)*(refBox$y - ylim[1])/diff(ylim)
      zcoords <- zlim.scaled[1] + diff(zlim.scaled)*(refBox$z - zlim[1])/diff(zlim)
      panel.3dscatter(x=xcoords, y=ycoords, z=zcoords, type='l',  rot.mat, distance, ...)

      # Reference line for minimum posterior expectation (Bayes estimate)
      zcoords <- zlim.scaled[1] + diff(zlim.scaled)*(refBox$minz - zlim[1])/diff(zlim)
      panel.3dscatter(x=xcoords, y=ycoords, z=zcoords, type='l', lty=2, rot.mat, distance, ...)

      # Reference line for maximum posterior expectation (Bayes estimate)
      zcoords <- zlim.scaled[1] + diff(zlim.scaled)*(refBox$maxz - zlim[1])/diff(zlim)
      panel.3dscatter(x=xcoords, y=ycoords, z=zcoords, type='l', lty=2, rot.mat, distance, ...)       
    }
    
    tobj <- wireframe(z ~ x*y, data=ds, pretty=TRUE,
      scales=list(arrows=FALSE, col="black"),  # show axese with ticks (z.ticks=2), 
      par.settings = list(axis.line = list(col = "transparent")), # remove the default outer frame lines
      par.box=list(col='black', lwd=0.25), 
      panel.3d.wireframe=panel.3d.wireframe, 
      ...)
    print(tobj, ...)
    rvals <- list(ds=ds, xtms=xtms, clevels=clevels)
  }
  
  if (xreg) {
		x$coords <- switch(m0shape, 
           "eqns2d"=fn.2dxreg(...),
           "circle2d"=fn.2dxreg(...), 
           "eqns3d"=fn.eqns3dxreg(...), 
           "sphere3d"=fn.sphere3dxreg(...)) 
  } else {
		x$coords <- fn.eqns2d(...)
	}
  invisible(x)
}
NULL