RMark vs MARK comparison

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.

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.