Fix some bugs and add tests

7 files changed, 55 insertions, 7 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index fe34bfc..972f834 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: unmarked
-Version: 1.2.5.9002
-Date: 2022-06-17
+Version: 1.2.5.9003
+Date: 2022-07-22
 Type: Package
 Title: Models for Data from Unmarked Animals
 Authors@R: c(
diff --git a/R/predict.R b/R/predict.R
index cd14637..68efbb7 100644
--- a/R/predict.R
+++ b/R/predict.R
@@ -132,9 +132,12 @@ setGeneric("get_formula", function(object, type, ...){
 })
 
 setMethod("get_formula", "unmarkedFit", function(object, type, ...){
-  if(type == "state") return(as.formula(paste("~", object@formula[3], sep="")))
-  if(type == "det") return(as.formula(object@formula[[2]]))
-  stop("Invalid type")
+  if(type == "state"){
+    return(as.formula(paste("~", object@formula[3], sep="")))
+  } else if(type == "det"){
+    return(as.formula(object@formula[[2]]))
+  }
+  NULL
 })
 
 # When newdata is data.frame/raster, get original dataset
@@ -574,6 +577,12 @@ setMethod("get_orig_data", "unmarkedFitGDR", function(object, type, ...){
 # bespoke predict method since it has numerious unusual options
 # and requires bootstrapping
 
+# This method is used by simulate but not by predict
+setMethod("get_formula", "unmarkedFitOccuMulti", function(object, type, ...){
+  switch(type, state=object@stateformulas,
+               det=object@detformulas)
+})
+
 setMethod("predict", "unmarkedFitOccuMulti",
      function(object, type, newdata,
               #backTransform = TRUE, na.rm = TRUE,
@@ -748,6 +757,12 @@ setMethod("predict", "unmarkedFitOccuMulti",
 
 # bespoke predict method since it requires bootstrapping
 
+# This method is used by simulate by not by predict
+setMethod("get_formula", "unmarkedFitOccuMS", function(object, type, ...){
+  switch(type, psi=object@psiformulas, phi=object@phiformulas,
+               det=object@detformulas)
+})
+
 setMethod("predict", "unmarkedFitOccuMS",
      function(object, type, newdata,
               #backTransform = TRUE, na.rm = TRUE,
diff --git a/R/simulate.R b/R/simulate.R
index db10552..5df7247 100644
--- a/R/simulate.R
+++ b/R/simulate.R
@@ -109,9 +109,9 @@ generate_random_effects <- function(coefs, fit){
         }
 
         if(!is.factor(lvldata)){
-          stop("Random effect covariates must be specified as factors with guide argument", call.=FALSe)
+          stop("Random effect covariates must be specified as factors with guide argument", call.=FALSE)
         }
-        b <- rnorm(length(levels(lvldata)), 0, old_coefs[signame])
+        b <- stats::rnorm(length(levels(lvldata)), 0, old_coefs[signame])
         names(b) <- rep(paste0("b_",i), length(b))
         new_coefs <- c(new_coefs, b)
         coefs[[i]] <- new_coefs
diff --git a/tests/testthat/test_distsamp.R b/tests/testthat/test_distsamp.R
index 9c42119..9d6fe04 100644
--- a/tests/testthat/test_distsamp.R
+++ b/tests/testthat/test_distsamp.R
@@ -265,6 +265,7 @@ test_that("getP works with distsamp",{
 
 test_that("distsamp works with random effects",{
 
+  set.seed(123)
   data(linetran)
   umf <- unmarkedFrameDS(y=as.matrix(linetran[,1:4]), siteCovs=linetran[,6:7],
                          survey="line", tlength=linetran$Length, unitsIn='m',
@@ -301,4 +302,10 @@ test_that("distsamp works with random effects",{
   pr <- lapply(mods,  function(x) predict(x, "state"))
   expect_true(all(sapply(pr, inherits, "data.frame")))
 
+  # Make sure simulate accounts for random effects
+  s <- simulate(hn, nsim=30)
+  avg <- apply(sapply(s, function(x) apply(x,1,sum)),1, mean)
+  # average first count and predicted abundance should be highly correlated
+  pr <- predict(hn, 'state')
+  expect_true(cor(avg, pr$Predicted) > 0.7)
 })
diff --git a/tests/testthat/test_multinomPois.R b/tests/testthat/test_multinomPois.R
index ee51335..6b4b693 100644
--- a/tests/testthat/test_multinomPois.R
+++ b/tests/testthat/test_multinomPois.R
@@ -237,6 +237,12 @@ test_that("multinomPois can fit models with random effects",{
   expect_equivalent(dim(pr), c(100, 4))
   expect_equivalent(dim(pr2), c(5,4))
 
+  # Make sure simulate accounts for random effects
+  s <- simulate(fm, nsim=30)
+  avg <- apply(sapply(s, function(x) x[,1]),1, mean)
+  # average first count and predicted abundance should be highly correlated
+  expect_true(cor(avg, pr$Predicted) > 0.7)
+
   umf2@y[1,1] <- NA
   umf2@y[2,] <- NA
   umf2@siteCovs$x1[3] <- NA
diff --git a/tests/testthat/test_powerAnalysis.R b/tests/testthat/test_powerAnalysis.R
index fbf75a4..d316a64 100644
--- a/tests/testthat/test_powerAnalysis.R
+++ b/tests/testthat/test_powerAnalysis.R
@@ -61,6 +61,17 @@ test_that("powerAnalysis method works",{
   pl <- unmarkedPowerList(template_model, effect_sizes, design=scenarios, nsim=10)
   expect_is(pl, "unmarkedPowerList")
 
+  # With random effect
+  set.seed(123)
+  rguide <- list(group=factor(levels=letters[1:20]))
+  rform <- list(state=~x+(1|group), det=~1)
+  rcf <- list(state=c(intercept=0, x=0.5, group=0.7), det=c(intercept=0))
+  umfr <- simulate("occu", formulas=rform, design=design, coefs=rcf, guide=rguide)
+  fm <- occu(~1~x+(1|group), umfr)
+  pa5 <- powerAnalysis(fm, rcf, nsim=10)
+  s <- summary(pa5)
+  expect_equal(nrow(s), 3)
+  expect_equal(s$Power[2], 1)
   })
 })
 
diff --git a/tests/testthat/test_simulate.R b/tests/testthat/test_simulate.R
index c1cda3e..04f1686 100644
--- a/tests/testthat/test_simulate.R
+++ b/tests/testthat/test_simulate.R
@@ -34,6 +34,15 @@ test_that("simulate can generate new datasets from scratch",{
   expect_true(is.factor(umf2@siteCovs$landcover))
   expect_equivalent(mean(umf2@siteCovs$elev), 2.01722, tol=1e-5)
 
+  # With random effect
+  set.seed(123)
+  rguide <- list(group=factor(levels=letters[1:20]))
+  rform <- list(state=~(1|group), det=~1)
+  rcf <- list(state=c(intercept=0, group=0.7), det=c(intercept=0))
+  umfr <- simulate("occu", formulas=rform, design=design, coefs=rcf, guide=rguide)
+  fm <- occu(~1~(1|group), umfr)
+  expect_equal(sigma(fm)$sigma, 0.6903913, tol=1e-5)
+
   # pcount
   set.seed(123)
   cf$alpha <- c(alpha=0.5)