`bootdht` and `dht2` incompatibility

[erex]

For analyses requiring dht2, e.g. non-geographic stratification or multipliers (not using the activity package, bootstrapping cannot be performed.

This is because bootdht only accepts as input dsmodel objects, not the sort of objects created by dht2

[LHMarshall]

@erex does this need to be resolved by a particular date? Is there an issue when using the activity package?

[erex]

No, no date. It should either be fixed (so the two functions are compatible) or there should be documentation to indicate (or test for) circumstances where a dsmodel object will cause dht2 to choke.

[LHMarshall]

Document this for the next milestone

[erex]

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.

[erex]

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.

[erex]

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|