For both model types the input data was structured as: ID, ch, Sex
(eg Fish1, 00001000000100000010000000000000..., M)
So far POPAN models seem to return fairly reasonable phi, p, and derived N estimates. However, the robust design models give fairly high p estimates and derived N values that are equal or slightly less than the number of fish captured. Similarly, the se for these derived N estimates are typically fewer than 1.
Do you have suggestions on where I could have gone wrong or how to extract population size? Below is the RMark code that I would use to generate the current top model. For model selection I also tried iterations where Gamma" and Gamma' were shared, as well as where p and c were shared.
Thanks!
Greg
- Code: Select all
RD_covars.proc <- process.data(data = robust_covars,
model = "Robust",
groups = "Sex",
time.intervals = time_intervals_dates)
RD.models = function (){
# S.dot = list(formula=~1)
# S.time = list(formula=~time)
# S.sex = list(formula=~group)
S.time.sex = list(formula=~time + group)
# GammaDoublePrime.Markov.dot = list(formula=~1, share=FALSE)
# GammaDoublePrime.Markov.time = list(formula=~time, share=FALSE)
# GammaDoublePrime.Markov.sex = list(formula=~group, share=FALSE)
GammaDoublePrime.Markov.time.sex = list(formula=~time + group)
# G' = 1/4
# GammaPrime.Markov.dot = list(formula=~1, share=FALSE)
GammaPrime.Markov.time = list(formula=~time)
# GammaPrime.Markov.sex = list(formula=~group, share=FALSE)
# GammaPrime.Markov.time.sex = list(formula=~time + group, share=FALSE)
# p = 7/7
# p.dot = list(formula=~1, share=FALSE)
# p.time = list(formula=~time, share=FALSE)
# p.sex = list(formula=~group, share=FALSE)
# p.time.sex = list(formula=~time+group, share=FALSE)
# p.session.time = list(formula=~time+session, share=FALSE)
# p.session.sex = list(formula=~session+group, share=FALSE)
p.session.time.sex = list(formula=~time+session+group)
# c = 7/7
# c.dot = list(formula=~1, share=FALSE)
# c.time = list(formula=~time, share=FALSE)
# c.sex = list(formula=~group, share=FALSE)
# c.time.sex = list(formula=~time+group, share=FALSE)
# c.session.time = list(formula=~time+session, share=FALSE)
# c.session.sex = list(formula=~session+group, share=FALSE)
c.session.time.sex = list(formula=~time+session+group)
cml=create.model.list("Robust")
results=mark.wrapper(cml, data=RD_covars.proc, output=FALSE)
return(results)
}
RD.models.list <- RD.models()
RD.models.list$S.time.sex.GammaDoublePrime.Markov.time.sex.GammaPrime.Markov.time.p.session.time.sex.c.session.time.sex$results$derived$`N Population Size`