Finish adding new tests

12 files changed, 208 insertions, 4 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 521adca..c8b43a3 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: unmarked
 Version: 1.1.1.9017
-Date: 2022-03-11
+Date: 2022-03-25
 Type: Package
 Title: Models for Data from Unmarked Animals
 Authors@R: c(
diff --git a/R/occuMulti.R b/R/occuMulti.R
index ad204e0..8b20dd1 100644
--- a/R/occuMulti.R
+++ b/R/occuMulti.R
@@ -58,7 +58,7 @@ occuMulti <- function(detformulas, stateformulas,  data, maxOrder,
       }
     }
     psi <- exp(f %*% t_dmF)
-    psi <- psi/rowSums(psi)
+    psi <- psi/rowSums(as.matrix(psi))
 
     #p
     p <- matrix(NA,nrow=nrow(y),ncol=S)
@@ -82,7 +82,7 @@ occuMulti <- function(detformulas, stateformulas,  data, maxOrder,
     pen <- penalty * 0.5 * sum(params^2)
 
     #neg log likelihood
-    -1 * (sum(log(rowSums(psi*prdProbY))) - pen)
+    -1 * (sum(log(rowSums(as.matrix(psi*prdProbY)))) - pen)
   }
   #----------------------------------------------------------------------------
 
