dht2 doubling Effort

[emreuter]

Hi, I think I may have uncovered a bug in dht2 -

I am currently working on some data to develop abundance estimates for penguins and seals. We have data taken from ship surveys and aerial surveys.

I have fitted two detection functions, one for each transect type, and want to combine these into one overall density and abundance estimate. I have used dht2 to specifiy the dataframes with observations and transect information, with ddf_id specifiying which detection function applies. There is no stratification necessary. However, when I run it, the Effort and therefore, Covered Area double. The number of sightings remains unchanged, and so the outcomes are roughly half of what I expect. I have checked my dataframes for repeat observations and that the correct lengths are provided to dht2.

When I run them individually using summary() their densities are 0.64 and 0.61, when run together in dht2, the density is 0.28.

Here is some code I wrote using minke to show what I mean, I have removed the stratification, and used 'north' and 'south' as how I treated my 'aerial' and 'ship' transects, producing individual detection functions and specifing that to dht2.

Please let me know if I have gone wrong somewhere, or you need more information. Thanks!

`library(Distance)
data(minke)

minke$Platform <- factor(minke$Region.Label, levels = c("North", "South"))
minke$ddf_id <- as.numeric(minke$Platform)
minke$Region.Label <- 'All'

minke.north <- subset(minke, Platform != 'South')
minke.north.ds <- ds(minke.north, truncation = 1)

minke.south <- subset(minke, Platform != 'North')
minke.south.ds <- ds(minke.south, truncation = 1)

minke.ddf <- list(minke.north.ds, minke.south.ds)

minke.obs <- unique(data.frame(
object = minke$object,
Sample.Label = minke$Sample.Label,
distance = minke$distance,
Region.Label = "All",
ddf_id = minke$ddf_id))

minke.effort <- unique(data.frame(
Sample.Label = minke$Sample.Label,
Effort = c(minke$Effort),
Region.Label = "All",
ddf_id = minke$ddf_id,
Platform = minke$Platform))

minke.strat <- data.frame(Region.Label = 'All', Area = sum(unique(minke$Area)))

minke.dht <- dht2(
observations = minke.obs,
transects = minke.effort,
geo_strat = minke.strat,
ddf = minke.ddf,
binomial_var = FALSE,
strat_formula = ~1)

minke.dht
sum(minke.effort$Effort)`

[LHMarshall]

@emreuter If I am understanding correctly then you have 2 surveys (one ship and one aerial) covering the same area. This is a single study region. You would like to combine these 2 surveys to produce a single estimate of abundance for each species. I was discussing this with @lenthomas recently and we felt that this was a scenario in which it is suitable to use the stratification = "effort_sum" option in dht2 (which is actually an effort weighted mean). However, I have been working on a vignette to explore these options and I have discovered that it is not currently permitted to have this stratification option alongside multiple detection functions. We will need to look at this further. I did however manage to simulate some data and as a workaround fit a single detection function but had survey as a covariate, meaning that the detection functions could vary in some way for the different platforms - not sure if this is an option for you?

I think we need to look at this more as it should make sense to be able to provide different detection functions in this scenario. I'm not sure what the gain is using dht2 when doing it this way, there is a slight difference in estimated variace but I will look to discuss this more with @lenthomas

Simulation example below

# Simulate data
library(dsims)
set.seed(747)

# Make region and density surface
region <- make.region()
density <- make.density()
density <- add.hotspot(density, c(1000,250), 200, 1)
plot(density)

#Define population
pop.desc <- make.population.description(region,
                                        density,
                                        N=500)

# Define detectability - different for each survey
detect.1 <- make.detectability(scale.param = 10)
detect.2 <- make.detectability(scale.param = 25)

# Parallel line design
design.1 <- make.design(design = "systematic",
                        spacing = 100)

# Zigzag design
design.2 <- make.design(design = "eszigzag",
                        design.angle = 90,
                        spacing = 100)

# Basic analysis
analyses <- make.ds.analysis()

# Create simulations
sim.1 <- make.simulation(reps = 1,
                         design = design.1,
                         population.description = pop.desc,
                         detectability = detect.1,
                         ds.analysis = analyses)

sim.2 <- make.simulation(reps = 1,
                         design = design.2,
                         population.description = pop.desc,
                         detectability = detect.2,
                         ds.analysis = analyses)


survey.1 <- run.survey(sim.1)
survey.2 <- run.survey(sim.2)

plot(survey.1)
plot(survey.2)

# Need to update the data before fitting the models
data.1 <- survey.1@dist.data
data.2 <- survey.2@dist.data
# Make data.2 unique in terms of transects and objects
data.2$object <- data.2$object + max(data.1$object)
data.2$Sample.Label <- data.2$Sample.Label + max(data.1$Sample.Label)

# Put modified data back
survey.2@dist.data <- data.2

mod.1 <- analyse.data(analyses, survey.1)
mod.2 <- analyse.data(analyses, survey.2)

plot(mod.1)
plot(mod.2)

# Add survey identifyer
data.1$survey <- "Plane"
data.2$survey <- "Ship"
# Create full data set
all.data <- rbind(data.1, data.2)

# But this fails!
dht2.results <- dht2(ddf = list(mod.1, mod.2),
                     flatfile = all.data,
                     stratification = "effort_sum",
                     binomial_var = FALSE,
                     total_area = 1e+06,
                     strat_formula = ~survey)

# Create an analysis with survey as a covariate
analyses.all <- make.ds.analysis(dfmodel = list(~survey),
                                 truncation = 50)
# Fit a single detection function
mod.all <- analyse.data(analyses.all, all.data)

plot(mod.all)

# Run dht2
dht2.results <- dht2(ddf = mod.all,
                     flatfile = all.data,
                     stratification = "effort_sum",
                     binomial_var = FALSE,
                     total_area = 1e+06,
                     strat_formula = ~survey)

# Compare results
dht2.results

mod.all$dht$individuals$N

[emreuter]

Thanks Laura, I had considered this myself - but they need different truncations. The aerial transects were taken in a helicopter where the 100m either side of the trackline were obscured by the fuselage, so there needed to be a left truncation that was not necessary in the ship transects. The tuncations were aerial (100-800m) & ship (0-800m). My understanding is that in order for them to be used as covariates in the same detection function, the truncations need to be the same. Is that right? Thanks for your thoughts!

[LHMarshall]

@emreuter yes that is correct. Ok, I am looking through this code at the moment any way so I will get back to you next week with code which will do the calculations manually and will later look to extending dht2 to include this functionality. Laura