coa_tagint-patch

.Rbuildignore

@@ -6,7 +6,9 @@
 ^\.github$
 ^Dockerfile$
 ^_pkgdown\.yml$
+^\.devcontainer$
 ^docs$
 ^pkgdown$
 ^CNAME$
 ^air\.toml$
+^src/StanExports_*

.gitignore

@@ -5,3 +5,7 @@ data-raw
 .DS_Store
 docs
 .vscode
+src/RcppExports.cpp
+src/stanExports_*
+*.o
+*.so

R/COA_Tag_Integrated.R

@@ -2,20 +2,22 @@
 
 #' Fits a test-tag integrated Bayesian Spatial Point Process model to estimate individual centers of activity from acoustic telemetry data using Stan
 #'
-#' @param nind   Number of tagged individuals
-#' @param nrec   Number of receivers
-#' @param ntime  Number of time steps
-#' @param ntest  Number of test tags
-#' @param ntrans Number of expected transmissions per tag per time interval
-#' @param y      Array of detection data, where row = individual, column = time step, and matrix = receiver
-#' @param test   Array of test tag detection data, where row = individual tag, column = time step, and matrix = receiver
-#' @param recX   Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency)
-#' @param recY   Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency)
-#' @param xlim   East-west boundaries of spatial extent (receiver array + buffer)
-#' @param ylim   North-south boundaries of spatial extent (receiver array + buffer)
-#' @param testX  Test tag coordinates in the east-west direction (should be projected and scaled for computational efficiency)
-#' @param testY  Test tag coordinates in the north-south direction (should be projected and scaled for computational efficiency)
-#' @param ndraws to be passed to `generated_quantities`. Changes the number of draws. Default is 10.
+#' @param n_ind   Number of tagged individuals
+#' @param n_rec   Number of receivers
+#' @param n_time  Number of time steps
+#' @param n_test  Number of test tags
+#' @param n_trans Number of expected transmissions per tag per time interval
+#' @param n_trans_test Number of expected transmissions per sentinel tag per time interval
+#' @param det      Array of detection data, where row = individual, column = receiver, and matrix = time step
+#' @param det_test   Array of test tag detection data, where row = individual tag, column = receiver, and matrix = time step
+#' @param rec_x   Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency)
+#' @param rec_y   Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency)
+#' @param x_lim   East-west boundaries of spatial extent (receiver array + buffer)
+#' @param y_lim   North-south boundaries of spatial extent (receiver array + buffer)
+#' @param test_x  Test tag coordinates in the east-west direction (should be projected and scaled for computational efficiency)
+#' @param test_y  Test tag coordinates in the north-south direction (should be projected and scaled for computational efficiency)
+#' @param decay  desired decay function. Currently one of "gaussian" or "logistic". Default is "gaussian".
+#' @param n_draws to be passed to `generated_quantities`. Changes the number of draws. Default is 10.
 #' @param ... Additional arguments passed to `sampling` from `rstan`.
 #' This can include setting `chains`, `iter`, `warmup`, and `control`. Please see
 #' `rstan::sampling()` for more info.
@@ -26,42 +28,46 @@
 #'
 #' @seealso [rstan::sampling()]
 #' @export
