Detection function fitting finding less optimal model also see mrds issue 108 https://github.com/DistanceDevelopment/mrds/issues/108

[Globicephala]

Good night, today when giving a class on DSMs, we found an issue with the ds function of the newest version of the Distance package. When estimating detection probabilities with a half-normal key function, the average probability didn't change regardless of the value we inserted for the max_adjustments argument in the ds function. This is something that is not happening when using the uniform or hazard-rate keys.

When using a version of Distance that I had installed in August, for the data we were working with, average p was about 0.50 without any adjustment, and 0.34 when using a maximum of two adjustments. However, with the newest version of the Distance package, the average p was always 0.50, despite any change on maximum adjustments.

This was an issue that was common to all my students, and once I installed the most recent version of Distance, it started happening to me. Looking forward for this issue to be solved.

Attached are screenshots showcasing this issue. First, with the previous version of Distance:

Now the most recent one:

[erex]

@Globicephala Thanks for raising the issue. Without your data and more complete information about your August analysis, I can't give a definitive answer to your question.

One distinction I note between what you provide from August and the most recent output is the appearance that the most recent analysis was conducted with distances binned; however it is not clear if the distances were binned or exact in your August analysis. There is no message "distbegin" and "distend" in data from your August result. Could it be that the difference between August and now is in the data rather than in the code?

In your first current result, where you set max_adjustments=2, there are two models fit to your data; one with no adjustments (returning an AIC of 1998) and a second with a single adjustment (returning an AIC of 2077). In that event, ds() will return the result for the model with no adjustments. This is the same result you got when specifying max_adjustments=0. Hence the later two analyses you present are consistent because the model with the smallest AIC is the model without adjustments; this situation is unaffected by the value you provide to max_adjustments.

Tell us about the potential that the August analysis was conducted on exact distance data as that might explain the situation.

[Globicephala]

For both cases the analysis was done using binned distances, the data is the exact same. Here's a screenshot of it.

Here are plots of the detection functions. Using the previous version of the package

And the new one:

To show you what happens if I use exact distance to fit the detection function:

I also double checked with the students if the data was the same, and it was. While I was using a version of the package I downloaded in August, I did the code last week.

Thank you for the explanation on the AICs from the models, I completely missed it.

Today I've tried with a dataset for another animal species and cosine adjustments were added. So I wonder what could be different in the way that the detection functions are fit between the two packages.

[erex]

I confess I don't have the answer to your question. Grasping at straws, I note obs_data appears to be a tibble, rather than a data frame. I would not expect unusual behaviour from a tibble, but I have not sent a tibble to ds. On the off chance this has an effect, could you convert your tibble to a data frame before passing it to ds. If that has no effect, I have no answer. We did alter the fitting algorithm for models with key functions and adjustment terms in the latest release (see news).

[Globicephala]

Thank you for the suggestion. Changing the data from a tibble to a data fram had no impact:

It's probably due to the new way that you're fitting the models with key functions, then. Nevertheless, this is a huge difference in estimated detectability, which impacts dramatically the abundance estimates for this species. Thank you for sending the news link.

I've now tested for other data that I had worked with and the detection functions, regardless of the key, give different average p's, and for the same data and cosine adjustments, the AIC scores are different. I'm discussing this with my advisor, Dr. Marques from St. Andrews, and if he gives me the good to go sign, which I think he will, I will paste the data here.

Thank you for all the help.

[erex]

Ahhh, you're working with Tiago, that's good to know. I'm wondering if the fitting algorithm struggles with binned data, with only four bins.

Not that this is related to the results you are obtaining, but why are you binning the data, rather than using exact distances?

[Globicephala]

Yeah I'm working with him, having quite some luck on the advisors department.

I'm binning the data because, in reality, I don't have exact distances. When providing you the detection function with exact distances, the distances I created were the mean between distbegin and distend. The truncation distance considered for data collection was 300m, any sighting beyond that distance was considered as beyond 300m, which I don't have a clue what that means. This is a dataset that was collected opportunistically and we are trying to see what is the best we can do with it.

