Merge pull request #146 from DistanceDevelopment/mcds-dot-exe

`MCDS.exe` support for `Distance`

DESCRIPTION

@@ -11,13 +11,13 @@ Description: A simple way of fitting detection functions to distance sampling
     Horvitz-Thompson-like estimator) if survey area information is provided. See
     Miller et al. (2019) <doi:10.18637/jss.v089.i01> for more information on
     methods and <https://examples.distancesampling.org/> for example analyses.
-Version: 1.0.7
+Version: 1.0.8
 URL: https://github.com/DistanceDevelopment/Distance/
 BugReports: https://github.com/DistanceDevelopment/Distance/issues
 Language: en-GB
 Depends:
     R (>= 3.5.0),
-    mrds (>= 2.2.1)
+    mrds (>= 2.2.9)
 Imports:
     dplyr,
     methods,
@@ -35,4 +35,4 @@ Suggests:
     readxl
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3

NAMESPACE

@@ -48,8 +48,8 @@ importFrom(dplyr,inner_join)
 importFrom(dplyr,left_join)
 importFrom(dplyr,mutate)
 importFrom(dplyr,mutate_if)
+importFrom(dplyr,reframe)
 importFrom(dplyr,select)
-importFrom(dplyr,summarize)
 importFrom(dplyr,summarize_all)
 importFrom(dplyr,ungroup)
 importFrom(dplyr,vars)

NEWS

@@ -1,3 +1,8 @@
+Distance 1.0.8
+----------------------
+
+* Support for using MCDS.exe from Distance for Windows to run analyses. You can now download MCDS.exe using mrds::download_MCDS_dot_exe run analyses using this engine, rather (or in tandem with) the usual R optimizers provided in mrds. ds will pick the best (by likelihood) and return that. See ?ds and ?"mcds-dot-exe" for more details.
+
 Distance 1.0.7
 ----------------------
 

R/dht2.R

@@ -54,7 +54,7 @@
 #' @importFrom stats qt na.omit predict terms var qnorm
 #' @importFrom dplyr group_by group_by_at mutate ungroup select distinct
 #' mutate_if if_else summarize_all "%>%" filter_at inner_join anti_join
-#' bind_rows left_join arrange vars summarize
+#' bind_rows left_join arrange vars
 #' @importFrom mrds DeltaMethod
 #' @section Data:
 #' The data format allows for complex stratification schemes to be set-up. Three

R/ds.R

@@ -100,8 +100,21 @@
 #' used.
 #' @param max_adjustments maximum number of adjustments to try (default 5) only
 #' used when `order=NULL`.
-#' @param er_method encounter rate variance calculation: default = 2 gives the method of Innes et al, using expected counts in the encounter rate. Setting to 1 gives observed counts (which matches Distance for Windows) and 0 uses binomial variance (only useful in the rare situation where study area = surveyed area). See [`dht.se`][mrds::dht.se] for more details.
-#' @param dht_se should uncertainty be calculated when using `dht`? Safe to leave as `TRUE`, used in `bootdht`.
+#' @param er_method encounter rate variance calculation: default = 2 gives the
+#' method of Innes et al, using expected counts in the encounter rate. Setting
+#' to 1 gives observed counts (which matches Distance for Windows) and 0 uses
+#' binomial variance (only useful in the rare situation where study area =
+#' surveyed area). See [`dht.se`][mrds::dht.se] for more details.
+#' @param dht_se should uncertainty be calculated when using `dht`? Safe to
+#' leave as `TRUE`, used in `bootdht`.
+#' @param optimizer By default this is set to 'both'. In this case 
+#'   the R optimizer will be used and if present the MCDS optimizer will also 
+#'   be used. The result with the best likelihood value will be selected. To 
+#'   run only a specified optimizer set this value to either 'R' or 'MCDS'. 
+#'   See [`mcds_dot_exe`][mrds::mcds_dot_exe] for setup instructions.
+#' @param winebin If you are trying to use our MCDS.exe optimizer on a
+#'   non-windows system then you may need to specify the winebin. Please 
+#'   see [`mcds_dot_exe`][mrds::mcds_dot_exe] for more details.
 #' @param dht.group deprecated, see same argument with underscore, above.
 #' @param region.table deprecated, see same argument with underscore, above.
 #' @param sample.table deprecated, see same argument with underscore, above.
