Calibration of Camera Densities for Moose to Aerial Surveys by WMU
Dave Huggard, Marcus Becker
July 4, 2019
Objective
To compare density estimates for moose from ABMI cameras to estimates from aerial surveys, using Wildlife Management Units (WMU’s). The results can be used to calibrate camera density estimates for this species to the aerial survey estimates, which might be considered more reliable. The analysis also looks at how much of a contribution different numbers of cameras in a WMU could make to a density estimate that combines aerial survey and camera results.
# Import data
df_comp <- read_csv("./aerial-camera-combined_full.csv")
df_2019 <- read_csv("./OSM_Ungulate_Survey_Data_2019.csv")
ab_wmu1 <- st_read("C:/Users/mabec/Documents/R/ABMI-Moose/Data/SpatialData/AB_WMU-OSR-Intersect.shp",
stringsAsFactors = FALSE, quiet = TRUE)
# Obtain confidence intervals for 2019 estimates
df_size <- ab_wmu1 %>%
separate(WMUNIT_COD, into = c("zeroes", "WMU"), sep = 2, remove = FALSE) %>%
select(WMU, kmsquareac) %>%
st_set_geometry(NULL)
df_2019_ci <- df_2019 %>%
separate(Range, into = c("Lower", "Upper"), sep = "-", remove = TRUE) %>%
mutate(Upper = str_replace_all(Upper, ",", "")) %>%
mutate(Lower = as.numeric(Lower),
Upper = as.numeric(Upper),
WMU= as.character(WMU)) %>%
left_join(df_size, by = "WMU") %>%
# Calculate confidence intervals based on size of WMU
mutate(CIL = Lower / kmsquareac,
CIU = Upper / kmsquareac) %>%
select(WMU, Survey_Year, CIL, CIU)
# Prepare data for modeling
df_comp1 <- df_comp %>%
mutate(WMU = as.character(WMU.num)) %>%
left_join(df_2019_ci, by = c("WMU", "Survey_Year")) %>%
mutate(CIL = ifelse(is.na(CIL.x), CIL.y, CIL.x),
CIU = ifelse(is.na(CIU.x), CIU.y, CIU.x)) %>%
select(-c(CIL.x, CIL.y, CIU.x, CIU.y)) %>%
rename(Aerial.mean = Density,
Aerial.lci = CIL,
Aerial.uci = CIU) %>%
setNames(gsub("Moose", "Camera", names(.))) %>%
filter(!is.na(Camera.mean) & !is.na(Aerial.mean)) %>%
filter(Camera.mean < 5) %>% # There was one WMU with an over-active moose.
# Calculate weights for regression
mutate(wt1 = sqrt(n))
# write_csv(df_comp1, "./aerial-camera-combined_2014-19.csv")Methods
Density estimates for moose were calculated for each ABMI camera, using the time-in-field-of-view method. This analysis method and its assumptions are described elsewhere. Camera estimates are made for summer and winter seasons, then these are averaged with equal weight to compensate for any differences in seasonal deployment times between cameras. Average density for each WMU is calculated from all the cameras in the WMU. Confidence intervals are based on a compound distribution of binomial presence/absence and log-normal abundance-given-presence.
Aerial survey estimates were provided by the Government of Alberta. Estimates were provided with 90% confidence intervals. We do not currently have information on the analysis methods, the population that the confidence intervals are for, or assumptions of the methods.
We fit models of camera density as a function of aerial survey density, including GAMs with smoothing splines using both normal and log-normal (log-link) error distributions, and a linear model both with and without an intercept. Points were weighted in inverse proportion to the width of the camera confidence intervals, which vary widely due to large differences in number of cameras per WMU and inherent variability of camera estimates. Confidence interval width for aerial estimates were a remarkably consistent proportion of the mean estimate, and so were not used for weighting points.
One outlying point with aerial density of 0.5 per km\(^2\) but camera density of 7.1 per km\(^2\) was omitted. The extreme camera estimate is from a WMU with only 4 cameras, and is largely due to a single camera with an extended visit from one moose. The 90% confidence intervals for that camera estimate are 1.4-35.3 per km\(^2\), showing an extremely uncertain estimate. One point with an outlying aerial estimate of 0.77 per km^\(2\) was included in the model fits.
# Smoothing splines - log-link and identity
# Dave's GAM
m.gam1 <- gam(data = df_comp1, Camera.mean ~ s(Aerial.mean, k=6, m=2), family = "gaussian", weights = wt1)
# Marcus' GAM
m.gam2 <- gam(data = df_comp1, Camera.mean ~ s(Aerial.mean), family = "gaussian", method = "REML")
# Marcus' GAM w/ weights suggested by Dave
m.gam3 <- gam(data = df_comp1, Camera.mean ~ s(Aerial.mean), family = "gaussian", method = "REML", weights = wt1)
# gam.check(m.gam1)
x<-seq(0,0.75,0.01)
predmgam1 <- predict(m.gam1, newdata = data.frame(Aerial.mean = x), se.fit = TRUE)
predline1 <- data.frame(x, predmgam1$fit)
# plot.gam1 <- plot(m.gam1, shade = TRUE, shade.col = "lightblue", shift = coef(m.gam1)[1], ylab = "Camera.mean")
# plot.gam2 <- plot(m.gam2, shade = TRUE, shade.col = "lightblue", shift = coef(m.gam2)[1], ylab = "Camera.mean")
# Linear model
m.lm <- lm(Camera.mean ~ Aerial.mean -1, data = df_comp1, weights = wt1)
tidymlm <- tidy(m.lm)
# summary(m.lm)
tidymlm <- tidymlm %>%
mutate(inv.c = 1 / estimate,
inv.c.lci = 1 / (estimate + 1.65 * std.error),
inv.c.uci = 1 / (estimate - 1.65 * std.error))Results
Calibrating camera estimates to aerial survey estimates for moose
There is a general positive relationship between camera estimates and aerial-survey estimates of moose across WMU’s, but wide scatter as densities increase (Figure 1).
The very wide confidence intervals on the GAM included the linear fit line. The linear models with and without intercepts were very similar (although of course the shape of the CI’s differs).
Because the models produced similar results and the linear model without intercept is the simplest for calibration, we use that model. The calibration coefficient, 1/slope of the linear model without intercept, is 0.406 (90% CI’s: 0.341-0.503). That means that the aerial estimate for moose density in a WMU is 0.406 times the camera estimate. Equivalently, the camera estimate is 2.46 times higher than the aerial estimate. If the aerial survey results are taken as the “gold standard” to which camera estimates need to be calibrated, then the error around the calibration value has to be added to any error in the camera estimate itself.
plot1 <- df_comp1 %>%
ggplot(mapping = aes(x = Aerial.mean, y = Camera.mean)) +
geom_smooth(method = "lm", se = FALSE, aes(color = "Linear"), size = 2) +
geom_line(data = predline1, mapping = aes(x = x, y = predmgam1.fit, color = "GAM"),
size = 2) +
geom_point(size = 3) +
theme_classic() +
labs(x = expression(Moose~Density~via~Aerial~(per~km^2)),
y = expression(Moose~Density~via~Camera~(per~km^2)),
color = "Model:",
title = "Comparison of moose density estimates between two monitoring methods",
subtitle = "Aerial ungulate surveys and remote cameras",
caption = "The linear model has an adjusted R-squared value of 0.68") +
scale_x_continuous(breaks = seq(0,0.8,0.1)) +
scale_y_continuous(breaks = seq(0,3,0.5)) +
theme(axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
legend.title = element_text(size = 13),
legend.text = element_text(size = 13),
legend.position = "bottom")
plot1
Figure 1
kable(tidymlm, digits = 2)| term | estimate | std.error | statistic | p.value | inv.c | inv.c.lci | inv.c.uci |
|---|---|---|---|---|---|---|---|
| Aerial.mean | 2.46 | 0.29 | 8.59 | 0 | 0.41 | 0.34 | 0.5 |
The substantial overestimation of moose densities by cameras (assuming the aerial densities are close to the true value) is expected. ABMI cameras are put in open micro-habitats so that vegetation doesn’t hide animals for at least 5m; moose prefer those open areas for foraging, particularly in summer. Additionally, moose are clearly attracted to the cameras themselves, often going right up to the camera to investigate it. This inflates densities estimates by increasing the time that moose spend in the camera’s field-of-view and also reducing the effective detection distance that shows what area the camera is surveying. Species less associated with open vegetation and less curious about the camera would probably show less inflated density estimates. Aerial surveys also have assumptions that may not be met, which probably includes perfect detectability right under the plane.
Contribution of cameras to a combined WMU estimate
When there are both cameras and aerial surveys in a WMU, one option is to use a combined estimate. The camera estimate would first be calibrated to the aerial survey estimate using the above calibration coefficient and its uncertainty, then combined with the aerial survey estimate using standard inverse-variance weighting. We looked at how different numbers of cameras per WMU would contribute to a combined density estimate.
The first step was using the current camera results to show how the error associated with the camera estimate decreases with increasing number of cameras in the WMU (Figure 2). We fit a GAM smoothing spline to the relationship between the SE (expressed as a proportion of the mean) and the square root of the sample size in the WMU’s, using a log-normal error distribution. The square-root transformation was used based on statistical principles (SE = 1/sqrt(n)), and the log-normal distribution to account for the necessarily positive value of SE’s. The fitted relationship is shown in Figure 2.
df_comp1 <- df_comp1 %>%
mutate(se_prop_mean = Camera.se / Camera.mean) %>% # SE expressed as a proportion of the mean
mutate(se_prop_mean = na_if(se_prop_mean, "NaN"))
mm.gam <- gam(se_prop_mean ~ s(I(sqrt(n)), k = 4, m = 2), family = gaussian(link = "log"), data = df_comp1)
# summary(mm.gam)
# plot(mm.gam)
pred2 <- predict(mm.gam, newdata = data.frame(n = 0:200), se.fit = TRUE)
pred2line <- data.frame(n = 0:200, fit = exp(pred2$fit))
pred2line_l <- data.frame(n = 0:200, fit.l = exp(pred2$fit - 1.65*pred2$se.fit))
pred2line_u <- data.frame(n = 0:200, fit.u = exp(pred2$fit + 1.65*pred2$se.fit))
plot2 <- df_comp1 %>%
ggplot(mapping = aes(x = n, y = se_prop_mean)) +
geom_point(size = 3) +
geom_line(data = pred2line, mapping = aes(x = n, y = fit),
color = "blue", size = 2) +
geom_line(data = pred2line_u, mapping = aes(x = n, y = fit.u),
color = "blue", linetype = 2) +
geom_line(data = pred2line_l, mapping = aes(x = n, y = fit.l),
color = "blue", linetype = 2) +
labs(x = "Cameras per WMU",
y = "SE/mean for camera density",
title = "Relationship between number of cameras and error associated with the estimate",
subtitle = "Error expressed as SE as a proportion of the mean density estimate",
caption = "Solid line is fitted relationship, with dotted 90% confidence intervals.") +
scale_x_continuous(breaks = seq(0, 200, 50)) +
scale_y_continuous(breaks = seq(0, 0.8, 0.1)) +
theme_classic() +
theme(axis.title = element_text(size = 14),
axis.text = element_text(size = 12))
plot2
Figure 2
We then use Monte Carlo simulations with a range from 3 to 1000 cameras per WMU. Density estimates from that sample of cameras were generated from a log-normal distribution with the SE predicted from the relationship in Figure 2. The mean camera estimate for the WMU was then converted to the equivalent aerial survey density using the calibration coefficient, and its uncertainty, reported above. The aerial surveys have a limited range of SE’s based on the available confidence intervals, all being fairly near to the average of 0.161 * the mean value for the WMU. Using the Monte-Carlo-simulated SE of the camera estimate and the 0.161 SE for the aerial surveys, we calculated the weight of the camera estimate in an inverse-variance-weighted combined estimate as 1/SEcamera\(^2\) / (1/SEcamera\(^2\) + 1/SEaerial\(^2\)).
Even with 1000 cameras in a WMU, the cameras were only expected to contribute 0.406 (about 40%) of the weight of the combined estimate. With small numbers of cameras, their contribution was negligible. Given the assumptions of this analysis, and the confidence intervals we were given for the aerial surveys, feasible numbers of cameras would make no real improvements to density estimates of WMU’s that had aerial surveys.
Alternative statistical methods of estimating confidence intervals on aerial surveys might give different results. There are also assumptions to the aerial surveys that may not be met, which probably include perfect detectability right under the plane and a particular pattern of detectability decay with distance. Cameras might make a more valuable contribution if there was less certainty in the aerial results, or questions about how well their assumptions are met in different WMUs.
# Monte Carlo Simulations
iter <- 100000
n <- c(3,5,10,30,50,100,300,500,1000) # Cameras in WMU
y1 <- exp(predict(mm.gam, newdata = data.frame(n=n))) # Expected SE/mean at each level of n
c1 <- 1/rnorm(iter,tidymlm$estimate,tidymlm$std.error) # Distribution of inverse calibration coefficients
sda<-NULL
for (i in 1:length(n)) {
# Distribution of expected means with sample size n (and grand mean = 1), calculated on log scale (by delta method, se of log(x) = se of x when mean = 1)
dm<-exp(rnorm(iter,0,y1[i]))
dm<-dm/mean(dm) # Correct to mean=1 for log effect
sda[i]<-sd(dm*c1)/mean(dm*c1) # And those density estimates calibrated to aerial - SE/mean
}
ci.range<-(df_comp1$Aerial.uci-df_comp1$Aerial.lci)/df_comp1$Aerial.mean # Aerial survey CI range as a percentage of the mean
sd.aerial<-ci.range/2/1.65 # Effective SE as a proportion of the mean, assuming 90% CI's
sd.aerial.mean<-mean(sd.aerial,na.rm=TRUE)
p.camera <- (1 / sda^2) / ((1 / sda^2) + (1 / sd.aerial.mean^2))
q <- data.frame(Cameras = n, sd.Cameras = sda, sd.Aerial = sd.aerial.mean, Cont.Cameras = p.camera)
q <- q %>%
mutate_at(1:4, round, 3)
kable(q)| Cameras | sd.Cameras | sd.Aerial | Cont.Cameras |
|---|---|---|---|
| 3 | 0.647 | 0.161 | 0.058 |
| 5 | 0.628 | 0.161 | 0.062 |
| 10 | 0.601 | 0.161 | 0.067 |
| 30 | 0.534 | 0.161 | 0.083 |
| 50 | 0.497 | 0.161 | 0.095 |
| 100 | 0.433 | 0.161 | 0.122 |
| 300 | 0.318 | 0.161 | 0.204 |
| 500 | 0.262 | 0.161 | 0.274 |
| 1000 | 0.196 | 0.161 | 0.404 |