secr is an R package written by Murray Efford designed to compute maximum likelihood and 'inverse prediction' estimates from spatially explicit CMR data. As with all R libraries, secr must be installed (downloaded from one of the CRAN mirror sites) and loaded before running. In lab we will generally have secr and other libraries installed but you will need to load the library in your session:
>library(secr)
If secr is not installed you will get an error message to that effect, in which case you need to go the the main menu "Packages/Install Package". You will be asked to choose a CRAN mirror site (pick one nearby [in your hemisphere at least] for speed) ; you of course need to be connected to the internet.
Help with secr can be gotten by typeing ?secr at the prompt:
>?secr
which will provide html help and also point you to several pdf manuals; 2 are provided here.
General introduction
Data formatting and input
Data input
There are several options in secr for inputting data. The most basic is to read from 2 text files, one containing the labels and x-y coordinates for the detector array (e.g., grid of camera traps), the other containing individual records for each animal capture (ID, capture occasions, and capture location). Other files, such as file containing the boundaries of the study area, or a habitat mask stipulating areas of habitat/ non-habitat within the study area, but are not required. I created 2 small files to illustrate this process, a detector file and a capture file. The detector file can look like this:
# Detector X Y
A1 -1500 -1500
A2 -1500 -1250
A3 -1500 -1000
A4 -1500 -750
where each line represents a detector (trap) location, the labels are any specified names for the detectors, and the numbers are the grid coordinates (in this case, all negative with respect to some specified origin). The capture data can look like this
# Session ID Occasion Detector
MatakitakiStoats 1 1 A1
MatakitakiStoats 1 2 A3
MatakitakiStoats 1 3 A1
MatakitakiStoats 2 1 A3
MatakitakiStoats 2 2 A2
MatakitakiStoats 9 1 A4
where the session s just a label specified by the investigator (e.g., species, location, year), ID is the individual animal ID (3 unique individuals in this case), Occasion is the sampling occasion (1, 2, or 3), and Detector is the label referencing the detector at which each detection ("capture") occurred; obviously these need to match the labels in the detector file, but not all detectors in the array will have a detection necessarily.
These files are then read into a secr data object like so
example<-read.capthist(captfile="captures.txt", trapfile="traps.txt", detector = "single", fmt = "trapID")
where in this example it is specified that the detector is "single", meaning that no more than one animal can be detected at time; in many cases this will be "multiple", for instance with camera traps where of course >1 animal can be detected in an occasion.
Building, running, and comparing models
Once the secr data object is defined we can perform other tasks including building and comparing models. I will use a built-in example for brushtail possums (Trichosurus vulpecula) in New Zealand, an exotic invasive from Australia (no this is not a mis-spelling of oppossum!). The data are loaded by the command
>data(possum)
and the command
>possumCH
displays the capture histories for each of the 5 sampling occasions; note that the entries 0 signify not captured and 9=80, 77, etc. correspond to the line in the trapping database.
1 2 3 4 5
1501 0 80 0 0 84
1502 0 78 0 76 77
1503 0 77 0 0 0
1504 0 0 80 0 0
1506 0 0 113 114 0
1509 0 0 0 0 163
...
...
..
The trap locations are displayed by
>traps(possumCH)
which is a dataframe with 180 lines (so the numbers in the capture history frame are between 0 and 180). We can do some useful preliminary analysis. For instance
> summary(possumCH)
provides summary statistics that if they have no other function may help diagnose data entry errors:
Object class capthist
Detector type single
Detector number 180
Average spacing 20 m
x-range 2697793 2699225 m
y-range 6077470 6078201 m
Counts by occasion
1 2 3 4 5 Total
n 61 51 56 66 73 307
u 61 19 9 14 9 112
f 32 21 22 18 19 112
M(t+1) 61 80 89 103 112 112
losses 0 0 0 0 0 0
detections 61 51 56 66 73 307
detectors visited 61 51 56 66 73 307
detectors used 180 180 180 180 180 900
The commands
>plot(possummask)
>plot(possumCH, tracks = TRUE, add = TRUE)
>plot(traps(possumCH), add = TRUE)
>lines(possumarea)
produce some overlaid plots of the study area, habitat mask, and trap grid (again, the study boundaries and habitat mask would have to be read in from text files, possibly created in GIS:
Note that the trapping array actually is composed of 5 grids spread throughout the study area. Not only is this allowed, it's probably a better design (in terms of representation not to mention logistics) than either spreading the traps out uniformly over the whole area, or clustering them in 1 huge grid (say in the center).
Models are created in secr using a syntax very similar to that which we have seen for previous examples in RMark, Unmarked, etc. There are 3 types of parameters in secr : density (D), home range center detection (g0), and movement scaling parameter (sigma), and models basically consist of combinations of factors including covariates that influence variation in these parameters. The simplest (default) model would be
model = list(D~1, g0~1, sigma~1)
Depending on the type of detector, there may be the possibility of only a single detection of an animal at a given trap-occasion, or multiple detections. In the first case, the detections are modeled as Bernoulli responses (1 detected or 0 not detected), in the second as Poisson counts. The possum example is based on binary detections. Often, even if the response theoretically involves counts (e.g., camera trapping), it is sufficient to reduce to data to binary (detected or not at each occasion), as in capture-recapture, since binomial models are faster and tend to produce very similar results to those from Poisson models, especially if there are few multiple detections.
Since this is the default model, it is the one created if no model statement is included. Putting this altogether for the possum example, I created and ran 3 models and performed AIC comparison. I focussed on variation in g0, since I had no a priori reason to think that either density or movement rates varied either over time (it's a closed model, dammit!) or with respect to covariates. Analysis is S-L-O-W; there's a lot of data information here and calculations to perform, even for 'simple' models.
>#Null model (default)- no variation in parameters
>possum.model.0.NPT <- secr.fit(subset(possumCH,NPT), mask =possummask, trace=F)
>#behavior model (1 time change in g0 after first capture)
>possum.model.b.NPT <- secr.fit(subset(possumCH,NPT),model=g0~b, mask = possummask, trace=F)
>#time model
>possum.model.t.NPT <- secr.fit(subset(possumCH,NPT),model=g0~t, mask = possummask, trace=F)
>AIC (possum.model.b.NPT,possum.model.0.NPT,possum.model.t.NPT)
model detectfn npar logLik AIC AICc dAICc AICwt
possum.model.t.NPT D~1 g0~t sigma~1 halfnormal 7 -1026.765 2067.530 2068.697 0.000 0.5389
possum.model.0.NPT D~1 g0~1 sigma~1 halfnormal 3 -1031.687 2069.373 2069.613 0.916 0.3409
possum.model.b.NPT D~1 g0~b sigma~1 halfnormal 4 -1031.647 2071.294 2071.698 3.001 0.1202
The results indicate best support (among these 3 candidates) for the time-specific model. The real parameter estimates for that moddel are below, indicating possum densities of about 1.5/ ha and movement scaling on the order of 52m.
Fitted (real) parameters evaluated at base levels of covariates
link estimate SE.estimate lcl ucl
D log 1.5385009 0.16487941 1.2477734 1.8969669
g0 logit 0.2160281 0.03638561 0.1531576 0.2956951
sigma log 52.0102277 2.71255275 46.9596821 57.6039629
The package also contains a fairly easy to implement model averaging function. Model averaged estimates for the real parameters are gotten by:
>#model averaging
>model.average(possum.model.0.NPT,possum.model.b.NPT,possum.model.t.NPT)
As an aside here, note that whereas possums moved in a 50-m or so radius, average trap spacing was about 20m (except obviously in the areas outside the individual grids). So individual possums had a good chance of being captured on >1 detector. This is crucial, because if animals do not move among detectors there is no basis for estimating sigma, the scaling parameter, and the method falls apart. I have experienced this for a white-tailed deer study, in which cameras were places 1.5- 2km apart, but deer moved far less than this (700 or so m for bucks, <500 m for does), resulting in virtually no recaptures between cameras. Thus a general recommendation for study design is to place at least 2 detectors per animal home range, with more being better.
The R script for the above examples is contained here.
Next:SCR using BUGS/ JAGS in R