RMark v MARK
Up to this point, we have used some fairly simple procedures based on Lincoln-Petersen that could be solved by hand, or the calculations easily coded in a spreadhseet or simple R program. From this point forward we will be using more complex models based on maximum likelihood that compute estimates and confidence intervals and perform certain additional functions such as model comparison, model averaging, and goodness of fit evaluations.
There are a number of choices available but the most popular and convenient toolbox for these analyses is based around the program MARK. MARK itself runs out of a graphical user interface (GUI) that while powerful can at times be cumbersome to use. For class we will use the R interface to MARK, package RMark. This package requires the installation of MARK (because MARK is still doing the computing), but models are created using the R scripting language within RMark. In my view this offers 2 advantages over native MARK. First, the models are created using a notational syntax that is a variant on what we have already used in unmarked() to conduct occupancy and distance analysis, so will appear somewhat familiar to you. Second, you will find that one of the more tedious chores of these analyses is getting the data in the right format , getting data into programs, and summarizing results, all of which can be much easier using the data handling power of R.
Basis of MARK models - PIM and Design matrices
Note: In the current versions of MARK and RMark N (abundance) is no longer a real parameter. Instead, f0 the number of unmarked animals is the real parameter and estimates of N are obtained as derived parameters by adding to f0 the number of marked animals in the study (which is assumed to be known and fixed in a closed study). I have tried to reflect this notational change in this content but in case I refer to N, I mean "f0".
Regardless of whether we use RMark or MARK, we need to be familiar with some basic concepts for structuring a model. There are 2 essential elements: the Parameter Index Matrix (PIM) and the Design Matrix (DM). The PIM is basically a list of parameters by different types, with an indexing number associated with each. The parameter types and their indices are determined by (1) the basic data type (e.g., closed vs . open CMR), and (2) the dimensions of the data (number of sampling occasions and groups). We will illustrate for the basic Closed CMR model structure and a single group (stratum) with 5 capture occasions. This model and data structure would have a PIM with 3 different parameter types: capture probabilities (p), recapture probabilties (c) and abundance (f0 = N- marked animals). A basic PIM allowing these 3 parameter types to have distinct values from each other and over occasions would number the parameters:
p1-p5 1:5
c1-c-4 6:9
f0 10
That is, there are 5 distinct capture probabilities allowed over the 5 occasions, followed by 4 recapture probabilities (note that while these are labeled c1-c4 they refer to probability of recapture of a previously marked animal on occasions 2-5). Finally, the last of these parameters is unmarked animals (f0), which because we are assuming closure obviously does not vary. Note that it is important that parameters have distinct indices, that is ordinarily we will not allow parameters of different types to "share" an index. An exception to this will be seen below, where we allow p and c to share indices (since under some models we want these parameters to be the same).
The DM operates on a basic model structure define by the PIM, and creates particular models under specific assumptions. For example
specified a model in which p and c are constant over time but different from each other (and of course from f0); this model has 3 parameters. By comparison the model
indicates that capture and recapture are the same (p=c) and constant, and has 2 parameters. Note that the columns of the design matrix wind up corresponding to the number or paramters we are trying to estimate in the model.
In practice we will create models using the scripting language in RMark, as will be seen in the examples below. RMark essentially combines the PIM and DM steps into one step, eliminating parameter indices if these are not ultimately used in the design. To illustrate, if wind up building a model with both time and behavior effects, we would use a script similar to this to define capture probability
(~time+c)
and the model would need to keep trakc of the 1-5 time indices for capture and recapture. If by contrast we were only interested in a behavior effect, a model like
(~1+c)
would be used, and there would be no need to keep track of time indices (since the effect do not involve time). This will be clearer with specific examples, below. The ~1 indicates that we are replacing the time effects for occasions 1-5 with an intercept.
Also note that we could in theory create a model in which the number of parameters estimated equaled the original list in the PIM (in the above example 10), but in fact that will almost never be the best model, and it will usually be better to construct a series of simpler, alternative models.
Setting up an analysis in RMark
The first step in a RMark (or MARK) analysis is getting the data in the right format, which for RMark will be in a special dataframe. For the moment we're going assume we have our data formatted, and use a canned example to illustrate model building and estimation. We'll then come back and show how to put data into RMark.
First we need to load the RMark library and fetch a data example. The following commands will fetch a classic example by Edwards and Eberhardt.
> library(RMark)
> data(edwards.eberhardt)
The command
> edwards.eberhardt
displays the data. You will see that the data is a matrix of 76 rows and 18 columns, containing 1's and 0's. The rows are the capture histories for individual animals (76 in all) and the 1/0 indicates capture (1) or not capture (0) in each of 18 occasions.
The basic syntax for reading in the data, creating a model, and getting estimates is fairly simple, and looks generically like this
model_object<-mark(data,model, model.parameters)
"Data" tells us where the data is, "model" specifies the model structure (in this case, closed CMR),and model.parameters will indicate the specific model we wish to build. The command
mb<-mark(edwards.eberhardt,model="Closed")
will actually create a simple (but not the simplest) model for this structure and store the model and any results in the object "mb". By default, the basic model has a single parameter for each parameter type (p, c, and N) , so no variation in capture over time. There are therefore 3 parameters. We can run this and see some results right away:
mb<-mark(edwards.eberhardt,model="Closed")
.........................
...............
> mb$results$real
estimate se lcl ucl fixed note
p g1 t1 0.0868518 0.0215415 0.0528969 0.1393947
c g1 t2 0.0809816 0.0095560 0.0641189 0.1017966
f0 g1 a0 t1 17.8918560 10.8161440 -3.3077870 39.0915000
mb$results$derived
estimate se lcl ucl
1 93.89185 10.81614 81.99214 129.4231
> mb$results$AICc
[1] 381.5498
Note by the way that the explicit formulation of this model gives exactly the same results
pdot.cdotshared=list(formula=~1+c,share=TRUE)
mb_e<-mark(edwards.eberhardt,model="Closed", model.parameters=list(p=pdot.cdotshared))
So, under this model, we get estimates of N of approximately 94, of capture of p =0.087, and of recapture of 0.081, so we might conclude that there is (slight) evidence of trap shyness. But do we really need this much complexity in the model? We do this by creating a model in which p=c, that is the probability of catching an animal that is already marked is no different than catching an unmarked animal. This requires 2 lines: 1 to specify a formula for how p is modeled, and the second that incorporates that formula in a list (actually a 'list of lists'). We use the 'share= TRUE' convention because this will allow p and c (technically different parameter types) to take on common values. The formula ~1 specifies an intercept only for the design matrix, so constant p=c.
>pdotshared=list(formula=~1,share=TRUE)
> m0<-mark(edwards.eberhardt,model="Closed", model.parameters=list(p=pdotshared))
STOP NORMAL EXIT
....
> m0$results$real
estimate se lcl ucl fixed note
p g1 t1 0.081955 0.008816 0.0662522 0.1009771
f0 g1 a0 t1 20.258773 6.878605 6.7767062 33.7408400
> m0$results$derived
estimate se lcl ucl
1 96.25877 6.878605 86.60328 114.7067
> m0$results$AICc
[1] 379.6029
Note that the estimate of N is a bit higher (about 96) but also that the AIC value is lower, 379.6 vs. 381.5. As we have seen in previous labs, we use AIC to guide model selection, or preferably, we can simply average across the alternative models, and take model uncertainty into account. We'll come back to this after we run a few more models.
Let's run an additional model, this time allowing p to vary among occasions (time) but p=c otherwise.
> ptimeshared=list(formula=~time,share=TRUE)
> mt<-mark(edwards.eberhardt,model="Closed",model.parameters=list(p=ptimeshared))
> mt$results$real
estimate se lcl ucl fixed note
p g1 t1 9.465900e-02 3.073600e-02 4.921970e-02 1.743549e-01
p g1 t2 8.414150e-02 2.906480e-02 4.202320e-02 1.613629e-01
p g1 t3 9.465910e-02 3.073600e-02 4.921980e-02 1.743550e-01
p g1 t4 1.472475e-01 3.775660e-02 8.740720e-02 2.373980e-01
p g1 t5 8.414140e-02 2.906480e-02 4.202320e-02 1.613628e-01
p g1 t6 5.258830e-02 2.318190e-02 2.181280e-02 1.213959e-01
p g1 t7 1.893181e-01 4.228040e-02 1.197926e-01 2.860802e-01
p g1 t8 1.156943e-01 3.377600e-02 6.410880e-02 1.999216e-01
p g1 t9 4.207080e-02 2.079510e-02 1.572320e-02 1.077369e-01
p g1 t10 3.155290e-02 1.806130e-02 1.012560e-02 9.401680e-02
p g1 t11 1.682827e-01 4.011370e-02 1.034390e-01 2.619014e-01
p g1 t12 5.258840e-02 2.318190e-02 2.181290e-02 1.213960e-01
p g1 t13 2.103530e-02 1.478950e-02 5.230800e-03 8.071700e-02
p g1 t14 7.362360e-02 2.726840e-02 3.502820e-02 1.482131e-01
p g1 t15 9.465910e-02 3.073600e-02 4.921970e-02 1.743550e-01
p g1 t16 1.536719e-09 3.998270e-06 -7.835073e-06 7.838146e-06
p g1 t17 4.207070e-02 2.079500e-02 1.572310e-02 1.077367e-01
p g1 t18 1.051767e-01 3.230150e-02 5.658910e-02 1.872037e-01
f0 g1 a0 t1 1.907807e+01 6.613441e+00 6.115724e+00 3.204041e+01
> mt$results$derived
estimate se lcl ucl
1 95.07807 6.613441 85.85801 112.9215
> mt$results$AICc
[1] 354.5968
Now we get a lower AIC (the lowest yet) and a value of N in between the last 2 (95 vs. 94 and 96). Maybe we're getting close.
So far we have run 3 of the "classic" models (or their equivalents) discussed by Otis et al. (1978):
The AIC results for these models are summarized below
model npar AICc DeltaAICc weight Deviance
3 p(~time)c()f0(~1) 19 354.5968 0.00000 9.999949e-01 305.2648
1 p(~1)c()f0(~1) 2 379.6029 25.00616 3.715168e-06 364.8259
2 p(~1)c(~1)f0(~1) 3 381.5498 26.95302 1.403539e-06 364.7640
and we can use the AIC weights to calculate a model-averaged estimated for N (code attached to do this):
However, we haven't considered the model that has been shown, time and again, to be the needed because of the presence of individual trap heterogeneity in capture probabilities. This is the heterogeneity model or Mh, and it was part of the original Otis et al. (1978) architecture. However, the original Mh was based on a jackknife procedure and does not involve a likelihood, so cannot be compared to other models via AIC.
Programs MARK and RMARK instead use a likelihood approach based on finite mixtures. This models specifiy that animals are composed of 2 or more (we will use 2) groups, one with a higher and the other a lower capture probability. An additional parameter pi specifies the probability that an animal is a member. These mixture models require invocation of a new class of models (Het or FullHet). We will use the Full Hetergeneity class, which allows for models that combine individual (h) , behavioral (b) and time (t) variation. We will create 3 mixture models that involve individual heterogeneity alone (h), with the presence of time (th) and with behavior (bh); for comparison we will also look at some models within the FullHet class that eliminate the effect of the mixture (so we can compare AIC) values. For instance the Mh model produces
>>pdotmixture=list(formula=~mixture,share=TRUE)
>>ee.closed.Mh=mark(edwards.eberhardt,model="FullHet",model.parameters=list(p=pdotmixture))
> ee.closed.Mh$results$real
closed.Mh$results$real
estimate se lcl ucl fixed note
pi g1 m1 0.1547608 0.0783461 0.0535732 0.3719563
p g1 t1 m1 0.1795149 0.0487418 0.1026379 0.2950415
p g1 t1 m2 0.0360236 0.0196284 0.0121914 0.1016498
f0 g1 a0 t1 59.4770160 36.6027910 -12.2644560 131.2184900
>
estimate se lcl ucl
1 135.477 36.60279 95.58636 256.6112
> ee.closed.Mh$results$AICc
[1] 369.6449
6 Models are considered below (the original 3 plus 3 heterogeneity models):
> ee.results
model npar AICc DeltaAICc weight
6 pi(~1)p(~time + mixture)c()f0(~1) 21 343.2755 0.00000 9.965298e-01
5 pi(~1)p(~time)c()f0(~1) 19 354.5968 11.32123 3.468299e-03
4 pi(~1)p(~mixture)c()f0(~1) 4 369.6449 26.36932 1.872688e-06
1 pi(~1)p(~1)c()f0(~1) 2 379.6029 36.32739 1.288538e-08
2 pi(~1)p(~1)c(~1)f0(~1) 3 381.5498 38.27425 4.867917e-09
3 pi(~1)p(~mixture)c(~1)f0(~1) 5 385.4574 42.18184 6.899525e-10
The top-ranked model now involves both time and heterogeneity; furthermore the model-averaged estimate is now 132.28; note though the fairly large unconditional SE; this reflects the fact that there is quite a bit of variation in the estimates of N among models (even though most of the weight is on the top model).
R code is attached that creates user-defined function to run the above examples and call a built-in model averaging procedure for the non-mixture and mixture set of models. The code also creates input files that can be read by MARK, for users who prefer to operate in the MARK GUI.
Finally, the mixture approach to dealing with heterogeneity has a number of pitfalls, and often does not work well with field data. Sometimes, heterogeneity in p can be described and predicted by measurable attributes of individuals, in which case individual covariate models are to be preferred. We will return to covariate models in some other labs, but I provide an example where a covariate is used to successfully model heterogeneity.
Additional notes
Note that the object ee.results object above contains all the individual model results, a general feathure of these kinds of 'collecting function'. The command summary(ee.results) gives a quick list of the models that have been run, and the command collect(ee.results) gives the AIC table for all the models run. You can find out more about the contents of this (or other R objects) by the command
>attributes(ee.results)
You can drill down farther into the object and get attributes of these sub-objects, e.g.,
> ee.results$ee.closed.Mt$results
will get you details for the Mt model run earlier,and
>ee.results$ee.closed.Mt$results$AICc
gives a particular result-- namely the corrected AIc value for the model.
You could write code in R to pull out AIC values, compute model weights, and perform tasks such as model averaging. However, our friend Jeff Laake has helpfully written some very useful functions to perform these sorts of tasks for you, which reads in all the results in the ee.results object and produces model averaged estimates (see the attached example).
Comparison of RMark and MARK for Edwards - Eberhardt example