Multistate survival estimates for an empty state

[aschmidt]

Hello,
I’m attempting to run a multistate model with seabird data in RMark. The states are non-breeder (stratum 1) and breeder (stratum 2). All individuals are marked as chicks so they are all known age and start in the non-breeder state. The earliest age of breeding is 2 so individuals can transition to breeding state between age1 and age2. I’ve binned ages so that the first ageclass includes age 0-1 and the second bin is age2+. I’ve been attempting to run models with time varying survival, recapture, and transition probabilities.

I noticed that the models are generating estimates for the survival of juveniles (age 0-1) in the breeding state. Since there are no juveniles in the breeding state in the data, I don’t see how these are being generated unless there is something wrong with how I have specified the model? If they are not estimable they should show up as a constant default value right? I tried fixing these survival values to 0 and it doesn’t appear to affect the other estimates.

I’ve included the code for the model specification below. I am running several combinations of different parameters so I’m using create.model.list and mark.wrapper. This code here just has the most parameterized model I was running. Since age0 individuals can’t change state, I have fixed Psi to equal 0 for those individuals transitioning to stratum1. I’m also uncertain about the coding for Psi. I essentially want Psi to include the additive effects of ageclass, stratum, and time (like for S and p). Is the right way to do it? If you have any insight into these issues, I would greatly appreciate it. Thanks!

--AS

Code: Select all

# Design data is called BC.ddl, processed data is called BC.process
S_ast.p_ast.Psi_ast=function() {
  # S
  S.ast = list(formula = ~ageclass+stratum+time)
 
  #p
  p.ast = list(formula = ~ageclass+stratum+time)
 
  # Psi
  # create index of rows containing transitions of age0 individuals to stratum2
  Psi0.indices = as.numeric(row.names(BC.ddl$Psi[BC.ddl$Psi$Age == 0 & BC.ddl$Psi$tostratum==2,])
 
  Psi.ast <- list(formula = ~-1+ageclass:stratum:tostratum:time,
                  fixed = list(index = Psi0.indices, value = 0))
 
 
  cml=create.model.list("Multistrata")
  results=mark.wrapper(cml,data=BC.process,ddl=BC.ddl,output=FALSE)
  return(results)
}


[jlaake]

Sorry for my delayed response. I was in the field and out of contact.

You are getting an estimate for the stratum because you are using an additive model for survival. If you fix. the transition probability to 0 then the value of survival will not matter.

Your model for Psi is fine as long as you want to estimate all of those parameters which would be specific to age and time and their interaction. Let me know if you need further clarification.

--jeff

[aschmidt]

OK, thanks Jeff. I think I understand the survival issue now. As for the Psi coding, I believe what I want is the interaction between stratum and tostratum but just additive effects of ageclass and time (no interactions). Hopefully this will get me that?

Code: Select all

Psi0.indices1 = as.numeric(row.names(BC.ddl$Psi[BC.ddl$Psi$Age == 0
                                                  & BC.ddl$Psi$tostratum==2,]))
  Psi.as.PDO_NB <- list(formula = ~stratum:tostratum+ageclass+time,remove.intercept=TRUE,
                        fixed = list(index = Psi0.indices1, value = 0))


--AS

[jlaake]

It should but with any formula (excluding individual covariates) you can check by using model.matrix with the formula on your design data. It will yield a design matrix with the intercept and possibly columns that sum to zero which you can ignore. If x is the dm that model.matrix returns then x[,colSums(x)>0] will give you the dm RMark will use except that it will also remove the first column intercept with remove.intercept=TRUE.

--jeff