Merge pull request #247 from kenkellner/plotMarginal

Add plotEffects function

7 files changed, 274 insertions, 4 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 5454b14..5a04732 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: unmarked
-Version: 1.2.5.9010
-Date: 2022-12-29
+Version: 1.2.5.9011
+Date: 2023-03-07
 Type: Package
 Title: Models for Data from Unmarked Animals
 Authors@R: c(
@@ -46,6 +46,7 @@ Collate: 'classes.R' 'unmarkedEstimate.R' 'mapInfo.R' 'unmarkedFrame.R'
     'unmarkedCrossVal.R' 'piFun.R' 'vif.R' 'makePiFun.R' 'posteriorSamples.R'
     'nmixTTD.R'
     'gdistremoval.R'
+    'plotEffects.R'
     'mixedModelTools.R'
     'power.R' 
     'simulate.R'
diff --git a/NAMESPACE b/NAMESPACE
index 448be4f..7273ce1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,7 +7,7 @@ importFrom(stats, confint, fitted, coef, vcov, predict, update, profile,
            pnorm, qchisq, qnorm, quantile, rbinom,
            reshape, rmultinom, rnbinom, rpois, runif, sd, uniroot,
            update.formula, sigma)
-importFrom(graphics, plot, hist, abline)
+importFrom(graphics, plot, hist, abline, axis, lines, points, polygon, segments)
 importFrom(utils, head, read.csv)
 importFrom(grDevices, devAskNewPage, dev.interactive, palette.colors)
 importFrom(MASS, mvrnorm)
@@ -51,7 +51,7 @@ exportMethods(backTransform, coef, confint, coordinates, fitted, getData,
               "siteCovs<-", summary, update, vcov, yearlySiteCovs,
               "yearlySiteCovs<-", "[", smoothed, projected, nonparboot, logLik,
               LRT, ranef, bup, crossVal, posteriorSamples, sigma, randomTerms, 
-              optimizePenalty, unmarkedPowerList)
+              optimizePenalty, unmarkedPowerList, plotEffectsData, plotEffects)
 
 S3method("print", "unmarkedPostSamples")
 
