I am almost at the end of my first ‘batch’ of analysis using RMARK and have hit a final hurdle, hoping that you can help again. I have tried to provide as much info as possible here, but can send the full script/ input file off-forum if required.
This is a 16-year data-set, the time interval is every 6-months, so 32 occasions. The initial age is 6 months of age for all birds in the study (the current analysis is only dealing with the adult population; in this study system the majority of birds are ringed on the nest so age is known).
In the current analysis I have 2 age structures:
• S1 with three age classes (young= [0.5,1), Middle=[1,2), Old=[2,16]
• S2 with two age classes (subadult = [0.5,1), Adult=[1,16] – i.e. the adult group here is the Middle and Old group from S1 lumped together.
I need to include time*age interactions, so I set up the dummy variables for S1 and S2 (code below)
- Code: Select all
#S1 - three age-class structure
#S1 age bins for age*time interaction using age class structure from S10 in previous step:
IAAAGE1.ddl=add.design.data(IAAAGE1.process,IAAAGE1.ddl,"Phi","age",bins=c(0.5,1,2,16),right=FALSE,name="S1")
#
#Add a field to the design data that is equivalent except it is a numeric dummy coding variable 1.
#Add field for young age group - this one is 0.5 only (note brackets)
#
IAAAGE1.ddl$Phi$S1Young=0
IAAAGE1.ddl$Phi$S1Young[IAAAGE1.ddl$Phi$S1=="[0.5,1)"]=1
#
#then add field for Middle age class - note brackets - this one is 1 to 1.5
IAAAGE1.ddl$Phi$S1Middle=0
IAAAGE1.ddl$Phi$S1Middle[IAAAGE1.ddl$Phi$S1=="[1,2)"]=1
#then add old age class - note brackets - this one is 2 to 16 including 16
IAAAGE1.ddl$Phi$S1Old=0
IAAAGE1.ddl$Phi$S1Old[IAAAGE1.ddl$Phi$S1=="[2,16]"]=1
#
#~~~~~~
# s2 - a two-age class age structure - comparing S1 and S2 will test for senescence in survival. In this age structure middle and old from s1 are merged into adult
#
#S2 age bins for age*time intercation subadult and adult
IAAAGE1.ddl=add.design.data(IAAAGE1.process,IAAAGE1.ddl,"Phi","age",bins=c(0.5,1,16),right=FALSE,name="S2")
#
#Add a field to the design data that is equivalent except it is a numeric dummy coding variable 1.
#Add field for young age group - this one is 0.5 only (note brackets)
#
IAAAGE1.ddl$Phi$S2Subadult=0
IAAAGE1.ddl$Phi$S2Subadult[IAAAGE1.ddl$Phi$S2=="[0.5,1)"]=1
#
#then add field for adult age class - as there are only 2 groups can do this one by subtraction, it is 1-S2Subadult
IAAAGE1.ddl$Phi$S2Adult=1-IAAAGE1.ddl$Phi$S2Subadult
I checked the ddl for Phi and all was as expected.
I set up the models using create.model.list and mark.wrapper (code below)
- Code: Select all
do.analysis=function()
{
#create forumalae for p first:
# P is constant
p.dot=list(formula=~1)
#
#Create forumulae for S:
#
# A fully age dependent but time invariant survival model
Phi.age=list(formula=~age)
#
#A fully time dependent survival model
Phi.time=list(formula=~time)
#
# A constant survival model
Phi.dot=list(formula=~1)
#
#using s1 age classes:
#Fully time*age structured Phi model. cannot use ymo*time because there are no olds for time 1. The minus 1 removes the intercept which is not needed (see C-15 in RMARK guidelines)
Phi.S1Young.time.S1Middle.time.S1Old.time=list(formula=~-1+S1Young:time+S1Middle:time+S1Old:time)
#Limited time-age interaction: young and middle age class vary with time but Old is constant
Phi.S1Young.time.S1Middle.time.S1Old.dot=list(formula=~S1Young:time+S1Middle:time)
#Limited time-age interaction: young and old age class vary with time but middle is constant
Phi.S1Young.time.S1Middle.dot.S1Old.time=list(formula=~S1Young:time+S1Old:time)
#Limited time-age interaction: middle & old age class vary with time but young is constant
Phi.S1Young.dot.S1Middle.time.S1Old.time=list(formula=~S1Middle:time+S1Old:time)
#Limited time-age interaction: old age class varies with time but young and middle constant
Phi.S1Young.dot.S1Middle.dot.S1Old.time=list(formula=~S1Old:time)
#Limited time-age interaction: middle age class varies with time but young and old constant
Phi.S1Young.dot.S1Middle.time.S1Old.dot=list(formula=~S1Middle:time)
#Limited time-age interaction: young age class varies with time but middle and old constant
Phi.S1Young.time.S1Middle.dot.S1Old.dot=list(formula=~S1Young:time)
#Age class only NO age*time interaction (YMO constant with respect to time but differ with age)
Phi.S1=list(formula=~S1)
#
#Using S2 age structure (where middle and old are lumped together)
#Fully time*age structured phi model. cannot use ymo*time because there are no olds for time 1. The minus 1 removes the intercept which is not needed (see C-15 in RMARK guidelines)
Phi.S2Subadult.time.S2Adult.time=list(formula=~-1+S2Subadult:time+S2Adult:time)
#limited time-age interaction: Subadult time dep but adult constant. the intercept is the adult value
Phi.S2Subadult.time.S2Adult.dot=list(formula=~S2Subadult:time)
#limited time-age interaction: Subadult constant, adult time dep, the intercept is the subadult value
Phi.S2Subadult.dot.S2Adult.time=list(formula=~S2Adult:time)
#Age class only NO time interaction (SA constant with respect to time, but differ with age)
Phi.S2=list(formula=~S2)
#
#
#use create model list (CML) to construct the combinations of models and store in cml object
cml=create.model.list("CJS")
#use mark.wrapper to fit each model in mark. adjust = FALSE is NOT included here because I do NOT want to accept MARK's parameter counts, I want to adjust them.
model.list=mark.wrapper(cml,data=IAAAGE1.process,ddl=IAAAGE1.ddl)
#Return the list of model results as the value of the function
return(model.list)
}
#
#Run analysis by assigning it to object:
IAAAGE1results=do.analysis()
The problem:
When I started to look at some of the real parameter estimates for the models for S1, I realised that for the models in which only one of the 3 age classes (young middle or old) had time interaction, and the other 2 age classes were constant, the output only gives real parameter estimates for the age class with time dependency, plus one intercept for phi, but the other (constant) ageclass parameter etsimate is missing from the output. I assume this is down to how I have set up the model structure in the function do.analysis, (the output for models with 2 ageclasses:time and 3ageclasses:time were fine), but having checked the source material on this and tinkered about with it for some time I can’t spot what I have done wrong here and hope someone can put me back on the right track.
Many thanks for all the help provided
Lianne