R/theta_y_plot_corr.R

In eehh-stanford/yada: Yet Another Demographic Analysis package

#' @title Plot the correlation vector in theta_y across Metropolis samples
#'
#' @description \code{theta_y_plot_corr} plots the correlation vector, s, in theta_y across Metropolis samples given a bundle of relevant variables generated by \code{bundle_theta_y_sample}.
#'
#' @param bundle The bundle of variables (see \code{bundle_theta_y_sample})
#' @param bundle For optional subsetting of variables; first index [optional]
#' @param bundle For optional subsetting of variables; second index [optional]
#' @param decim An integer for sampling (decimating) data to reduce the number of points for very large samples [optional; default = 1 for all samples]
#' @param fileName Filename for pdf save file [optional]
#'
# @keywords
#' @export
#'
#' @author Michael Holton Price <MichaelHoltonPrice@gmail.com>

theta_y_plot_corr <- function(bundle,sub1=NA,sub2=NA,decim=1,fileName=NA) {

  if(!bundle$haveNames) {
    varNames <- sapply(as.character(1:bundle$hp$J),function(s){paste('ord',s,sep='')})
    varNames <- c(varNames,sapply(as.character(1:bundle$hp$K),function(s){paste('cont',s,sep='')}))
  } else {
    varNames <- bundle$varNames
  }

  numVar <- bundle$hp$J+bundle$hp$K
  numCorr <- choose(numVar,2)
  combMat <- combn(numVar,2)
  if(any(is.na(sub1)) || any(is.na(sub2))) {
    sub1 <- combMat[1,]
    sub2 <- combMat[2,]
  }

  legNames <- rep('',length(sub1))
  comb <- rep(NA,length(sub1))
  for(z in 1:length(sub1)) {
    comb[z] <- which(sub1[z] == combMat[1,] & sub2[z] == combMat[2,])
    legNames[z] <- paste(varNames[sub1[z]],varNames[sub2[z]],sep='')
  }

  ind <- seq(1,length(bundle$theta_yList),by=decim)
  toFile <- !is.na(fileName)
  options(warn=-1) # in case length(sub1) < 3
  colVect <- RColorBrewer::brewer.pal(n=length(comb),name='Set1')
  options(warn=0)

  if(length(colVect) < bundle$hp$M[j]) {
    colVect <- colVect[0:(length(comb)-1) %% length(colVect) + 1]
  }

  valList <- lapply(bundle$theta_yList,function(x){covMat2CorrVect(x$Sigma)[comb]})

  yMin <- Inf
  yMax <- -Inf
  yMin <- min(yMin,unlist(valList))
  yMax <- max(yMax,unlist(valList))

  if(bundle$haveKnown) {
    v <- covMat2CorrVect(bundle$known$theta_y$Sigma)[comb]
    yMin <- min(yMin,v)
    yMax <- max(yMax,v)
  }
  if(toFile) {
    pdf(fileName)
  }
  plot(1,type="n",xlab="Sample", ylab='correlation',xlim=c(0,max(ind)),ylim=c(yMin,yMax))
  for(z in 1:length(comb)) {
    points(ind,sapply(valList,function(x){x[z]})[ind],col=colVect[z],pch=20)
    if(bundle$haveKnown) {
      lines(c(0,max(ind)),c(1,1)*v[z],col=colVect[z],lwd=3)
    }
  }

  legend('left',legend=legNames,col=colVect,lty=1,lwd=3)
  if(toFile) {
    dev.off()
  }
}