R/plot.scalescape.R

In benjaminiuliano/scalescape: Distance-Weighted Landscape Analysis

#plot.scalescape
#' Plotting scalescape objects
#'
#' \code{plot} method for classes \code{scalescape} and \code{scalescape.boot}. For objects of
#' class \code{scalescape}, \code{summary.scalescape} returns two plots:
#' \enumerate{
#'   \item the log-likelihood values for the range parameter up to the distance (in meters),
#'    specified as the \code{max.radius} in \code{landscape_matrix}.
#'    \item the shape of the weighting function, i.e.how the weight of  the landscape predictor
#'    decreases with distance (in meters), given the maximum likelihood range value.
#'  }
#' For objects of class \code{scalescape.boot}, \code{summary.scalescape} returns two histograms:
#' \enumerate{
#'   \item the frequency with which different range parameter values were identified as optimal
#'   \item the distribution of deviance values generated by the bootstraps, from which the p-value is derived
#' }
#'
#' @param x an object of class \code{scalescape} or \code{scalescape.boot}
#' @param ... arguments to be passed to methods
#'
#' @export

plot.scalescape <- function(x, ...) {
  if(!is.element("scalescape.boot", class(x))) {
    sol <- weighting_funct(x$opt.range, mod0=x$mod0, landscape.formula = x$landscape.formula,
                           data = x$data, max.Dist = x$max.Dist, landscape.vars=x$landscape.vars, weight.fn = x$weight.fn, return.coef = T)
    par(mfrow=c(length(x$opt.range),2), mai=c(1,1,.3,.3))
    for(i.range in 1:length(x$opt.range)){
      w <- data.frame(range=x$max.Dist[i.range] * .005*(1:199), logLik = NA, b = NA)
      for(i in 1:199) {
        opt.par <- x$opt.range
        opt.par[i.range] <- .005*i
        z <- weighting_funct(par=opt.par, mod0=x$mod0, landscape.formula = x$landscape.formula,
                             data = x$data, max.Dist = x$max.Dist, landscape.vars=x$landscape.vars, weight.fn = x$weight.fn, return.coef = T)
        if(length(z)>1){
          w$logLik[i] <- z$logLik
          w$b[i] <- z$coef[3]
        }
      }

      #1
      plot(logLik ~ range, data=w, typ="l", xlab="Distance", ylab="logLik")
      points(x$opt.range[i.range], x$logLik[i.range], col="red")
      mtext(side=3, paste0("Range for ",names(x$opt.range)," = ",round(x$opt.range[i.range], digits=3)))

      #2
      # Gaussian weightings
      if(x$weight.fn=="Gaussian") curve(exp(-0.5*(x/(x$opt.range[i.range]))^2),
                                        from = 0, to = x$max.Dist[i.range], ylim=c(0,1), xlab="Distance", ylab="Weighting")
      abline(v=x$opt.range[i.range], col="red", lty=2)

      # exponential weightings
      if(x$weight.fn=="exponential") curve(exp(-x/(x$opt.range[i.range])),
                                           from = 0, to = x$max.Dist[i.range], ylim=c(0,1), xlab="Distance", ylab="Weighting")
      abline(v=x$opt.range[i.range], col="red", lty=2)
    }
  }else{
      par(mfcol=c(1,1 + length(x$mod.full$opt.range)), mai=c(1,1,.3,.3))

      #1
      if(length(x$mod.full$opt.range) == 1){
        hist(x$mod.full$max.Dist * x$logLik.values$opt.range.boot,
             xlab="Range", main=paste0("Bootstrap for ", names(x$mod.full$opt.range)))
      }else{
        for(i.range in 1:length(x$mod.full$opt.range))
          hist(x$mod.full$max.Dist[i.range] * x$logLik.values[,paste0("opt.range.boot.",i.range)],
               xlab="Range", main=paste0("Bootstrap for ", names(x$mod.full$opt.range)[i.range]))
      }

      #2
      hist(x$logLik.values$dev, xlab="Deviance", main=paste0("Bootstrap P = ", round(x$P, digits=4)),
           xlim=c(-.01,max(c(x$logLik.values$dev,x$logLik.values$dev))))
      lines(c(x$logLik.values$dev,x$logLik.values$dev), c(0,nrow(x$logLik.values)), col="red")
    }
  }