+
 COA_TagInt <- function(
-  nind,
-  nrec,
-  ntime,
-  ntest,
-  ntrans,
-  y,
-  test,
-  recX,
-  recY,
-  xlim,
-  ylim,
-  testX,
-  testY,
-  ndraws = NULL,
+  n_ind,
+  n_rec,
+  n_time,
+  n_test,
+  n_trans,
+  n_trans_test,
+  det,
+  det_test,
+  rec_x,
+  rec_y,
+  x_lim,
+  y_lim,
+  test_x,
+  test_y,
+  decay = "gaussian",
+  n_draws = NULL,
   ...
 ) {
   # First move everything into a list
   standata <- list(
-    nind = nind,
-    nrec = nrec,
-    ntime = ntime,
-    ntest = ntest,
-    ntrans = ntrans,
-    y = y,
-    test = test,
-    recX = recX,
-    recY = recY,
-    xlim = xlim,
-    ylim = ylim,
-    testX = testX,
-    testY = testY
+    n_ind = n_ind,
+    n_rec = n_rec,
+    n_time = n_time,
+    n_test = n_test,
+    n_trans = n_trans,
+    n_trans_test = n_trans_test,
+    det = det,
+    det_test = det_test,
+    rec_x = rec_x,
+    rec_y = rec_y,
+    x_lim = x_lim,
+    y_lim = y_lim,
+    test_x = test_x,
+    test_y = test_y
   )
 
   # validate this list prior to sending it to the model
-  exp_len <- expected_lengths(recX = recX, recY = recY, ntest_len = ntest)
+  exp_len <- expected_lengths(rec_x = rec_x, rec_y = rec_y, n_test_len = n_test)
 
   validate_standata(standata, exp_len)
   # set rstan options
@@ -70,53 +76,68 @@ COA_TagInt <- function(
   options(mc.cores = parallel::detectCores())
 
   # fit model
-  fit_model <- rstan::sampling(
-    stanmodels$COA_Tag_Integrated,
-    data = standata,
-    ...
-  )
+  if (decay == "gaussian") {
+    fit_model <- rstan::sampling(
+      stanmodels$COA_Tag_Integrated_gaussian,
+      data = standata,
+      ...
+    )
+  } else if (decay == "logistic") {
+    fit_model <- rstan::sampling(
+      stanmodels$COA_Tag_Integrated_logistic,
+      data = standata,
+      ...
+    )
+  } else {
+    stop("decay parameter must be one of \"gaussian\" or \"logistic\".")
+  }
 
   # Save chains after discarding warmup
-  fit_estimates <- as.data.frame(fit_model) # Note this returns parameters and latent states/derived values
+  fit_estimates <- as.data.frame(fit_model)
+  # Note this returns parameters and latent states/derived values
 
   # Summary statistics and convergence diagnostics
-  fit_summary <- rstan::summary(fit_model, pars = c("p0", "sigma"))$summary
-  #fit_summary <- fit_sum$summary
+  if (decay == "gaussian") {
+    fit_summary <- rstan::summary(fit_model, pars = c("p0", "sigma"))$summary
+  } else if (decay == "logistic") {
+    fit_summary <- rstan::summary(fit_model, pars = c("p0"))$summary
+  }
 
   # How much time did fitting take (in minutes)?
   fit_time <- sum(print(rstan::get_elapsed_time(fit_model))) / 60
-  # # calculate generated quantities
+
+  # calculate generated quantities
   fit_generated_quantities <- generated_quantities(
     model = fit_model,
     standata = standata,
-    ndraws = ndraws
+    n_draws = n_draws
   )
   # transform gq into matrix
   tran_fit_gq <- transform_gq(fit_generated_quantities)
 
   # Extract COA estimates
-  coas <- array(NA, dim = c(ntime, 7, nind))
+  coas <- array(NA, dim = c(n_time, 7, n_ind))
   dimnames(coas)[[2]] <- c(
-    'time',
-    'x',
-    'y',
-    'x_lower',
-    'x_upper',
-    'y_lower',
-    'y_upper'
+    "time",
+    "x",
+    "y",
+    "x_lower",
+    "x_upper",
+    "y_lower",
+    "y_upper"
   )
   ew <- NULL
   ns <- NULL
 
-  for (i in 1:nind) {
-    coas[, 1, i] <- seq(1, ntime, 1)
+  for (i in 1:n_ind) {
+    coas[, 1, i] <- seq(1, n_time, 1)
     ew <- dplyr::select(
       fit_estimates,
-      dplyr::starts_with(paste("sx[", i, ",", sep = ''))
+      dplyr::starts_with(paste("x[", i, ",", sep = ""))
     )
     ns <- dplyr::select(
       fit_estimates,
-      dplyr::starts_with(paste("sy[", i, ",", sep = ''))
+      dplyr::starts_with(paste("y[", i, ",", sep = ""))
     )
     coas[, 2, i] <- apply(ew, 2, stats::median)
     coas[, 3, i] <- apply(ns, 2, stats::median)
@@ -128,15 +149,15 @@ COA_TagInt <- function(
 
   coas <- as.data.frame(coas[,, 1])
   # Extract time-varying detection probability estimates
-  d_probs <- array(NA, dim = c(nrec, ntime))
+  d_probs <- array(NA, dim = c(n_rec, n_time))
   p0_est <- NULL
 
-  for (i in 1:ntime) {
+  for (i in 1:n_time) {
     p0_est <- dplyr::select(
       fit_estimates,
-      dplyr::starts_with(paste("p0[", i, ",", sep = ''))
+      dplyr::starts_with(paste("p0[", i, ",", sep = ""))
     )
-    for (j in 1:nrec) {
+    for (j in 1:n_rec) {
       d_probs[j, i] <- stats::median(p0_est[, j])
     }
   }
@@ -152,13 +173,13 @@ COA_TagInt <- function(
     tran_fit_gq
   )
   names(model_results) <- c(
-    'model',
-    'summary',
-    'time',
-    'coas',
-    'detection_probs',
-    'all_estimates',
-    'generated_quantities'
+    "model",
+    "summary",
+    "time",
+    "coas",
+    "detection_probs",
+    "all_estimates",
+    "generated_quantities"
   )
   return(model_results)
 }