bootdht giving output inconsistent with analytical results

[LHMarshall]

I'm not convinced everything is fine with bootdht, even when merely dealing with geographic stratification. Here is an example with the sikadeer data set, indirect survey with creation and decay rate plus strata.

library(Distance)
data(sikadeer)
cu <- convert_units("centimeter", "kilometer", "square kilometer")
easy <- ds(sikadeer, key="hn", truncation="10%", convert_units = cu)
# Create list of multipliers
mult <- list(creation = data.frame(rate=25, SE=5),
 decay = data.frame(rate=163, SE=14.2))
print(mult)
deer.ests <- dht2(easy, flatfile=sikadeer, strat_formula=~Region.Label,
 convert_units=cu, multipliers=mult, 
 stratification="geographical")
print(deer.ests, report="density")
plain.estimates <- attr(deer.ests, which = "density")
analy.cis <- plain.estimates[, c("Region.Label", "LCI", "UCI")]

How does the bootstrap handle this?

bootnosf <- bootdht(easy, flatfile = sikadeer, resample_transects = TRUE,
 convert_units=cu, summary_fun=bootdht_Dhat_summarize,
 multipliers = mult, cores=10, nboot=200)
library(dplyr)
boot.cis <- bootnosf %>%
 group_by(Label) %>%
 summarize(quant025 = quantile(Dhat, probs = 0.025), 
 quant975 = quantile(Dhat, probs = 0.975))
print(boot.cis)

Output

Analytical

analy.cis
 Region.Label LCI UCI
1 A 43.006276 127.032567
2 B 19.374792 71.306900
3 C 1.905592 7.983742
4 E 1.101750 12.052725
5 F 9.747757 22.992961
6 G 1.573342 43.363109
7 H 1.008389 2.378582
8 J 4.705814 11.100050
9 Total 12.643932 34.370969

Bootstrap result

as.data.frame(boot.cis)
 Label quant025 quant975
1 A 51.3049436 98.679598
2 B 19.2213441 55.138872
3 C 2.7838666 5.298189
4 E 0.7472652 6.551853
5 F 13.2874121 16.994035
6 G 4.2230542 27.781516
7 H 1.3745599 1.758004
8 J 6.4146128 8.204017
9 Total 118.9993812 188.353451

Stratum-specific confidence interval bounds are roughly comparable, but the interval for the entire study area are quite disparate. I don't suspect bootdht is calculating a weighted average, weighted by stratum areas. It cannot be taking a simple average because thse values are larger than any stratum-specific estimates. I don't know what calculation is being made by bootdht

Whatever the situation, the documentation for bootdht that states
Note that only geographical stratification as supported in dht is allowed.

is not accurate.

Originally posted by @erex in #155

[LHMarshall]

Eric's investigations continued...

I continue my investigation of bootdht and the way it deals with "unusual" situations. Here my investigation focuses upon multipliers.

It is indeed the case that bootdht has a multipliers argument, in line 194, bootdht calls dht2_multipliers to add the variance of multipliers to the flatfile.

Line 297 of bootdht calls bootit which does the real work of the function, the multipliers are carried into bootit via an argument.

Line 47 of bootit multiplies together all multipliers and carries out the multiplier adjustments to the estimates (line 61).

In this workflow, I cannot find where the uncertainty in the multipliers is brought into the calculation. I am of the opinion that the bootstrap does incorporate multiplier adjustments to estimates, but does not propagate the uncertainty in the multipliers into the abundance or density estimates.

Minimum reproducible example: multiplier uncertainty exclusion
Altered the sikadeer data set to perform a bootstrap with 30%CVs in creation and decay rate and similar analysis where the multiplier CVs were set to 0. The reported stratum-specific confidence intervals are virtually identical after 1000 bootstrap replicates. Evidence to me that multiplier uncertainty is not propogated in bootdht.

library(Distance)
data(sikadeer)
cu <- convert_units("centimeter", "kilometer", "square kilometer")
easy <- ds(sikadeer, key="hn", truncation="10%", convert_units = cu)
mult.se <- list(creation = data.frame(rate=25, SE=25*0.3),
 decay = data.frame(rate=163, SE=163*0.3))
mult.nose <- list(creation = data.frame(rate=25, SE=0),
 decay = data.frame(rate=163, SE=0))
deer.ests <- dht2(easy, flatfile=sikadeer, strat_formula=~Region.Label,
 convert_units=cu, multipliers=mult.se, 
 stratification="geographical")
print(deer.ests, report="density")

boot.with.se <- bootdht(easy, flatfile = sikadeer, resample_transects = TRUE,
 convert_units=cu, summary_fun=bootdht_Dhat_summarize,
 multipliers = mult.se, 
 cores=10, nboot=1000)
boot.no.se <- bootdht(easy, flatfile = sikadeer, resample_transects = TRUE,
 convert_units=cu, summary_fun=bootdht_Dhat_summarize,
 multipliers = mult.nose, 
 cores=10, nboot=1000)

plain.estimates <- attr(deer.ests, which = "density")
analy.cis <- plain.estimates[, c("Region.Label", "LCI", "UCI")]

library(dplyr)
boot.se <- boot.with.se %>%
 group_by(Label) %>%
 summarize(LowerBound = quantile(Dhat, probs = 0.025), 
 UpperBound = quantile(Dhat, probs = 0.975))

boot.without <- boot.no.se %>%
 group_by(Label) %>%
 summarize(LowerBound = quantile(Dhat, probs = 0.025), 
 UpperBound = quantile(Dhat, probs = 0.975))

knitr::kable(as.data.frame(boot.se), caption="Multipliers with 30%CV", 
 row.names = FALSE, digits=3)
knitr::kable(as.data.frame(boot.without), caption="Multipliers without uncertainty", 
 row.names = FALSE, digits=3)

Output

Table: Multipliers with 30%CV

|Label | LowerBound| UpperBound|
|:-----|----------:|----------:|
|A | 51.246| 100.438|
|B | 19.011| 55.353|
|C | 2.684| 5.371|
|E | 0.839| 7.135|
|F | 13.201| 16.930|
|G | 4.243| 28.815|
|H | 1.366| 1.751|
|J | 6.373| 8.173|
|Total | 117.509| 194.044|

Table: Multipliers without uncertainty

|Label | LowerBound| UpperBound|
|:-----|----------:|----------:|
|A | 50.911| 98.927|
|B | 19.176| 52.884|
|C | 2.647| 5.381|
|E | 0.860| 7.180|
|F | 13.149| 16.847|
|G | 4.189| 28.461|
|H | 1.360| 1.743|
|J | 6.348| 8.133|
|Total | 117.460| 186.972|