I hope this answers your question.

[erex]

@Globicephala You may have seen the comment on the related issue I raised under mrds (it is actually mrds that performs the detection function model fitting).

There are two arguments in the ds function you could investigate. See if setting mono_method="solnp" to use the optimiser that was used in previous versions of mrds.

You could also weaken the enforcement of monotonicity of the fitted model with monotonicity="weak"

Let us know what results from either of those two experiments.

[Globicephala]

@erex Thank you for the help. I've been following the other issue that you raised. Here are the results.

mono_method:

With this change in the argument, the results approach what was done with the previous versions of the packages.

Monotonicity:

As you can see, if you compare with the previous functions I've pasted, it gives a different result from either.

Regarding sharing data, Tiago told me to wait a bit untill he returns from the US, so he can also take a look into this properly. Now, considering all these results, what would you recommend me to o for now? Obviously this will most likely be easier if/when you have the dataset, but I'm all ears.

[erex]

Thanks for taking the time to work through the experiments. I think progress has been made. Sifting through all the permutation of model fitting you've done, here is where I think we are:

Your model from last August produced:
preference for 2 adjustment terms, AIC 1844, with $\hat{P}$=0.341

Your model today, setting mono_method="solnp" produced:
preference for 2 adjustment terms, AIC 1936, with $\hat{P}$=0.333

That is reassuring, in that using the optimiser used in the previous version of the package produces results biologically indistinguishable to the detection function model results from last August. That should give you a recipie to go forward for the moment.

After Tiago returns, we can revisit the subject.

[Globicephala]

@erex Tiago has returned and here is the dataset. If you can figure out more stuff about this issue please let me know.
Especially which version of mrds makes sense to use.

obs_data.csv

Thank you in advance.

[erex]

Thanks for the data; I don't know when time will be available to dig deeply into the situation.

[erex]

Summary of issue (with data provided)

mrds 3.0.0 optimiser, using slsqp produces substantially different fitted model compared to old optimisation methods solnp and MCDS

library(Distance)
glob <- read.csv(file="https://github.com/user-attachments/files/18163069/obs_data.csv")
newoptR <- ds(data=glob, max_adjustments = 2, mono_method = "slsqp", optimizer = "R")
oldoptR <- ds(data=glob, max_adjustments = 2, mono_method = "solnp", optimizer = "R")
optMCDS <- ds(data=glob, max_adjustments = 2, optimizer = "MCDS")

> summary(newoptR)$ds

Summary for ds object
Number of observations :  737 
Distance range         :  0  -  300 
AIC                    :  1998.288 
Optimisation           :  mrds (nlminb) 

Detection function:
 Half-normal key function 

Detection function parameters 
Scale coefficient(s): 
            estimate         se
(Intercept) 4.800132 0.02780724

                        Estimate          SE         CV
Average p              0.5008159  0.01260542 0.02516976
N in covered region 1471.5985190 53.27991039 0.03620547
> summary(oldoptR)$ds

Summary for ds object
Number of observations :  737 
Distance range         :  0  -  300 
AIC                    :  1836.716 
Optimisation           :  mrds (solnp) 

Detection function:
 Half-normal key function with cosine adjustment terms of order 2, 3 

Detection function parameters 
Scale coefficient(s): 
            estimate         se
(Intercept) 4.883668 0.03300186

Adjustment term coefficient(s):  
              estimate         se
cos, order 2 0.5133234 0.06902655
cos, order 3 0.1135018 0.11234024


Strict monotonicity constraints were enforced.
                        Estimate           SE         CV
Average p              0.3334279   0.01870346 0.05609446
N in covered region 2210.3728542 140.68512708 0.06364769

Strict monotonicity constraints were enforced.
> summary(optMCDS)$ds

Summary for ds object
Number of observations :  737 
Distance range         :  0  -  300 
AIC                    :  1836.716 
Optimisation           :  MCDS.exe 

Detection function:
 Half-normal key function with cosine adjustment terms of order 2, 3 

