Encounter rate variance of Buckland 2001 cannot be requested by `ds`

[erex]

er_method argument in ds is ignored by dht

Either observed counts (er_method=1) or expected counts (er_method=2) can be requested. Default is "2". When user request er_method=1 results are still using er_method=2

I believe the problem lies within dht, however, users calling ds from the Distance package will not get the results they expect if they specify er_method=1

MRE

library(Distance)
data(Savannah_sparrow_1981)
cu <- convert_units("meter", NULL, "hectare")
pooler1 <- ds(Savannah_sparrow_1981, transect="point", key="hn", adjustment = "cos",
           convert_units = cu, er_method = 1)  # encounter var uses observed counts
pooler1$dht$individuals$summary

pooler2 <- ds(Savannah_sparrow_1981, transect="point", key="hn", adjustment = "cos",
              convert_units = cu, er_method = 2)  # encounter var uses expected counts
pooler2$dht$individuals$summary

se.ER and cv.ER are identical regardless of value of er_method.

Debug trail

• Inside ds
  • dht_option object varflag is set to value of `er_option
• Inside dht
  • options$varflag is given default value of 2
  • dht.se is called with options as an argument
• Inside dht.se
  • encounter rate variance computed depending upon value of options$varflag

[erex]

Supplement to this issue

I checked whether the lack of difference in ER variance estimates was caused by either a) miscommunication between ds and dht and b) lack of covariate in the Savannah sparrow detection function model.

I turned my attention to the amakihi data set to introduce a covariate in the detection function and I ran the analysis in ddf rather than through ds. Result below: my interpretation is that I still can't induce differences in ER variance estimation using varflag:

data("amakihi")  # multiple years as strata
bust <- unflatten(amakihi)
detsize <- ddf(dsmodel=~mcds(key="hn", formula=~MAS),
               data=bust$data, method="ds", meta.data = list(point=TRUE, width=82.5))
summary(detsize)
meth1 <- dht(detsize, region.table = bust$region.table, sample.table = bust$sample.table,
             obs.table = bust$obs.table, options = list(varflag=1))
meth1$individuals$summary
meth2 <- dht(detsize, region.table = bust$region.table, sample.table = bust$sample.table,
             obs.table = bust$obs.table, options = list(varflag=2))
meth2$individuals$summary

Result

> meth1 <- dht(detsize, region.table = bust$region.table, sample.table = bust$sample.table,
+              obs.table = bust$obs.table, options = list(varflag=1))
> meth1$individuals$summary
  Region      Area CoveredArea Effort    n   k       ER     se.ER      cv.ER
1   1292  855298.6    855298.6     40  140  40 3.500000 0.3121472 0.08918492
2    194  876681.1    876681.1     41  172  41 4.195122 0.3088521 0.07362173
3    493  876681.1    876681.1     41  231  41 5.634146 0.3289972 0.05839344
4    494  513179.2    513179.2     24  141  24 5.875000 0.2712859 0.04617632
5    495  855298.6    855298.6     40  212  40 5.300000 0.4938078 0.09317129
6    792  876681.1    876681.1     41  146  41 3.560976 0.1945877 0.05464449
7    793  855298.6    855298.6     40  201  40 5.025000 0.2984737 0.05939775
8  Total 5709118.2   5709118.2    267 1243 267 4.655431 0.1262862 0.02712664
> meth2 <- dht(detsize, region.table = bust$region.table, sample.table = bust$sample.table,
+              obs.table = bust$obs.table, options = list(varflag=2))
> meth2$individuals$summary
  Region      Area CoveredArea Effort    n   k       ER     se.ER      cv.ER
1   1292  855298.6    855298.6     40  140  40 3.500000 0.3121472 0.08918492
2    194  876681.1    876681.1     41  172  41 4.195122 0.3088521 0.07362173
3    493  876681.1    876681.1     41  231  41 5.634146 0.3289972 0.05839344
4    494  513179.2    513179.2     24  141  24 5.875000 0.2712859 0.04617632
5    495  855298.6    855298.6     40  212  40 5.300000 0.4938078 0.09317129
6    792  876681.1    876681.1     41  146  41 3.560976 0.1945877 0.05464449
7    793  855298.6    855298.6     40  201  40 5.025000 0.2984737 0.05939775
8  Total 5709118.2   5709118.2    267 1243 267 4.655431 0.1262862 0.02712664

[LHMarshall]

Added debug lines at 254 in dht.se to print out vc2 values which should be calculated differently based on varflag - these values differ as expected.

meth1 <- dht(detsize, region.table = bust$region.table, sample.table = bust$sample.table, obs.table = bust$obs.table, options = list(varflag=1))
#VC2 for varflag =  1
#[1] 1160.1559 1185.8751 1500.0899  385.9715 3269.8088  516.9155 1229.0232

meth2 <- dht(detsize, region.table = bust$region.table, sample.table = bust$sample.table, obs.table = bust$obs.table, options = list(varflag=2))
#VC2 for varflag =  2
#[1] 1270.3755 1206.1956 1441.8924  459.4968 3809.8326  563.3517 1428.4070

Check if other tables in the output have different variances - yes they do

 meth1$individuals$D
#  Label     Estimate           se         cv          lcl          ucl        df
#1  1292 0.0004465283 4.274438e-05 0.09572603 0.0003686428 0.0005408692  51.72535
#2   194 0.0005335459 4.357829e-05 0.08167673 0.0004532601 0.0006280528  60.49008
#3   493 0.0007565767 5.046891e-05 0.06670693 0.0006623785 0.0008641710  67.91837
#4   494 0.0008290664 4.895083e-05 0.05904331 0.0007368170 0.0009328655  61.02400
#5   495 0.0007175651 7.082074e-05 0.09869591 0.0005887583 0.0008745519  49.08281
#6   792 0.0004745935 3.011270e-05 0.06344944 0.0004182646 0.0005385084  72.42554
#7   793 0.0006900687 4.707747e-05 0.06822142 0.0006023283 0.0007905900  67.65184
#8 Total 0.0006232863 2.619752e-05 0.04203128 0.0005739514 0.0006768619 847.12622
 meth2$individuals$D

#  Label     Estimate           se         cv          lcl          ucl        df
#1  1292 0.0004465283 4.447191e-05 0.09959483 0.0003657713 0.0005451153  50.55352
#2   194 0.0005335459 4.388059e-05 0.08224332 0.0004527398 0.0006287746  60.11195
#3   493 0.0007565767 4.971306e-05 0.06570790 0.0006637249 0.0008624180  69.18879
#4   494 0.0008290664 5.172401e-05 0.06238825 0.0007316702 0.0009394276  53.79278
#5   495 0.0007175651 7.585370e-05 0.10570984 0.0005804841 0.0008870178  47.58620
#6   792 0.0004745935 3.109973e-05 0.06552919 0.0004164996 0.0005407905  69.41720
#7   793 0.0006900687 4.988830e-05 0.07229469 0.0005973597 0.0007971659  63.21186
#8 Total 0.0006232863 2.673753e-05 0.04289767 0.0005729701 0.0006780211 764.63549

Where do ER variances come from in the summary table?

The summary.table already exists before dht.se is called (line 380 dht.R) and these values are not used in the subsequent calculations.

The ER values in the summary.table are obtained using code in lines 307 - 351 in dht.R. The ER variances are calculated using the varn function which applies the specified variance estimator (R2, O2, etc.) but is not affected by varflag. This looks like the values in the summary table are applicable when using varflag = 1 but they are not what is used when varflag = 2 (which is the default). Having the summary table different to the default is somewhat confusing. Looks like I could change the summary table, currently this code is difficult to debug as tables.dht extists inside dht (advise making a function that exists outside that function) would also have to change an sapply call to a for loop over strata I think.

    # ER variance for each stratum, total below
    var.er <- sapply(split(Nhat.by.sample, Nhat.by.sample$Region.Label),
                     function(x) varn(x$Effort.x, x$n, type=options$ervar))

   ...


    # calculate encounter rate
    er <- summary.table$n/summary.table$Effort

    # get totals of easy stuff
    if(numRegions > 1){
      summary.table <- data.frame(Region=c(levels(summary.table$Region),
                                           "Total"),
                                  rbind(summary.table[, -1],
                                        apply(summary.table[, -1], 2, sum)))
    }

    # get total encounter rate and its variance
    if(numRegions > 1){
      er <- c(er, sum(er * summary.table$Area[1:length(er)])/
                    sum(summary.table$Area[1:length(er)]))
      total.var.er <- sum(var.er * summary.table$Area[1:length(var.er)]^2)/
                        sum(summary.table$Area[1:length(var.er)])^2
      var.er <- c(var.er, total.var.er)
    }

    summary.table$ER <- er
    summary.table$se.ER <- sqrt(var.er)
    summary.table$cv.ER <- summary.table$se.ER/summary.table$ER
    summary.table$cv.ER[summary.table$ER==0] <- 0

    # set missing values to 0
    summary.table$ER[is.nan(summary.table$ER)] <- 0
    summary.table$se.ER[is.nan(summary.table$se.ER)] <- 0
    summary.table$cv.ER[is.nan(summary.table$cv.ER)] <- 0

[LHMarshall]

Fix this in the summary table so that it uses the varflag

[LHMarshall]

Check that summary function is just printing a table and that no additional calculations are carried out

[erex]

Further comparison of er_method=1 and er_method=2 when covariate is in detection function. Mock example with ducknests

ducknests.pdf

[LHMarshall]

@erex thanks for those checks I've created a separate issue (see above)