I have been working with multistate models in RMark, and I have learned some nice R tricks from Jeff Laake through email. Upon his suggestion, I will try to summarize some of these, so others can have access to the information until Appendix C is updated.
1. How to constrain some of the transition rates to be the same (e.g., B to C and C to C)?
Jeff's suggestion:
This is accomplished in RMark by creating a dummy variable 0/1 in the design data (e.g., "bc"):
ddl$Psi$bc=0
ddl$Psi$bc[ddl$Psi$stratum%in%c("B","C")&ddl$Psi$tostratum=="C"]=1
Then one can use "bc" in the model formula. E.g., if you were to use the formula ~bc+time then the intercept would be for the other strata and bc would be for C and D and the time effect would be additive. (Note that in this example other strata are constrained to be the same. One can also create separate design data columns for other parameters.)
2. How to fix some of the model parameters to zero (or another value)?
R help file for make.mark.model indicates that there are a number of options for fixing parameters in RMark: one can fix certain indices, times, ages, cohorts, or groups. Fixing indices is the most flexible one, and if you want to fix, for example, some of the transition rates to zero, it is the only option. However, note that these index numbers are the design data indices (row numbers), and not the PIM indices (it was not very clear to me reading the help file). PIM indices are stacked on top of each other, but design data indices are not, which makes index identification much easier. To see this more clearly, run the following code:
data(mstrata)
mproc=process.data(mstrata,model="Multistrata")
m.ddl=make.design.data(mproc)
test=mark(mproc,m.ddl,model.parameters=list(S=list(fixed=list(index=c(1,4,6),value=0)),p=list(fixed=list(index=c(3,8),value=1))))
m.ddl$S
m.ddl$p
Here's Jeff's suggestion on how to identify index numbers in a sleek way and fix the corresponding parameters in a model:
> One of the ways you can extract them is to construct a selection on the design data. What if you could only start in stratum C but couldn't transition there. Then all the Psi for tostratum="C" from A or B would be 0. Now also consider a strange situation in which it stays in B for only one time step and automatically transitions to A so Psi from B to A =1. Below in a series of steps I show how to construct the two sets of indices and values and then create a list for the fixed Psi parameters which is in turn used in the call to mark. Note that this is all built on the previous example which created mproc and m.ddl for the mstrata data. The way I've done it below is the best way to do things because it will extract the correct indices even if you were to change the data to create more capture occasions or even strata. If you were to specify the indices specifically by number then you'ld have to regenerate the set of indices. The code below is invariant to that kind of change. Also by creating a separate myPsifixed object it can be used with any of the formula and models for that data. This type of code can be used with any of the parameters and fields in the design data.
toC.zeroes=as.numeric(row.names(m.ddl$Psi[m.ddl$Psi$tostratum=="C",]))
toC.zeroes.values=rep(0,length(toC.zeroes))
fromBtoA=as.numeric(row.names(m.ddl$Psi[m.ddl$Psi$stratum=="B"&m.ddl$Psi$tostratum=="A",]))
fromBtoA.values=rep(1,length(fromBtoA))
myPsifixed=list(index=c(toC.zeroes,fromBtoA),value=c(toC.zeroes.values, fromBtoA.values))
test=mark(mproc,m.ddl,model.parameters=list(Psi=list(fixed=myPsifixed)) )
Two useful R functions used in this code: (1) "row.names" extracts the design data indices, and (2) "as.numeric" classifies them as numeric value.
Thanks to Jeff Laake for being very kind and generous in replying to my questions and providing the example code.
Arpat