standard errors for Phi equal zero in multistate model

[brouwern]

I am using multistate models set up via RMark to study perennial plants. I am modeling the effects of an experimental treatment on life history transitions as well as the effects of herbivory, a binary time-varying individual covariate. The models seems to work with either the treatment effect OR the herbivory covariate, but including both in the model frequently gives me standard errors that equal zero. This occurs for both additive (treatment + herbivory) and multiplicative models (treatment x herbivory) with some variation on how I specify the model in RMark.

My focal transitions are coded "N" and "D," with "D" being a hidden state; I am therefore applying treatment and covariate effects to the term in my model [N:toD]. My understanding from general R use is that

Code: Select all

 N:toD:time:treatment:eaten

and

Code: Select all

N:toD:time + N:toD:time:eaten + N:toD:time:treatment + N:toD:time:eaten:treatment

should result in similar if not identical models. However, the number of parameters with SE = 0 seems to depend on how I specify my models. Example RMark code is given below.

Are there any inherent issues with having multiple factors affect a transition within multistate model or is this likely an estimability issue related to my dataset? Ideally I'd like to build a single model that compares the effects of ambient levels of herbivory (which is fairly common) with our experimental treatment.

Some more details:
-As noted above, "D" is a hidden state with p = 0. Above-ground states ("Y" and "N") have p = 1.
-Because I am interested in the dynamics of the hidden state I constrain survival to be constant over time.
-N is about 200 individuals.

Treatment x Herbivory effect, version 1

Code: Select all

Psi.1  <-  list(formula = ~-1 + 
                          N:toY:time  +                 
                          N:toD:time:treatment:eaten +
                          Y:toY:time  +
                          Y:toN:time  +
                          D:toN:time  +
                          D:toD:time)   


  Treatment x Herbivory effect, version 1
  Psi.2  <-  list(formula = ~-1 +   #Remove intercept (?)
                            N:toY:time  +                 
                            N:toD:time + N:toD:time:eaten + N:toD:time:treatment + N:toD:time:eaten:treatment  +
                            Y:toY:time  +
                            Y:toN:time  +
                            D:toN:time  +
                            D:toD:time)

[egc]

This is an RMark question -- so is being moved to that subforum.

[jlaake]

Hard to say exactly. Where the parameters with se=0 at boundaries? From looking at the formulas with N=200 I'd say that it was likely over-parameterized. If these are plants I'm wondering why you are using MARK? Don't you always know where they are at?

Was your subject header supposed to be Psi instead of Phi?

Purely on a formula basis in R, it is x*y*z=x+y+z+x:y+x:z+y:z+x:y:z. Strictly : specifies an interaction and * includes main effects and interactions. You seem to be using : with two meanings. Now if all of the variables are factors then it is somewhat interchangeable. x*y=-1+x:y because x:y will create each combination.

In the models you are describing you are mixing numeric with factor variables and in that case the : is just multiplicative for the numeric variable. I assume eaten=herbivory which is 0/1 but you didn't say what treatment was - 0/1 numeric or factor. Regardless the 2 formula you specified are not the same because at least one is numeric.

I suggest you add a dummy variable called eaten with each value =1 in the design data. Then use model.matrix with your formula and the design data frame and look at the resulting design matrix so you can see and understand what it is doing. For individual covariates, that is what RMark does. It adds a dummy variable with all 1s for the individual covariate in the design data, constructs the model.matrix and then replaces it with the covariate name.

--jeff