mark model function

[megpetrie]

Hello,
I'm using the model function presented in C.13 of the RMark appendix. It works great except when I try to input formulas for more than one interaction. (year, hatchery, and reartype are factor variables). Here is my code:

sh.models=function()
{
Phi.timeplusyear=list(formula=~time+year)
Phi.timeplusyearplushatchery=list(formula=~time+year+hatchery)
Phi.timeplusyearpluspop=list(formula=~time+year+pop)
Phi.timeplusyearplusreartype=list(formula=~time+year+reartype)
Phi.yearpluspop=list(formula=~year+pop)
Phi.timetimesreartypeplusyear=list(formula=~time*reartype+year)
Phi.timetimeshatcheryplusyear=list(formula=~time*hatchery+year)
Phi.timetimeshatchery=list(formula=~time*hatchery)
Phi.timetimeshatcheryplustimetimesyear=list(formula~time*hatchery+time*year)
Phi.timetimeshatcherytimesyear=list(formula~time*hatchery *year)
Phi.timetimeshatcheryplushatcherytimesyear=list(formula~time*hatchery+hatchery*year)
p.timeplusyear=list(formula=~time+year)
cml=create.model.list("CJS")
results=mark.wrapper(cml,data=sh.process,ddl=sh.ddl,adjust=FALSE)
return(results)
}

When I call the model results the only models I get results for are the models with no or only one interaction. This is the model table I get:

Phi p model npar AICc DeltaAICc weight Deviance
6 ~time * hatchery + year ~time + year Phi(~time * hatchery + year)p(~time + year) 18 1854.127 0.00000 9.980680e-01 513.0779
5 ~time * hatchery ~time + year Phi(~time * hatchery)p(~time + year) 16 1866.622 12.49448 1.932047e-03 529.6859
3 ~time + year + pop ~time + year Phi(~time + year + pop)p(~time + year) 17 1897.057 42.93029 4.752992e-10 558.0666
2 ~time + year + hatchery ~time + year Phi(~time + year + hatchery)p(~time + year) 14 1906.047 51.91958 5.308433e-12 573.2112
8 ~year + pop ~time + year Phi(~year + pop)p(~time + year) 16 1911.873 57.74568 2.882999e-13 574.9371
7 ~time * reartype + year ~time + year Phi(~time * reartype + year)p(~time + year) 16 1924.208 70.08038 0.000000e+00 587.2717
4 ~time + year + reartype ~time + year Phi(~time + year + reartype)p(~time + year) 13 1961.761 107.63377 0.000000e+00 630.9707
1 ~time + year ~time + year Phi(~time + year)p(~time + year) 13 1977.360 123.23257 0.000000e+00 646.5695

When I run these >1-interaction models outside the function they return results. Is there something I need to add to the function so that it recognizes these formulae?
As a side note, I also don't get results for formulas with fixed parameters.

Thanks and let me know if you need more information.

Megan

[bacollier]

Megan,
Not sure on this being the issue ,but do you intend your models to be the interaction between just those 2 variables? e.g., in R speak time*reartype actually means time + reartime + time*reartype, or did you just want time*reartype? I wonder if part of the problem is that, for instance you have time*hatchery + hatchery*year which in R would (I think be interpreted as): time + hatchery + time*hatchery + hatchery + year + year*hatchery which is obviously not right ( I have not tested this in RMark but I cannot see how it would work)?

So, try replacing the " * " with ":" for the interactions of interest where you don't want the main effects as well (I suspect the interaction is what you are interested in)?

Also, and I am not sure on this, you might try not using create.model.list() and mark.wrapper() in your function and just run them straight through and see if that works.

Bret

[jlaake]

Megan-

I've been in the field away from email. There is a small typo in Bret's email.

time*reartype = time + reartime + time:reartype

Now the confusing part is that if you have x:y and both are factor variables then you can also use ~-1+ x:y to get the full interaction model but the difference is that with this latter method, each parameter is specified separately whereas with time*reartype the model is specified as an intercept + main effects + interactions. The -1+x:y approach is needed if you want to use an identity DM with say the the sin link. An identity DM can often have better convergence but then you can interpret the parameters as additive amounts. That was just for clarification and I don't think it has anything to do with your problem. When you run a sequence of models like this sometimes the error messages are lost in the output stream on the screen. I suggest you add the run=FALSE in the mark.wrapper. This will run through the list of models and will construct them but not run them. That will help to identify whether it is constructing each model as Bret suggested and if there was any problem in constructing it. What puzzles me is that you said that you ran these models separately and they ran but I'm not sure what that means. If it runs a model and it doesn't converge then it will not be added to the results list so this could also be a convergence issue. Dig around some more and let me know what you find. If you can't figure it out, send me the code and data and I'll work it out.

--jeff