diff --git a/R/plotEffects.R b/R/plotEffects.R
new file mode 100644
index 0000000..c77f5aa
--- /dev/null
+++ b/R/plotEffects.R
@@ -0,0 +1,109 @@
+#type_to_covs <- function(umf, type){
+#  if(type %in% c("state","psi","lam","lambda","sigma","dist")){
+#    return(methods::slot(umf, "siteCovs"))
+#  } else if(type %in% c("det","rem","fp","b")){
+#    return(methods::slot(umf, "obsCovs"))
+#  } else if(type %in% c("phi","transition","col","ext","gamma","omega","iota")){
+#    return(methods::slot(umf, "yearlySiteCovs"))
+#  }
+#  return(NULL)
+#}
+
+get_base_newdata <- function(object, type){
+  covs <- get_orig_data(object, type)
+  out <- lapply(covs, function(x){
+    if(is.numeric(x)){
+      return(median(x, na.rm=TRUE))
+    } else if(is.factor(x)){
+      return(factor(levels(x)[1], levels=levels(x)))
+    }  else {
+      stop("Unknown column type")
+    }
+  })
+  as.data.frame(out)
+}
+
+get_cov_seq <- function(covariate, object, type){
+  cov_values <- get_orig_data(object, type)[[covariate]]
+  if(is.numeric(cov_values)){
+    rng <- range(cov_values, na.rm=TRUE)
+    return(seq(rng[1], rng[2], length.out=100))
+  } else if(is.factor(cov_values)){
+    return(factor(levels(cov_values), levels=levels(cov_values)))
+  } else {
+    stop("Unknown covariate type")
+  }
+}
+
+setGeneric("plotEffectsData", function(object, ...) standardGeneric("plotEffectsData"))
+
+setMethod("plotEffectsData", "unmarkedFit",
+  function(object, type, covariate, level=0.95, ...){
+
+  #umf <- umf_to_factor(object@data)
+  nd <- get_base_newdata(object, type)
+  if(! covariate %in% names(nd)){
+    stop("Covariate not in this submodel", call.=FALSE)
+  }
+  values <- get_cov_seq(covariate, object, type)
+  nd <- nd[rep(1, length(values)),,drop=FALSE]
+  nd[[covariate]] <- values
+
+  pr <- predict(object, type=type, newdata=nd, level=level, ...)
+  pr$covariate <- covariate
+  pr$covariateValue <- values
+  pr
+})
+
+setGeneric("plotEffects", function(object, ...) standardGeneric("plotEffects"))
+
+setMethod("plotEffects", "unmarkedFit",
+  function(object, type, covariate, level=0.95, ...){
+
+  # Get data for plot
+  plot_data <- plotEffectsData(object, type, covariate, level, ...)
+
+  # Is covariate a factor?
+  is_factor <- is.factor(plot_data$covariateValue)
+
+  # If not a factor (i.e., numeric)
+  if(!is_factor){
+    # Setup basic plot structure
+    plot(plot_data$covariateValue, plot_data$Predicted, type='l',
+         xlab=covariate, ylab=paste("Predicted", tolower(object[type]@name)),
+         ylim=c(min(plot_data$lower, na.rm=T), max(plot_data$upper, na.rm=T)))
+    #Draw error ribbon
+    polygon(x=c(plot_data$covariateValue, rev(plot_data$covariateValue)),
+            y=c(plot_data$upper, rev(plot_data$lower)),
+            border=NA, col='gray90')
+    #Draw line on top of ribbon
+    lines(plot_data$covariateValue, plot_data$Predicted)
+
+  # If covariate is a factor
+  } else {
+    # Get number of unique factor levels (groups)
+    ngroup <- nrow(plot_data)
+    # Calculate width of error bar tips based on # of factor levels
+    bw <- ngroup/10
+
+    # Add basic plot of points
+    plot(1:ngroup, plot_data$Predicted,
+         xlab=covariate, ylab=paste("Predicted", tolower(object[type]@name)),
+         xlim=c(1-bw/2,ngroup+bw/2),
+         xaxt='n',
+         ylim=c(min(plot_data$lower, na.rm=T), max(plot_data$upper, na.rm=T)))
+
+    # Draw new x-axis with factor labels
+    axis(1, at=1:ngroup, labels=plot_data$covariateValue)
+
+    # Draw error bars
+    for (i in 1:ngroup){
+      segments(i, plot_data$lower[i], i, plot_data$upper[i])
+      segments(i-bw/2, plot_data$lower[i], i+bw/2, plot_data$lower[i])
+      segments(i-bw/2, plot_data$upper[i], i+bw/2, plot_data$upper[i])
+    }
+    #Draw points on top of bars
+    points(1:ngroup, plot_data$Predicted, pch=19)
+  }
+
+})
diff --git a/man/plotEffects.Rd b/man/plotEffects.Rd
new file mode 100644
index 0000000..d176d43
--- /dev/null
+++ b/man/plotEffects.Rd
@@ -0,0 +1,100 @@
+\name{plotEffects}
+\alias{plotEffects}
+\alias{plotEffects-methods}
+\alias{plotEffects,unmarkedFit-method}
+\alias{plotEffectsData}
+\alias{plotEffectsData-methods}
+\alias{plotEffectsData,unmarkedFit-method}
+
+\title{Plot marginal effects of covariates in unmarked models}
+
+\description{This function generates a plot visualizing the effects
+ of a single covariate on a parameter (e.g. occupancy, abundance) in an unmarked 
+ model. If the covariate is numeric, the result is a line plot with an error
+ ribbon where the x-axis is the range of the covariate and the y-axis is the 
+ predicted parameter value. If the covariate is an R factor (i.e., categorical),
+ the x-axis instead contains each unique value of the covariate. 
+
+ All covariates in the model besides the one being plotted are held either at
+ their median value (if they are numeric) or at their reference level (if they
+ are factors).
+
+ Some types of unmarked models may require additional arguments, which are passed
+ to the matching \code{predict} method. For example, \code{unmarkedFitOccuMulti} 
+ models require the \code{species} argument to be included in the function 
+ call in order to work properly.
+
+ If you want to customize a plot, the easiest approach is to get data 
+ formatted for plotting using \code{plotEffectsData}, and use that. If
+ you want to see and/or modify the code used by \code{plotEffects} to generate
+ the default plots, run \code{getMethod("plotEffects", "unmarkedFit")} in
+ the R console.
+}
+
+\usage{
+\S4method{plotEffects}{unmarkedFit}(object, type, covariate, level=0.95, ...)
+\S4method{plotEffectsData}{unmarkedFit}(object, type, covariate, level=0.95, ...)
+}
+
+\arguments{
+  \item{object}{A fitted model inheriting class \code{unmarkedFit}}
+  \item{type}{Submodel in which the covariate of interest can be found, for
+    example \code{"state"} or \code{"det"}. This will depend on the fitted model}
+  \item{covariate}{The name of the covariate to be plotted, as a character string}
+  \item{level}{Confidence level for the error ribbons or bars}
+  \item{...}{Other arguments passed to the \code{predict} function, required
+    for some \code{unmarkedFit} types such as \code{unmarkedFitOccuMulti}}
+} 
+
+\value{A plot (\code{plotEffects} or a data frame (\code{plotEffectsData})
+  containing values to be used in a plot.
+}
+
+\author{Ken Kellner \email{contact@kenkellner.com}}
+
+\examples{
+
+\dontrun{
+
+# Simulate data and build an unmarked frame
+set.seed(123)
+dat_occ <- data.frame(x1=rnorm(500))
+dat_p <- data.frame(x2=rnorm(500*5))
+
+y <- matrix(NA, 500, 5)
+z <- rep(NA, 500)
+
+b <- c(0.4, -0.5, 0.3, 0.5)
+
+re_fac <- factor(sample(letters[1:5], 500, replace=T))
+dat_occ$group <- re_fac
+re <- rnorm(5, 0, 1.2)
+re_idx <- as.numeric(re_fac)
+
+idx <- 1
+for (i in 1:500){
+  z[i] <- rbinom(1,1, plogis(b[1] + b[2]*dat_occ$x1[i] + re[re_idx[i]]))
+  for (j in 1:5){
+    y[i,j] <- z[i]*rbinom(1,1,
+                    plogis(b[3] + b[4]*dat_p$x2[idx]))
+    idx <- idx + 1
+  }
+}
+
+umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p)
+
+# Fit model
+(fm <- occu(~x2 ~x1 + group, umf))
+
+# Plot marginal effects of various covariates
+plotEffects(fm, "state", "x1")
+plotEffects(fm, "state", "group")
+plotEffects(fm, "det", "x2")
+
+# Get raw data used for a plot
+plotEffectsData(fm, "state", "group")
+
+# See code used by plotEffects so you can edit it yourself and customize the plot
+methods::getMethod("plotEffects", "unmarkedFit")
+}
+}
diff --git a/tests/testthat/test_gdistremoval.R b/tests/testthat/test_gdistremoval.R
index 8f453a5..3e205fb 100644
--- a/tests/testthat/test_gdistremoval.R
+++ b/tests/testthat/test_gdistremoval.R
@@ -343,6 +343,12 @@ test_that("gdistremoval can fit models",{
   expect_is(fl, "unmarkedFitList")
   ms <- modSel(fl)
   expect_is(ms, "unmarkedModSel")
+
+  #plotEffects
+  pdf(NULL)
+  plotEffects(fit, 'lambda', "sc1")
+  plotEffects(fit, 'rem', 'oc1')
+  dev.off()
 })
 
 test_that("gdistremoval predict method works",{
diff --git a/tests/testthat/test_occuTTD.R b/tests/testthat/test_occuTTD.R
index f152cde..0fb7d81 100644
--- a/tests/testthat/test_occuTTD.R
+++ b/tests/testthat/test_occuTTD.R
@@ -382,6 +382,12 @@ test_that("occuTTD can fit a dynamic model",{
   expect_equivalent(dim(s[[1]]), c(100,4))
   r <- residuals(fit)
   expect_equivalent(dim(r), c(100,4))
+  pdf(NULL)
+  plotEffects(fit, 'psi', 'elev')
+  plotEffects(fit, 'col', 'forest')
+  plotEffects(fit, 'det', 'wind')
+  plotEffects(fit, 'det', 'obs')
+  dev.off()
 
   #Check ranef
   r <- ranef(fit)
diff --git a/tests/testthat/test_plotEffects.R b/tests/testthat/test_plotEffects.R
new file mode 100644
index 0000000..3aea3ec
--- /dev/null
+++ b/tests/testthat/test_plotEffects.R
@@ -0,0 +1,48 @@
+context("plotEffects")
+
+skip_on_cran()
+
+set.seed(123)
+dat_occ <- data.frame(x1=rnorm(500))
+dat_p <- data.frame(x2=rnorm(500*5))
+
+y <- matrix(NA, 500, 5)
+z <- rep(NA, 500)
+
+b <- c(0.4, -0.5, 0.3, 0.5)
+
+re_fac <- factor(sample(letters[1:5], 500, replace=T))
+dat_occ$group <- re_fac
+re <- rnorm(5, 0, 1.2)
+re_idx <- as.numeric(re_fac)
+
+idx <- 1
+for (i in 1:500){
+  z[i] <- rbinom(1,1, plogis(b[1] + b[2]*dat_occ$x1[i] + re[re_idx[i]]))
+  for (j in 1:5){
+    y[i,j] <- z[i]*rbinom(1,1,
+                    plogis(b[3] + b[4]*dat_p$x2[idx]))
+    idx <- idx + 1
+  }
+}
+
+umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p)
+
+fm <- occu(~x2 ~x1 + group, umf)
+
+test_that("plotEffects works", {
+  
+  pdf(NULL)
+  plotEffects(fm, "state", "x1")
+  plotEffects(fm, "state", "group")
+  plotEffects(fm, "det", "x2")
+  dev.off()
+
+  expect_error(plotEffects(fm, "state", "x2"))
+
+  dat <- plotEffectsData(fm, "state", "group")
+
+  expect_true(inherits(dat, "data.frame"))
+  expect_equal(nrow(dat), 5)
+
+})