CITE THIS NOTEBOOK: Estimating Species Richness with Camera Traps by Andrew Mashintonio, Wolfram Community, OCT 4 2022.
ORIGINAL ARTICLE: Andrew F. Mashintonio, Grant M. Harris, David R. Stewart, Matthew J. Butler, Jim Sanderson, Gareth Russell, (2022). Estimating Species Richness with Camera Traps: Modeling the Effects of Delay Period, Deployment Length, Number of Sites, and Interference Imagery. Wildlife Society Bulletin. https://doi.org/10.1002/wsb.1357
CITE THIS NOTEBOOK: Estimating Species Richness with Camera Traps by Andrew Mashintonio, Wolfram Community, OCT 4 2022.
ORIGINAL ARTICLE: Andrew F. Mashintonio, Grant M. Harris, David R. Stewart, Matthew J. Butler, Jim Sanderson, Gareth Russell, (2022). Estimating Species Richness with Camera Traps: Modeling the Effects of Delay Period, Deployment Length, Number of Sites, and Interference Imagery. Wildlife Society Bulletin. https://doi.org/10.1002/wsb.1357
ORIGINAL ARTICLE: Andrew F. Mashintonio, Grant M. Harris, David R. Stewart, Matthew J. Butler, Jim Sanderson, Gareth Russell, (2022). Estimating Species Richness with Camera Traps: Modeling the Effects of Delay Period, Deployment Length, Number of Sites, and Interference Imagery. Wildlife Society Bulletin. https://doi.org/10.1002/wsb.1357
Introduction
Introduction
Camera traps used to estimate species richness are commonly set to have as small a delay period as possible to maximize the chances of observing each species, particularly rare ones. However, this leads to an overabundance of imagery to process and manage. The presence of empty images, or images taken without any target species in them, can further complicate the effect of the delay period setting on species capture. Increasing the delay period can reduce the number of images recorded while still capturing nearly all of the species present at a site, even when the number of empty images in a dataset is high.
This notebook first contains code for subsampling camera trap data at various delay periods to simulate the effects of using a longer delay on species richness estimates, the number of images recorded, and the probability of capturing a rare species. In addition, it provides code for simulating a camera trapping dataset with the inclusion of various levels of blank images (a.k.a. interference imagery). Lastly, it provides code for modeling the effect of persistent interference on species richness estimates under various other conditions such as return time and persistence time of the target organism. Each section below contains the code for the functions and sample outputs and models. It does not contain the data used to generate the outputs; that is available through Zenodo.
This notebook first contains code for subsampling camera trap data at various delay periods to simulate the effects of using a longer delay on species richness estimates, the number of images recorded, and the probability of capturing a rare species. In addition, it provides code for simulating a camera trapping dataset with the inclusion of various levels of blank images (a.k.a. interference imagery). Lastly, it provides code for modeling the effect of persistent interference on species richness estimates under various other conditions such as return time and persistence time of the target organism. Each section below contains the code for the functions and sample outputs and models. It does not contain the data used to generate the outputs; that is available through Zenodo.
Subsample Data
Subsample Data
Functions
Functions
The following four functions are called within the ‘accumulationData’ function, defined later, and do not need to be set up to run manually.
In[]:=
findLatestStart[trapData_,interval_,dateRange_]:=Module[{lastDay,lastStartTime},(*Determinethelastdayinthedataset*)lastDay=trapData[[All,All,3]][[-1]];(*Determinethelatestsubsamplingcanstarttoendbythelastday*)lastStartTime=lastDay-dateRange;(*Findtheimageinthedatasetthatisclosest(butnotpast)thelaststarttime*)SelectFirst[trapData,#[[1,3]]>lastStartTime[[1]]&][[1,3]]]
subSample[data_,delay_,dateRange_,startTime_]:=Module[{simpleTimes,start,end,numDropped,keepQ},(*Pullallofthetimestampsinthedataset*)simpleTimes=data[[All,1,3]];(*Statethestarttimeandgeneratetheendtimebasedonthedeploymentlength*)start=Round[AbsoluteTime[startTime]];end=start+dateRange;(*Selecttheimagesthatareafterthestarttime*)simpleTimes=Select[simpleTimes,#≥start&];(*Determinehowmanyimageswereexcludedpriortothestarttime*)numDropped=Length[data]-Length[simpleTimes];(*Selecttheimagesthatarepriortotheendtime*)simpleTimes=Select[simpleTimes,#<end&];(*Rescalethedata*)simpleTimes=If[Length[simpleTimes]>0,simpleTimes-simpleTimes[[1]],simpleTimes];(*Subsamplethedataaccordingtothedelayperiod*)keepQ=If[Length[simpleTimes]>0,subSampleC[Round[simpleTimes],Round[delay]],simpleTimes];Prepend[Drop[data,numDropped][[Flatten[Position[keepQ,1]]]],{{data[[1,1,1]],"Init",start,99}}]]
In[]:=
subSampleC=Compile[{{times,_Integer,1},{delay,_Integer}},Module[{n,keepQ,lastTrigger,nextTrigger},n=Length[times];(*Initializelistofsamplestokeep*)keepQ=Table[0,{n}];keepQ[[1]]=1;(*Initializecounters*)lastTrigger=1;nextTrigger=2;(*Dothefollowinguntilwehittheendofthelist*)While[nextTrigger≤n,(*Checktheinterval,andifitisgreaterthanthedelaykeepthatnextcaptureandresetittobethelastcapture*)If[(times[[nextTrigger]]-times[[lastTrigger]])≥delay,keepQ[[nextTrigger]]=1;lastTrigger=nextTrigger;nextTrigger++,(*...otherwiselookatthenextcapture*)nextTrigger++]];(*Returntheindexlistofsamplestokeep*)keepQ]];
In[]:=
accumulation[sampleData_,bycatch_]:=Module[{acc,spp,localAcc,localSpp,i,sample,j,candidate},(*Beginwithempylists,addtothemwhenimagesaresampled*)acc={{0,0,{}}};spp={};For[i=2,i≤Length[sampleData],i++,localAcc=0;localSpp={};sample=sampleData[[i]];For[j=1,j≤Length[sample],j++,candidate=sample[[j,2]];If[!(MemberQ[spp,candidate]||MemberQ[bycatch,candidate]),localAcc++;localSpp=Append[localSpp,candidate];spp=Append[spp,candidate];];];If[localAcc>0,acc=Append[acc,{QuantityMagnitude[DateDifference[sampleData[[1,1,3]],sampleData[[i,1,3]]]],acc[[-1,2]]+localAcc,localSpp}]];];acc]
This function takes as input the data to be subsampled, the number of sites (cameras) to subsample, the number of different combinations of sites to compare, a list of delay period intervals (e.g., {10, 60, 600}), the deployment length of the subsample, the sample size (n), and a list of bycatch objects to ignore. The output is a list of the number of target species for each sample at each delay period.
accumulationData[data_,sites_,comparisonN_,delayList_,deployment_,n_,bycatch_]:=Module[{traps,IDcodes,cameras,codes,codesCheck,accumulations,sampleData,years,yearsCommon,positions,lastStarts1,lastStarts2,lastStarts,lastStart,firstStarts,firstStart,startDates,startDate,data2,subSamplesCombined,accumulations1,sample,trap},traps=sites;If[traps>1,(*Ifmorethanonetrap,firstgeneratealistoftrapstosubsample*)IDcodes=If[Length[data]traps,{Range[traps]},(*Subsamplefromalltraps;otherwisegeneratealistofallcameras*)cameras=Range[Length[data]];(*Ensurethatthelengthoftimethecameraisinoperationisgreaterthanthedeploymentinputvalue;otherwise,removethecamerafromconsideration*)cameras=Table[If[DateDifference[data[[camera]][[1,1,3]],data[[camera]][[-1,1,3]]][[1]]>Max[deployment/3600/24],camera],{camera,cameras}]/.NullSequence[];If[Length[cameras]≥traps,(*Iftherearestillmorecamerasthanthenumbertosubsamplefrom,randomlyselectasetofcamerasequaltothesitesinput,andrepeatforthenumberofdifferentcomparisonstobemade*)codes=Union[Table[Sort[Union[RandomSample[cameras,traps]]],{comparisonN}]];(*Thebelowfunctioncheckstomakesurethatthecodesgeneratedabovearecorrect*)codesCheck[input_]:=Module[{output},output=Union[Sort[Join[input,{RandomSample[cameras,traps]}]]];output=Map[If[Length[#]traps,#]&,output]//.NullSequence[]];NestWhile[codesCheck,codes,Length[#]<comparisonN&,1,15000(*Settingamaxof15000loopsisacatchagainstthesitutaionwherethetotaltrapcombinationsislessthanthecomparisonN,inwhichcasethefunctionwouldrunforever.Thenumber1isaplaceholder*)],Print[Length[cameras]];(*Ifafterremovingsomecameras,thenumberofremainingcamerasisequaltothesamplesize:*)traps=Length[cameras];{cameras}]];(*Summarizethelistofencountersforeachtargetspeciesundereachsamplingscheme.*)accumulations=Table[(*Loopoverdifferentdelays*)Print["Delay = "<>ToString[j]<>" seconds"];ParallelTable[(*Loopoverdifferentcombinationsoftraps*)Print["Traps "<>ToString[IDcodes[[i]]]];sampleData=If[Length[data]traps,data,(*Usealldataorloopthroughdifferenttrapcombinations*)Flatten[Table[Take[data,{IDcodes[[i,c]]}],{c,traps}],1]];If[n>1,(*Iftherearemultiplesamples*)(*Findthelastdateatwhichthesamplingcanstartandincludethefulldeploymentlength*)lastStart=Min[Map[findLatestStart[#,Max[delayList],deployment]&,sampleData]];(*Findthefirststartdateforeachtrap*)firstStart=Max[Table[sampleData[[i,1,1,3]],{i,1,traps}]];(*Generatealistofdatestostartsubsamplesfrom*)startDates=DateRange[firstStart,lastStart,DateDifference[firstStart,lastStart]/(n-1)];(*Subsampledatafromthelistforeachofthestartdates*)subSamplesCombined=Table[Drop[SortBy[Flatten[Map[subSample[#,j,deployment,s]&,sampleData],1],#[[1,3]]&],(traps-1)],{s,startDates}];(*Accumulatealistoftargetspeciesforeachsubsample*)accumulations1=DeleteCases[Map[accumulation[#,bycatch]&,subSamplesCombined],{{0,0,{}}}];accumulations1,(*Ifthereisonlyonesubsample,dothesameasabove,butwithdifferentcode*)lastStart=Min[Map[findLatestStart[#,Max[delayList],deployment]&,sampleData]];firstStart=Max[Table[sampleData[[i,1,1,3]],{i,1,traps}]];startDates=DateRange[firstStart,lastStart,DateDifference[firstStart,lastStart]/100];startDate=RandomChoice[startDates];subSamplesCombined=Drop[SortBy[Flatten[Map[subSample[#,j,deployment,startDate]&,sampleData],1],#[[1,3]]&],(traps-1)];accumulations1=DeleteCases[accumulation[subSamplesCombined,bycatch],{{0,0,{}}}];accumulations1],{i,1,Length[IDcodes]}],{j,delayList}],(*Ifonly1cameralocation(site),n=1,dothesameasabovebutwithdifferentcode*)trap=RandomChoice[data];accumulations=(lastStart=findLatestStart[trap,Max[delayList],deployment];startDates=DateRange[0,lastStart,DateDifference[0,lastStart]/100];startDate=RandomChoice[startDates];Table[Print["Delay = "<>ToString[j]<>" seconds"];sample=Rest[subSample[trap,j,deployment,startDate]];DeleteCases[accumulation[sample,bycatch],{0,0,{}}],{j,delayList}])]]
Output
Output
Each study area is subsampled across varying number of sites, deployment length, and delay period. Each combination of these variables will be sampled n = 100 times. For each sample n, up to 25 different combinations of sites will be sampled, resulting in a total number of subsamples of up to 2500.
Set parameters.
Check cameras.
One camera missing for the 3-year subsample.
Species Richness
Species Richness
To compare richness estimates across study areas, we divided the richness found in each subsample by the total richness as measured by the full dataset (proportion of species richness, Rp).
Functions
Functions
This function determines the proportion of the total number of species for one subsample.
This function determines the proportion of the total number of species across the whole dataset.
Output
Output
First check to make sure that none of the species are only found at the site that was not sampled during the 3-year deployment.
Two fewer species.
Badger and jackrabbit are not found in any of the traps used during the 3-year subsampling period. Therefore, reduce total species richness for these subsamples.
Use the subsampled data to calculate the mean and SD of Rp for each combination of deployment length, number of sites, and delay.
Models
Models
There are four models (10s, 60s, 300s, 3600s) for each study area.
Number of Images
Number of Images
Functions
Functions
The function takes as input the data to be subsampled, the number of sites (cameras) to subsample, the number of different combinations of sites to compare, a list of delay period intervals (e.g., {10, 60, 600}), the deployment length of the subsample, the sample size (n), and a list of bycatch objects to ignore. The output is a list of the number of images for each sample at each delay period.
Output
Output
Set the delay periods:
Count the number of images across different subsamples for deployment length, number of sites, and delay.
Take the mean and SD of the data.
Models
Models
There are four models (10s, 60s, 300s, and 3600s) for each study area.
Check parameters of model and fit to data.
Subsampling Ratio
Subsampling Ratio
Model the proportion of subsamples containing a rare species (subsampling ratio) under different combinations of deployment length, number of sites, and delay.
Functions
Functions
This function takes as input the subsampled data and the target species. The output is the proportion of subsamples (under a given combination of deployment length, number of sites, and delay) that resulted in at least one capture of the target species.
Output
Output
Use the subsampled data to determine the proportion of subsamples that resulted in at least one capture of the target species under each combination of deployment length, number of sites, and delay.
Models
Models
Interference Imagery
Interference Imagery
This analysis will incorporate varying levels of interference imagery into simulated camera trap data. The levels are 0%, 5%, 10%, 25%, 50%, and 75% of all images. Once these datasets are generated, they will be analyzed to determine the effect of interference level on Rp under varying levels of deployment length, number of sites, and delay.
Functions
Functions
This function takes as input the total of the capture frequencies from the target species and the desired proportion of interference images. The output is the “capture” frequency of interference events necessary to achieve the desired proportion of interference images.
This function generates a trap history based on species encounter rates:
This function generates encounter times (in seconds) based on the desired distributions and number of images; the distributions in the function are based off of actual data.
This function takes as input the list of capture frequencies per species, total deployment length of the simulation (in days), the total number of sites, the desired interference weights (output from above), and the default delay period. The output is a simulated dataset formatted in the same way as the real data.
This function does the same as ‘propSpecies’, except it removes the interference images from the count.
Data Generation
Data Generation
Use the same set of 20 species and capture frequencies for all 1000 simulations.
First, set the capture frequencies of target species.
Add what type of daily distribution each should have. The rare species can have a random distribution, since they will not appear much in the dataset. For all other species, they should be given a normal distribution around noon or midnight.
The common species should all have yearly distributions around June. The moderate species can have any type of distribution, and rare species will be random (label doesn’t matter).
Finally, determine the number of sites each species should appear at.
Set interference levels. First determine the total capture frequencies of the target species for reference, then set the desired proportions of interference images. Use the function to generate the necessary capture frequency for interference images to obtain the desired proportions.
Set other parameters.
Run the simulations (Note: this is a computation-intensive step and may require saving parts of the output, clearing the storage, and repeating until all simulations are complete. We performed 10 runs of 100 simulations each).
Output
Output
Subsample the simulated data in the same way as the real data.
Loop across interference level (within data), deployment length (L), number of sites (no), and delay period (delay) for each site to calculate species counts. Each combination of variables will be sampled n=1 times. 1 combination of sites will be sampled for each n sample. This entire process will be repeated for all 1000 simulated datasets.
Loop across interference level (within data), deployment length (L), number of sites (no), and delay period (delay) for each site to calculate species counts. Each combination of variables will be sampled n=1 times. 1 combination of sites will be sampled for each n sample. This entire process will be repeated for all 1000 simulated datasets.
Take the mean and s.d. Rp (of 1000 samples) for each combination of deployment length, number of sites, and delay, loop this over 6 levels of interference imagery.
Determine the effect of interference images on Rp.
Take the proportion (of 1000 samples) for each species captures for each combination of deployment length, number of sites, and delay, table this over 6 levels of interference.
Determine the effect of interference images on rare species subsampling ratio (SR).
Trigger Frequency
Trigger Frequency
The idea is to simulate a sequence of identical interference events (meaning anything other than a target species we are interested in: other animals, tree branches, etc.) and see how much of the time the camera is idle as result. From that we can calculate the effect on capturing the target animal.
The data structure of a simulation is a sequence of triples: {“appearance time”,”capture type”, “reactivation time”}. Appearance time is simply when a (simulated) animal makes an appearance in front of the camera (whether captured by that camera or not). Capture type is a string with three possible values: “None”, “Delayed”, “Live”. “None” means no capture because the camera was idle, “Delayed” means no immediate capture because the camera was idle, but the animal was still there when it became active and is captured at that time. “Live” means that the animal appeared when the camera is active and is captured right away. For a live capture, the reactivation time is the capture (appearance) time plus the delay period. For a delayed capture, the reactivation time is the previous reactivation time plus a(nother) delay period.
Functions
Functions
Housekeeping function.
This function introduces a new animal after a random (Exponential) interval and adds the appropriate record to the timeseries.
This code makes a graphic from a simulation output. Blue points are live captures, Green points are delayed captures (the point shows the appearance time, not the capture time), black points are non-captures, and the red lines show when the camera is idle.
As well a making graphics, we can calculate different metrics from the simulation output…
…or from a new, long simulation.
The reason we might want to know this metric is that for a rare species, we can calculate pCap, the probability that any one appearance will result in a capture give the specified level of interference. It will be the non-idle fraction (live captures) + the idle fraction * the persistence time of the rare species / the delay period (delayed captures). Knowing that, we can calculate the expected number of successful captures for any occurrence frequency and window of camera activity (see end below).
The below function takes as input the fraction of time the camera is live (from above), the delay period, the persistence time of the target species (how long it stays in front of the camera, in seconds), the frequency in which the target species returns to the camera (calculated from return time), and the length of time the camera is operational. Optionally, you can add the number of independent sites (second function below). The output is the absolute probability of at least one capture of the target species with those parameters.
The below function takes as input the fraction of time the camera is live (from above), the delay period, the persistence time of the target species (how long it stays in front of the camera, in seconds), the frequency in which the target species returns to the camera (calculated from return time), and the length of time the camera is operational. Optionally, you can add the number of independent sites (second function below). The output is the relative probability of at least one capture of the target species with those parameters.
Output
Output
First, vary the interference frequency.
Obviously, the more interference, the more time the camera is idle, but it is non-linear because empty images interfere with themselves when they are common. Happily, the relationship is very well fit by a simple, two-parameter nonlinear function:
The reason this is good is that we can store just this parameter instead of the entire set of results while we vary other things like delay and persistence.
Next, determine the capture probability for species of varying return times across combinations of delay, interference frequency, and interference persistence.
Return time 7 days (1 week).
Return time 30 days (1 month).
Return time 120 days (4 months).
Sample graphs: