In timcdlucas/RandEM: Population density estimates via Random Encounter Models in R

Read the packaged data file. These data are from a study of hares in Northern Ireland, published as Caravaggi et al. (2016). The first 5 columns of any data set must follow the same format in terms of column contents.

library(RandEM)

data(hDat)

Alternatively, generate relevant simulated data.

remDat <- data.frame(cbind(site <- sample(1:4, 200, replace=TRUE),
                           cam <- sample(1:20, 200, replace = TRUE),
                           count <- sample(1:3, 200, replace = TRUE),
                           dist <- runif(200, min=0, max=0.006),
                           theta <- runif(200, min=0, max=0.6)))

colnames(remDat) <- c("site","cam", "count", "dist", "theta")

Now we'll check the first few rows of data, check the sample size for each site, and check data consistency.

head(hDat) 
tapply(hDat[,3], hDat[,1], length) 
sum(tapply(hDat[,3], hDat[,1], length)) == nrow(hDat)

We'll use two functions to split the data according to survey site (split_dat) and store the number of rows. These data will be used for iterative resampling during bootstrapping.

grpDat <- split_dat(hDat)

Check that the code has done what we wanted

lapply(grpDat, head)

lapply(grpDat, nrow)

The REM requires that the user provide survey effort (i.e. camera hours; tm) and the distance travelled by the focal species in 24 hours (given in km; v). Here, we assign values to tm and v.

tm <- 1880
v <- 0.89

We then use the REM function (rem) to calculate densities for each site. Change the grpDat value to specify a different site.

rem(dat = grpDat[[1]], tm, v)

If tm and v differ for each survey site, we can specify them alongside the REM function, as below. Note that if the focal species is a constant, v should not change.

rem(dat = grpDat[[1]], tm = 3600, v = 1.4) 
rem(dat = grpDat[[2]], tm = 3360, v = 1.4)

If there is only one survey site, we can calulate its density, by specifying the data frame rather than the group. Note that the two examples below differ as one calls the stored tm and v, while the other defines them.

rem(hDat, tm, v) 

rem(hDat, tm = 3360, v = 1.4)

Now we'll calculate variance for each study site. First, define the number of bootstrapping iterations:

nboots <- 1000

Now we use the bootstrapping function boot_sd on each group dataframe n (i.e.nboots) times. Note that tm and v must be the same across all sites. If this is not the case, split your data into separate dataframes before bootstrapping and run the function on each group individually.

remsD <- lapply(grpDat, boot_sd, tm, v, nboots)

Now calculate the standard deviation for each site.

remsSD <- lapply(remsD, sd)
remsSD

Alternatively, we can forgo all the above and use the remBoot function to split and group data and calculate variance (SD and/or 95% confidence intervals) in one go. remBoot requires that all sites have the same values for tm and v. As before, we'll define tm, v and nboots a priori:

tm <- 1880
v <- 0.89
nboots <- 1000

output <- remBoot(hDat, tm, v, nboots, error_stat = c("sd"))
output

Alternatively, we could pass the values to the function directly.

output <- remBoot(hDat, tm = 2870, v = 1.8, nboots, error_stat = c("ci"))
output

We can then collate our data into a single data frame and generate a simple bar plot with error bars to lok at relative differences between sites.

require(ggplot2)

site <- seq(1:5)
den <- seq(1:5)
estDat <- data.frame(site, den)
estDat[1,2] <- rem(dat = grpDat[[1]], tm, v) 
estDat[2,2] <- rem(dat = grpDat[[2]], tm, v) 
estDat[3,2] <- rem(dat = grpDat[[3]], tm, v)
estDat[4,2] <- rem(dat = grpDat[[4]], tm, v)
estDat[5,2] <- rem(dat = grpDat[[5]], tm, v)
estDat$sd <- as.numeric(remsSD)

ggplot(estDat, aes(x = site, y = estDat$den)) +  
  geom_bar(position = position_dodge(), stat="identity", fill="light blue") + 
  geom_errorbar(aes(ymin=den-sd, ymax=den+sd)) +
  scale_x_continuous(name="Survey site") +
  scale_y_continuous(name="Animals.km2") +
  ggtitle("Animal density estimates calculated with REM and remBoot")