`dht2()` incorrectly handling binned distances with truncation

[erex]

Email from user:

I'm having a few problems calculating density/abundance using dht2 and flatfiles that I have never had before.
dht2 doesn't seem to be choosing the right observations based on truncation distance. In the screenshot attached, it shows that there are six observations being used, but there are 7 in the file and only 5 are within the truncation distance (2km).
Additionally, in some other cases, I get this error:
Error in $<-.data.frame(*tmp*, "distance", value = 0.0709) : replacement has 1 row, data has 0
The distance value 0.0709 is specified as the left truncation distance in the detection function - but there are never detections at exactly that distance so I don't know why this error would show up in only some instances.

I'm not sure how much testing dht2 has received when using binned data with truncation.

[erex]

Follow up with original user still shows challenges when taking exact distance data, specifying both cutpoints and truncation. Nothing to do with dht2, but rather an issue with ds.

Demonstration of problem:

bins <- c(0.0709, .2315, .463, .926, 1.852)
AllHWhr <- ds(data=All_HW, key="hr", adjustment=NULL, cutpoint = bins,
              truncation=list(left=0.0709, right = 3.704),
              optimizer = "MCDS")

results in

Summary for ds object
Number of observations :  64 
Distance range         :  0.0709  -  3.704 
AIC                    :  212.8108 

In contrast, the following call to ds

bins2 <- c(0.0709, .2315, .463, .926, 1.852, 3.704)
AllHWhr.bins2 <- ds(data=All_HW, key="hr", adjustment=NULL, cutpoint = bins2,
                    truncation=list(left=0.0709, right = 3.704),
                    optimizer = "MCDS")

results in

Summary for ds object
Number of observations :  81 
Distance range         :  0.0709  -  3.704 
AIC                    :  255.4857 

I presume the values labelled Distance range are taken from truncation and in the first call to ds the reported Distance range is incorrect.

[LHMarshall]

I cannot replicate the above issue, I did not have access to the original data but I reconstructed what I could from the E-mail to the list which seemed to contain all the relevant distances.

The ds documentation states the following:
"When the data is binned, the right truncation is the largest bin end point."

This means that the right truncation is ignored when the cutpoint values are supplied. It is strange that it is then reported in the summary output (rather than the largest value in the binning cutpoints), this should be changed to the this value instead.

I never saw 6 observations I only saw either 5 or 7 as expected.

eg.data <- data.frame(
  Region.Label = rep("MyStratum", 15),
  Area = rep(9002, 15),
  Sample.Label = c("-7051", "-7100","-7042","-7042" ,"-7042" , "-7015", "-7006", "-7123", "-7114", "-7105", "-7056", "-7038", "-7029", "-7020", "-7011"),
  Effort = c(173.33232, 180.20727, 137.14122, 137.14122, 137.14122, 237.22610, 137.02171, 46.50701, 46.45782, 71.88510, 40.01047, 39.42089, 93.37925, 68.58188, 184.06937),
  Aircraft = rep(649, 15),
  distance = c(2.77800, 0.34725, 1.38900, 2.77800, 0.69450,0.69450, 0.34735, rep(NA,8)),
  size = c(3,1,1,1,2,1,1, rep(NA,8))
)

bins <- c(0.0709, .2315, .463, .926, 1.852)
AllHWhr <- ds(data=eg.data, key="hr", 
              adjustment=NULL, cutpoint = bins,
              truncation=list(left=0.0709, right = 3.704),
              optimizer = "MCDS")

summary(AllHWhr$ddf)

#Summary for ds object
#Number of observations :  5 
#Distance range         :  0.0709  -  3.704 
#AIC                    :  17.03973 
#Optimisation           :  MCDS.exe

#Summary statistics:
#  Region Area CoveredArea  Effort n  k          ER       se.ER     cv.ER
#1 MyStratum 9002    10574.07 1455.24 5 13 0.003435858 0.001405006 0.4089244


bins2 <- c(0.0709, .2315, .463, .926, 1.852, 3.704)
AllHWhr.bins2 <- ds(data=eg.data, key="hr", 
                    adjustment=NULL, cutpoint = bins2,
                    truncation=list(left=0.0709, right = 3.704),
                    optimizer = "MCDS")

summary(AllHWhr.bins2$ddf)
#Summary for ds object
#Number of observations :  7 
#Distance range         :  0.0709  -  3.704 
#AIC                    :  25.69781 
#Optimisation           :  MCDS.exe

