# \donttest{
# This example is excluded from testing to reduce package check time
data(robust)
run.robust=function()
{
#
# data from Robust.dbf with MARK
# 5 primary sessions with secondary sessions of length 2,2,4,5,2
#
time.intervals=c(0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0)
#
# Random emigration, p=c varies by time and session, S by time
#
S.time=list(formula=~time)
p.time.session=list(formula=~-1+session:time,share=TRUE)
GammaDoublePrime.random=list(formula=~time,share=TRUE)
model.1=mark(data = robust, model = "Robust",
time.intervals=time.intervals,
model.parameters=list(S=S.time,
GammaDoublePrime=GammaDoublePrime.random,p=p.time.session),threads=2,delete=TRUE)
#
# Random emigration, p varies by session, uses Mh but pi fixed to 1,
# S by time.This model is in the example Robust with MARK but it is
# a silly example because it uses the heterogeneity model but then fixes
# pi=1 which means there is no heterogeneity.Probably the data were
# not generated under Mh. See results of model.2.b
#
pi.fixed=list(formula=~1,fixed=1)
p.session=list(formula=~-1+session,share=TRUE)
model.2.a=mark(data = robust, model = "RDHet",
time.intervals=time.intervals,
model.parameters=list(S=S.time,
GammaDoublePrime=GammaDoublePrime.random,
p=p.session,pi=pi.fixed),threads=2,delete=TRUE)
#
# Random emigration, p varies by session, uses Mh and in this
# case pi varies and so does p across
# mixtures with an additive session effect.
#
pi.dot=list(formula=~1)
p.session.mixture=list(formula=~session+mixture,share=TRUE)
model.2.b=mark(data = robust, model = "RDHet",
time.intervals=time.intervals,
model.parameters=list(S=S.time,
GammaDoublePrime=GammaDoublePrime.random,
p=p.session.mixture,pi=pi.dot),threads=2,delete=TRUE)
#
# Markov constant emigration rates, pi varies by session,
# p=c varies by session, S constant
# This model is in the example Robust with MARK
# but it is a silly example because it
# uses the heterogeneity model but then fixes pi=1
# which means there is no heterogeneity.
# Probably the data were not generated under Mh.
# See results of model.3.b
#
S.dot=list(formula=~1)
pi.session=list(formula=~session)
p.session=list(formula=~-1+session,share=TRUE)
GammaDoublePrime.dot=list(formula=~1)
GammaPrime.dot=list(formula=~1)
model.3.a=mark(data = robust, model = "RDHet",
time.intervals=time.intervals,
model.parameters=list(S=S.dot,
GammaPrime=GammaPrime.dot,
GammaDoublePrime=GammaDoublePrime.dot,
p=p.session,pi=pi.session),threads=2,delete=TRUE)
#
# Markov constant emigration rates, pi varies by session,
# p=c varies by session+mixture, S constant. This is model.3.a
# but allows pi into the model by varying p/c by mixture.
#
S.dot=list(formula=~1)
pi.session=list(formula=~session)
GammaDoublePrime.dot=list(formula=~1)
GammaPrime.dot=list(formula=~1)
model.3.b=mark(data = robust, model = "RDHet",
time.intervals=time.intervals,
model.parameters=list(S=S.dot,
GammaPrime=GammaPrime.dot,
GammaDoublePrime=GammaDoublePrime.dot,
p=p.session.mixture,pi=pi.session),threads=2,delete=TRUE)
#
# Huggins Random emigration, p=c varies by time and session,
# S by time
# Beware that this model is not quite the same
# as the others above that say random emigration because
# the rates have been fixed for the last 2 occasions.
# That was done with PIMS in the MARK example and
# here it is done by binning the times so that times 3 and 4
# are in the same bin, so the time model
# has 3 levels (1,2, and 3-4). By doing so the parameters
# become identifiable but this may not be
# reasonable depending on the particulars of the data.
# Note that the same time binning must be done both for
# GammaPrime and GammaDoublePrime because the parameters are
# the same in the random emigration model. If you
# forget to bin one of the parameters across time it will fit
# a model but it won't be what you expect as it will
# not share parameters. Note the use of the argument "right".
# This controls whether binning is inclusive on the right (right=TRUE)
# or on the left (right=FALSE). Using "right" nested in the list
# of design parameters is equivalent to using it as a calling
# argument to make.design.data or add.design.data.
#
S.time=list(formula=~time)
p.time.session=list(formula=~-1+session:time,share=TRUE)
GammaDoublePrime.random=list(formula=~time,share=TRUE)
model.4=mark(data = robust, model = "RDHuggins",
time.intervals=time.intervals,design.parameters=
list(GammaDoublePrime=list(time.bins=c(1,2,5))),
right=FALSE, model.parameters=
list(S=S.time,GammaDoublePrime=GammaDoublePrime.random,
p=p.time.session),threads=2,delete=TRUE)
return(collect.models())
}
robust.results=run.robust()
#
# You will receive a warning message that the model list
# includes models of different types which are not compatible
# for comparisons of AIC. That is because
# the runs include closed models which include N
# in the likelihood and Huggins models which don't include
# N in the likelihood. That can be avoided by running
# the two types of models in different sets.
#
robust.results
# }
Run the code above in your browser using DataLab