Merge remote-tracking branch 'upstream/master' into add_terra_support

14 files changed, 334 insertions, 35 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 39e668c..62da3ef 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/R/unmarkedFitList.R b/R/unmarkedFitList.R
index 2974678..67e94ae 100644
--- a/R/unmarkedFitList.R
+++ b/R/unmarkedFitList.R
@@ -42,32 +42,52 @@ setClass("unmarkedFitList",
 
 
 # constructor of unmarkedFitList objects
-fitList <- function(..., fits) {
-    if(length(list(...)) > 0 & !missing(fits))
-        stop("Do not use both the '...' and 'fits' arguments")
-    if(missing(fits)) {
-        fits <- list(...)
-        isList <- sapply(fits, function(x) is.list(x))
-        if(sum(isList) > 1)
-            stop("Specify models as common-seperated objects, or use fits = 'mylist'")
-        if(isList[1L]) {
-            warning("If supplying a list of fits, use fits = 'mylist'")
-            fits <- fits[[1L]] 	# This is allowed for back-compatability.
-            }
-        if(is.null(names(fits))) {
-            c <- match.call(expand.dots = FALSE)
-            names(fits) <- as.character(c[[2]])
-            warning("Your list was unnamed, so model names were added as object names")
-            }
-        }
+fitList <- function(..., fits, autoNames=c("object","formula")) {
+  autoNames <- match.arg(autoNames)
+  if(length(list(...)) > 0 & !missing(fits))
+    stop("Do not use both the '...' and 'fits' arguments")
+  if(missing(fits)) {
+    fits <- list(...)
+    isList <- sapply(fits, function(x) is.list(x))
+    if(sum(isList) > 1)
+      stop("Specify models as common-seperated objects, or use fits = 'mylist'")
+    if(isList[1L]) {
+      warning("If supplying a list of fits, use fits = 'mylist'")
+      fits <- fits[[1L]] 	# This is allowed for back-compatability.
+    }
     if(is.null(names(fits))) {
-        names(fits) <- as.character(1:(length(fits)))
-        warning("Your list was unnamed, so model names were added as c('1','2',...)")
+      message("Your list was unnamed, so model names were added automatically")
+      if(autoNames=="formula"){
+        if(inherits(fits[[1]], c("unmarkedFitOccuMulti", "unmarkedFitOccuMS"))){
+          warning("simple formula not available, naming based on object instead", call.=FALSE)
+          autoNames <- "object"
+        } else {
+          names(fits) <- lapply(fits, function(x) gsub(" ", "", as.character(deparse(x@formula))))
         }
-    umfl <- new("unmarkedFitList", fits=fits)
-    return(umfl)
+      }
+      if(autoNames=="object"){
+        c <- match.call(expand.dots = FALSE)
+        names(fits) <- as.character(c[[2]])
+      }
     }
-
+  }
+  if(is.null(names(fits))) {
+    message("Your list was unnamed, so model names were added automatically")
+    if(autoNames=="formula"){
+        if(inherits(fits[[1]], c("unmarkedFitOccuMulti", "unmarkedFitOccuMS"))){
+        warning("simple formula not available, naming as numbers instead", call.=FALSE)
+        autoNames <- "object"
+      } else {
+        names(fits) <- lapply(fits, function(x) gsub(" ", "", as.character(deparse(x@formula))))
+      }
+    }
+    if(autoNames=="object"){
+      names(fits) <- as.character(1:(length(fits)))
+    }
+  }
+  umfl <- new("unmarkedFitList", fits=fits)
+  return(umfl)
+}
 
 setMethod("summary", "unmarkedFitList", function(object) {
     fits <- object@fits
diff --git a/man/fitList.Rd b/man/fitList.Rd
index e5047bd..e01e315 100644
--- a/man/fitList.Rd
+++ b/man/fitList.Rd
@@ -1,11 +1,12 @@
 \name{fitList}
 \alias{fitList}
 \title{constructor of unmarkedFitList objects}
-\usage{fitList(..., fits)}
+\usage{fitList(..., fits, autoNames=c("object", "formula"))}
 \description{Organize models for model selection or model-averaged prediction.}
 \arguments{
 \item{...}{Fitted models. Preferrably named.}
 \item{fits}{An alternative way of providing the models. A (preferrably named) list of fitted models.}
+\item{autoNames}{Option to change the names \code{unmarked} assigns to models if you don't name them yourself. If \code{autoNames="object"}, models in the \code{fitList} will be named based on their R object names.  If \code{autoNames="formula"}, the models will instead be named based on their formulas. This is not possible for some model types.}
 }
 \note{Two requirements exist to conduct AIC-based model-selection and model-averaging in unmarked. First, the data objects (ie, unmarkedFrames) must be identical among fitted models. Second, the response matrix must be identical among fitted models after missing values have been removed. This means that if a response value was removed in one model due to missingness, it needs to be removed from all models.
 } 
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_fitList.R b/tests/testthat/test_fitList.R
index c33c03d..444a3a0 100644
--- a/tests/testthat/test_fitList.R
+++ b/tests/testthat/test_fitList.R
@@ -37,5 +37,9 @@ test_that("fitList operations work",{
 
   se <- SE(mt)
   expect_equal(dim(se), c(2,5))
-
+  
+  fl <- expect_message(fitList(fm, fm2, autoNames='formula'))
+  expect_equal(names(fl@fits), c("~o1+o2~x", "~1~x"))
+  fl <- expect_message(fitList(fits=list(fm, fm2), autoNames='formula'))
+  expect_equal(names(fl@fits), c("~o1+o2~x", "~1~x"))
 })
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_modSel.R b/tests/testthat/test_modSel.R
index 14fffc2..b6d1485 100644
--- a/tests/testthat/test_modSel.R
+++ b/tests/testthat/test_modSel.R
@@ -12,11 +12,11 @@ test_that("fitLists can be constructed",{
 
   fits1.1 <- fitList(m1=fm1, m2=fm2)
   expect_equal(names(fits1.1@fits), c("m1","m2"))
-  expect_warning(fits1.2 <- fitList(fm1, fm2))
+  expect_message(fits1.2 <- fitList(fm1, fm2))
   expect_equal(names(fits1.2@fits), c("fm1","fm2"))
   fits2.1 <- fitList(fits = list(m1=fm1, m2=fm2))
   expect_equal(names(fits2.1@fits), c("m1","m2"))
-  expect_warning(fits2.2 <- fitList(fits = list(fm1, fm2)))
+  expect_message(fits2.2 <- fitList(fits = list(fm1, fm2)))
   expect_equal(names(fits2.2@fits), c("1","2"))
 
   expect_equal(fits1.1, fits2.1)
diff --git a/tests/testthat/test_occuFP.R b/tests/testthat/test_occuFP.R
index f1381bb..fc61461 100644
--- a/tests/testthat/test_occuFP.R
+++ b/tests/testthat/test_occuFP.R
@@ -22,7 +22,7 @@ test_that("occuFP model can be fit",{
   m1 <- occuFP(detformula = ~ METH, FPformula = ~1,
                stateformula = ~ habitat, data = umf1)
   expect_equal(names(m1), c("state","det","fp"))
-  expect_warning(fl <- fitList(m1,m1))
+  expect_message(fl <- fitList(m1,m1))
   expect_is(fl,"unmarkedFitList")
   expect_equal(length(fl@fits), 2)
 
diff --git a/tests/testthat/test_occuMS.R b/tests/testthat/test_occuMS.R
index aa0d024..7c25acf 100644
--- a/tests/testthat/test_occuMS.R
+++ b/tests/testthat/test_occuMS.R
@@ -227,7 +227,9 @@ test_that("occuMS can fit the multinomial model",{
   expect_equivalent(r@post[1,,1], c(0,0.5222,0.4778), tol=1e-4)
 
   #Check fitList
-  expect_warning(fl <- fitList(fit_C, fit_C))
+  expect_message(fl <- fitList(fit_C, fit_C))
+  expect_message(expect_warning(fl <- fitList(fit_C, fit_C, autoNames='formula')))
+  expect_message(expect_warning(fl <- fitList(fits=list(fit_C, fit_C), autoNames='formula')))
   expect_is(fl,"unmarkedFitList")
   expect_equivalent(length(fl@fits), 2)
 
diff --git a/tests/testthat/test_occuMulti.R b/tests/testthat/test_occuMulti.R
index fbc264a..d37cbc5 100644
--- a/tests/testthat/test_occuMulti.R
+++ b/tests/testthat/test_occuMulti.R
@@ -42,7 +42,9 @@ test_that("occuMulti can fit simple models",{
   expect_equivalent(det, rep(1,length(detlist)), tolerance= 1e-4)
 
   #Check fitList
-  expect_warning(fl <- fitList(fm, fm))
+  expect_message(fl <- fitList(fm, fm))
+  expect_message(expect_warning(fl <- fitList(fm, fm, autoNames='formula')))
+  expect_message(expect_warning(fl <- fitList(fits=list(fm, fm), autoNames='formula')))
   expect_is(fl,"unmarkedFitList")
   expect_equivalent(length(fl@fits), 2)
 
@@ -335,7 +337,7 @@ test_that("occuMulti predict method works",{
 
   #fitList with maxOrder set
   fm2 <- occuMulti(c("~1","~1"), c("~1","~1"), umf, maxOrder=1)
-  expect_warning(fl2 <- fitList(fm, fm2))
+  expect_message(fl2 <- fitList(fm, fm2))
   expect_is(fl2, "unmarkedFitList")
   ms <- modSel(fl2)
   expect_is(ms, "unmarkedModSel")
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)
+
+})