Here we demonstrate the use of the alternative optimization engine mcds.exe in the Distance and mrds packages.
This vignette requires the packages Distance version 1.0.8 or later and mrds version 2.2.9 or later. It is designed for use in the Microsoft Windows operating system – the mcds.exe engine currently only has experimental support for MacOS or Linux.
Objectives
- Download the mcds.exe optimization engine
- Demonstrate its use in a simple line transect example (golf tee dataset) via the
Distancepackage - Demonstrate the same example via the
mrdspackage - Demonstrate its use in a point transect example (wren data) where one of the optimizers does not work well (gives a negative estimated detection probability)
- Demonstrate its use to speed up an analysis of camera trap distance sampling data (duiker data) via the
Distancepackage - Discuss when using the alternative optimization engine may be useful.
Introduction
The Distance package is designed to provide a simple way to fit detection functions and estimate abundance using conventional distance sampling methodology (i.e., single observer distance sampling, possibly with covariates, as described by Buckland et al. (2015)). The main function is ds. Underlying Distance is the package mrds – when the function ds is called it does some pre-processing and then calls the function ddf in the mrds package to do the work of detection function fitting. mrds uses maximum likelihood to fit the specified detection function model to the distance data using a built-in algorithm written in R.
An alternative method for analyzing distance sampling data is using the Distance for Windows software (Thomas et al., 2010). This software also uses maximum liklihood to fit the detection function models, and relies on software written in the programming language FORTRAN to do the fitting. The filename of this software is MCDS.exe.
In a perfect world, both methods would produce identical results given the same data and model specification, since the likelihood has only one maximum. However, the likelihood surface is sometimes complex, especially when monotonicity constraints are used (which ensures the estimated detection probability is flat or decreasing with increasing distance when adjustment terms are used) or with “overdispersed” or “spiked” data (see Figure 2 in Thomas et al. (2010)), and so in some (rare) cases one or other piece of software fails to find the maximum. To counteract this, it is possible to run both the R-based optimizer and MCDS.exe from the ds function within the Distance package or the ddf function within mrds package.
Another motivation for using the MCDS.exe software from within R is that the R-based optimizer is sometimes slow to converge and so using MCDS.exe in place of the R-based optimizer can save significant time, particularly when doing a nonparametric bootstrap for large datasets.
This vignette demonstrates how to download and then use the MCDS.exe sofware from within the Distance and mrds packages. For more information, see the MCDS.exe help page within the mrds package.
Downloading and verifying MCDS.exe
The program MCDS.exe does not come automatically with the Distance or mrds packages, to avoid violating CRAN rules, so you must first download it from the distance sampling website. You can check whether MCDS.exe is installed already or not by loading the Distance library:
Loading required package: mrdsThis is mrds 2.2.9
Built: R 4.2.3; ; 2023-07-20 01:36:15 UTC; windows
MCDS.exe not detected, single observer analyses will only be run using optimiser in mrds R library. See ?MCDS for details.
Attaching package: 'Distance'The following object is masked from 'package:mrds':
create.binsIf MCDS.exe is not installed, then you will receive the message
MCDS.exe not detected, single observer analyses will only be run using optimiser in mrds R library. See ?MCDS for details.
In this case, you need to download it from the Distancesampling.org web site:
download.file("http://distancesampling.org/R/MCDS.exe", paste0(system.file(package="mrds"),"/MCDS.exe"), mode = "wb")Now if you reload the Distance package, the MCDS.exe not detected message should not be shown:
Attaching package: 'Distance'The following object is masked from 'package:mrds':
create.binsNow that this software is available, both it and the R optimizer will be used by default for each analysis; you can also choose to use just one or the other, as shown below.
Example with Golf Tee data
Both MCDS.exe and the R-based optimizer
This example (of golf tee data, using only observer 1) is taken from the R help for the ds function: (There is a warning about cluster sizes being coded as -1 that can be ignored.)
#Load data
data(book.tee.data)
tee.data <- subset(book.tee.data$book.tee.dataframe, observer==1)
#Fit detection function - default is half-normal with cosine adjustments
ds.model <- ds(tee.data, truncation = 4)Starting AIC adjustment term selection.Fitting half-normal key functionAIC= 311.138Fitting half-normal key function with cosine(2) adjustmentsAIC= 313.124
Half-normal key function selected.No survey area information supplied, only estimating detection function.summary(ds.model)
Summary for distance analysis
Number of observations : 124
Distance range : 0 - 4
Model : Half-normal key function
AIC : 311.1385
Optimisation: mrds (nlminb)
Detection function parameters
Scale coefficient(s):
estimate se
(Intercept) 0.6632435 0.09981249
Estimate SE CV
Average p 0.5842744 0.04637627 0.07937412
N in covered region 212.2290462 20.85130344 0.09824906Assuming you have MCDS.exe installed, the default is that both it and the R-based optimizer are run. Both give the same result in this example, and when this happens the result from the R-based optimizer is used. You can see this from the line of summary output:
Optimisation: mrds (nlminb)
where mrds is the R package that the Distance package relies on, and nlminb is the R-based optimizer.
You can see the process of both optimizers being used by setting the debug_level argument of the ds function to a value larger than the default of 0 and then examining the output:
ds.model <- ds(tee.data, truncation = 4, debug_level = 1)Starting AIC adjustment term selection.Fitting half-normal key functionDEBUG: initial values = -0.1031529 Running MCDS.exe...Command file written to C:\Users\lt5\AppData\Local\Temp\RtmpCSltvm\cmdtmp53e8649a53d1.txtStats file written to C:\Users\lt5\AppData\Local\Temp\RtmpCSltvm\stat53e832ca2d2c.txtDEBUG: initial values = 0.6632378
DEBUG: Convergence!
Iteration 0.0
Converge = 0
nll = 154.5692
parameters = 0.6632378 MCDS.exe log likehood: -154.5697MCDS.exe pars: 1.941067mrds refitted log likehood: -154.5692276mrds refitted pars: 0.6632378
DEBUG: Convergence!
Iteration 0.0
Converge = 0
nll = 154.5692
parameters = 0.6632435 AIC= 311.138Fitting half-normal key function with cosine(2) adjustmentsDEBUG: initial values = -0.1031529 0 Running MCDS.exe...Command file written to C:\Users\lt5\AppData\Local\Temp\RtmpCSltvm\cmdtmp53e85790220e.txtStats file written to C:\Users\lt5\AppData\Local\Temp\RtmpCSltvm\stat53e85b436126.txtDEBUG: initial values = 0.6606793 -0.0159333
DEBUG: Convergence!
Iteration 0.0
Converge = 0
nll = 154.5619
parameters = 0.6606793, -0.0159333 MCDS.exe log likehood: -154.5624MCDS.exe pars: 1.936107, -0.0159333mrds refitted log likehood: -154.5619307mrds refitted pars: 0.6606793, -0.0159333
Iter: 1 fn: 154.5619 Pars: 0.66068 -0.01591
Iter: 2 fn: 154.5619 Pars: 0.66069 -0.01592
solnp--> Completed in 2 iterations
DEBUG: Convergence!
Iteration 0.0
Converge = 0
nll = 154.5619
parameters = 0.6606853, -0.0159233 AIC= 313.124
Half-normal key function selected.No survey area information supplied, only estimating detection function.First the half-normal with no adjustments is run; for this model the MCDS.exe sofware is run first, followed by the R-based (mrds) optimizer. Both converge and both give the same nll (negative log-likelihood) or 154.5692, giving an AIC of 311.138. The model with half-normal and a cosine adjustment of order 2 is then fitted to the data, with first the MCDS.exe optimizer and then the R-based optimizer. Again both give the same result of nll 154.5619 and an AIC of 313.124. This is higher than the AIC with no adjustments so half-normal with no adjustments is chosen.
In this case, both optimizers produced the same result, so there is no benefit to run MCDS.exe.
Specifying which optimzier to run
As we said earlier, the default behaviour when MCDS.exe has been downloaded is to run both MCDS.exe and the R-based optimizer. However, the optimizer argument can be used to specify which to use – either both, R or MCDS. Here is an example with just the MCDS.exe optimizer:
ds.model <- ds(tee.data, truncation = 4, optimizer = "MCDS")Starting AIC adjustment term selection.Fitting half-normal key functionAIC= 311.138Fitting half-normal key function with cosine(2) adjustmentsAIC= 313.124
Half-normal key function selected.No survey area information supplied, only estimating detection function.Demonstration using ddf in mrds package
Here we demonstrate using both optimizers in the ddf function, rather than via ds.
#Half normal detection function
ddf.model <- ddf(dsmodel = ~mcds(key = "hn", formula = ~1), data = tee.data, method = "ds",
meta.data = list(width = 4))
#Half normal with cos(2) adjustment
ddf.model.cos2 <- ddf(dsmodel = ~mcds(key = "hn", adj.series = "cos", adj.order = 2, formula = ~1),
data = tee.data, method = "ds", meta.data = list(width = 4))
#Compare with AIC
AIC(ddf.model, ddf.model.cos2) df AIC
ddf.model 1 311.1385
ddf.model.cos2 2 313.1239#Model with no adjustment term has lower AIC; show summary of this model
summary(ddf.model)
Summary for ds object
Number of observations : 124
Distance range : 0 - 4
AIC : 311.1385
Optimisation : mrds (nlminb)
Detection function:
Half-normal key function
Detection function parameters
Scale coefficient(s):
estimate se
(Intercept) 0.6632435 0.09981249
Estimate SE CV
Average p 0.5842744 0.04637627 0.07937412
N in covered region 212.2290462 20.85130344 0.09824906As an exercise, fit using just the MCDS.exe optimizer:
ddf.model <- ddf(dsmodel = ~mcds(key = "hn", adj.series = "cos", adj.order = 2,
formula = ~1), data = tee.data, method = "ds",
meta.data = list(width = 4),
control = list(optimizer = "MCDS"))
summary(ddf.model)
Summary for ds object
Number of observations : 124
Distance range : 0 - 4
AIC : 313.1239
Optimisation : MCDS.exe
Detection function:
Half-normal key function with cosine adjustment term of order 2
Detection function parameters
Scale coefficient(s):
estimate se
(Intercept) 0.6606782 0.1043327
Adjustment term coefficient(s):
estimate se
cos, order 2 -0.01593274 0.1351281
Estimate SE CV
Average p 0.5925856 0.08165144 0.1377884
N in covered region 209.2524623 31.22790760 0.1492356Point transect example - wren data
This is an example of point transect data for a bird (wren), from Buckland (2006). In this case one of the optimizers fails correctly to constrain the detection function so the probability of detection is more than zero at all distances, and so we use the other optimizer for inference.
We load the wren 5 minute example dataset and define cutpoints for the distances (they were collected in intervals).
The following call to ds gives several warnings. Some warnings are about the detection function being less than zero at some distances. There is also a warning about the Hessian (which is used for variance estimation), but this relates to the Hermite(4, 6) model (i.e., two Hermite adjustment terms of order 4 and 6) which is not chosen using AIC and so this warning can be ignored.
wren5min.hn.herm.t100 <- ds(data=wren_5min, key="hn", adjustment="herm",
transect="point", cutpoints=bin.cutpoints.100m)Warning in create_bins(data, cutpoints): Some distances were outside
bins and have been removed.Starting AIC adjustment term selection.Fitting half-normal key functionAIC= 427.471Fitting half-normal key function with Hermite(4) adjustmentsWarning in check.mono(result, n.pts = control$mono.points): Detection
function is less than 0 at some distancesWarning in check.mono(result, n.pts = control$mono.points): Detection
function is less than 0 at some distancesAIC= 422.228Fitting half-normal key function with Hermite(4,6) adjustmentsWarning: First partial hessian is singular and second-partial hessian is NULL, no hessian
Warning: Detection function is less than 0 at some distances
Warning: Detection function is less than 0 at some distancesAIC= 423.255
Half-normal key function with Hermite(4) adjustments selected.Warning in mrds::check.mono(model, n.pts = 20): Detection function is
less than 0 at some distancessummary(wren5min.hn.herm.t100)
Summary for distance analysis
Number of observations : 132
Distance range : 0 - 100
Model : Half-normal key function with Hermite polynomial adjustment term of order 4
Strict monotonicity constraints were enforced.
AIC : 422.2284
Optimisation: MCDS.exe
Detection function parameters
Scale coefficient(s):
estimate se
(Intercept) 12.08697 1e+05
Adjustment term coefficient(s):
estimate se
herm, order 4 0.5723854 0.07889437
Estimate SE CV
Average p 0.4399177 0.0253497 0.05762374
N in covered region 300.0561563 26.0954740 0.08696863
Summary statistics:
Region Area CoveredArea Effort n k ER se.ER cv.ER
1 Montrave 33.2 2010619 64 132 32 2.0625 0.1901692 0.09220324
Abundance:
Label Estimate se cv lcl ucl
1 Total 0.004954625 0.00053871 0.1087287 0.003988055 0.006155458
df
1 57.84101
Density:
Label Estimate se cv lcl ucl
1 Total 0.0001492357 1.62262e-05 0.1087287 0.0001201222 0.0001854054
df
1 57.84101The MCDS.exe optimizer is the chosen one (see the `Optimisation’ line of output).
The warnings persist if only the MCDS.exe optimizer is used:
wren5min.hn.herm.t100.mcds <- ds(data=wren_5min, key="hn", adjustment="herm",
transect="point", cutpoints=bin.cutpoints.100m,
optimizer = "MCDS")Warning in create_bins(data, cutpoints): Some distances were outside
bins and have been removed.Starting AIC adjustment term selection.Fitting half-normal key functionAIC= 427.471Fitting half-normal key function with Hermite(4) adjustmentsWarning in check.mono(result, n.pts = control$mono.points): Detection
function is less than 0 at some distancesWarning in check.mono(result, n.pts = control$mono.points): Detection
function is less than 0 at some distancesAIC= 422.228Fitting half-normal key function with Hermite(4,6) adjustmentsWarning: First partial hessian is singular and second-partial hessian is NULL, no hessian
Warning: Detection function is less than 0 at some distances
Warning: Detection function is less than 0 at some distancesAIC= 423.255
Half-normal key function with Hermite(4) adjustments selected.Warning in mrds::check.mono(model, n.pts = 20): Detection function is
less than 0 at some distancesLooking at a plot of the fitted object, it seems that the evaluated pdf is less than 0 at distances close to the truncation point (approx. 95m and greater):
plot(wren5min.hn.herm.t100.mcds, pdf = TRUE)
What appears to be happening here is a failure of the optimization routine to appropriately constrain the model parameters so that the detection function is valid. This happens on occasion (the routines aren’t perfect!) and where it does we recommend trying the other optimization routine. Here we use the R-based optimizer:
wren5min.hn.herm.t100.r <- ds(data=wren_5min, key="hn", adjustment="herm",
transect="point", cutpoints=bin.cutpoints.100m,
optimizer = "R")Warning in create_bins(data, cutpoints): Some distances were outside
bins and have been removed.Starting AIC adjustment term selection.Fitting half-normal key functionAIC= 427.471Fitting half-normal key function with Hermite(4) adjustmentsAIC= 422.743Fitting half-normal key function with Hermite(4,6) adjustmentsAIC= 424.52
Half-normal key function with Hermite(4) adjustments selected.Here the fitted AIC for the chosen model (half normal with one Hermite adjustment of order 4) is 422.74, higher than that with the MCDS.exe optimizer (which was 422.23), which explains why the MCDS.exe optimizer fit was chosen when we allowed ds to choose freely. However, the detection function fit from MCDS.exe was invalid, because it went lower than 0 at about 95m, while the fit with the R-based optimizer looks valid:
plot(wren5min.hn.herm.t100.r, pdf = TRUE)
Hence in this case, we would use the R-based optimizer’s fit.
Camera trap example
For this example, it helps if you are familiar with the Analysis of camera trapping data vignette on the distanceexamples web site.
We first read in the Duiker data.
#Read in data and set up data for analysis
DuikerCameraTraps <- read.csv(file="https://datadryad.org/stash/downloads/file_stream/73221",
header=TRUE, sep="\t")
DuikerCameraTraps$Area <- DuikerCameraTraps$Area / (1000*1000)
DuikerCameraTraps$object <- NA
DuikerCameraTraps$object[!is.na(DuikerCameraTraps$distance)] <- 1:sum(!is.na(DuikerCameraTraps$distance))
#Specify breakpoints and truncation
trunc.list <- list(left=2, right=15)
mybreaks <- c(seq(2,8,1), 10, 12, 15)Then we fit the detection function selected in the camera trap vignette, uniform plus 3 cosine adjustment terms, and time how long the fitting takes:
start.time <- Sys.time()
uni3.r <- ds(DuikerCameraTraps, transect = "point", key="unif", adjustment = "cos",
nadj=3, cutpoints = mybreaks, truncation = trunc.list, optimizer = "R")Fitting uniform key function with cosine(1,2,3) adjustmentsAIC= 44012.394
Summary for distance analysis
Number of observations : 10284
Distance range : 2 - 15
Model : Uniform key function with cosine adjustment terms of order 1,2,3
Strict monotonicity constraints were enforced.
AIC : 44012.39
Optimisation: mrds (nlminb)
Detection function parameters
Scale coefficient(s):
NULL
Adjustment term coefficient(s):
estimate se
cos, order 1 0.93541299 0.01504249
cos, order 2 -0.05304180 0.02437441
cos, order 3 -0.08043313 0.01557445
Estimate SE CV
Average p 3.288267e-01 1.348162e-02 0.04099917
N in covered region 3.127483e+04 1.306897e+03 0.04178751
Summary statistics:
Region Area CoveredArea Effort n k ER
1 Tai 40.37 21858518573 31483179 10284 21 0.0003266506
se.ER cv.ER
1 8.970466e-05 0.2746196
Abundance:
Label Estimate se cv lcl ucl
1 Total 5.776078e-05 1.603804e-05 0.2776632 3.276637e-05 0.0001018211
df
1 20.90147
Density:
Label Estimate se cv lcl ucl
1 Total 1.430785e-06 3.972762e-07 0.2776632 8.116514e-07 2.522197e-06
df
1 20.90147Fitting takes quite a while! - 2 mins. Here we try the MCDS.exe optimizer:
start.time <- Sys.time()
uni3.mcds <- ds(DuikerCameraTraps, transect = "point", key="unif", adjustment = "cos",
nadj=3, cutpoints = mybreaks, truncation = trunc.list, optimizer = "MCDS")Fitting uniform key function with cosine(1,2,3) adjustmentsAIC= 44012.211
Summary for distance analysis
Number of observations : 10284
Distance range : 2 - 15
Model : Uniform key function with cosine adjustment terms of order 1,2,3
Strict monotonicity constraints were enforced.
AIC : 44012.21
Optimisation: MCDS.exe
Detection function parameters
Scale coefficient(s):
NULL
Adjustment term coefficient(s):
estimate se
cos, order 1 0.93518220 0.01504583
cos, order 2 -0.05345965 0.02438049
cos, order 3 -0.08073799 0.01557817
Estimate SE CV
Average p 3.290679e-01 1.349917e-02 0.04102246
N in covered region 3.125191e+04 1.306645e+03 0.04181008
Summary statistics:
Region Area CoveredArea Effort n k ER
1 Tai 40.37 21858518573 31483179 10284 21 0.0003266506
se.ER cv.ER
1 8.970466e-05 0.2746196
Abundance:
Label Estimate se cv lcl ucl
1 Total 5.771844e-05 1.602649e-05 0.2776666 3.274219e-05 0.000101747
df
1 20.9025
Density:
Label Estimate se cv lcl ucl
1 Total 1.429736e-06 3.9699e-07 0.2776666 8.110524e-07 2.520361e-06
df
1 20.9025This took only 8 secs. Hence, for some datasets, it may be quicker to use the MCDS.exe optimizer. This makes a particularly big difference if using the nonparametric bootstrap to estimate variance.
Discussion
We have shown how to fit distance sampling detection functions (for single platform data) using either the R-based optimizer built into the ddf function (via calling ddf or, more likely, calling the ds function in the Distance package) or the MCDS.exe analysis engine used by Distance for Windows. In the vast majority of cases both fitting methods give the same result, and so there is no need to use both. However, the only downside is that fitting takes longer, as each is called in turn. If you have downloaded the MCDS.exe file and want to speed things up, you can use just the R-based optimizer by specifying optimizer = "R" in the call to ds or ddf, or just the MCDS.exe optimizer with optimizer = "MCDS".
Some situations where the two may produce different results are given below.
Detection functions that are close to non-monotonic or close to zero at some distances. When adjustment terms are used in the detection function, then constraints are required to prevent the fitted function from having “bumps” where detection probability increases with increasing distance and also to prevent detection probability from becoming less than zero. The former are called monotonicity constraints and are set using the
monotonicityargument indsor in themeta.dataargument inddf; monotonicity is set on by default. In practice, monotonicity and values less than zero are monitored at a finite set of distances between the 0 and the right truncation point, and (for historical reasons) this set of distances is different for theR-based andMCDS.exeoptimizers. This typically makes no difference to the optimization, but particularly in borderline cases it can result in different fitted functions. Plotting the fitted functions (as we did in the wren example above) can reveal when there is an issue with a fitted function, and if this occurs the associated optimizer should not be used. In the future we plan to bring the two into line so they use the same distances for checking.Detection functions with many adjustment terms. The two optimizers use different algorithms for optimization: the
R-based optimizer uses a routine callednlminbwhileMCDS.exeuses a nonlinear constrained optimizer routine produced by the IMSL group. In cases where there are multiple adjustment terms, and hence several parameters to estimate (that are often correlated) the likelihood maximization is harder, and one or other routine can sometimes fail to find the maximum. In this case, choosing the routine with the higher likelihood (i.e., lower negative log-likelihod, or equivalently lower AIC) is the right thing to do, and this is the default behaviour of the software.Detection functions that are “overdispersed” or with a “spike” in the detection function close to zero distance. Similarly to the above, the detection function can then be hard to maximize and hence on or other optimizer can fail to find the maximum. Solution is as above. Overdispersed data is common in camera trap distance sampling because many detections can be generated by the same individual crossing in front of the camera.
If you are interested in seeing more comparisons of the optimizers on various datasets, we maintain a test suite of both straightforward and challenging datasets together with test code to run and compare the two optimizers – this is available at the MCDS_mrds_compare repository.
If you encounter difficulties when using both optimizers, one possible troubleshooting step is to run the analysis first choosing one optimizer (e.g., specifing the argument optimizer = "MCDS") and then choosing the other (optimizer = "R"). This allows you clearly to see what the output of each optimizer is (including any error messages) and facilitates their comparison.
One other criterion to favour one optimizer over the other is speed, and we have found that for large datasets the MCDS.exe optimizer is quicker.
One thing to note is that the MCDS.exe file will get deleted each time you update the mrds package, so you’ll need to re-download the file if you want to continue using the MCDS.exe optimizer. As shown above, this only requires running one line of code.
Going forward, we plan to work on making further improvements to the R-based optimizer and it is possible that at some point in the future we are confident this optimizer is uniformly better (in terms of better fit and speed) than the MCDS.exe one. If this happens, we will update this vignette and also make announcements on the distance sampling list.