diff --git a/tests/testthat/test_gdistsamp.R b/tests/testthat/test_gdistsamp.R
index 76c5bae..0cbfefe 100644
--- a/tests/testthat/test_gdistsamp.R
+++ b/tests/testthat/test_gdistsamp.R
@@ -88,6 +88,13 @@ test_that("gdistsamp with halfnorm keyfunction works",{
                             dist.breaks=breaks,
                             tlength=rep(transect.length, R), numPrimary=T)
 
+    # R and C give same result
+    fm_R <- gdistsamp(~1, ~1, ~1, umf, output="density", se=FALSE, engine="R",
+                      control=list(maxit=1))
+    fm_C <- gdistsamp(~1, ~1, ~1, umf, output="density", se=FALSE, engine="C",
+                      control=list(maxit=1))
+    expect_equal(coef(fm_R), coef(fm_C))
+
     #When output = density
     #fm_R <- gdistsamp(~1, ~1, ~1, umf, output="density", se=FALSE, engine="R")
     fm_C <- gdistsamp(~1, ~1, ~1, umf, output="density", se=FALSE, engine="C")
@@ -199,6 +206,14 @@ test_that("gdistsamp with halfnorm keyfunction works",{
     # Check error when formula has random effects
     expect_error(gdistsamp(~(1|dummy), ~1, ~1, umf))
 
+    # R and C engines return same result
+    fm_C <- gdistsamp(~elevation, ~1, ~chaparral, jayumf, output='density',
+                      engine="C", se=F, control=list(maxit=1))
+    fm_R <- gdistsamp(~elevation, ~1, ~chaparral, jayumf, output='density',
+                      engine="R", se=F, control=list(maxit=1))
+    expect_equal(coef(fm_C), coef(fm_R))
+
+
 })
 
 test_that("gdistsamp with uniform keyfunction works",{
@@ -239,6 +254,13 @@ test_that("gdistsamp with uniform keyfunction works",{
                             dist.breaks=breaks,
                             tlength=rep(transect.length, R), numPrimary=T)
 
+    # R and C engines return same result
+    fm_R <- gdistsamp(~par1, ~par2, ~1, umf, output="density",
+                      keyfun="uniform", se=FALSE, engine="C", control=list(maxit=1))
+    fm_C <- gdistsamp(~par1, ~par2, ~1, umf, output="density",
+                      keyfun="uniform", se=FALSE, engine="R", control=list(maxit=1))
+    expect_equal(coef(fm_R), coef(fm_C))
+
     #fm_R <- gdistsamp(~par1, ~par2, ~1, umf, output="density",
     #                  keyfun="uniform", se=FALSE, engine="R")
     fm_C <- gdistsamp(~par1, ~par2, ~1, umf, output="density",
@@ -310,6 +332,13 @@ test_that("gdistsamp with exp keyfunction works",{
                             dist.breaks=breaks,
                             tlength=rep(transect.length, R), numPrimary=T)
 
+    # R and C engines return same result
+    fm_C <- gdistsamp(~par1, ~par2, ~par3, umf, output="density", se=FALSE,
+                      keyfun="exp",engine="C", control=list(maxit=1))
+    fm_R <- gdistsamp(~par1, ~par2, ~par3, umf, output="density", se=FALSE,
+                      keyfun="exp",engine="R", control=list(maxit=1))
+    expect_equal(fm_C@AIC, fm_R@AIC)
+
     #fm_R <- gdistsamp(~par1, ~par2, ~par3, umf, output="density", se=FALSE,
     #                  keyfun="exp",engine="R")
     fm_C <- gdistsamp(~par1, ~par2, ~par3, umf, output="density", se=FALSE,
@@ -342,6 +371,13 @@ test_that("gdistsamp with exp keyfunction works",{
 
     expect_equivalent(coef(fm_C),c(-3.1952,-0.1244,3.7986,3.6574),tol=1e-4)
     #expect_equivalent(coef(fm_R),coef(fm_C),tol=1e-4)
+
+    fm_C <- gdistsamp(~elevation, ~1, ~chaparral, jayumf, output='density',
+                      keyfun="exp",engine="C", se=F,control=list(maxit=1))
+    fm_R <- gdistsamp(~elevation, ~1, ~chaparral, jayumf, output='density',
+                      keyfun="exp",engine="R", se=F, control=list(maxit=1))
+    expect_equal(fm_C@AIC, fm_R@AIC, tol=1e-5)
+
 })
 
 test_that("gdistsamp with hazard keyfunction works",{
@@ -383,6 +419,13 @@ test_that("gdistsamp with hazard keyfunction works",{
                             dist.breaks=breaks,
                             tlength=rep(transect.length, R), numPrimary=T)
 
+    # R and C engines give same result
+    fm_C <- gdistsamp(~par1, ~par2, ~par3, umf, output="density", se=FALSE,
+                      keyfun="hazard",engine="C", control=list(maxit=1))
+    fm_R <- gdistsamp(~par1, ~par2, ~par3, umf, output="density", se=FALSE,
+                      keyfun="hazard",engine="R", control=list(maxit=1))
+    expect_equal(fm_C@AIC, fm_R@AIC, tol=1e-5)
+
     #fm_R <- gdistsamp(~par1, ~par2, ~par3, umf, output="density", se=FALSE,
     #                  keyfun="hazard",engine="R")
     fm_C <- gdistsamp(~par1, ~par2, ~par3, umf, output="density", se=FALSE,
@@ -415,6 +458,12 @@ test_that("gdistsamp with hazard keyfunction works",{
 
     expect_equivalent(coef(fm_C),c(0.8584,-0.04336,-0.70738,3.2762,0.1807),tol=1e-4)
     #expect_equivalent(coef(fm_R),coef(fm_C),tol=1e-3)
+
+    fm_C <- gdistsamp(~elevation, ~1, ~chaparral, jayumf, output='density',
+                      keyfun="hazard",engine="C", se=F,control=list(maxit=1))
+    fm_R <- gdistsamp(~elevation, ~1, ~chaparral, jayumf, output='density',
+                      keyfun="hazard",engine="R", se=F, control=list(maxit=1))
+    expect_equal(fm_C@AIC, fm_R@AIC, tol=1e-3)
 })
 
 test_that("predict works for gdistsamp",{
diff --git a/tests/testthat/test_gmultmix.R b/tests/testthat/test_gmultmix.R
index d0f37d7..79c99c7 100644
--- a/tests/testthat/test_gmultmix.R
+++ b/tests/testthat/test_gmultmix.R
@@ -247,3 +247,21 @@ test_that("MRR custom piFun works",{
   expect_equal(dim(umf.cr1@obsToY)[1] , 6)
   expect_equal(dim(umf.cr1@obsToY)[2] , 14)
 })
+
+test_that("R and C++ engines give identical results",{
+  y <- matrix(0:3, 5, 4)
+  siteCovs <- data.frame(x = c(0,2,3,4,1))
+  siteCovs[3,1] <- NA
+  obsCovs <- data.frame(o1 = 1:20, o2 = exp(-5:4)/20)
+  yrSiteCovs <- data.frame(yr=factor(rep(1:2, 5)))
+
+  umf <- unmarkedFrameGMM(y = y, siteCovs = siteCovs, obsCovs = obsCovs,
+        yearlySiteCovs = yrSiteCovs, type="removal", numPrimary=2)
+  expect_warning(fm_R <- gmultmix(~x, ~yr, ~o1 + o2, data = umf, K=23,
+                                  engine="R", control=list(maxit=1)))
+  expect_warning(fm_C <- gmultmix(~x, ~yr, ~o1 + o2, data = umf, K=23,
+                                  engine="C", control=list(maxit=1)))
+  expect_equal(coef(fm_R), coef(fm_C))
+
+
+})
diff --git a/tests/testthat/test_gpcount.R b/tests/testthat/test_gpcount.R
index 41dbefb..3d980a6 100644
--- a/tests/testthat/test_gpcount.R
+++ b/tests/testthat/test_gpcount.R
@@ -96,3 +96,23 @@ test_that("gpcount function works", {
   expect_error(gpcount(~(1|dummy),~1,~1,umf))
 
 })
+
+test_that("gpcount R and C++ engines give same results",{
+
+  y <- matrix(c(0,0,0, 1,0,1, 2,2,2,
+                3,2,3, 2,2,2, 1,1,1,
+                NA,0,0, 0,0,0, 0,0,0,
+                3,3,3, 3,1,3, 2,2,1,
+                0,0,0, 0,0,0, 0,0,0), 5, 9, byrow=TRUE)
+  siteCovs <- data.frame(x = c(0,2,-1,4,-1))
+  obsCovs <- list(o1 = matrix(seq(-3, 3, length=length(y)), 5, 9))
+  yrSiteCovs <- list(yr=matrix(c('1','2','2'), 5, 3, byrow=TRUE))
+
+
+  expect_warning(umf <- unmarkedFrameGPC(y = y, siteCovs = siteCovs, obsCovs = obsCovs,
+        yearlySiteCovs = yrSiteCovs, numPrimary=3))
+
+  fm <- gpcount(~x, ~yr, ~o1, data = umf, K=23, control=list(maxit=1))
+  fmR <- gpcount(~x, ~yr, ~o1, data = umf, K=23, engine="R", control=list(maxit=1))
+  expect_equal(coef(fm), coef(fmR))
+})
diff --git a/tests/testthat/test_nmixTTD.R b/tests/testthat/test_nmixTTD.R
index 473f92b..759d0dc 100644
--- a/tests/testthat/test_nmixTTD.R
+++ b/tests/testthat/test_nmixTTD.R
@@ -174,3 +174,29 @@ test_that("nmixTTD can fit a NB/weib model",{
   sim <- simulate(fit, 2)
   r <- ranef(fit)
 })
+
+test_that("R and C++ engines give identical results",{
+
+  set.seed(123)
+  dens <- exp(log(mu.dens) + beta1 * covDens)
+  N <- rpois(M, dens) # Realized density per site
+  lambda <- exp(log(mu.lambda) + alpha1 * covDet) # per-individual detection rate
+  ttd <- NULL
+  for(i in 1:nrep){
+    expect_warning(ttd <- cbind(ttd,rexp(M, N*lambda[,i])))
+  }
+  ttd[N == 0,] <- 5 # Not observed where N = 0; ttd set to Tmax
+  ttd[ttd >= Tmax] <- 5 # Crop at Tmax
+  umf <- unmarkedFrameOccuTTD(y = ttd, surveyLength=5,
+                            siteCovs = data.frame(covDens=covDens,
+                                                  cdens2=rnorm(length(covDens)),
+                                                  cdens3=rnorm(length(covDens))),
+                            obsCovs = data.frame(covDet=as.vector(t(covDet)),
+                                                 cdet2=rnorm(length(covDet)),
+                                                 cdet3=rnorm(length(covDet))))
+
+  fit <- nmixTTD(~covDens, ~covDet, data=umf, K=max(N)+10, control=list(maxit=2))
+  fitR <- nmixTTD(~covDens, ~covDet, data=umf, K=max(N)+10,
+                 engine="R", control=list(maxit=2))
+  expect_equal(coef(fit), coef(fitR))
+})
diff --git a/tests/testthat/test_occuMulti.R b/tests/testthat/test_occuMulti.R
index 6d37027..cb899cd 100644
--- a/tests/testthat/test_occuMulti.R
+++ b/tests/testthat/test_occuMulti.R
@@ -560,3 +560,29 @@ test_that("Mismatched NAs are identified in unmarkedFrameOccuMulti",{
   expect_true(any(pre_na[1] != pre_na))
   expect_true(!any(post_na[1] != post_na))
 })
+
+test_that("R and C++ engines give same results",{
+
+  y <- list(matrix(rep(0:1,10)[1:10],5,2),
+            matrix(rep(0:1,10)[1:10],5,2))
+
+  set.seed(123)
+  N <- dim(y[[1]])[1]
+  J <- dim(y[[1]])[2]
+  occ_covs <- as.data.frame(matrix(rnorm(N * 3),ncol=3))
+  names(occ_covs) <- paste('occ_cov',1:3,sep='')
+
+  det_covs <- as.data.frame(matrix(rnorm(N*J*2),ncol=2))
+  names(det_covs) <- paste('det_cov',1:2,sep='')
+
+  umf <- unmarkedFrameOccuMulti(y = y, siteCovs = occ_covs, obsCovs = det_covs)
+  stateformulas <- c('~occ_cov1','~occ_cov2','~occ_cov3')
+  detformulas <- c('~det_cov1','~det_cov2')
+
+  fm <- occuMulti(detformulas, stateformulas, data = umf, se=FALSE,
+                  control=list(maxit=1))
+  fmR <- occuMulti(detformulas, stateformulas, data = umf, se=FALSE,
+                  engine="R",control=list(maxit=1))
+  expect_equal(coef(fm), coef(fmR))
+
+})
diff --git a/tests/testthat/test_occuPEN.R b/tests/testthat/test_occuPEN.R
index 6c9a5a9..ba13fde 100644
--- a/tests/testthat/test_occuPEN.R
+++ b/tests/testthat/test_occuPEN.R
@@ -117,6 +117,13 @@ test_that("occuPEN can fit models with covariates",{
   expect_error(fm <- occuPEN(~ 1 ~ 1, data = umf,pen.type="Ridge"))
   expect_error(fm <- occuPEN_CV(~ o1 + o2 ~ x, data = umf,lambda=c(0)))
   expect_error(fm <- occuPEN_CV(~ o1 + o2 ~ x, data = umf,foldAssignments=c(1,2,3,4,5),k=6))
+
+  # nonparboot
+  nbp <- nonparboot(fm, B=2)
+  expect_is(nbp@covMatBS, "matrix")
+  nbp_cv <- nonparboot(fmCV, B=2)
+  expect_is(nbp_cv@covMatBS, "matrix")
+
 })
 
 test_that("occuPEN can handle NAs",{
diff --git a/tests/testthat/test_occuRN.R b/tests/testthat/test_occuRN.R
index a596afa..2b59383 100644
--- a/tests/testthat/test_occuRN.R
+++ b/tests/testthat/test_occuRN.R
@@ -5,6 +5,11 @@ test_that("occuRN can fit models",{
   data(birds)
   woodthrushUMF <- unmarkedFrameOccu(woodthrush.bin)
 
+  # R and C engines give same result
+  fm_R <- occuRN(~ obsNum ~ 1, woodthrushUMF, engine="R", K=5, control=list(maxit=1))
+  fm_C <- occuRN(~ obsNum ~ 1, woodthrushUMF, engine="C", K=5, control=list(maxit=1))
+  expect_equal(fm_R@AIC, fm_C@AIC)
+
   # survey occasion-specific detection probabilities
   fm_C <- occuRN(~ obsNum ~ 1, woodthrushUMF, engine="C", K=10)
   #fm_R <- occuRN(~ obsNum ~ 1, woodthrushUMF, engine="R")
diff --git a/tests/testthat/test_pcount.R b/tests/testthat/test_pcount.R
index 49dc572..59daa1d 100644
--- a/tests/testthat/test_pcount.R
+++ b/tests/testthat/test_pcount.R
@@ -191,3 +191,21 @@ expect_true(all(names(fmnb@estimates@estimates)==c("state","det","alpha")))
 fm <- pcount(~(1|group)~1, umf2, K=19)
 expect_true(sigma(fm)$Model[1]=="p")
 })
+
+test_that("pcount R, C++ and TMB engines give same results",{
+
+  y <- matrix(c(
+      8,7,7,8,
+      6,7,7,5,
+      8,8,7,8,
+      4,5,5,5,
+      4,4,3,3), nrow=5, ncol=4, byrow=TRUE)
+  siteCovs <- data.frame(x = c(0,2,3,4,1))
+  obsCovs <- data.frame(o1 = seq(-1, 1, length=length(y)))
+  umf <- unmarkedFramePCount(y = y, siteCovs = siteCovs, obsCovs = obsCovs)
+  fmC <- pcount(~ o1 ~ x, data = umf, K=30, control=list(maxit=1))
+  fmT <- pcount(~ o1 ~ x, data = umf, K=30, control=list(maxit=1), engine="TMB")
+  fmR <- pcount(~ o1 ~ x, data = umf, K=30, control=list(maxit=1), engine="R")
+  expect_equal(coef(fmC), coef(fmR))
+  expect_equal(coef(fmC), coef(fmT), tol=1e-7)
+})
diff --git a/tests/testthat/test_powerAnalysis.R b/tests/testthat/test_powerAnalysis.R
index eb9ded9..e377aec 100644
--- a/tests/testthat/test_powerAnalysis.R
+++ b/tests/testthat/test_powerAnalysis.R
@@ -19,6 +19,14 @@ test_that("powerAnalysis method works",{
   s <- summary(pa)$Power
   expect_true(s[2]>0.7)
 
+  # output printout
+  out <- capture.output(pa)
+  expect_equal(out[5], "Power Statistics:")
+
+  # update
+  pa_up <- update(pa, alpha=0.5)
+  expect_is(pa_up, "unmarkedPower")
+
   # fewer sites
   set.seed(123)
   pa2 <- powerAnalysis(template_model, effect_sizes, design=list(M=50, J=3), nsim=10)
@@ -39,6 +47,13 @@ test_that("powerAnalysis method works",{
   # list
   pl <- unmarkedPowerList(list(pa, pa2, pa3, pa4))
   expect_is(pl, "unmarkedPowerList")
+  s <- summary(pl)
+  expect_is(s, "data.frame")
+
+  pdf(NULL)
+  pl_plot <- plot(pl)
+  expect_is(pl_plot,"list")
+  dev.off()
 
   # generate list
   scenarios <- expand.grid(M=c(50,100), J=c(2,3))
@@ -78,3 +93,23 @@ test_that("custom datasets can be passed to powerAnalysis",{
   expect_error(powerAnalysis(fit, coefs=coefs, datalist=pco_umf))
   expect_error(powerAnalysis(fit, coefs=coefs, datalist=conv_umf, nsim=20))
 })
+
+test_that("powerAnalysis can be run in parallel",{
+  skip_on_cran()
+  skip_on_ci()
+  forms <- list(state=~elev, det=~1)
+  coefs <- list(state=c(intercept=0, elev=-0.4), det=c(intercept=0))
+  design <- list(M=50, J=3) # 300 sites, 8 occasions per site
+  occu_umf <- simulate("occu", formulas=forms, coefs=coefs, design=design)
+
+  template_model <- occu(~1~elev, occu_umf)
+  nul <- capture.output(expect_error(powerAnalysis(template_model)))
+
+  effect_sizes <- list(state=c(intercept=0, elev=-0.4), det=c(intercept=0))
+  set.seed(123)
+  pa <- powerAnalysis(template_model, coefs=effect_sizes, alpha=0.05, nsim=3,
+                      parallel=TRUE)
+  expect_is(pa, "unmarkedPower")
+
+
+})
diff --git a/vignettes/occuMulti.Rnw b/vignettes/occuMulti.Rnw
index 03ed0ad..21afa3b 100644
--- a/vignettes/occuMulti.Rnw
+++ b/vignettes/occuMulti.Rnw
@@ -105,7 +105,7 @@ names(MesoCarnivores)
 Presence/absence matrices for the three species are in list elements \code{bobcat}, \code{coyote}, and \code{redfox}, and \code{sitecovs} contains the site-level covariate data.
 
 Using this information, we will create an \code{unmarkedFrameOccuMulti} object. 
-You can information by looking at the help file for \code{unmarkedFrameOccuMulti}:
+You can get more information by looking at the help file for \code{unmarkedFrameOccuMulti}:
 
 <<eval=FALSE>>=
 ?unmarkedFrameOccuMulti