DistanceDevelopment / Distance

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

Camera trapping example produces seemingly non-monotonic solution #133

Closed lenthomas closed 2 years ago

lenthomas commented 2 years ago

Analysis of the camera trapping example with the QAIC-best model, unif + 3 cosine adjustments, seems to produce a solution in which the detection probability increases beyond about 13 meters. Distance 1.0.6 mrds 2.2.7

DuikerCameraTraps <- read.csv(file="https://datadryad.org/stash/downloads/file_stream/73221", 
                              header=TRUE, sep="\t")

library(Distance)
conversion <- convert_units("meter", NULL, "square kilometer")
trunc.list <- list(left=2, right=15)
mybreaks <- c(seq(2, 8, 1), 10, 12, 15)

#Set start points and then it converges
uni2 <- ds(DuikerCameraTraps, transect = "point", key="unif", adjustment = "cos",
           nadj=2,
           cutpoints = mybreaks, truncation = trunc.list)
uni3 <- ds(DuikerCameraTraps, transect = "point", key="unif", adjustment = "cos",
           nadj=3,
           cutpoints = mybreaks, truncation = trunc.list, 
           initial_values = list(adjustment = c(as.numeric(uni2$ddf$par), 0)))

summary(uni3)
plot(uni3, showpoints = FALSE)
Screenshot 2022-10-09 212657 Screenshot 2022-10-09 212736

Note that Distance for Windows produces a fit with a slightly higher AIC (so lower negative log-likelihood) but without the non-monotonicity. I would have thought that both would check monotonicity at the same evaluation points, so I struggle to see how this happens (below is from the DuikerDaytime sample project in DistWin, analysis 6).

Screenshot 2022-10-09 212446 Screenshot 2022-10-09 212513 Screenshot 2022-10-09 212600
dill commented 2 years ago

mrds uses roughly the same rule as Distance for Windows to create the monotonicity constraint grid. Specifically, the CDS engine uses this code:

      NCNLIN=MAXNLC/2
      XLAT=WIDTH/FLOAT(NCNLIN)**1.5
c      XLAT=WIDTH/FLOAT(NCNLIN)
      DO I=1,NCNLIN
          XGRID(I)=FLOAT(I)**1.5*XLAT
c          XGRID(I)=FLOAT(I)*XLAT
      END DO
      NCNLIN=MAXNLC

where MAXNLC takes a value of 20.

As it stands, mrds uses a similar rule, but modified to take into account left truncation:

getRefPoints<- function(no_d, int.range){

  # previous versions just used width and assumed that left truncation
  # was at zero, now using int.range to tell us the interval

  # note this currently doesn't work with multiple integration ranges
  # i.e. when int.range is a matrix w. > 1 row. This should have been caught
  # and an error throw before now though.

  xlat <- (int.range[2]-int.range[1])/(no_d^1.5)
  ref_points <- double(no_d)
  for(i in 1:no_d){
    ref_points[i] <- (i^1.5) * xlat
  }

  ref_points <- ref_points+int.range[1]

  return(ref_points)
}

The default value of no_d is 20. So we're using twice as many grid points to detect monotonicity.

Fiddling with the mrds code above to make it more like that in Distance for Windows yielded worse results.

dill commented 2 years ago

Using the Distance for Windows parameters in Distance/mrds yields the same plot, so I'm inclined to suspect this is a plotting issue rather than a fitting one...

control <- uni3$ddf$control
control$initial$adjustment <- c(0.9352, -0.05346, -0.08074)
control$nofit <- TRUE
rr <- ddf(dsmodel = ~cds(key = "unif", formula = ~1, adj.series = "cos", 
          adj.order = c(1, 2, 3), adj.scale = "width"), data = uni3$ddf$data, 
          method = "ds", meta.data =uni3$ddf$meta.data, control = control)
plot(rr)

Screenshot 2022-10-10 at 11 23 55

dill commented 2 years ago

It was the implementation of the scaling for adjustments. The code used width-left rather than just width when the scaling option was set to width. Fix in 9945e71.

lenthomas commented 2 years ago

Great that you found the scaling issue @dill ! 🥇 Regarding the points at which to check for constraints, do you think I should change mcds.exe so that it also takes account of left truncation? It makes sense to me, as g(x) for 0 >= x >= [left trucation point] should be 0 when it comes to things like calculating average p so no point in checking for non-monotonicity there.

lenthomas commented 2 years ago

Another thought: I just got the following warning message during fitting in an analysis (I can provide the analysis if required): Warning in mrds::check.mono(model, n.pts = 20) : Detection function is greater than 1 at some distances From this, it appears that there is a monotonicity check at 20 evenly spaced points between the left and right truncation points. This seems inconsistent with the constrained optimization that is done, where monotonicity is checked at 10 non-evenly spaced points between left and right truncation. Shouldn't these be consistent? Just a thought (and I'm happy to spin out into a separate issue as a feature enhancement request.)

dill commented 2 years ago

Moving these to a separate issue about consistency between Distance for Windows and mrds here.