R session_REM
Paulina Pawlikowska, Dr Antonio Uzal
2024-11-15
In this session you will learn how to calculate the density of animals using Random Encounter Model (Rowcliffe et al., 2008). In this exercise we will calculate the density of hedgehogs in urban and rural habitat using data from Schaus et al. (2020). The code was also adapted from this publication.
Create environment
Set working directory and download the necessary packages
There are a few packages available for calculating density with REM. Today we will be using remBoot (Caravaggi, 2017).
#install.packages("devtools")
#devtools::install_github("arcaravaggi/remBoot")
#install.packages("activity")
#install.packages("dplyr")
library(remBoot)
library(readxl)
library(dplyr)Let’s see what format of the data is required for the REM analysis
#?remBootAs we can see in the Help section, to be able to run the model using remBoot, the data has to be structured in a certain way:
each row will represent one camera trap
there will be 4 columns: survey location, trapping rate (number of individuals captured), radial* distance to the detected animal (in km), and angle to the detected animal (in radians**)
*radial-> meaning that it’s a radius of the circular space around the camera trap
**radian->unit of angular measure, a radian is approx. 57 degrees
Additional variables required by the model are:
tm- The total number of hours all cameras were left in-situ at a focal site
v- The distance traveled by the focal species in 24 hours, in kilometers
Let’s upload the data and change it to a data frame
camDat <- read_excel("camDat.xlsx")
camDat<-as.data.frame(camDat)
head(camDat)## Site count dist theta
## 1 Southwell 1 0.0023500 0.6981317
## 2 Southwell 1 0.0015200 0.6283185
## 3 Southwell 3 0.0005800 0.7446738
## 4 Southwell 4 0.0008075 0.7243116
## 5 Southwell 2 0.0009950 0.4886922
## 6 Southwell 1 0.0007000 0.3839724When you take a look at the camDat dataset, you will notice that there are 2 site locations- urban (Southwell) and rural (Brackenhurst). We will be calculating the density for each location separately, so we need to split them into two separate groups. We first assign them numerical values
camDat$Site[camDat$Site=="Southwell"] <- 1
camDat$Site[camDat$Site=="Brackenhurst2017"] <- 2Make sure that your data has the correct format. Here we will change Site and count to an integer
camDat$Site<-as.integer(camDat$Site)
camDat$count<-as.integer(camDat$count)and split by survey
grpDat <- split_dat(camDat)Day range
One of the variables we need to calculate is day range (v), so distance traveled by an animal in a day (24h) in kilometers. Day range is a product of animal speed and time spent active.
Speed
Speed is extracted from camera trap footage. It can be calculated by measuring the distance that the animal traveled in a sequence of pictures that camera trap had captured and dividing it by the time that the animal was present in the camera trap’s detection zone. More automated approaches are being developed, where it will possible to extract the speed from images captured by cameras that were previously calibrated- see e.g. Miles et al., 2024.
In our example the mean speed was previously calculated:
0.77km/h for Southwell
0.50km/h for Brackenhurst
Activity
To calculate the time spent active, we will use package activity
library(activity)We will calculate the activity from camera trap footage. We will use the time space between the earliest possible recording of a hedgehog and the latest one as bounds and fit activity model to it. To identify the earliest and latest recording of hedgehogs in Southwell, we we will need a dataset that includes the time of each capture.
cap_time1 <- read_excel("activity.xlsx")
head(cap_time1)## # A tibble: 6 × 3
## Site Date Time
## <chr> <dttm> <dttm>
## 1 Southwell 2016-05-23 00:00:00 1900-01-01 23:48:00
## 2 Southwell 2016-05-23 00:00:00 1900-01-01 21:32:00
## 3 Southwell 2016-05-26 00:00:00 1900-01-01 23:45:00
## 4 Southwell 2016-05-27 00:00:00 1900-01-01 00:00:00
## 5 Southwell 2016-05-29 00:00:00 1900-01-01 00:13:00
## 6 Southwell 2016-05-24 00:00:00 1900-01-01 22:04:00We will be using a function fitact from activity package. Let’s take a look at what are the requirements of the function:
#?fitactThe function requires that time of capture is in a format of radian time of day.
Radians in that context are a way of expressing time. Imagine an analog clock, however instead of 12 hours, it represents the whole day (24 hours). Each hour can be represented as an angle between that hour and the starting point- in hour case 00:00.
However, angles expressed in degrees may cause problems in mathematical operations, therefore it is better to express this in radians. To make a full cycle on our “Radian clock” would be equal to 2pi. Half of the clock, so 12:00 would be represented as pi, because we divide 2pi/ 2 (because it’s half of the whole circle). Through this operation we can express each time of day in that form.
This will help us in analysing activity of animals, because it lets us express time in a circular form. We are not interested in the timeline of hedgehog captures, for example that a hedgehog was seen on the 5th of July at 23:44 and then on the 8th of July at 2:44. We want to pool together all hedgehog captures from one period and see if there were significantly more records around certain hours of day, and that is what expressing time in radians allows us to do.
rad_time <- gettime(cap_time1$Time, tryFormats = "%Y-%m-%d %H:%M:%S", scale = "radian")The gettime function converts the hours in POSIXct or character format to radians, which will be recognized by R.
To plot activity model we will have to choose the bounds of activity of hedgehogs. This is because they are nocturnal species and we will only be including night time in the analysis. To choose the bounds, we will need to identify the earliest and the latest time of capturing a hedgehog. As they are nocturnal species, the beginning of their activity will be the start of the night. Therefore, the first video after noon will be the beginning of their activity and the last before noon will be the end of their activity.
We need to create an element defining noon time. Remember that we are now operating in radians- the full cycle on hour radian clock is 2*pi; therefore noon is 2*pi divided by 2, so =pi.
#identifying the first video
noon_time <- pi
#filter the times after noon
times_after_noon <- rad_time[rad_time > noon_time]
#choose the lowest- beginning of activity
first_video <- min(times_after_noon)
first_video## [1] 5.598144#identifying the last video
#filter the times before noon
times_before_noon <- rad_time[rad_time < noon_time]
#choose the highest
last_video <- max(times_before_noon)
last_video## [1] 1.125737The first video is 5.598144, and the last video is 1.125737. We will use astroFns package and rad2hms function to identify the time of the first and the last video in the time format.
#install.packages("astroFns")
library(astroFns)
first<-rad2hms(rad = 5.598144 ,places = 0)
last<-rad2hms(rad = 1.125737 ,places = 0)
first## [1] "21:23:00"last## [1] "04:18:00"The first video was taken at 21:23 and the last video was taken at 4:18. We will use that information later.
Now we can fit activity model, using fitact function and choosing bounds as the argument
act1<-fitact(rad_time,bounds=c(first_video,last_video), sample = "data")act1Results of the model:
act=0.7879
se= 0.0726
To interpret this result we should multiply the activity rate (act=0.79) by the time between the two bounds.
Let’s recall the time for the first and last video of hedgehog recording
first 21:23
last 4:18
The spread between these two hours is approx. 7 hours. Therefore 0.79 multiplied by 7h will give 5,53, meaning that hedgehogs are active on average 5,53 hours of the night.
We then can plot it:
plot(act1,centre="night")
In this figure you can see the activity of hedgehogs. We can see that their activity starts after dawn and stops before dusk. We can also observe fluctuations in their activity levels with the highest peak after midnight.
Now we finally can calculate the day range. To do this we multiply the average speed by activity level:
s_Southwell<- 0.77 #speed for Southwell in km/h
v <- 0.79 * s_Southwell
v## [1] 0.6083Results are in km/day, so according to our calculations, hedgehogs in this site travel 0.61 km/day on average.
Trapping effort
The last element necessary for the density calculation is the trapping effort (tm), so the total number of hours that cameras were working. To do this we will sum all days, on which the cameras were active and multiply them by the number of hours, when they were recording. As hedgehogs were captured only at night, again we will use the bounds between the first and the last capture of hedgehogs.
Summing all days
We will need a dataset containing information on camera deployments.
deployment_cam <- read_excel("Southwell_cam_deployment.xlsx")
head(deployment_cam)## # A tibble: 6 × 6
## Habitat Site Camera_ID Nights_active X Y
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 Urban Southwell2016 S102 5 469007 353782
## 2 Urban Southwell2016 S103 4 469218 353874
## 3 Urban Southwell2016 S104 6 469631 353355
## 4 Urban Southwell2016 S105 8 469729 353988
## 5 Urban Southwell2016 S106 5 469919 353727
## 6 Urban Southwell2016 S107 8 469639 353781The file contains a column “Nights_active” with a number of nights that each camera was active
cam_days <- sum(as.numeric(deployment_cam$Nights_active))
cam_days## [1] 746Cameras were active for 746 days in total. We now multiply the number of camera days by number of hours they were active (7 hours).
tm<- cam_days*7
tm## [1] 5222Our trapping effort is 5222 hours.
Density estimation
Now we have all the metrics to calculate density using REM.
We will use function rem. We define our arguments tm and v
rem(grpDat[[1]], tm = 5222, v = 0.61)## [1] 29.16472And now standard deviation using remBoot function.
output <- remBoot(grpDat[[1]], tm = 5222, v = 0.61,nboots = 1000, error_stat = c("sd"))
output## $sd
## $sd[[1]]
## [1] 3.777056The density of hedgehogs in Southwell is approx. 29 hedgehogs/ km2.
Exercise 1
Try to calculate density for Brackenhurst 2017 using the same method.
Remember that you already have data frame for rem model at Brack saved as:
grpDat[[2]]## Site count dist theta
## 37 2 1 0.002410000 0.3141593
## 38 2 2 0.000900000 0.4537856
## 39 2 1 0.002200000 0.6632251
## 40 2 3 0.001490000 0.5235988
## 41 2 1 0.000800000 0.2094395
## 42 2 1 0.000920000 0.0000000
## 43 2 1 0.000600000 0.4188790
## 44 2 1 0.002200000 0.1745329
## 45 2 1 0.000890000 0.7330383
## 46 2 9 0.001821111 0.4033650Activity dataset for Brackenhurst is available as:
acti_Brack <- read_excel("activity_Brack.xlsx")and camera deployment file is called:
deployment_Brack <- read_excel("Brack_cam_deployment.xlsx")In which habitat the density of hedgehogs was higher? Is it consistent with other available data on hedgehogs?
Additional possibilities of the model
Now we will look at other possibilities of calculating the variables of the model. As the interest in REM is growing, new functions and packages are being developed. Now we will look at functions described in Palencia et al., (2022). We will be using data adapted from this publication as well.
In this exercise we will assume that the probability of detecting an animal in the camera trap’s field of view is uniform across the whole area. In reality, the probability of detecting an animal decreases with distance and with the angle of detection, so for future applications of the model, consider fitting appropriate models for both distance and angle (see more: Palencia et al., 2022)
Please note, this exercise uses only a part of the code discussed in the publication and a modified dataset. If you are interested in the full analysis and results, please take a look at the original publication and code description available at this link:
#link to detailed code: https://github.com/PabloPalencia/CameraTrappingAnalysis/blob/main/REM/REM%20_vignette.pdf#Install the necessary packages:
#devtools::install_github("PabloPalencia/trappingmotion")
#And load them
library(trappingmotion)
library(dplyr)
library(readxl)The species we will look at is the Iberian hare Lepus granatensis.
#Load the data
iberian_dat <- read_excel("Iberian_hare_all.xlsx")
#Take a look at the data structure
glimpse(iberian_dat)## Rows: 90
## Columns: 8
## $ Point_ID <dbl> 102, 102, 102, 102, 102, 102, 102, 102, 103, 103, 104, 10…
## $ Sp <chr> "iberianhare", "iberianhare", "iberianhare", "iberianhare…
## $ Date <dttm> 2020-04-28, 2020-04-28, 2020-04-29, 2020-04-30, 2020-05-…
## $ H_first <dttm> 1899-12-31 00:25:41, 1899-12-31 00:27:38, 1899-12-31 21:…
## $ Speed.m.s <chr> "0.17499999999999999", "0.9", "0.74285714285714299", "1.4…
## $ Interval.min <dbl> 2, 3, 3, 3, 2, 3, 2, 1, 3, 1, 1, 4, 2, 2, 3, 2, 3, 2, 2, …
## $ Dist_det <chr> "6.4999999999999997E-3", "NA", "6.0999999999999995E-3", "…
## $ Ang_det <chr> "0.24434609499999999", "0.38397243499999995", "0.94247779…#We will first focus on the speed estimation, so the column "Speed.m.s"
#Change data to numeric
iberian_dat$Speed.m.s<- as.numeric(iberian_dat$Speed.m.s)
#check for NA's
is.na(iberian_dat)
#NA's are present
#removing NA's
iberian_dat<-na.omit(iberian_dat)#check
which(is.na(iberian_dat$Speed.m.s))## integer(0)#should be 0As we can see in the dataset, speeds for each sequence are provided in meters per second.
Speed estimation
First we will calculate the average speed of iberian hares using meanspeed function. The function calculates the harmonic mean of provided speeds. In this function, the authors give you the option to account for the behaviors of animals in camera trap’s field of view. First different speeds are clustered using identbhvs function and the output is used as an argument in meanspeed function. The speed measurements should be in m/s.
identbhvs(iberian_dat$Speed.m.s)
Four different speed clusters were identified. Now we calculate speed for each of them.
speed_data<- meanspeed()
speed_data## speeds speed_se n_seq
## 1 2.127993308 0.107416691 5
## 2 0.003338448 0.002809752 41
## 3 0.519895515 0.043938325 12
## 4 1.273111983 0.050503985 12The output is a table consisting of mean speed for each behavioral state, speed standard error and number of sequences withing each behavioral category.
Day range
The package also includes a dedicated function for calculating day range, called dayrange. It requires as arguments: activity rate (values between 0 and 1), activity standard error, and the output of the meanspeed function.
Let’s start from calculating the activity. As iberian hares are not nocturnal, we will not be calculating the activity within bounds, but throughout the whole day.
#load the package
library(activity)
#change to data frame
iberian_acti <- as.data.frame(iberian_dat)
# Change time to radian time of day as in the previous example
rad_time_ib <- gettime(iberian_acti$H_first, tryFormats = "%Y-%m-%d %H:%M:%S", scale = "radian")
#fit activity model, set bounds as NULL
ib_act<-fitact(rad_time_ib,bounds= NULL, sample = "data")
ib_act
#plot activity pattern
plot(ib_act,centre="night")
The model calculated the activity rate as 0.3414 (SE=0.03974), which can be interpreted as 0.3414*24h=8,19h, so hares in that site are active approx. 8h a day.
Now we can calculate the day range. The output of this function is in km/day as in our previous example.
dayrange(act=ib_act@act[1], act_se=ib_act@act[2], speed_data)## Day range (Km/day) 0.1676522
## Day range SE (Km/day) 0.08469679Now the only missing variables for REM calculation will be trapping effort and the number of captures.
Trapping effort
Let’s start with the trapping effort.
#Load the data
trapping_ib <- read_excel("Trapping_iberian.xlsx")
trapping_ib <- as.data.frame(trapping_ib)This time, the data has a bit different structure than in our previous example.
head(trapping_ib)## CAM 43948 43949 43950 43951 43952 43953 43954 43955 43956 43957 43958 43959
## 1 101 0 1 1 1 1 1 1 1 0 0 0 0
## 2 102 0 1 1 1 1 1 1 1 1 1 1 1
## 3 103 0 1 1 1 1 1 1 1 1 1 1 1
## 4 104 0 1 1 1 1 1 1 1 1 1 1 1
## 5 105 0 1 1 1 1 1 1 1 1 1 1 1
## 6 106 0 1 1 1 1 1 1 0 0 0 0 0
## 43960 43961 43962 43963 43964 43965 43966 43967 43968 43969 43970 43971 43972
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 1 1 1 1 1 1 1 0 0 0 0 0 0
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 1 1 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0
## 43973 43974 43975 43976 43977 43978 43979 43980 43981 43982 43983 43984 43985
## 1 1 0 0 0 0 0 0 0 0 0 0 0 0
## 2 1 1 1 1 1 1 1 1 1 1 1 1 0
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1 1 1 1 1 1 1
## 6 1 1 1 1 1 1 1 1 1 1 0 0 0
## 43986 43987 43988 43989 43990 43991 43992 43993 43994 43995 43996 43997 43998
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0 0 0 0
## 3 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1 1 1 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0
## 43999 44000 44001 44002 44003 44004 44005 44006 44007 44008 44009 44010 44011
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0 0 0 0
## 3 1 1 1 1 1 1 1 1 0 0 0 0 0
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0
## 44012 44013 44014 44015 44016 44017 44018 44019 44020 44021 44022 44023 44024
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0
## 44025 44026 44027 44028 44029 44030 44031 44032 44033 44034 44035 44036 44037
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 1 1 1 1 1 1 1 1 1 1 1 1 1
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0
## 44038 44039 44040
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 1 1 1
## 5 0 0 0
## 6 0 0 0Each row represents the operativity of one camera trap, it gives 1 when the camera trap was working properly on a given day and 0 when it wasn’t.
To calculate the survey effort we will first summarize values from all the rows
#counting the number of trapping days- summarizing all the columns, except from the CAM column
trap_days <- sum(rowSums(trapping_ib, na.rm = FALSE, dims = 1))- sum(trapping_ib$CAM)
trap_days## [1] 821#821 trap daysNow we will calculate the trapping hours. We were not using bounds of the activity this time, as we did with hedgehogs, so this time we will be multiplying camera days by 24 hours.
#counting trap hours
t_hours<-trap_days*24
t_hours## [1] 19704#trapping hours= 19704hOur trapping effort is 19704 hours.
Remember that for the rem function we need a dataframe containing specific variables. We will need a column with the number of captures of hares in each camera. In iberian_dat we have one row for each capture of a hare. Therefore there can be multiple rows for one camera (Point_ID). To get the frequency of captures of each camera we will use the table function to create a new element with that data.
#adding a column with a number of days that each camera was operating (except CAM column)
trapping_ib$enc_days<-rowSums(trapping_ib, dims = 1)-trapping_ib$CAM
#creating new element with data on sequences and frequencies of each
seq_freq<-data.frame(table(iberian_dat$Point_ID))
head(seq_freq)## Var1 Freq
## 1 102 6
## 2 104 1
## 3 105 3
## 4 106 1
## 5 107 1
## 6 108 2And now we will merge the data about the frequency to our operativity matrix. For merging we will use the CAM column in our matrix and Var1 in our new element.
#add a column to the trapping data frame, connecting information about how many captures each camera had;assign NA's as 0
trapping_ib <- merge(trapping_ib, seq_freq, by.x= "CAM", by.y = "Var1", all.x = TRUE); trapping_ib[is.na(trapping_ib)] <- 0Now we will have to calculate the mean speed and mean angle for each camera.
iberian_dat<-as.data.frame(iberian_dat)
iberian_dat$Dist_det<-as.numeric(iberian_dat$Dist_det)
iberian_dat$Ang_det<-as.numeric(iberian_dat$Ang_det)
dist_by_cam<-iberian_dat %>% group_by(Point_ID) %>% summarise(mean_dist=mean(Dist_det, na.rm=TRUE))
head(dist_by_cam)## # A tibble: 6 × 2
## Point_ID mean_dist
## <dbl> <dbl>
## 1 102 0.00544
## 2 104 0.0012
## 3 105 0.00513
## 4 106 0.0041
## 5 107 0.0035
## 6 108 0.0061ang_by_cam<-iberian_dat %>% group_by(Point_ID) %>% summarise(mean_ang=mean(Ang_det, na.rm=TRUE))
head(ang_by_cam)## # A tibble: 6 × 2
## Point_ID mean_ang
## <dbl> <dbl>
## 1 102 0.529
## 2 104 0.593
## 3 105 0.396
## 4 106 0.175
## 5 107 0.105
## 6 108 0.663In trapping_ib we have some rows with 0, meaning there were 0 captures of hares in that camera. We will remove the rows now, as we don’t need them in our data frame for the density calculation (note that the cameras were still included in the trapping effort).
#remove the zeros
trapping_clean <- trapping_ib[(!trapping_ib$Freq==0), ]As we will later be merging a few datasets through the Point_ID column, we have to create a column with such a name in the trapping_clean dataset
trapping_clean<- dplyr::rename(trapping_clean, Point_ID=CAM)Now from the trapping_clean data frame we will create a new element with the columns that we will use in the subsequent analysis.
#creating a new element with Freq and Point_ID
simple_freq<-trapping_clean[,c("Freq","Point_ID")]
head(simple_freq)## Freq Point_ID
## 2 6 102
## 4 1 104
## 5 3 105
## 6 1 106
## 7 1 107
## 8 2 108And now we merge all the columns of interest into one dataset.
#joining captures with mean distance and mean angle by Point_ID
rem_dats <- simple_freq %>%
left_join(dist_by_cam, by = "Point_ID") %>%
left_join(ang_by_cam, by = "Point_ID")
head(rem_dats)## Freq Point_ID mean_dist mean_ang
## 1 6 102 0.005440000 0.5294165
## 2 1 104 0.001200000 0.5934119
## 3 3 105 0.005133333 0.3956080
## 4 1 106 0.004100000 0.1745329
## 5 1 107 0.003500000 0.1047198
## 6 2 108 0.006100000 0.6632251We will not need the information on Point_ID anymore, so let’s change the values to 1 for each row, as required by remBoot.
rem_dats <- rem_dats %>%
mutate(Point_ID = 1)
head(rem_dats)## Freq Point_ID mean_dist mean_ang
## 1 6 1 0.005440000 0.5294165
## 2 1 1 0.001200000 0.5934119
## 3 3 1 0.005133333 0.3956080
## 4 1 1 0.004100000 0.1745329
## 5 1 1 0.003500000 0.1047198
## 6 2 1 0.006100000 0.6632251#change to numeric
rem_dats$Point_ID<-as.numeric(rem_dats$Point_ID)As we saw in our previous analysis, rem function requires certain organisation of data. Let’s reorder the columns as required by the function.
#reorder columns as required by rem function
rem_dats <- rem_dats %>%
select(Point_ID,Freq,mean_dist,mean_ang)Density estimation
Now we can calculate the density.
#recall the other values that we calculated earlier
t_hours<-19704
v<-0.16
#calculate the density
rem(rem_dats,t_hours,v)## [1] 6.503561#and standard deviation
density_hare <- remBoot(rem_dats, t_hours, v,nboots = 1000, error_stat = c("sd"))
density_hare## $sd
## $sd[[1]]
## [1] 1.43176References:
Caravaggi, A. (2017). remBoot: an R package for random encounter modelling. Journal of Open Source Software, 2(10), 176.
Miles, V., Woodroffe, R., Donnelly, C. A., Brotherton, P. N., Ham, C., Astley, K., … & Rowcliffe, M. (2024). Evaluating camera‐based methods for estimating badger (Meles meles) density: Implications for wildlife management. Ecological Solutions and Evidence, 5(3), e12378.
Palencia, P., Barroso, P., Vicente, J., Hofmeester, T. R., Ferreres, J., & Acevedo, P. (2022). Random encounter model is a reliable method for estimating population density of multiple species using camera traps. Remote Sensing in Ecology and Conservation, 8(5), 670-682.
Schaus, J., Uzal, A., Gentle, L. K., Baker, P. J., Bearman‐Brown, L., Bullion, S., … & Yarnell, R. W. (2020). Application of the Random Encounter Model in citizen science projects to monitor animal densities. Remote Sensing in Ecology and Conservation, 6(4), 514-528.
https://github.com/PabloPalencia/CameraTrappingAnalysis/blob/main/REM/REM%20_vignette.pdf