www.phidot.org

My closed robust design multistrata model has 5 strata, including "U", the unobserved stratum. I'd like one of my models to have transitions to and from "U" be constant, while having another set of transitions (ones that make no sense biologically) fixed at zero, and have all remaining transitions be estimated. I can't think of a way to have these three different definitions of Psi in the same model.

Here's where I am so far:

I have been trying to create grouping indices for the U transitions (object name: PsiU), and for the nonsensical transitions (object name: Psi0), below is the code i've been trying to use to create Psi0:

Code: Select all

Psi0=as.numeric(row.names(crdms2.ddl$Psi[crdms2.ddl$Psi$stratum=="A"&crdms2.ddl$Psi$tostratum=="H",]),          row.names(crdms2.ddl$Psi[crdms2.ddl$Psi$stratum=="A"&crdms2.ddl$Psi$tostratum=="J",]),          row.names(crdms2.ddl$Psi[crdms2.ddl$Psi$stratum=="A"&crdms2.ddl$Psi$tostratum=="S",]),             row.names(crdms2.ddl$Psi[crdms2.ddl$Psi$stratum=="H"&crdms2.ddl$Psi$tostratum=="A",]),             row.names(crdms2.ddl$Psi[crdms2.ddl$Psi$stratum=="H"&crdms2.ddl$Psi$tostratum=="S",]),             row.names(crdms2.ddl$Psi[crdms2.ddl$Psi$stratum=="J"&crdms2.ddl$Psi$tostratum=="A",]),             row.names(crdms2.ddl$Psi[crdms2.ddl$Psi$stratum=="J"&crdms2.ddl$Psi$tostratum=="H",]),             row.names(crdms2.ddl$Psi[crdms2.ddl$Psi$stratum=="S"&crdms2.ddl$Psi$tostratum=="H",]),             row.names(crdms2.ddl$Psi[crdms2.ddl$Psi$stratum=="S"&crdms2.ddl$Psi$tostratum=="J",]))

the code is taken from someone else's script (see https://github.com/jlaake/RMark/blob/master/RMark/R/RMark-package.R and scroll down to the part about mulitistrata models) and they only used it to index a single transition. when i add more transitions, the code still runs but it only actually includes the first transition listed in the grouping index. So if you have a fix for that, i'd love to see it. I suppose i could create a set of indices and just list them in the parameter definition statement below, but there must be a simpler way.

I am able to define Psi so that the Psi0 transitions (or at least the one of them that is actually included in Psi0) are set to zero while estimating all other transitions as follows:

Code: Select all

Psistos.specfix=list(formula=~1+stratum:tostratum,fixed=list(index=Psi0,value=0))


It seems like it shouldn't be possible to add in the PsiU transitions held constant, but I am hoping I'm wrong and someone else knows how to do this.

Thanks for any help you can give!

It is quite possible but it has nothing to do with fixing the real parameters. Copying others work can be useful but it doesn't lead to your full understanding. Also while there is nothing wrong with the way you fixed real parameters, there is a better way that was created a few versions ago.

See http://www.phidot.org/forum/viewtopic.php?f=21&t=2887 which is the first fixed post on this sub-forum.

Now to your question. You have specified the formula as ~1+stratum:tostratum. First of all it should have been ~-1+stratum:tostratum to remove the intercept. stratum:tostratum specifies that the probability of moving from one stratum to another are all different for both the stratum and the tostratum. Thus it would not allow you to fix movement to and from U to be constant. One might think you should use ~stratum*tostratum; however, with the default subtract.stratum values the parameters differ across strata. In other words A can go to B and C, but B can only go to A and C because AtoA and BtoB are computed by subtraction. Thus, stratum and tostratum are not fully crossed and using ~stratum*tostratum will create too many parameters. But ~1+stratum:tostratum will create one too many as well. Whereas, ~-1+stratum:tostratum will remove the intercept and create a separate parameter that is specific to the stratum (e.g, there will be a AtoB, AtoC, BtoA, BtoC etc). You can get a better understanding of this if you use the R function model.matrix with your formula and design data for Psi. This is all fully explained in the 2015 RMark workshop notes and in other places in the documentation.

So how do you do it. I haven't tried this but the following should work. Create a variable in ddl$Psi called notU and assign it the value 0 unless stratum==U or tostratum==U with something like

Code: Select all

ddl$Psi$notU=ifelse(ddl$Psi$stratum=="U" | ddl$Psi$tostratum=="U",1,0)

Then use the formula ~notU:stratum:tostratum. The intercept should be the constant probability of moving from or to U and the remaining parameters are relative to the intercept. If you don't want the other parameters to be relative (treatment contrast) to the intercept, create another variable called tofromU as

Code: Select all

ddl$Psi$tofromU=1-ddl$Psi$notU

then use the formula ~-1+tofromU+notU:stratum:tostratum.

Try all of this out with the molde.matrix function to look at the results to help your understanding and to make sure my logic is correct.