Computation of CVs for individual and group density

[erex]

Because of a question to the Distance list, I did some manual calculation of CVs for the ClusterExercise sample data set. I found disagreement in the manner in which CVs were reported.

Zip file of .Rmd and .html below

CV comparison.zip

[erex]

Similar computations carried out with output (not raw data) provided by Arthur (user referenced in Distance list question).
arthurs cvs.pdf

[dill]

@lenthomas and I had a long chat about this. There are quite a few things going on here and we ended-up with more questions than answers... Some thoughts and reminders to self:

1. the comparison you made in the .Rmd is not quite right as the default estimator in mrds is the Innes et al (2002) so your use of CV(ER) to compare with the total CV shouldn't work (the missing component is given in equation 11 of that paper)
2. we think that the cluster abundance (i.e., setting size=1) should give the same answers for Innes and Buckland et al (2001) variance estimators when there are no covariates in the detection function
3. setting the example to use the Buckland et al (2001) estimator (setting er_method=1) still leads to the wrong answer, so there's still something going wrong here (more on this later)
4. there is a weird rescaling of the variance which we currently don't understand for the Buckland et al (2001) case, which we assume is the source of the issue
5. when cluster size is not a covariate (but there are other covariates in the detection function) then Distance for Windows uses Buckland et al (2001) not Innes et al (as you might assume from what's in the Advanced Distance Sampling book or Marques et al (2003)).
6. Distance for Windows doesn't do analytic variance estimation when cluster size is included in the detection function (bootstrap only)
7. comment in the MCDS source suggests we should look into this c. 2005

[dill]

I followed-up on your calculations @erex by configuring the same analysis in Distance for Windows. The pooled estimates are as follows:

Pooled Estimates:
                         Estimate      %CV     df     95% Confidence Interval
                        ------------------------------------------------------
                 DS     0.21934E-01   31.02    14.17 0.11459E-01  0.41985E-01
                 D      0.51634E-01   33.75    17.95 0.25895E-01  0.10296    
                 N       36935.       33.75    17.95  18523.       73647.    

Note the CV here is very close to that reported by mrds/Distance when the Buckland et al (2001) method is used:

Density:
  Label   Estimate          se        cv         lcl        ucl       df
1 North 0.01931297 0.007467003 0.3866315 0.008591546 0.04341371 12.44828
2 South 0.04141444 0.010105480 0.2440086 0.024909968 0.06885418 16.58302
3 Total 0.02193104 0.006801881 0.3101486 0.011457426 0.04197893 14.16840

(Comparing the Total and DS lines.)

I think the difference here comes from the fact that internally the R (maybe also Fortran) code doesn't simply add the CVs, but instead calculates the variances (actually covariance matrices) of the components on the abundance/density scale rather than on the their own scale (i.e., uncertainty about abundance/density with respect to p or n/L). These are then summed (rather than using the CV-squared trick).

[erex]

Thanks @dill for continuing to work through this. I presume the results you present in comment uses the ClusterExercise data set and a hazard rate key function without adjustments and er_method=1. Am I right?

That results in agreement for DistWin and ds when group size is not considered.

That gets us to point 2 in the comment above, correct? I'm assuming there is still detective work to do regarding uncertainty calculations with group size?

In the end, I don't yet know how to get back to Arthur with his question about how to perform variance components assessment for his situation. The last paragraph

R (maybe also Fortran) code doesn't simply add the CVs, but instead calculates the variances (actually covariance matrices) of the components on the abundance/density scale rather than on the their own scale (i.e., uncertainty about abundance/density with respect to p or n/L). These are then summed (rather than using the CV-squared trick).

is something I didn't know and I'm unsure how I would explain it to a user who wanted to get their hands dirty. Is this written down somewhere?

[lenthomas]

Another progress report (me and Dave are working on this at present):

We have worked through the calculations in the original question and have found an issue with the way that CV in total encounter rate is calculated in post-processing. The CV in total density is calculated correctly, but CV in ER (which is calculated separately, after CV density is calculated - i.e., not used in calculation of CV density) is not. This is why the manual calculations of overall CV were not working. Here's a demonstration that when we calculate CV in ER manually and correctly, we get the same answer as ds on CV in density:

data(ClusterExercise)
groups <- ds(ClusterExercise, key="hr", er_method = 1)
s.groups <- summary(groups)
s.groups2 <- s.groups$dht$clusters$summary
#Show CV on ER
s.groups2

#Manual calculation of Total cv.ER

#Just keep strata estimates
s.groups3 <- s.groups2[1:2, ]
#Total ER is area-weighted average of ER in each stratum
ER.tot <- with(s.groups3, sum(Area * ER) / sum(Area))
#Total var is area^2-weighted average of var in each stratum
var.ER.tot <- with(s.groups3, sum(Area ^ 2 * se.ER ^ 2) / sum(Area) ^ 2)
#CV is sqrt(var / ER^2)
cv.ER.tot <- sqrt(var.ER.tot / ER.tot ^ 2)

#Note this is different from what s.groups2 reports
print(cv.ER.tot)

#Work out total CV on density
cv.p <- s.groups$ds$average.p.se[1,1] / s.groups$ds$average.p
cv.tot <- sqrt(cv.p ^ 2 + cv.ER.tot ^ 2)
print(cv.tot)

#Note the above cv.tot *is* what s.groups$dht$clusters$D
s.groups$dht$clusters$D

#So, ds seems to be getting the total CV on encounter rate wrong, but
# the total CV on density right.  This appears then to be "just" a 
# post-processing issue in summary

Dave has isolated the bug - it occurs in mrds in dht.R line 310.

(Unfortunately, total CV in ER is not reported by DistWin which made it a bit harder to double check, but manual calculation did the trick.)

[lenthomas]

This bug also affected the calculated total ER, which was also reported incorrectly.

[lenthomas]

@erex could you please send (via email so private) Arthur's data to @dill and me, so we can check fixing this also resolves his issue, mentioned above.

[lenthomas]

BTW @erex, we are also discussing how to enhance the documentation around DistWin and the R packages to make the methods for point and variance estimation with and without covariates more transparent. Will be happy to jointly discuss this at a DistDev meeting.

[lenthomas]

Another BTW for the record: the CV on density of clusters is not affected by whether we use the Buckland et al. (2001) method or the Innes et al. (2002) method:

data(ClusterExercise)
groups <- ds(ClusterExercise, key="hr", er_method = 1)
tmp <- summary(groups)
tmp$dht$clusters$D

groups <- ds(ClusterExercise, key="hr", er_method = 2)
tmp <- summary(groups)
tmp$dht$clusters$D

Both the same:

This is because things cancel out so the two variance estimators are the same mathematically.

However, the CV on density of individuals does depend on the variance estimation method.

data(ClusterExercise)
groups <- ds(ClusterExercise, key="hr", er_method = 1)
tmp <- summary(groups)
tmp$dht$individuals$D

groups <- ds(ClusterExercise, key="hr", er_method = 2)
tmp <- summary(groups)
tmp$dht$individuals$D

Gives

Different, with the CV calculated using Innes et al. (method = 2) smaller than Buckland et al (method = 1). This is because the former now includes an additional component for CV(E(s)) while the latter just uses the line-to-line CV in variance of estimated individuals, which is the same as the line-to-line CV in variance of estimated clusters since there are no covariates. In this sense, method = 2 is ignoring cluster size variance, and so appears to be underestimating overall variance. (This would not happen if there is a covariate in the df model that is correlated with cluster size, or even cluster size is in the model.)

This topic needs more thinking regarding what is the best default, and then some documentation!