www.phidot.org

Hello,
I am running a RDOccupEG model. I have 70 paired sites that were surveyed in the early 1900s and again in the early 2000s. I would like to test for an "era" effect in my p models because detection probably differs between these time periods. However, I have not been able to specify and test the effect of the primary session (i.e. era) in my models. Neither time nor Time test differences between the primary sessions.
The following code is from the RDOccupancy example in R, see help(RDOccupancy), however I modified it to be a RDOccupEG model. This example data has 35 independent 'sites' on which presence/absence of a breeding bird is recorded 3 times annually for 3 years, thus 3 primary sessions.
Can anyone please provide the necessary code to test the effect of primary session in a p model, assuming that this is mathematically possible.
Thanks so much!
Pete

Code: Select all


data(RDOccupancy)         
run.RDExample=function()
{
Psi.area=list(formula=~samplearea)
Psi.cover=list(formula=~cover)
Psi.areabycover=list(formula=~samplearea*cover)
Psi.dot=list(formula=~1)
Psi.time=list(formula=~time)
p.Time=list(formula=~Time)
p.dot=list(formula=~1)
p.occ=list(formula=~occ)
p.area=list(formula=~sample.area)
p.coverbyocc=list(formula=~occ*cover)
gam.dot=list(formula=~1)
epsilon.dot=list(formula=~1)
gam.area=list(formula=~samplearea)
epsilon.area=list(formula=~samplearea)
#  setting time intervals for 3 primary sessions with secondary session length of 3,3,3
time_intervals=c(0,0,1,0,0,1,0,0)
#  Initial data processing for RMARK RDOccupEG
RD_process=process.data(RDOccupancy, model="RDOccupEG", time.intervals=time_intervals)
RD_ddl=make.design.data(RD_process)
#  Candidate model list
 #  1.  all parameters are constant
model.p.dot.Psi.dot.gam.dot.epsilon.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=TRUE)
#   2.  detection varies by Time. Occupancy, colonization, extinction are constant
model.p.Time.Psi.dot.gam.dot.epsilon.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.Time, Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=TRUE)
#   3.  Occupancy varies by time, detection is constant, colonization is constant, extinction is constant
#   Psi.time model doesn't work.
#model.p.dot.Psi.time.gam.dot.epsilon.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.time, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=TRUE)
 #  4.  Occupancy varies by area, detection is constant, colonization varies by area, extinction is constant
model.p.dot.Psi.area.gam.area.epsilon.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.area, Gamma=gam.area, Epsilon=epsilon.dot), invisible=TRUE)
 #  5.  Occupancy varies by cover, detection is constant, colonization varies by area, extinction is constant
model.p.dot.Psi.cover.gam.area.epsilon.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.dot, Psi=Psi.cover, Gamma=gam.area, Epsilon=epsilon.dot), invisible=TRUE)
 #  6.  Occupancy is constant, detection is session dependent, colonization is constant, extinction is constant
model.p.occ.Psi.dot.gam.dot.epsilon.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.occ, Psi=Psi.dot, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=TRUE)
 #  7.  Occupancy varied by area, detection is session dependent, colonization is constant, extinction is constant
model.p.occ.Psi.area.gam.dot.epsilon.dot<-mark(RD_process, RD_ddl, model.parameters=list(p=p.occ, Psi=Psi.area, Gamma=gam.dot, Epsilon=epsilon.dot), invisible=TRUE)
#
#  Return model table and list of models
#
return(collect.models())
}
#  This runs the 7 models above-Note that if you use invisible=FALSE in the above model calls
#  then the mark.exe prompt screen will show as each model is run.
robustexample<-run.RDExample()  #This runs the 6 models above
#  Outputting model selection results
robustexample                                                #  This will print the model selection results on the R Screen
options(width=150)                                           #  Sets page width to 100 characters
sink("results.table.txt")                                    #  Captures screen output to file
print.marklist(robustexample)                                #  Sends output to file
sink()                                                       #  Returns output to screen
system("notepad results.table.txt", invisible=FALSE, wait=FALSE)         #  Allows you to view results in notepad
#  Examine the output for Model 1:  Psi(.), p(.), Gamma(.)
robustexample$model.p.dot.Psi.dot.gam.dot                    #  Opens MARK results file in text editor
robustexample$model.p.dot.Psi.dot.gam.dot$results$beta       #  View beta estimates for specified model in R
robustexample$model.p.dot.Psi.dot.gam.dot$results$real       #  View real estimates for specified model in R
#  Examine the output for the best fitting model (Model 5: Psi(.), p(occ), Gamma(.))
robustexample$model.p.occ.Psi.dot.gam.dot                   #  Opens MARK results file in text editor
robustexample$model.p.occ.Psi.dot.gam.dot$results$beta        #  View beta estimates for specified model in R
robustexample$model.p.occ.Psi.dot.gam.dot$results$real        #  View real estimates for specified model in R
robustexample$model.p.occ.Psi.dot.gam.dot$results$beta.vcv  #  View estimated variance/covariance matrix in R
#  Examine model average estimates for the parameters of interest
model.average(robustexample, "p", vcv=TRUE)                    #  View model averages estimates for session-dependent detection probabilities 
model.average(robustexample, "Psi", vcv=TRUE)                  #  View model averaged estimate for Psi (Occupancy)
model.average(robustexample, "Gamma", vcv=TRUE)               #  View model averaged estimate for Gamma (Colonization)
#
#  Compute real estimates across the range of covariates for a specific model paramter
#  using Model 6
#
#  Identify indices we are interested in predicting
#  see covariate.predictions for information on index relationship to real parameters
summary.mark(robustexample$model.p.occ.Psi.area.gam.dot, se=TRUE) 
#  Define data frame of covariates to be used for analysis
ha<-sort(RDOccupancy$samplearea)
#  Predict parameter of interest (Psi) across the range of covariate data of interest
Psi.by.Area<-covariate.predictions(robustexample, data=data.frame(samplearea=ha), indices=c(1))
Psi.by.Area[1]     # View dataframe of real parameter estimates without var-cov matrix printing (use str(Psi.by.Area) to evaluate structure))


You want to be using session rather than time. "time" or "Time" refers to occasions within a session. To do what you want, define a column in your ddl$p that you can call era. Assign era to be 1 for sessions in the later era and 0 for older era (or vice versa). For example,

myddl$p$era=0
myddl$p$era[myddl$p$session%in%20:35]=1

Then you can use era in your model for p. Don't expect one of these templates to work straight away for a data set other than the example. They are only examples and the code the script contains may not match the type of covariates that have been setup for your data. -jeff

Jeff,
Thanks a lot for your response to my post.
I used your suggested advice on how to define sessions.
I assume that the 20:35 in your code referred to the last 15 rows of data in the example I posted. However, given that each row has capture history data for all of the sessions, it seems that era (or session) needs to be defined by session number, not row. Thus for the example data, one can specify the 3 sessions (in this case testing the 1st session against the 2nd & 3rd) using
RD_ddl$p$era=0
RD_ddl$p$era[RD_ddl$p$session%in%2:3]=1

or in my case, with two eras (sessions):
RD_ddl$p$era=0
RD_ddl$p$era[RD_ddl$p$session%in%2]=1

This worked well when I included era in my p models.
Thanks again!
-Pete