#Summary statistics:
#  Region Area CoveredArea  Effort n  k          ER       se.ER     cv.ER
#1 MyStratum 9002    10574.07 1455.24 7 13 0.004810202 0.001898489 0.3946796

[LHMarshall]

I tried adding the following checks in get_truncation but this is not the only function that performs actions and checks on the truncation distances / cutpoins - further checks and possible manipulation are done in ds and possibly other places. It would take some in depth exploration to safely fix the output so that it matches the values that are used in the data filtering process.

# Final check that truncation is not greater than largest cutpoint
  if(!is.null(cutpoints)){
    if(width > cutpoints[length(cutpoints)]){
      # Make truncation largest cutpoint
      width <- cutpoints[length(cutpoints)]
    }
    # And that the left truncation is not smaller than the first cutpoint
    if(left < cutpoints[1]){
      # Make left truncation smallest cutpoint
      left <- cutpoints[1]
    }
  }

[LHMarshall]

Put in a warning at the minimum for the next release but ideally sort this out and tidy up!

[LHMarshall]

This looks more problematic than first thought...

First model: right truncation is larger than the largest cutpoint

bins <- c(0.0709, .2315, .463, .926, 1.852)
# Model where right truncation and largest cutpoint don't match
AllHWhr <- ds(data=eg.data, key="hr", 
              adjustment=NULL, cutpoint = bins,
              truncation=list(left=0.0709, right = 3.704),
              optimizer = "R")

summary(AllHWhr$ddf)

# Summary for ds object
# Number of observations :  5 
# Distance range         :  0.0709  -  3.704 
# AIC                    :  16.68454 
# Optimisation           :  mrds (nlminb) 
# 
# Detection function:
#   Hazard-rate key function 
# 
# Detection function parameters 
# Scale coefficient(s): 
#   estimate       se
# (Intercept) 0.1099626 2187.166
# 
# Shape coefficient(s):  
#   estimate       se
# (Intercept) 3.368696 99976.08
# 
# Estimate         SE        CV
# Average p            0.2942046 0.06674907 0.2268797
# N in covered region 16.9949738 7.45909927 0.4389003

plot(AllHWhr$ddf)
#Last interval endpoint did not include 3.704. It was reset.

The plot looks to indicate two points that are equal height in the first two bins and the thirds looks to be zero.

Note the abundance estimates of the cluster for this first model

AllHWhr$dht$clusters$N
#  Label Estimate       se        cv      lcl      ucl       df
#1 Total  14.4683 6.766054 0.4676469 5.603587 37.35672 14.88361

Second model: right truncation is the same as the largest cutpoint

# In this model the right truncation distance and the largest bin width are the same
bins <- c(0.0709, .2315, .463, .926, 1.852)
AllHWhr2 <- ds(data=eg.data, key="hr", 
              adjustment=NULL, cutpoint = bins,
              truncation=list(left=0.0709, right = 1.852),
              optimizer = "R")

summary(AllHWhr2$ddf)

# Summary for ds object
# Number of observations :  5 
# Distance range         :  0.0709  -  1.852 
# AIC                    :  16.65696 
# Optimisation           :  mrds (nlminb) 
# 
# Detection function:
#   Hazard-rate key function 
# 
# Detection function parameters 
# Scale coefficient(s): 
#   estimate        se
# (Intercept) -0.1416633 0.8976892
# 
# Shape coefficient(s):  
#   estimate       se
# (Intercept) 1.281357 2.226187
# 
# Estimate        SE        CV
# Average p           0.5504213 0.3016126 0.5479668
# N in covered region 9.0839503 5.6742586 0.6246466

This plot looks different to the first

The abundance estimates for the clusters also differ

AllHWhr2$dht$clusters$N
#Label Estimate       se        cv      lcl      ucl      df
#1 Total 15.77469 10.78563 0.6837301 3.606816 68.99179 6.74857

[LHMarshall]

ds() lines 406 - 412

    if(!is.null(left)){
      if(cutpoints[1]!=left){
        stop("The first cutpoint must be 0 or the left truncation distance!")
      }
    }else if(cutpoints[1]!=0){
      stop("The first cutpoint must be 0 or the left truncation distance!")
    }

Error message says "The first cutpoint must be 0 or the left truncation distance!" but code is triggered in the following situation where there was no value supplied for a left truncation distance. If you supply cutpoints then the user has to define the left truncation distance.

# data setup
data(book.tee.data)
egdata <- book.tee.data$book.tee.dataframe
egdata <- egdata[!duplicated(egdata$object),]

ds(egdata,4,cutpoints=c(2,3,4))