Robust model

[Tara]

Hello

We have three years of prairie dog data that we are looking to analyze using a Robust model. We have version 1.9.7 of RMark and version 2.10.1 for R. The data consists of three primary periods (2007 – 2009) and 9 secondary periods. We have information on age (that we have binned into young and adult), sex, and location (5 Colonies). There was a crash (2007-2008) and then a rebound (2008-2009) in the population. Furthermore, very few animals reproduced in 2008 (only 3 pups were captured).

Below is a script we have used to try model survival and abundance however despite trying many different "s" formulas we are consistently getting either very small (ie.0.0000000) or very large(ie.232.00) SE for our beta coefficients which seems suspicious. Furthermore, pup survival for 2008-2009 is often high even on colonies where pups were never caught. As you can see in the script below we have tried to bin and share gamma prime and gamma double prime, in order to try and solve the identity problem for parameters of survival with a time model. However, these don’t seem to be working. We are not sure if the problem is with our model formulas or if our data is just possibiliy too sparse for some of these estimates. If you could provide any insight into what we are potentially doing wrong that would be fantastic. We can certainly send you our data file if that makes it easier to assess what is going on.

Thank-you for your time and any advice,

Here is our script:

AllColoniesNoWeightOne=import.chdata(“C:/Documents and Settings/taras/Desktop/RMarkFiles/
AllColoniesNoWeightAJ.txt", field.types=c(“f”,”f”,”f”))

AllColoniesNoWeightOne.process=process.data(AllColoniesNoWeightOne,begin.time = 2007,groups=c("age","Colony","Sex"),initial.ages=c(1,0),model="Robust",time.intervals=c(0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0))

AllColoniesNoWeightOne.ddl=make.design.data(AllColoniesNoWeightOne.process)

AllColoniesNoWeightOne.ddl=add.design.data(AllColoniesNoWeightOne.process,AllColoniesNoWeightOne.ddl,"S","age",bins=c(0,1,4),right=FALSE,name="YA", replace=T)

AllColoniesNoWeightOne.ddl$S$young=0
AllColoniesNoWeightOne.ddl$S$young[AllColoniesNoWeightOne.ddl$S$YA=="[0,1)"]=1
AllColoniesNoWeightOne.ddl$S$adult=1-AllColoniesNoWeightOne.ddl$S$young

AllColoniesNoWeightOne.ddl=add.design.data(AllColoniesNoWeightOne.process,AllColoniesNoWeightOne.ddl,"GammaDoublePrime","time",bins=c(2007,2008,2009),right=FALSE,name="Gamma", replace=T)
AllColoniesNoWeightOne.ddl
AllColoniesNoWeightOne.ddl=add.design.data(AllColoniesNoWeightOne.process,AllColoniesNoWeightOne.ddl,"GammaPrime","time",bins=c(2007,2008,2009),right=FALSE,name="Gamma", replace=T)

AllColoniesNoWeightGammasPCTimeColony.models=function()
{
S.Global=list(formula=~time*Sex*YA*Colony)
S.Model1=list(formula=~Colony + time + young:time + adult:time:sex)
Gamma.dot=list(formula=~1,share=T)
p.timesession=list(formula=~-1+session:time)
c.session=list(formula=~-1+session)
N.time=list(formula=~-1+session)

cml=create.model.list("Robust")
cml
results=mark.wrapper(cml,data=AllColoniesNoWeightOne.process,ddl=AllColoniesNoWeightOne.ddl,adjust=FALSE)
return(results)
}

AllColoniesNoWeightGammasPCTimeColony.results=AllColoniesNoWeightGammasPCTimeColony.models()
AllColoniesNoWeightGammasPCTimeColony.results
AllColoniesNoWeightGammasPCTimeColony.results$results$real

Thanks Tara

[jlaake]

I'm working with Tara offlist and will post a summary once it is complete. For starters, the model she specified had completely separate parameters for capture and recapture which reduces the information in the data. In general her models were far too general and the reason for her results. --jeff

[jlaake]

Tara-

I've not modelled any robust data in practice so some others with experience may want to chime in. I definitely recommend that you read the Robust chapter in Cooch and White. Below is a list of possible models to consider but I really don't have sufficient knowledge about the situation to devise a complete set to consider. One other model that just occurred to me was ~-1+Colony:time:session+c for because if there are in different locales then the temporal pattern of capture probabilities probably will vary across colonies. A few problems you had with your original specification. First as I mentioned you had completely separate models for c (recapture probability) and while that may be necessary I'd first start off with a +c and share=TRUE for p. This will add a single parameter and allow recapture probability to be different but follow the same temporal pattern. Second, you had used Gamma instead of GammaPrime or GammaDoublePrime. It didn't recognize Gamma and simply ignored it. Also, the share=TRUE can only be on one specific parameter. For p/c that are shared it must be with the p model specificaion and for GammaPrime/GammaDoublePrime it must be with GammaDoublePrime as shown below. The GammaDoublePrime specification below creates a single parameter that is equal for both it and GammaPrime. You'll have to decide what formulation is best for your situation. For N it really should be ~-1+group:session in this case. I'll have to patch the code because the default is ~group and should be ~-1+group:session. It is important to understand that the model for N is really a model for f0 which is the number that were not seen. In your case, f0=0 for the most part which occurs when the beta is a large negative number and so most of your betas are at a boundary and can have large se because they are at a boundary. To me it makes no sense to have ~session for N because that would mean that the number that was missed in each group was the same which might only make sense if the number seen and abundance was the same in each group. You had also set adjust=FALSE which uses the parameter counts from MARK which I don't think is a good idea in this case because you have many estimated parameters that are at a boundary. Those are properly estimated parameters but MARK will not count them. While the count from RMark is not perfect either because it will include any parameters that are truly confounded, I think it is better to count too many parameters than too few. Also, you can count parameters yourself and adjust them with adjust.npar in RMark.

Hopefully this has been helpful. --jeff
AllColoniesNoWeightGammasPCTimeColony.models=function()
{
S.Global=list(formula=~-1+time:Sex:YA:Colony)
S.Model1=list(formula=~Colony + time + young:time + adult:time:Sex)
S.time=list(formula=~time)
S.time.age=list(formula=~time+YA)
S.time.age.colony=list(formula=~time+YA+Colony)
S.time.age.colony.sex=list(formula=~time+YA+Colony+Sex)
GammaDoublePrime.dot=list(formula=~1,share=TRUE)
p.timesession=list(formula=~session:time+Colony,share=TRUE,remove.intercept=TRUE)
p.timesession.c=list(formula=~session:time+Colony+c,share=TRUE,remove.intercept=TRUE)
p.timesession.sex=list(formula=~session:time+Colony+Sex,share=TRUE,remove.intercept=TRUE)
p.timesession.c.sex=list(formula=~session:time+Colony+c+Sex,share=TRUE,remove.intercept=TRUE)
N.group.session=list(formula=~-1+group:session)
cml=create.model.list("Robust")
results=mark.wrapper(cml,data=AllColoniesNoWeightOne.process,ddl=AllColoniesNoWeightOne.ddl)
return(results)
}