Multi-state models with constraints on transitions

[Andre23]

I am running multi-state analysis on my data with 3 states: Pre-breeder (coded as state "1") - Breeder ("2") and Post-breeder ("3").
Individuals' age are separated into 3 classes for Survival: φ age(0, 1-3, 4+)

Transition matrix:
- Individuals transition from 1 to 2 only from age 2 minimum
- Individuals transition from 2 to 3 only from age 3 minimum
- Individuals transition back to 2 from stratum 3 only at age 4 minimum (has to go to post-breeder state first at age 3, then can go back to breeding state at age 4)
- Individuals can't go back to state "1" once they reached state "2" or "3"

Fixing detection:
- Individuals can't be seen at breeder state (2) before the age of 2

You can see in the code below all the fixed impossible transitions to 0, and the fixed detection at 0:

Code: Select all

# 2. Process the data----
ms_proc <- process.data(rmark_data, model = "Multistrata")

# 3. Create design data----
ms_ddl <-  make.design.data(ms_proc);str(ms_ddl)
{
# Create age classes for survival "S"
ms_ddl$S$ageclass3=cut(ms_ddl$S$Age,c(0, 1, 4, Inf), right =FALSE)
levels(ms_ddl$S$ageclass3)=c("Juveniles","Adults(1-3)","Adults(4+)")
# Create an indicator for juveniles (age0)
ms_ddl$S$is_juvenile <- ifelse(ms_ddl$S$ageclass3 == "Juveniles", 1, 0)
# Create an indicator for young adults (age1-3)
ms_ddl$S$youngA <- ifelse(ms_ddl$S$ageclass3 == "Adults(1-3)", 1, 0)
# Create an indicator for older adults (age4+)
ms_ddl$S$olderA <- ifelse(ms_ddl$S$ageclass3 == "Adults(4+)", 1, 0)

# Create age classes for re-capture "p"
ms_ddl$p$ageclass6=cut(ms_ddl$p$Age,c(-Inf,2,3,4,5,6,Inf),right=FALSE)
levels(ms_ddl$p$ageclass6) <- c("1","2","3","4","5","6+") # 6 age classes
#Fix p = 0 at Age 1 for breeding state
ms_ddl$p$fix[ms_ddl$p$stratum == "2" &
                 (ms_ddl$p$Age == 1)] <- 0
ms_ddl$p$fix[ms_ddl$p$stratum == "3" &
               (ms_ddl$p$Age < 3)] <- 0

# Create constraint matrix for transitions "Psi"
ms_ddl$Psi$fix <- NA
# Set impossible transitions to 0
ms_ddl$Psi$fix[ms_ddl$Psi$stratum == "1" & ms_ddl$Psi$tostratum == "3"] <- 0  # PB to PO
ms_ddl$Psi$fix[ms_ddl$Psi$stratum == "2" & ms_ddl$Psi$tostratum == "1"] <- 0  # B to PB
ms_ddl$Psi$fix[ms_ddl$Psi$stratum == "3" & ms_ddl$Psi$tostratum == "1"] <- 0  # PO to PB
#PB to B transition for age < 2 = 0
ms_ddl$Psi$fix[ms_ddl$Psi$stratum == "1" &
                 ms_ddl$Psi$tostratum == "2" &
                 (ms_ddl$Psi$Age < 2)] <- 0
#B to PO transition for age < 3 = 0
ms_ddl$Psi$fix[ms_ddl$Psi$stratum == "2" &
                 ms_ddl$Psi$tostratum == "3" &
                 (ms_ddl$Psi$Age < 3)] <- 0
#PO to B transition for age < 4 = 0
ms_ddl$Psi$fix[ms_ddl$Psi$stratum == "3" &
                 ms_ddl$Psi$tostratum == "2" &
                 (ms_ddl$Psi$Age < 4)] <- 0

# 4. Run model
ms_model <- mark(ms_proc, ms_ddl,
                 model.parameters = list(
                   S = list(formula = ~-1 + stratum:youngA + stratum:olderA + is_juvenile:time),
                   p = list(formula = ~stratum:ageclass6),
                   Psi = list(formula = ~-1+stratum:tostratum:age))
)

I would like to add constraints on transitions. Transition probability likely depends on the age of individuals but could also be impacted by resource availability ("SF") by changing their probability of going to a breeder state. Also, density could have a negative impact on recruitment, so it could impact transition probability from state 1 to state 2.

1. Is my way of fixing impossible transitions to 0 the correct way of doing it in Rmark?
2. When running the model, how should I code to test for the impact of SF or density on a specific transition (1 to 2)?

Thank you for your help

[jlaake]

1. Is my way of fixing impossible transitions to 0 the correct way of doing it in Rmark?

That is how it is done but you should look at the ddl$Psi to make sure you did what you wanted correctly. I didn't check. You can use R functions like apply or table to summarize fix by stratum and tostratum.

2. When running the model, how should I code to test for the impact of SF or density on a specific transition (1 to 2)?

You'll find that the ddl$Psi has a column for each stratum (1,2,3) and tostratum (to1,to2,to3). You can create a new column called 1to2 by multiplying columns 1 and to2. Then use 1to2:SF in your formula and it will limit the effect of SF to transitions from 1 to 2. Using a 0/1 field is useful way to limit the effects of a variable to certain subset of the parameters. If 71 year old brain's memory is correct, this is described somewhere in the documentation.

Now you should realize that Psi parameters sum to 1. Thus transitions from 1 to 1 (staying in pre-breeder) and 1 to 2 will be affected by SF. Now in this case this is fine and you want that. But because they have to sum to 1 in other circumstances the covariate will affect other real values. In your case it is only affecting the probability of staying a pre-breeder which you want.

[Andre23]

Thank you @jlaake !

I created specific transition columns in the ddl as you suggested:

Code: Select all

ms_ddl$Psi$trans_1to2 <- ms_ddl$Psi$`1` * ms_ddl$Psi$to2
ms_ddl$Psi$trans_2to3 <- ms_ddl$Psi$`2` * ms_ddl$Psi$to3
ms_ddl$Psi$trans_3to2 <- ms_ddl$Psi$`3` * ms_ddl$Psi$to2

Code: Select all

ms_model <- mark(ms_proc, ms_ddl,
                 model.parameters = list(
                   S = list(formula = ~-1 +
                              # Juvenile survival only in state 1
                              is_juvenile:time +
                              # Age 1 only in state 1
                              age1 +
                              # Age 2 in states 1 and 2
                              age2:stratum +
                              # Age 3 in all states
                              age3:stratum +
                              # Older adults (4+) in all states
                              olderA:stratum),
                   p = list(formula = ~stratum:ageclass6),
                   Psi = list(formula =  ~-1 + stratum:tostratum +trans_1to2:SF))

)

The model runs and the output shows the SF effect only for "trans_1to2" (beta estimates). It seems to work.

1. However, I noticed that my transitions probabilities do not sum to 1 (Psi s1 to2 = 0.17, Psi s2 to3 = 0.23, Psi s3 to2 = 0.37 / sum = 0.77)

2. I also can't see the transitions to "remain in the same state" in the output (e.g. 1 to 1, or 2 to 2?)

Am I missing something here? Could it be because of the way I fixed my transitions (above)?

Thanks

[jlaake]

MARK doesn't provide the probability that is computed by subtraction for the transition that is left out. This is typically the probability of remaining but can be changed to any stratum by using the subtract.stratum argument in make.design.data for Psi parameter. Use the function TransitionMatrix to compute all of the transitions.

[Andre23]

I tried using the subtract.stratum argument in make.design.data:

Code: Select all

ms_ddl <-  make.design.data(ms_proc,parameters=list(Psi=list(subtract.stratum=c("1","2","3"))))


Then using the TransitionMatrix function to compute all of the transitions:

Code: Select all

Psilist=get.real(mstrata.results[[4]],"Psi",vcv=T)
Psivalues=Psilist$estimates
TransitionMatrix(Psivalues[Psivalues$time==1,],vcv.real=Psilist$vcv.real)


The TransitionMatrix function gives me an empty output like this:

$TransitionMat
1 2 3
1 1 0 0
2 0 1 0
3 0 0 1

$se.TransitionMat
1 2 3
1 0 0 0
2 0 0 0
3 0 0 0

$lcl.TransitionMat
1 2 3
1 NaN NaN NaN
2 NaN NaN NaN
3 NaN NaN NaN

$ucl.TransitionMat
1 2 3
1 NaN NaN NaN
2 NaN NaN NaN
3 NaN NaN NaN

[jlaake]

Look at the intermediary results that you are getting. All of your previous code didn't indicate that you were using mark.wrapper with a marklist name of mstrata.results. It used a single call to mark with a resulting object ms_model. I can only guess that is the problem. To work out problems ALWAYS look at each step in the code to see where things went wrong.

Also subtract.stratum is only needed when you want to use something other than the default. I'm assuming your strata are "1","2","3" in the capture history in which case what you did was not needed because those are the default values.

[jhines]

1. However, I noticed that my transitions probabilities do not sum to 1 (Psi s1 to2 = 0.17, Psi s2 to3 = 0.23, Psi s3 to2 = 0.37 / sum = 0.77)

Transitions from each state should sum to 1. So, Psi s1 to1+Psi s1 to2+Psi s1 to3 should be 1. Same with Psi s2 to1+Psi s2 to2+Psi s2 to3 and Psi s3 to1+Psi s3 to2+Psi s3 to3. That is why Psi s1 to1 is not an estimated parameter. It can be computed as 1 - Psi s1 to2-Psi s1 to3.

The sum over different "from" strata do not need to sum to 1.

[Andre23]

Thank you all for the replies. Very helpful!