R/iplot.R

In mpplot: Plots for the R Package MCMCpack

iplot = function(model, focal, fixed, at=NULL, probs=c(0.025, 0.975), ylim=NULL, xlim=NULL, rug=TRUE, pcol=NULL, xlab=NULL, ylab=NULL, pch=NULL, las=NULL, ...) {
  # (1) identify variables in model
  mform = formula(attr(model, "call"))
  mterm = terms(mform)
  mvars = as.character(attr(mterm, "variables"))
  nvars = length(mvars)
  mvars = mvars[3:nvars]
  nvars = nvars - 2
  mdata = as.character(attr(model, "call"))[3] # dataset used
  # data type of focal and fixed variable:
  tyfoc = class(get(focal, pos=get(mdata)))
  tyfix = class(get(fixed, pos=get(mdata)))
  if(tyfoc=="factor") {tyfoc="character"}
  if(tyfix=="factor") {tyfix="character"}
  # switch if necessary
  if(tyfoc=="character" & tyfix=="numeric") {
    warning("Note: The provided focal variable is categorial and the fixed variable is continuous. Changing focal variable and fixed variable.")
    focal2 = focal
    focal  = fixed
    fixed  = focal2
    }
  # (2) mean (or first) values in data for all variables
  #mvalu = as.list(rep(NA, nvars)) # will be replaced
  mvalu = rep(NA, nvars) # will be replaced
  dtype = rep(NA, nvars)
  for(i in 1:nvars) {
    curvar = get(mvars[i], pos=get(mdata))
    if(is.factor(curvar) | is.character(curvar)) {
      mvalu[i] = as.character(curvar[1]) # first if categorical
      dtype[i] = "character" # explicit so that we can send it to mpredict()
      } else {
      mvalu[i] = as.numeric(mean(curvar, na.rm=TRUE))
      dtype[i] = "numeric"  # explicit so that we can send it to mpredict()
      } # mean if continuous
    }
  newd = data.frame(matrix(mvalu, nrow=1, byrow=FALSE)) #, col.names=mvars)
  colnames(newd) = mvars
  for(i in 1:nvars) {
    if(dtype[i]=="numeric") {newd[i] = as.numeric(newd[i])}
    if(dtype[i]=="character") {newd[i] = as.character(newd[i])}
    }
    #sapply(newd, class)
  # (3) predict
  pvlo = rep(NA, 3)
  pvhi = rep(NA, 3)
  ndmod = newd
  if(tyfoc=="character" & tyfix=="character") {
      # focal and fixed are characters
      if(is.null(at)) {
        allfix = unique(get(fixed, pos=get(mdata)))
        lenfix = length(allfix)
        xrange = c(allfix[1], allfix[lenfix])
        } else {
        xrange = at
        }
      ndmod[fixed] = xrange[1] # at (high)
      pvlo = mpredict(model, newdata=ndmod, probs=probs)
      ndmod[fixed] = xrange[2] # at (high)
      pvhi = mpredict(model, newdata=ndmod, probs=probs)
      # set range if not specified
      if(is.null(ylim)) {
        plim = range(c(pvlo, pvhi))
          } else {
          plim=ylim
        }
      } else {
      # focal is numeric, fixed is numeric or character
      if(is.null(xlim)) {
        frange = range(get(focal, pos=get(mdata)), na.rm=TRUE)
        } else{
        frange = xlim
      }
      if(tyfix=="character"){
        if(is.null(at)){
          ufixed = unique(get(fixed, pos=get(mdata)))
          ulengt = length(ufixed)
          xrange = c(ufixed[1], ufixed[ulengt])
          } else {
          xrange = at
          }
        } else {
        xrange = range(get(fixed, pos=get(mdata)), na.rm=TRUE)
        }
      fsteps = seq(frange[1], frange[2], length.out=50) ## make 50 variable XXX
      if(is.null(at)) {at=xrange}
      for(k in 1:50){
        ndmod[focal] = fsteps[k]
        ndmod[fixed] = xrange[1] # at (low)
        pvlo = rbind(pvlo, mpredict(model, newdata=ndmod, probs=probs))
        ndmod[fixed] = xrange[2] # at (high)
        pvhi = rbind(pvhi, mpredict(model, newdata=ndmod, probs=probs))
        }
        pvlo = pvlo[2:51,]
        pvhi = pvhi[2:51,]
        # set range if not specified
        if(is.null(ylim)) {
          plim = range(c(range(pvlo[,"Mean"]), range(pvhi[,"Mean"])))
          } else {
          plim=ylim
        }
      }
  # set colours if not specified
  if(is.null(pcol)) {pcol=c("#004488", "#DDAA33")} # Tol 2021, High-contrast colours
  ptempl = col2rgb(pcol[1], alpha=TRUE)
  pbackl = rgb(ptempl[1]/256, ptempl[2]/256, ptempl[3]/256, alpha=0.5)
  ptemph = col2rgb(pcol[2], alpha=TRUE)
  pbackh = rgb(ptemph[1]/256, ptemph[2]/256, ptemph[3]/256, alpha=0.5)
  #
  if(is.null(xlab)) {xlab=focal}
  if(is.null(ylab)) {ylab="predicted outcome"}
  if(tyfoc=="character" & tyfix=="character") {
    if(is.null(pch)) {pch=19}
    if(is.null(xlim)) {xlim=c(0.9, 2.1)}
    plot(c(1:2), c(pvlo["Mean"], pvhi["Mean"]), ylim=plim, xlim=xlim, bty="n", ylab=ylab, xlab=xlab, type="p", col=pcol, pch=pch, axes=FALSE)
    segments(1,pvlo[2], 1,pvlo[3], col=pcol[1])
    segments(2,pvlo[2], 2,pvlo[3], col=pcol[2])
    axis(1, at=c(1,2), labels=xrange)
    if(is.null(las)) {las=2}
    axis(2, las=2)
    } else {
    plot(fsteps[1:50], pvlo[,"Mean"], ylim=plim, bty="n", ylab=ylab, xlab=xlab,type="n", col=pcol[1], axes=FALSE, ...); axis(1, ...); axis(2, las=2, ...)
    polygon(fsteps[c(1:50, 50:1)], c(pvlo[,2], rev(pvlo[,3])), col=pbackl, border=NA)
    polygon(fsteps[c(1:50, 50:1)], c(pvhi[,2], rev(pvhi[,3])), col=pbackh, border=NA)
    #
    points(fsteps[1:50], pvlo[,"Mean"], col=pcol[1], type="l")
    points(fsteps[1:50], pvhi[,"Mean"], col=pcol[2], type="l")
    if(rug==TRUE){rug(get(focal, pos=get(mdata)))}
    }
 }