Hello,
I'm building a multistate model in RMark where some states are **transient** — they can only be occupied for a single occasion (the first month after capture). The model works correctly in Nimble (Bayesian), but I'm getting biased survival estimates in RMark and I believe the issue is related to the mlogit parameterisation of Psi.
## What I'm trying to do
I'm analysing GPS telemetry data (180 individuals, 77 deaths, 103 right-censored, ~80 monthly occasions, perfect detection via GPS + mortality sensors). I want to estimate monthly survival while separating a first-month post-capture effect from baseline survival, and then test whether processing time and relative tag weight affect survival (first month, subsequent, or both).
I use 6 states:
- **A (Y1)**: Juvenile, 1st month post-capture → **transient, cannot be occupied twice**
- **B (Y)**: Juvenile, subsequent months
- **C (A1)**: Adult, 1st month post-capture → **transient, cannot be occupied twice**
- **D (A)**: Adult, subsequent months
- **E (JD)**: Just dead (detected)
- **F (D)**: Long dead (absorbing, p = 0)
For state A (Y1), if the bird survives, it **must** transition to either B (stay juvenile) or D (mature to adult). It can never remain in A. Similarly, state C (A1) must transition to D. The self-transitions A→A and C→C have probability 0 by design.
## The mlogit problem
In MARK's multistate model, I understand that Psi is parameterised with a multinomial logit where "remain in the current state" is the reference category:
```
P(A→B) = exp(b1) / (1 + exp(b1) + exp(b2))
P(A→D) = exp(b2) / (1 + exp(b1) + exp(b2))
P(A→A) = 1 / (1 + exp(b1) + exp(b2)) ← reference, must be 0
```
Since P(A→A) must equal 0, both b1 and b2 are pushed to +∞. Any pair of large values with the right ratio gives the same likelihood, creating a ridge. This results in:
```
estimate se
Psi:stratumA:tostratumB 14.11 664.96
Psi:stratumA:tostratumD 15.03 472.18
```
The `$real` output shows both Psi ≈ 1.0 (summing to ~2.0):
```
Psi sA toB: 0.9999993
Psi sA toD: 0.9999997
```
I tried `ddl$Psi$fix[ddl$Psi$stratum == "A" & ddl$Psi$tostratum == "A"] <- 0` but this doesn't change the mlogit reference — it just forces the reference probability to 0, which still requires b1, b2 → +∞.
**I believe this numerical instability is corrupting the survival estimates.** With `S ~ stratum` (one free S per state), RMark gives S(A) = 0.952 instead of the expected 0.854 from raw data (6 deaths / 41 at risk). The same model in Nimble — where the transition matrix is specified directly as probabilities summing to 1 — gives S(A) ≈ 0.85, matching the raw data.
## What I need
I need a way to tell RMark/MARK that for states A and C, the self-transition probability is structurally 0, and the remaining transitions should sum to 1 without using the self-transition as the mlogit reference. In other words:
- For state A: `P(A→B) + P(A→D) = 1` (estimated, no reference to A→A)
- For state C: `P(C→D) = 1` (already fixed)
- For state B: `P(B→B) + P(B→D) = 1` (normal mlogit, this one works fine)
## Questions
1. **Is there a way in RMark to change the mlogit reference category for Psi**, so that for state A, one of the real transitions (e.g., A→B) is the reference instead of the impossible A→A?
2. **Is there another way to parameterise Psi in RMark for transient states** — for example, using a logit link instead of mlogit for states with only two possible transitions?
4. Alternatively, **is there a completely different model structure in MARK** that would be more appropriate for estimating a first-month capture effect with known-fate telemetry data?
## My current RMark code
```r
dp <- process.data(KF_ms, model = "Multistrata",
strata.labels = c("A","B","C","D","E","F"),
groups = c("Species"))
ddl <- make.design.data(dp)
# S = 1 for dead states
ddl$S$fix <- NA
ddl$S$fix[ddl$S$stratum %in% c("E", "F")] <- 1
# Perfect detection
ddl$p$fix <- NA
ddl$p$fix[ddl$p$stratum %in% c("A", "B", "C", "D")] <- 1
ddl$p$fix[ddl$p$stratum == "E"] <- 1
ddl$p$fix[ddl$p$stratum == "F"] <- 0
# Transitions
ddl$Psi$fix <- 0
ddl$Psi$fix[ddl$Psi$stratum == "A" & ddl$Psi$tostratum == "A"] <- 0 # impossible self-transition
ddl$Psi$fix[ddl$Psi$stratum == "A" & ddl$Psi$tostratum == "B"] <- NA # estimated
ddl$Psi$fix[ddl$Psi$stratum == "A" & ddl$Psi$tostratum == "D"] <- NA # estimated
ddl$Psi$fix[ddl$Psi$stratum == "B" & ddl$Psi$tostratum == "D"] <- NA # estimated
ddl$Psi$fix[ddl$Psi$stratum == "C" & ddl$Psi$tostratum == "D"] <- 1 # deterministic
ddl$Psi$fix[ddl$Psi$stratum == "E" & ddl$Psi$tostratum == "F"] <- 1 # deterministic
ddl$S$cap_eff <- ifelse(ddl$S$stratum %in% c("A","C"), 1, 0)
mod <- mark(dp, ddl,
model.parameters = list(
S = list(formula = ~ stratum),
p = list(formula = ~ 1),
Psi = list(formula = ~ -1 + stratum:tostratum, link = "mlogit")
))
```
## Nimble version (works correctly)
For reference, here is how the same transitions are specified in Nimble, where the issue does not arise because probabilities are defined directly:
```r
# For state A (Y1): if survives, go to B or D
ps[1, 1, i, t] <- 0 # can't stay in A
ps[1, 2, i, t] <- S[i, t] * (1 - psi_mature) # A → B
ps[1, 4, i, t] <- S[i, t] * psi_mature # A → D
ps[1, 5, i, t] <- 1 - S[i, t] # A → JD (death)
# For state B (Y): stay juvenile or mature
ps[2, 2, i, t] <- S[i, t] * (1 - psi_mature) # B → B
ps[2, 4, i, t] <- S[i, t] * psi_mature # B → D
ps[2, 5, i, t] <- 1 - S[i, t] # B → JD
# Nimble result: S(A) ≈ 0.85, consistent with raw data (6/41 deaths)
# RMark result: S(A) = 0.952, inconsistent
```
Any help or pointers would be greatly appreciated. Thanks!
