Predict with occuCOP: tests and documentation

3 files changed, 114 insertions, 5 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index d85d28a..762d025 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -12,7 +12,7 @@ Authors@R: c(
     person("Jeff", "Hostetler", role="aut"),
     person("Rebecca", "Hutchinson", role="aut"),
     person("Adam", "Smith", role="aut"),
-    person("Léa", "Pautrel", role="ctb"),
+    person("Lea", "Pautrel", role="aut"),
     person("Marc", "Kery", role="ctb"),
     person("Mike", "Meredith", role="ctb"),
     person("Auriel", "Fournier", role="ctb"),
diff --git a/man/predict-methods.Rd b/man/predict-methods.Rd
index 908b3a3..4a5478f 100644
--- a/man/predict-methods.Rd
+++ b/man/predict-methods.Rd
@@ -19,14 +19,46 @@
 \alias{predict,unmarkedFitList-method}
 \alias{predict,unmarkedRanef-method}
 \title{ Methods for Function predict in Package `unmarked' }
+
+\usage{
+predict(object, type, newdata, 
+        backTransform = TRUE, na.rm = TRUE,
+        appendData = FALSE, level = 0.95, 
+        re.form = NULL, ...)
+}
+
 \description{
 These methods return predicted values from fitted model objects.
 }
+
+\arguments{
+
+  \item{object}{An \linkS4class{unmarkedFit} object or an \linkS4class{unmarkedFitList} object. See methods.}
+  
+  \item{type}{A character string that depends upon the fitted models. You can see the possible types with \code{names(umfit)}, with \code{umfit} your \code{unmarkedFit} object.}
+  
+  \item{newdata}{Facultative: if no newdata, uses the actual data (\code{\link{getData}(umfit)}). When specified, \code{newdata} must be \code{\link{unmarkedFrame}}, \code{\link{data.frame}}, \code{\link{RasterLayer}}, \code{\link{RasterStack}}, or \code{\link{SpatRaster}}. See example.}
+  
+  \item{backTransform}{Logical specifying whether to back-transform parameters (\code{backTransform=TRUE}) or not (\code{backTransform=FALSE}).}
+
+  \item{na.rm}{Logical specifying whether to remove NAs (\code{na.rm=TRUE}) or not (\code{na.rm=FALSE}).}
+
+  \item{appendData}{Logical specifying whether to append the data (covariates) used for the prediction (\code{appendData=TRUE}) or not (\code{appendData=FALSE}).}
+
+  \item{level}{Numeric specifying the confidence interval level.}
+  
+  \item{re.form}{}
+  
+  \item{\dots}{Unused.}
+  
+}
+
+
 \section{Methods}{
 \describe{
 
 \item{\code{signature(object = "unmarkedFit")}}{
-"type" depends upon the fitted models. You can see the possible types with \code{names(fm), with \code{fm} your \code{unmarkedFit} object.}
+"type" depends upon the fitted models. You can see the possible types with \code{names(umfit), with \code{umfit} your \code{unmarkedFit} object.}
 }
 \item{\code{signature(object = "unmarkedFitList")}}{
 "type" depends upon the fitted models. You can see the possible types with \code{names(fl), with \code{fl} your \code{unmarkedFitList} object.}
@@ -73,5 +105,70 @@ predictive distribution.
 
 }}
 
+\examples{
+set.seed(123)
+options(max.print = 100)
+
+# Simulate an unmarkedFrame with count data in 100 sites for 3 sampling occasions
+(umf <- simulate(
+  "COP",
+  formulas = list(psi =  ~ habitat, lambda =  ~ rain),
+  coefs = list(
+    psi = c(intercept = qlogis(.2), habitatForest = .5, habitatGrassland = .8),
+    lambda = c(intercept = log(1), rain = -1)
+  ),
+  design = list(M = 100, J = 3),
+  guide = list(habitat = factor(levels = c("City", "Forest", "Grassland")))
+))
+
+# Fit the model
+(umfit <- occuCOP(umf, psiformula = ~ habitat, lambdaformula =  ~ rain, L1 = TRUE))
+
+# See the possible types for this unmarkedFit object
+names(umfit)
+
+# Back-transform the occupancy probability
+# = predict with the data furnished to the model
+predict(object = umfit, type = "psi", appendData = TRUE)
+
+# Because occupancy only depends on the habitat, 
+# we have the same information in a more concise form by specifying newdata
+# as the 3 habitat types in our covariates data.
+predict(
+  object = umfit,
+  type = "psi",
+  newdata = data.frame("habitat" = c("Forest", "Grassland", "City")),
+  appendData = TRUE
+)
+
+# We can look at the 90% confidence interval (instead of the default 95% CI)
+predict(
+  object = umfit,
+  type = "psi",
+  newdata = data.frame("habitat" = c("Forest", "Grassland", "City")),
+  level = 0.90,
+  appendData = TRUE
+)
+
+# Back-transform the detection rate
+predict(object = umfit, type = "lambda", appendData = TRUE)
+
+## This is not easily readable. We can show the results in a clearer way, by:
+##  - adding the site and observation
+##  - printing only the wind covariate used to get the predicted lambda
+cbind(
+  data.frame(
+    "site" = rep(1:numSites(umf), each = obsNum(umf)),
+    "obs" = rep(1:obsNum(umf), times = numSites(umf)),
+    "rain" = getData(umfit)@obsCovs
+  ),
+  predict(umfit, "lambda", appendData = FALSE)
+)
+
+# We can predict the detection rate with new detection covariates
+predict(umfit, "lambda", 
+        newdata = data.frame("rain" = c(0, 1, 2)),
+        appendData = TRUE)
+}
 
 \keyword{methods}
diff --git a/tests/testthat/test_occuCOP.R b/tests/testthat/test_occuCOP.R
index d19f390..e698709 100644
--- a/tests/testthat/test_occuCOP.R
+++ b/tests/testthat/test_occuCOP.R
@@ -414,7 +414,7 @@ test_that("We can simulate COP data", {
   expect_no_error(simulate(umfit))
 })
 
-test_that("occuCOP can fit models with covariates", {
+test_that("occuCOP can fit and predict models with covariates", {
   # Simulate data with covariates ----
   expect_no_error(umf <- simulate(
     "COP",
@@ -432,7 +432,7 @@ test_that("occuCOP can fit models with covariates", {
     guide = list(habitat = factor(levels = c("A", "B", "C")))
   ))
   
-  # Fit
+  # Fit ----
   expect_no_error(umfit <- occuCOP(
     umf,
     psiformula =  ~ habitat + elev,
@@ -450,4 +450,16 @@ test_that("occuCOP can fit models with covariates", {
     psistarts = c(0),
     lambdastarts = c(0,0)
   ))
-})
\ No newline at end of file
+  
+  # Predict ----
+  expect_no_error(predict(umfit, type = "psi"))
+  expect_no_error(predict(umfit, type = "lambda", appendData = TRUE))
+  expect_no_error(predict(umfit, type = "lambda", level = 0.5))
+  expect_no_error(predict(
+    umfit,
+    type = "psi",
+    newdata = data.frame("habitat" = c("A", "B", "C"), "elev" = c(0, 0, 0)),
+    appendData = TRUE
+  ))
+})
+