@@ -314,6 +327,8 @@ ds <- function(data, truncation=ifelse(is.null(cutpoints),
              convert_units=1, er_var=ifelse(transect=="line", "R2", "P3"),
              method="nlminb", quiet=FALSE, debug_level=0,
              initial_values=NULL, max_adjustments=5, er_method=2, dht_se=TRUE,
+             optimizer = "both",
+             winebin = NULL,
              # deprecated below here:
              dht.group,
              region.table,
@@ -496,7 +511,8 @@ ds <- function(data, truncation=ifelse(is.null(cutpoints),
   }
 
   # set up the control options
-  control <- list(optimx.method=method, showit=debug_level)
+  control <- list(optimx.method=method, showit=debug_level,
+                  optimizer = optimizer, winebin = winebin)
 
   # if initial values were supplied, pass them on
   if(!is.null(initial_values) & !aic.search){
@@ -685,7 +701,7 @@ ds <- function(data, truncation=ifelse(is.null(cutpoints),
   } # end for() over adjustments
 
   if(is.null(model)){
-    stop("No models could be fitted.")
+    stop("No models could be fitted.", call. = FALSE)
   }
 
   # check that hazard models have a reasonable scale parameter
@@ -714,6 +730,10 @@ ds <- function(data, truncation=ifelse(is.null(cutpoints),
         if(any(is.na(model$hessian))){
           message("Some variance-covariance matrix elements were NA, possible numerical problems; only estimating detection function.\n")
           dht.res <- NULL
+          # If hessian is NULL - with the exception of the unif with no adj
+        }else if(is.null(model$hessian) && !(model$ds$aux$ddfobj$type == "unif" && is.null(model$ds$aux$ddfobj$adjustment))){
+          message("No hessian, possible numerical problems; only estimating detection function.\n")
+          dht.res <- NULL
         }else{
           dht.res <- dht(model, region_table, sample_table,
                          options=dht_options, se=dht_se)
@@ -733,6 +753,10 @@ ds <- function(data, truncation=ifelse(is.null(cutpoints),
       if(any(is.na(model$hessian))){
         message("Some variance-covariance matrix elements were NA, possible numerical problems; only estimating detection function.\n")
         dht.res <- NULL
+        # If hessian is NULL - with the exception of the unif with no adj
+      }else if(is.null(model$hessian) && !(model$ds$aux$ddfobj$type == "unif" && is.null(model$ds$aux$ddfobj$adjustment))){
+        message("No hessian, possible numerical problems; only estimating detection function.\n")
+        dht.res <- NULL
       }else{
         dht.res <- dht(model, region_table, sample_table, obs_table,
                        options=dht_options, se=dht_se)

R/dummy_ddf.R

@@ -57,7 +57,7 @@ dummy_ddf <- function(data, width, left=0, transect="line"){
 #' @param object model object
 #' @param newdata how many 1s should we return?
 #' @param compute unused, compatibility with [`mrds::predict`][mrds::predict]
-#' @param int.range unused, compatability with [`mrds::predict`][mrds::predict]
+#' @param int.range unused, compatibility with [`mrds::predict`][mrds::predict]
 #' @param esw should the strip width be returned?
 #' @param \dots for S3 consistency
 #' @author David L Miller

R/make_activity_fn.R

@@ -9,7 +9,6 @@
 #' passed to `bootdht`. It is recommended that you try out the function before
 #' passing it to [`bootdht`]. See examples for a template for use.
 #'
-#' @inheritParams activity::fitact
 #' @param \dots parameters specified by activity::fitact
 #' @param detector_daily_duration by default we assume that detectors were able to detect animals for 24 hours, if they were only able to do this for some proportion of the day (say daylight hours), then adjust this argument accordingly
 #' @return a function which generates a single bootstrap estimate of

R/print.summary.dsmodel.R

@@ -42,14 +42,15 @@ print.summary.dsmodel <- function (x,...){
   cat("Number of observations : ", x$n, "\n")
   cat("Distance range         : ", x$left, " - ", x$width, "\n")
 
-  cat("\nModel :", model.description(model), "\n")
+  cat("\nModel       :", model.description(model), "\n")
   # Remind the user that monotonicity constraints were enforced
   if(x$mono & x$mono.strict){
     cat("\nStrict monotonicity constraints were enforced.\n")
   }else if(x$mono){
     cat("\nMonotonicity constraints were enforced.\n")
   }
-  cat("AIC   :", x$aic, "\n")
+  cat("AIC         : ", x$aic, "\n")
+  cat("Optimisation: ", x$optimise, "\n")
 
   # parameter summaries
   cat("\nDetection function parameters\n")

R/varNhat.R

@@ -1,5 +1,6 @@
 # Calculate the variance contribution of the detection function
 #  to abundance estimates
+#' @importFrom dplyr reframe
 varNhat <- function(data, model){
 
   # format the data
@@ -18,7 +19,7 @@ varNhat <- function(data, model){
   grp_dat <- data %>%
     select("Covered_area", "Area", "Sample.Label", !!strat_vars) %>%
     distinct() %>%
-    summarize(Covered_area = sum(.data$Covered_area),
+    reframe(Covered_area = sum(.data$Covered_area),
               Area         = .data$Area) %>%
     distinct()
 
@@ -38,7 +39,7 @@ varNhat <- function(data, model){
 
     res <- data %>%
       mutate(Nc = (.data$size/.data$p)/.data$rate) %>%
-      summarize(N = (.data$Area/.data$Covered_area) *
+      reframe(N = (.data$Area/.data$Covered_area) *
                      sum(.data$Nc, na.rm=TRUE)) %>%
       distinct()
 

tests/testthat/test_ds.R

@@ -187,7 +187,9 @@ test_that("Percentage truncation works when distances are missing",{
 
   data(minke)
 
-  expect_equal(ds(minke, truncation="15%", adjustment=NULL)$ddf$criterion,
+  expect_equal(ds(minke, truncation="15%", 
+                  adjustment=NULL, 
+                  optimizer = "R")$ddf$criterion,
                -8.1705496, tol=1e-5)
 })
 

tests/testthat/test_summarize.R

@@ -14,9 +14,9 @@ test_that("Error on different truncation distance", {
   skip_on_cran()
   set.seed(100)
   # some models
-  t4 <- ds(egdata, 4)
+  t4 <- ds(egdata, 4, optimizer = "R")
   t42 <- t4
-  t4hr <- suppressWarnings(ds(egdata, 4, key="hr"))
+  t4hr <- suppressWarnings(ds(egdata, 4, key="hr", optimizer = "R"))
   t3 <- suppressWarnings(ds(egdata, 3))
   t14 <- suppressWarnings(ds(egdata, list(left=1, right=4)))
 

tests/testthat/test_wrens.R

@@ -9,7 +9,8 @@ test_that("wren 5 minute counts works",{
   cu_wren_5min <- 1/sqrt(10000)
 
   w1_df_unif <- ds(wren_5min, transect="point", truncation=110, key="unif",
-                   convert_units=cu_wren_5min, order=c(1,2), er_var="P3")
+                   convert_units=cu_wren_5min, order=c(1,2), er_var="P3",
+                   optimizer = "R")
 
 
   # do the same thing with dht2