RMark and MARK equivalency for EE example

In an effort to somewhat de-mystify how RMark works in comparison to MARK, I've set up several closed CMR models for the Edwards-Eberhardt example in both frameworks, and show how they should produce virtually identical results, but there are some subtle distinctions along the way.  

A key point is that the 'unsimplified' PIM (parameter index matrix) structure for both approaches is the same, being determined by factors such as basis model structure and number of parameter types, number of groups, and (depending) number of sample occasions.  For example , for non-mixture closed CMR problem below we have 3 parameter types (p =capture, c= recapture, and N=abundance). There is one group (no strata), and there are 18 capture occasions.  The unsimplified PIM will have 18 rows for p, followed by 17 for c, and 1 for N = 36 rows.  

First, we will read in the data and run some non-mixture models in RMark, and then do the same thing using the built in model or  PIMs in the MARK interface.  

>#data_dir<-"C:/documents and settings/conroy/my documents/Dropbox/teaching/WILD8390/spring2014/week6/closed_cmr"

>setwd(data_dir)

>require(RMark)

>data(edwards.eberhardt)

We then create some shortcuts for parameter formulas that will later be used to create specific models 

>#model formulas

>pdotshared=list(formula=~1,share=TRUE)

>ptimeshared=list(formula=~time,share=TRUE)

>ptimeshared.c=list(formula=~time+c,share=TRUE)

The "share =TRUE" option allows the p and c parameters to overlap in the model, which we need to be able to do to create all the models we need (and which is ok because these are not really different parameter types like p and N are).  We then use these to create 3 models that correspond to the null model, behavior only model, and time only model). 

m0<-mark(edwards.eberhardt,model="Closed",model.parameters=list(p=pdotshared))

mb<-mark(edwards.eberhardt,model="Closed")

mt<-mark(edwards.eberhardt,model="Closed",model.parameters=list(p=ptimeshared))

We can pick any one of these to see the 'unsimplified' PIM structure, because that does not depend on the model

> ###unsimplified PIM structure models

> PIMS(m0,"p",simplified=F)

group = Group 1 

   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18

1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18

> PIMS(m0,"c",simplified=F)

group = Group 1 

   2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18

1 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

> PIMS(m0,"N",simplified=F)

group = Group 1 

  [,1]

1   36

This simply confirms that the first 18 parameters are capture probabilities, the next 17 are recapture, and the last is abundance. Depending on the specific model, the estimates may be the same across some of these dimensions, because the model is simpler than the PIM structure. E.g., 

>summary(m0)

Output summary for Closed model

Name : p(~1)c()N(~1) 

Npar :  2

-2lnL:  375.5941

AICc :  379.6029

Beta

               estimate        se       lcl       ucl

p:(Intercept) -2.416076 0.1171742 -2.645737 -2.186414

N:(Intercept)  3.008588 0.3395372  2.343095  3.674081

Real Parameter p

        1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18

 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955

Real Parameter c

        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18

 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955 0.081955

Real Parameter N

     [,1]

 96.25877

By comparison the simplified PIM structure for this model reveals that there are only 2 "real"  indices under this model: 1 for p and 1 for N, and 

group = Group 1 

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

1 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1

> PIMS(m0,"c")

group = Group 1 

  2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1

> PIMS(m0,"N")

group = Group 1 

  [,1]

1    2

> m0$results$real

         estimate       se        lcl         ucl fixed    note

p g1 t1  0.081955 0.008816  0.0662522   0.1009771              

N g1    96.258771 6.878605 86.6032790 114.7066900          

You can see a similar pattern for the other  models, e.g., Mb now has 3 real parameters (p, c, and N).  

You can collect these models into a results object and get real parameter model averaged estimates. e..g

ee.results<-collect.models()

model.average(ee.results,"N")

 The R script to run these examples is here. (Note: this script has the same name as the earlier one we used for this example but it is a different file!)

To run the same models in MARK we would read the data in using MARK format (slightly different with all the capture histories requiring a frequency (1 in this case) and a trailing ; , e.g., 111100100010000000 1;). The data and MARK project files are here.  If you open up the MARK project (which has pre-specified that these are closed capture models with no mixtures) , create the most general model p(t) c(t) N(), and look at the PIM matrices you see

 which is identical the simplifed PIM  provided by RMark .

You can run the equivalent MARK models, compare real parameter estimates, and confirm that the results are identical.

Mixture models

Mixture models add a couple of wrinkles to the above procedure but they are fairly straightforward. First, the number of mixtures must be specified when the data are models are first created (in the case of MARK, when the data are read into the model).  We then need to add an additional parameter type pi for the mixture probability.   We will use 2 mixtures, so there will be a single mixture probability (p= probability of being in group 1 and 1-p = probability of being in group 2).  This in turn means that in the general PIM structure there number of p and c parameter indices doubles (since we now have p's and c's for each group).    The RMark statements to create models now involve pi, p, and c . For example

>pmixture=list(formula=~mixture)

.cmixture=list(formula=~mixture)

>ee.closed.Mbh=mark(edwards.eberhardt,model="FullHet",model.parameters=list(p=pmixture,c=cmixture))

creates an "Mbh" model with both behavior and heterogeneity (mixture) effects.  The simplified PIM structure is.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

1 2 2 2 2 2 2 2 2 2  2  2  2  2  2  2  2  2  2

2 3 3 3 3 3 3 3 3 3  3  3  3  3  3  3  3  3  3

> PIMS(ee.closed.Mbh,"c")

group = Group 1 

  2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

1 4 4 4 4 4 4 4 4  4  4  4  4  4  4  4  4  4

2 5 5 5 5 5 5 5 5  5  5  5  5  5  5  5  5  5

> PIMS(ee.closed.Mbh,"N")

group = Group 1 

  [,1]

1    6

so there are 6 parameters: the mixture probability pi, p for mixture  1 and mixture 2,  c for mixture 1 and mixture 2, and N. .

We can create additional models:

>#Mh

>pdotmixture=list(formula=~mixture,share=TRUE)

>ee.closed.Mh=mark(edwards.eberhardt,model="FullHet",model.parameters=list(p=pdotmixture))

and

># Closed Mth - p different for time; mixture additive

>ptimemixtureshared=list(formula=~time+mixture,share=TRUE)

>ee.closed.Mth.additive=mark(edwards.eberhardt,model="FullHet",model.parameters=list(p=ptimemixtureshared),adjust=TRUE)

The last of these 'shares' the p an c parameters so that they can be allowed to be time varying but equal (no behavioral effects) under an additive model

I created equivalent full heterogeneity models in MARK, which you can view here.   Some of these were fairly easy (e.g., built in or via PIM manipulation ) but others like the Mth additive model required a lot of manipulation with the design matrix in MARK.  This is a good example where RMark provides a much easier flow in model construction (albeit with some needed understanding of syntax). 

The R script for both RMark examples is here.(Note: this script has the same name as the earlier one we used for this example but it is a different file!)

NOTES ON PARAMETER COUNTING IN RMARK VS MARK.

• By default, RMark assumes that the number of model parameters is equal to the number of columns of the design matrix

• MARK starts out with this assumption, but will adjust downward the number of estimable parameters if the data are sparse relative to model complexity

• This sets up a situation where AIC values for the same model and ordering for the same set of models potentially can differ between RMark and MARK, since AIC depends on the parameter count 

• The default computation of parameter count can be changed in either RMark or MARK to whatever is desired. In my view the RMark convention of counting based on the specified model is more sensible .