Unfortunately, the non-mixture models ignore 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 (HetClosed or FullHet). We will use the Full Hetergeneity class, which allows for models that combine individual (h) , behavioral (b) and time (t) variation. The HetClosed models, by contrast, combing individual (h) and time (t) variation, but not behavioral. See the List of MARK and RMark Models for details on models and parameter types.
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
estimate se lcl ucl fixed note
pi g1 m1 0.1547609 0.0783461 0.0535734 0.3719561
p g1 t1 m1 0.1795148 0.0487417 0.1026380 0.2950412
p g1 t1 m2 0.0360236 0.0196284 0.0121914 0.1016496
N g1 135.4770100 36.6027490 95.5863800 256.6109600
>
> ee.closed.Mh$results$AICc
[1] 369.6449
6 Models are considered below (the original 3 plus 3 heterogeneity models):
Otis notation RMark model notation in FullHet Parameter variation
M0 pi(~1)p(~1)c()N(~1) No variation in p, p=c
Mb pi(~1)p(~1)c(~1)N(~1) Behavior variation (p !=c) but constant in time
Mt pi(~1)p(~time)c()N(~1) Time variation in p, p=c
Mh pi(~1)p(~mixture)c()N(~1) Individual variation in p
Mth pi(~1)p(~time + mixture)c()N(~1) Individual and time variation in p
Mbh pi(~1)p(~mixture)c(~1)N(~1) Individual and behavior variation in p
> ee.results
model npar AICc DeltaAICc weight
6 pi(~1)p(~time + mixture)c()N(~1) 21 343.2755 0.00000 9.965298e-01
5 pi(~1)p(~time)c()N(~1) 19 354.5968 11.32123 3.468299e-03
4 pi(~1)p(~mixture)c()N(~1) 4 369.6449 26.36932 1.872688e-06
1 pi(~1)p(~1)c()N(~1) 2 379.6029 36.32739 1.288538e-08
2 pi(~1)p(~1)c(~1)N(~1) 3 381.5498 38.27425 4.867917e-09
3 pi(~1)p(~mixture)c(~1)N(~1) 5 385.4574 42.18184 6.899525e-10
> N<-subset(avg,par.index==72)
> Nhat<-N$estimate+nmarked
> seN<-N$se
> print("Model averaged N and se")
[1] "Model averaged N and se"
> Nhat
[1] 132.2824
> seN
[1] 33.53412
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 a user-defined function to run the above examples and call a built-in model averaging procedure for the 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 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.
>U<-model.average(ee.results,"N")
>#this is actually the estimate of the number of UNMARKED animals in the pop!
>#Add known number of marked animals to this to get Nhat
>nmarked<-nrow(edwards.eberhardt)
>Nhat<-U$estimate+nmarked
>seN<-N$se
print("Model averaged N and se")
Nhat
seN
[1] 132.2824
> seN
[1] 33.53412
Comparison of RMark and MARK for Edwards - Eberhardt example
See also Appendix C of the online MARK book
Next: Huggins model for heterogeneity and individual covariates