R/plot.sim.R

In lava: Latent Variable Models

##' @export
##' @export density.sim
density.sim <- function(x, ..., plot.type="single") {
    plot.sim(x, ..., scatter.plot=FALSE, plot.type=plot.type)
}

##' Plot method for simulation 'sim' objects
##'
##' Density and scatter plots
##' @examples
##' n <- 1000
##' val <- cbind(est1=rnorm(n,sd=1),est2=rnorm(n,sd=0.2),est3=rnorm(n,1,sd=0.5),
##'              sd1=runif(n,0.8,1.2),sd2=runif(n,0.1,0.3),sd3=runif(n,0.25,0.75))
##' 
##' plot.sim(val,estimate=c(1,2),true=c(0,0),se=c(4,5),equal=TRUE,scatter.plot=TRUE)
##' plot.sim(val,estimate=c(1,3),true=c(0,1),se=c(4,6),xlim=c(-3,3),
##' 	scatter.ylim=c(-3,3),scatter.plot=TRUE)
##' plot.sim(val,estimate=c(1,2),true=c(0,0),se=c(4,5),equal=TRUE,
##' 	plot.type="single",scatter.plot=TRUE)
##' plot.sim(val,estimate=c(1),se=c(4,5,6),plot.type="single",scatter.plot=TRUE)
##' plot.sim(val,estimate=c(1,2,3),equal=TRUE,scatter.plot=TRUE)
##' plot.sim(val,estimate=c(1,2,3),equal=TRUE,byrow=TRUE,scatter.plot=TRUE)
##' plot.sim(val,estimate=c(1,2,3),plot.type="single",scatter.plot=TRUE)
##' plot.sim(val,estimate=1,se=c(3,4,5),plot.type="single",scatter.plot=TRUE)
##' 
##' density.sim(val,estimate=c(1,2,3),density=c(0,10,10), lwd=2, angle=c(0,45,-45),cex.legend=1.3)
##' @aliases density.sim plot.sim
##' @export
##' @export plot.sim
##' @param x sim object
##' @param estimate columns with estimates
##' @param se columns with standard error estimates
##' @param true (optional) vector of true parameter values
##' @param names (optional) names of estimates
##' @param auto.layout Auto layout (default TRUE)
##' @param byrow Add new plots to layout by row
##' @param type plot type
##' @param ask if TRUE user is asked for input, before a new figure is drawn
##' @param col colour (for each estimate)
##' @param pch plot symbol
##' @param cex point size
##' @param lty line type
##' @param lwd line width
##' @param legend legend 
##' @param legendpos legend position
##' @param cex.legend size of legend text
##' @param plot.type 'single' or 'multiple' (default)
##' @param polygon if TRUE fill the density estimates with colour
##' @param density if non-zero add shading lines to polygon
##' @param angle shading lines angle of polygon
##' @param cex.axis Font size on axis
##' @param alpha Semi-transparent level (1: non-transparent, 0: full)
##' @param main Main title
##' @param cex.main Size of title font
##' @param equal Same x-axis and y-axis for all plots
##' @param delta Controls the amount of space around axis limits
##' @param ylim y-axis limits
##' @param xlim x-axis limits
##' @param ylab y axis label
##' @param xlab x axis label
##' @param rug if TRUE add rug representation of data to x-axis
##' @param rug.alpha rug semi-transparency level
##' @param line.col line colour (running mean, only for scatter plots)
##' @param line.lwd line width (running mean, only for scatter plots)
##' @param line.lty line type (running mean, only for scatter plots)
##' @param line.alpha line transparency 
##' @param scatter.ylab y label for density plots
##' @param scatter.ylim y-axis limits for density plots
##' @param scatter.xlim x-axis limits for density plots
##' @param scatter.alpha semi-transparency of scatter plot
##' @param scatter.col scatter plot colour
##' @param border border colour of density estimates
##' @param true.lty true parameter estimate line type
##' @param true.col true parameter colour
##' @param true.lwd true parameter line width
##' @param density.plot if TRUE add density plot
##' @param scatter.plot if TRUE add scatter plot 
##' @param running.mean if TRUE add running average estimate to scatter plot
##' @param ... additional arguments to lower level functions
plot.sim <- function(x,estimate,se=NULL,true=NULL,
             names=NULL,
             auto.layout=TRUE,
             byrow=FALSE,
             type="p",
             ask=grDevices::dev.interactive(),             
             col=c("gray60","orange","darkblue","seagreen","darkred"),
             pch=16,cex=0.5,lty=1,lwd=0.3,
             legend,
             legendpos="topleft",
             cex.legend=0.8,
             plot.type=c("multiple","single"),
             polygon=TRUE,
             density=0,
             angle=-45,
             cex.axis=0.8,
             alpha=0.2,
             main,
             cex.main=1,
             equal=FALSE,
             delta=1.15,
             ylim=NULL,
             xlim=NULL,
             ylab="",
             xlab="",
             rug=FALSE,
             rug.alpha=0.5,
             line.col=scatter.col,             
             line.lwd=1,
             line.lty=1,
             line.alpha=1,
             scatter.ylab="Estimate",
             scatter.ylim=NULL,
             scatter.xlim=NULL,
             scatter.alpha=0.5,
             scatter.col=col,
             border=col,
             true.lty=2,true.col="gray70",true.lwd=1.2,
             density.plot=TRUE,
             scatter.plot=FALSE,
             running.mean=scatter.plot,
             ...) {

    if (missing(estimate)) {
        estimate <- seq(ncol(x))
    }
    if (is.null(estimate)) {
        av <- apply(x[,drop=FALSE],2,function(z) cumsum(z)/seq(length(z)))
        graphics::matplot(x,type="p",pch=pch, cex=cex, col=col,...)
        graphics::matlines(av,type="l",col=col,lty=lty,...)
        if (!is.null(true)) abline(h=true,lty=true.lty,...)
        if (missing(legend)) legend <- colnames(x)
        if (!is.null(legend))
            graphics::legend(legendpos,legend=legend,bg="white",
                             col=scatter.col,lty=lty,pch=pch,...)
        return(invisible(NULL))
    }
    if (is.character(estimate)) {
        estimate <- match(estimate,colnames(x))
    }

    K <- length(estimate)
    est <- tru <- c()
    if (length(se)>0) {
        if (K==1 && !is.list(se))
            se <- list(se)
        else se <- as.list(se)
    } else {
        est <- estimate; tru <- true
    }

    for (i in seq_along(estimate)) {
        est <- c(est,list(rep(estimate[i],length(se[[i]]))))
        if (!is.null(true)) tru <- c(tru,list(rep(true[i],length(se[[i]]))))
    }

    if (length(se)>0) {
      for (i in seq_along(se)) {
            if (is.character(se[[i]])) se[[i]] <- match(se[[i]],colnames(x))
        }

    }
  ss <- summary.sim(x,estimate=unlist(est),se=unlist(se),true=unlist(tru),names=names)

    oldpar <- NULL
    on.exit({
        par(oldpar)
        return(invisible(ss))
    })

    single <- tolower(plot.type[1])=="single"

    nc <- (scatter.plot || running.mean) + density.plot
    nr <- min(6,K)
    if (single) nr <- 1
    
    if (auto.layout && (nc>1 || nr>1)) {
        oma.multi = c(2, 0, 2, 0)
        mar.multi = c(1.5, 4.1, 1, 1)
        opt <- list(mar=mar.multi, oma=oma.multi,
                      cex.axis=cex.axis,las=1,
                      ask=FALSE)
        if (byrow) {
            opt <- c(opt, list(mfrow=c(nr,nc)))
        } else {
            opt <- c(opt, list(mfcol=c(nc,nr)))
        }
        oldpar <- do.call(par, opt)
    }
  dys <- c()
  maxdy <- 0

  if (density.plot)
    for (i in seq(K)) {
            ii <- estimate[i]
            y <- as.vector(x[,ii,drop=TRUE])
            dy <- stats::density(y,...)
            dys <- c(dys,list(dy))
            maxdy <- max(maxdy,dy$y)
        }

    if (equal || single) {
        if (is.null(scatter.ylim)) {
            rg <- range(x[,estimate])
            rg <- rg+c(-1,1)*abs(diff(rg)*(delta-1))
            scatter.ylim <-  rep(list(rg),K)
        }
        if (density.plot) {
            if (is.null(ylim)) ylim <- rep(list(c(0,maxdy*delta)),K)
            if (is.null(xlim)) xlim <- scatter.ylim
        }
    }

    if (!is.null(ylim)) {
        if (!is.list(ylim)) ylim <- list(ylim)
        ylim <- rep(ylim,length.out=K)
    }
    ylab <- rep(ylab,length.out=K)
    if (!is.null(scatter.ylim)) {
        if (!is.list(scatter.ylim)) scatter.ylim <- list(scatter.ylim)
        scatter.ylim <- rep(scatter.ylim,length.out=K)
    }
    if (!is.null(xlim)) {
        if (!is.list(xlim)) xlim <- list(xlim)
        xlim <- rep(xlim,length.out=K)
    }
    if (missing(main)) {
        main <- NULL
        if (!missing(names)) main <- names
        else if (K>1 && !single) main <- colnames(ss)
    }
    if (!is.null(main)) main <- rep(main,length.out=K)
    if (missing(lty)) {
        lty <- rep(1,K)
        if (single || !polygon) {
            lty <- 1:20
        }
    }

    my.scatter.sim <- function(i,add=FALSE,colors,...) {
        ii <- estimate[i]
        if (!missing(colors)) {
            scatter.col <- line.col <- true.col <- colors[1]
        }
        y <- as.vector(x[,ii,drop=TRUE])
        args <- list(y,ylab=scatter.ylab[i],col=Col(scatter.col[1],scatter.alpha),cex=cex,pch=pch,type=type)
        if (!is.null(scatter.ylim)) args <- c(args,list(ylim=scatter.ylim[[i]]))
        if (scatter.plot) {
            if (!add) {
                do.call(graphics::plot,args)
            } else {
                do.call(graphics::points,args)
            }
        }
        if (running.mean) {
            lines(cumsum(y)/seq_along(y),col=line.col[1],lwd=line.lwd,lty=line.lty)
            if (!is.null(true))
                abline(h=true[i],lty=true.lty,col=true.col[1],lwd=true.lwd)
        }
    }

    my.density.sim <- function(i,add=FALSE,colors,
                               alphas=alpha,
                               auto.legend=TRUE,
                               densities=NULL,
                               angles=angle,
                               ...) {
        ii <- estimate[i]
        y <- as.vector(x[,ii,drop=TRUE])
        if (!missing(colors)) {
            col <- border <- colors
            col <- true.col <- colors
        }
        if (density.plot) {
            dy <- stats::density(y)
            if (is.null(ylim)) {
                density.ylim0 <- c(0,max(dy$y)*delta)
            } else {
                density.ylim0 <- ylim[[i]]
            }
            if (is.null(xlim)) {
                density.xlim0 <- range(dy$x)
            } else {
                density.xlim0 <- xlim[[i]]
            }
            if (!add) graphics::plot(0, 0, type="n",
                                     main="", ylab=ylab, xlab=xlab,
                                     ylim=density.ylim0, xlim=density.xlim0)
            if (polygon) {
              with(dy, graphics::polygon(c(x,rev(x)),
                                         c(y,rep(0,length(y))),
                                         col=Col(col[1], alpha=alphas[1]),
                                         border=NA,density=densities[1],angle=angles[1]))
                if (!is.null(border)) with(dy, lines(x,y,col=border[1],lty=lty[1],lwd=lwd[1]))
            } else {
                graphics::lines(dy,main="",lty=lty[1],col=col[1],lwd=lwd[1])
            }
            if (rug) graphics::rug(y,col=Col(col[1],rug.alpha[1]))
            if (!is.null(main) && !(running.mean || scatter.plot)) {
                title(main[i],cex.main=cex.main)
            }
            if (!is.null(true)) {
                abline(v=true[i],lty=true.lty,col=true.col,lwd=true.lwd)
            }
            if (!is.null(se) && !any(is.na(se[[i]]))) {
                se.pos <- match(se[[i]],unlist(se))
                ns <- length(se.pos)+1
                se.alpha <- rep(alphas,length.out=ns)[-1]
                se.border <- rep(border,length.out=ns)[-1]
                se.col <- rep(col,length.out=ns)[-1]
                se.lty <- rep(lty,length.out=ns)[-1]
                se.lwd <- rep(lwd,length.out=ns)[-1]
                xx <- dy$x
                for (j in seq_along(se.pos)) {
                    if (polygon) {
                        yy <- dnorm(xx,mean=ss["Mean",se.pos[j]],sd=ss["SE",se.pos[j]])
                        if (se.alpha[j]>0) graphics::polygon(c(xx,rev(xx)),c(yy,rep(0,length(yy))),col=Col(se.col[j],alpha=se.alpha[j]),border=NA,density=densities[j],angle=angles[j])
                        if (!is.null(border)) lines(xx,yy,col=se.border[j],lty=se.lty[j],lwd=se.lwd[j])
                    } else {
                        graphics::curve(dnorm(x,mean=ss["Mean",se.pos[j]],sd=ss["SE",se.pos[j]]),lwd=se.lwd[j],lty=se.lty[j],col=se.col[j],add=TRUE)
                    }
                 }
                if (auto.legend) legend <- c("Kernel",colnames(ss)[se.pos])
                if (!is.null(legend)) {
                    if (polygon) {
                        dcol <- c(col[1],se.col)
                        fill <- Col(dcol,alpha)
                        fill[which(alpha==0)] <- NA
                        border[which(alpha==0)] <- NA
                        dcol[which(alpha!=0)] <- NA
                        graphics::legend(legendpos,legend,
                                         fill=fill,
                                         lty=1,
                                         col=dcol,
                                         border=border,
                                         cex=cex.legend)
                    } else {
                        graphics::legend(legendpos,legend,
                                         col=c(col[1],se.col),
                                         lty=c(lty[1],se.lty),
                                         lwd=c(lwd[1],se.lwd),
                                         cex=cex.legend)
                    }
                }
            }

        }
    }

    if (single) {
        nk <- unlist(lapply(se,length))
        col <- rep(col,length.out=K)
        for (i in seq(K)) {
            my.scatter.sim(i,add=(i>1),colors=col[i])
        }
        if (!is.null(main) && !byrow) {
            title(main[1],cex.main=cex.main)
        }
        if (missing(legend)) legend <- colnames(x)[estimate]
        legendold <- legend
        legend <- NULL
        alpha <- rep(alpha,length.out=K)
        density <- rep(density,length.out=K)
        angle <- rep(angle,length.out=K)
        for (i in seq_len(K)) {
            alphas <- if (K==1L) alpha else alpha[i]
            densities <- density[i]
            if (!is.null(densities) && densities<1) densities <- NULL
            if (length(se)>0) alphas <- c(alphas,rep(0,nk[i]))
            my.density.sim(i, add=(i>1),colors=col[i],alphas=alphas,
                           densities=densities,
                           angles=angle[i],
                           auto.legend=FALSE)
        }
        if (!is.null(legendold)) {
            legend <- rep(legendold,length.out=K)
            graphics::legend(legendpos,legend,
                             fill=Col(col,alpha),border=col,cex=cex.legend)
        }

    } else {
        for (i in seq(K)) {
            my.scatter.sim(i)
            if (!is.null(main) && !byrow && scatter.plot) {
                title(main[i],cex.main=cex.main)
            }
            my.density.sim(i,auto.legend=missing(legend))
            if (i==1 && ask && K>1) par(ask=ask)
        }
    }

}