Detection function parameters 
Scale coefficient(s): 
            estimate        se
(Intercept)  4.88371 0.0330043

Adjustment term coefficient(s):  
              estimate         se
cos, order 2 0.5133043 0.06902667
cos, order 3 0.1134970 0.11233955


Strict monotonicity constraints were enforced.
                        Estimate           SE         CV
Average p              0.3334432   0.01870459 0.05609528
N in covered region 2210.2714564 140.67989613 0.06364824

Strict monotonicity constraints were enforced.

[Globicephala]

Hi again,

Thank you for sending the message, and for looking into the data with more detail. I was discussing this with Tiago just now, who might comment in this issue soon enough.

What your output shows, which is something that we had already discussed here, is that the new optimiser, for this dataset, is not going beyond the first cosinde adjustment, as the AIC goes up when compared to the half-normal without cosine adjustments. This, in my opinion, doesn't necessarily mean that if we added another cosine adjustment, that the AIC wouldn't drop again. This leads me to believe that there might be an issue with the new optimiser, as it doesn't lead to the models with lowest AIC.

It will probably make the most sense to fit the detection function using the older optimisation methods.

[LHMarshall]

@Globicephala thanks for bring this to attention.

@lenthomas @erex I've run the 3 optimisers forcing 0, 1 and 2 adjustments - there seems to be an issue with the optimisation with the slsqp for one adjustment - this also means that with the slsqp optimiser and ds for automatic selection of the number of adjustments it never tries 2 adjustments. This is concerning as the models with 2 adjustments have the lowest AIC values. See table of AIC values and plots below.

glob <- read.csv(file="obs_data.csv")

# New optimiser
newoptR_0 <- ds(data=glob, nadj = 0, mono_method = "slsqp", optimizer = "R")
newoptR_1 <- ds(data=glob, nadj = 1, mono_method = "slsqp", optimizer = "R")
newoptR_2 <- ds(data=glob, nadj = 2, mono_method = "slsqp", optimizer = "R")

# Old optimiser
oldoptR_0 <- ds(data=glob, nadj = 0, mono_method = "solnp", optimizer = "R")
oldoptR_1 <- ds(data=glob, nadj = 1, mono_method = "solnp", optimizer = "R")
oldoptR_2 <- ds(data=glob, nadj = 2, mono_method = "solnp", optimizer = "R")

# MCDS optimiser
optMCDS_0 <- ds(data=glob, nadj = 0, optimizer = "MCDS")
optMCDS_1 <- ds(data=glob, nadj = 1, optimizer = "MCDS")
optMCDS_2 <- ds(data=glob, nadj = 2, optimizer = "MCDS")


results <- data.frame(optimiser = c("slsqp", "solnp", "MCDS"), 
                      adj0 = c(AIC(newoptR_0)$AIC, AIC(oldoptR_0)$AIC, AIC(optMCDS_0)$AIC), 
                      adj1 = c(AIC(newoptR_1)$AIC, AIC(oldoptR_1)$AIC, AIC(optMCDS_1)$AIC), 
                      adj2 = c(AIC(newoptR_2)$AIC, AIC(oldoptR_2)$AIC, AIC(optMCDS_2)$AIC))

results
#  optimiser     adj0     adj1     adj2
#1     slsqp 1998.288 2077.666 1836.717
#2     solnp 1998.288 1855.506 1836.716
#3      MCDS 1998.288 1855.506 1836.716


par(mfrow = c(3,3))
plot(newoptR_0, main = "slsqp no adj")
plot(newoptR_1, main = "slsqp 1 adj")
plot(newoptR_2, main = "slsqp 2 adj")
plot(oldoptR_0, main = "solnp no adj")
plot(oldoptR_1, main = "solnp 1 adj")
plot(oldoptR_2, main = "solnp 2 adj")
plot(optMCDS_0, main = "MCDS no adj")
plot(optMCDS_1, main = "MCDS 1 adj")
plot(optMCDS_2, main = "MCDS 2 adj")