DistanceDevelopment / Distance

Simple distance sampling analysis
GNU General Public License v3.0
9 stars 8 forks source link

Case where mcds.exe fits with negative pdf #160

Open lenthomas opened 1 year ago

lenthomas commented 1 year ago

R 4.2.3; mcds 2.2.8.9005; Distance 1.0.6.9003

Documenting this case where mcds.exe fits a pdf with some negative values as the possible basis for one example in a DistanceExamples page showing differences between mcds.exe and the R optimizer. Could possibly use this as part of a story of what to do when you get a warning message. In this case the warning is of 0 pdf, and on investigation it turns out to be mcds.exe that's generating it so a possible solution is just to use the R optimizer in this case.

(Could also potentially investigate what's going on in the mcds optimizer at some point - but I'm not raising this as a DistWin issue for now.)

library(Distance)
data("wren_5min")
conversion.factor <- convert_units("meter", NULL, "hectare")
bin.cutpoints.100m <- bin.cutpoints <- c(0, 10, 20, 30, 40, 60, 80, 100)
#Gives warnings that detection function is less than 0 in some places
wren5min.hn.herm.t100 <- ds(data=wren_5min, key="hn", adjustment="herm", 
                            transect="point", cutpoints=bin.cutpoints.100m,
                            convert_units=conversion.factor)
#warning present in mcds.exe fit
wren5min.hn.herm.t100.mcds <- ds(data=wren_5min, key="hn", adjustment="herm", 
                            transect="point", cutpoints=bin.cutpoints.100m,
                            convert_units=conversion.factor,
                            skip_R = TRUE)
#warning absent in R optimizer fit
wren5min.hn.herm.t100.r <- ds(data=wren_5min, key="hn", adjustment="herm", 
                                 transect="point", cutpoints=bin.cutpoints.100m,
                                 convert_units=conversion.factor,
                                 skip_mcds = TRUE)
#both select the same model (adj(4)) but have different parameter estiataes
summary(wren5min.hn.herm.t100.mcds)
summary(wren5min.hn.herm.t100.r)
#lnl higher for mcds.exe fit
wren5min.hn.herm.t100.mcds$ddf$lnl #-209.1142
wren5min.hn.herm.t100.r$ddf$lnl #-209.3717
#.... but plotting clearly shows fitted pdf is negative around 100m
plot(wren5min.hn.herm.t100.mcds, pdf = TRUE)
plot(wren5min.hn.herm.t100.r, pdf = TRUE)

mcds fit: Screenshot 2023-03-21 070737

R optimizer fit: Screenshot 2023-03-21 070822