Interaction terms for nest survival in RMark

[Disco]

Greetings all,

I have re-read relevant chapters in the MARK book, RMark appendix, RMark guide, postings in the forums, and relevant literature but I am not finding the clarification I need about using interaction terms specifically in RMark nest survival models.

I wish to test how DSR varies with several covariates: Site, Habitat, TimeTrend, NestAge, and Year. Site and Habitat should autocorrelate as there are only 4 sites among 3 habitat types. I also suspect a relationship between Site & Habitat with Year because I did not sample all Sites/Habitats in each year of the study. To account for this, it is my understanding that I should include partially-crossed interaction terms for Site, Habitat, & Year (all are factor variables) to deal with the missing data.

I can get the fully crossed functions to work fine, for example:

Code: Select all

S=list(formula=~Site*Year)

However, the partial crosses yields unusually large delta AICs, for example:

Code: Select all

S=list(formula=~-1+Site:Year)

I saw that a suggested solution for MARK users was to use dummy codes, but it was my understanding that this wasn't necessary for RMark. If so, how do I proceed?

Any help would be greatly appreciated!
Deb

[jlaake]

We discussed this offlist. The problem was that Year and Habitat etc had been coded numerically and were being treated as a numeric individual covariates rather than as factor variables and used to define groups. Beware of this problem. RMark will tell you if you try to use a factor variable as a numeric covariate but if it is numeric in the data even if it is a factor variable it will proceed and use it as an numeric covariate. The difference is equivalent to a p(T) model (trend) versus a p(t) (varying by occasion) model. For something like time this can make sense but for something like habitat where there is no natural ordering it makes no sense but RMark cannot know.

--jeff

[Disco]

Many thanks for your help Jeff!

I have corrected the factor variable/grouping problem, but now I'm getting some strange results. In one case, a model including an interaction term is producing identical AIC values as one of the factors alone:

Code: Select all

 Site=mark(proc.myrexs,ddl,nocc=144,model="Nest",model.parameters=list(S=list(formula=~Site)), groups=c("Site","Habitat","Year"))
SiteHab=mark(proc.myrexs,ddl,nocc=144,model="Nest",model.parameters=list(S=list(formula=~-1+Site:Habitat)), groups=c("Site","Habitat","Year"))
 

Similarly, a three-way interaction is producing the exact same AIC values as a two-way interaction:

Code: Select all

SiteHabYear=mark(proc.myrexs,ddl,nocc=144,model="Nest",model.parameters=list(S=list(formula=~-1+Site:Habitat:Year)), groups=c("Site","Habitat","Year"))
SiteYear=mark(proc.myrexs,ddl,nocc=144,model="Nest",model.parameters=list(S=list(formula=~-1+Site:Year)), groups=c("Site","Habitat","Year"))
 

I've checked the manual, and it is my understanding that these aren't situations where I would use the remove.intercept argument, but perhaps I have missed something. It has also been very difficult to find information on how to handle three-way interactions, so I may have done something wrong there. Note that I am not using the full-cross operator * because I do not have all levels of Year represented at each Site, Habitat, and Year combination.

Can anyone shed some light onto the problem?
Thanks,
Deb

[jlaake]

How many unique site-habitat-year combinations do you have? Certainly site:habitat will produce the same as site because each site has only one type of habitat unless site changes habitat over time. Also site:habitat:year will produce the same thing as site, if each site was only sampled in one year.

You need to sit down and think about what these model formula do and what you are trying to ask with your data. This is where Evan has a very good point about RMark. When you use MARK interface it forces you to create the design matrix by hand which makes you fully understand the relationships you are proposing for the parameters. You can still do that with RMark by looking at the design matrix it creates or use model.matrix(~formula,ddl$parameter) to look at it. I use an analogy when I teach a workshop. If you were nailing pieces of wood, MARK is a hammer and if you are precise you can do a good job albeit slowly. RMark is a nail gun. It will get the job done quickly, but if you aren't careful you can create quite a mess quickly and hurt yourself. So make sure to read and understand the operating instructions.

Please don't take this post the personally. You are certainly not the first nor the last to make this mistake. If I'm incorrect about the assumptions with regard to site and year please do correct me.

regards --jeff

[Disco]

Jeff,

No personal offensive taken - although I did have to lick my nail gun wounds!

I indeed was caught up in thinking I made a coding mistake rather than a theoretical one. I have 9 unique combinations of site and year, and habitat does not add any more parameters.

Thanks for the push in the right direction, and the reminder that its time to re-RTFM.

Best,
Deb

[jlaake]

I hope nail gun wounds were not to deep.

So for clarity for me and other readers, the problem was inclusion of habitat. With multiple years for the same site, please confirm that you were getting differences between site and site:year models but no difference between site:year and site:habitat:year. That would make sense to me because each site doesn't change habitats.

Using model.matrix with your formula and design data is a good way to learn and understand what it is doing. However, note that model.matrix can create columns of all zeros. They are not used and RMark excludes them.

regards --jeff

[Disco]

That is correct. Site and Site:Year were different. The addition of habitat (Site:Habitat:Year) didn't change the results from the Site:Year model because there are 4 Sites and only 4 possible Site/Habitat combinations. Because the Habitat:Year model had more weight than the Site:Year alone, I ended up using Habitat:Year to build subsequent models.

Thanks again for your help,
Deb