I have a rather complicated question, however this is my first time building a multi-strata model in RMark.
I have spatially (as opposed to temporally) designed mark recapture experiment which looks like this:
The individuals are marked and released "R", survive a "reach" to "site" 2 ("W"), then can survive another reach to site 3 where they choose a stratum: "G" or "S" and then survive another reach to site 4 ("B")
First question: can I set up a ch where the stratum states are only available at time 3? In my head this looks like this: 11G1 or 10S1 or similar?
I have assumed that this is not possible, given that I could not find a single example which attempted this, therefore started to design a model wherein at site 1 you can only be stratum R (p = 1 for all other stratums but they won't be detected), then if you survive the reach to site 2 you automatically transition to W (Psi RtoW = 1, p = 1 for all stratums but "W"), again surviving the trip, you can choose "G" or "S" stratum (p = 1 for all stratums except "G" or "S"), before the final survival and forced transition (Psi GtoB = 1, StoB = 1, p = 1 for all other stratums). This is what that looks like:
- Code: Select all
> ch
# A tibble: 8 × 2
ch freq
<chr> <int>
1 R000 10053
2 R00B 3
3 R0S0 1
4 R0SB 1
5 RW00 534
6 RW0B 90
7 RWS0 157
8 RWSB 99
# Process the data
pr.dat = process.data(ch, model = "Multistrata")
# Create default design data
# pim.type = "time" keeps ddl from creating "marking cohorts"
# for each group and time.
gs.ddl = make.design.data(pr.dat,
parameters = list(S = list(pim.type = "time"),
p = list(pim.type = "time"),
Psi = list(pim.type = "time")))
# Since we are swapping space for time, "time" in ddl$S and ddl$Psi really
# means "reach" in our mark-recapture model. And "time" in ddl$p really
# means "detection site".
gs.ddl$S$reach = gs.ddl$S$time
gs.ddl$Psi$reach = gs.ddl$Psi$time
gs.ddl$p$site = gs.ddl$p$time
# Survival for a given Stratum can only be at the site that stratum occurs,
# Set all others to 0?
S_stratum_indices = c(which(gs.ddl$S$reach == 1 & gs.ddl$S$stratum %in% c("W","G","S","B")),
which(gs.ddl$S$reach == 2 & gs.ddl$S$stratum %in% c("R","G","S","B")),
which(gs.ddl$S$reach == 3 & gs.ddl$S$stratum %in% c("R","W","B"))) %>%
unique() %>% sort()
S_stratum_values = rep(0, length(S_stratum_indices))
# Survival and detection are confounded in the last reach.
# Let's set site 4 to 1. p4 would then be interpreted as Lambda, the
# joint probability of surviving the last reach and being detected at the last site.
site4_indices = which(gs.ddl$p$site == 4)
site4_values = rep(1, length(site4_indices))
# Since Detection Probability is dependent on site, but certain Strata are not
# possible to be detected at certain sites, lets set those to 1 (perfect detection, does not exist)
p_stratum_indices = c(which(gs.ddl$p$site == 1 & gs.ddl$p$stratum %in% c("W","G","S","B")),
which(gs.ddl$p$site == 2 & gs.ddl$p$stratum %in% c("R","G","S","B")),
which(gs.ddl$p$site == 3 & gs.ddl$p$stratum %in% c("R","W","B")),
which(gs.ddl$p$site == 4 & gs.ddl$p$stratum %in% c("R","W","G","S")))
p_stratum_values = rep(1, length(p_stratum_indices))
# Since site 2, stratum W is not online, detection is by definition 0
# p_stratum_indices2 = c(which(gs.ddl$p$site == 2 & gs.ddl$p$stratum %in% c("W")))
#
# p_stratum_values2 = rep(0, length(p_stratum_indices2))
# Combine the p indices and values.
p_fixed_indices = c(site4_indices, p_stratum_indices#, p_stratum_indices2
)
p_fixed_values = c(site4_values, p_stratum_values#, p_stratum_values2
)
# Since stratum is confounded with site set stratum for fixed sites
# At These sites you can only be one stratum and only transition
# to specific stratum: (e.g. reach 1 can only be from "R" to "W") however,
# at reach 3, one can transition from "R" to "B" if one was never detected at
# site 2 ("W") (Reach 1) or site 3 ("G" or "S") (Reach 2)
psi_fixed_0indices = c(which(gs.ddl$Psi$stratum != "R" &
gs.ddl$Psi$tostratum %!in% c("W") &
gs.ddl$Psi$reach == 1), # At reach 1 only R -> W
which(gs.ddl$Psi$stratum != "W" &
gs.ddl$Psi$tostratum %!in% c("G","S") &
gs.ddl$Psi$reach == 2), # At reach 2 only W -> G or S
which(gs.ddl$Psi$stratum %!in% c("G","S") &
gs.ddl$Psi$tostratum %!in% c("B") &
gs.ddl$Psi$reach == 3), # At reach 3 only G or S to B
which(gs.ddl$Psi$stratum %in% c("B") &
gs.ddl$Psi$tostratum %in% c("R","W","G","S")), # B can't return upstream
which(gs.ddl$Psi$stratum %in% c("G","S") &
gs.ddl$Psi$tostratum %in% c("R","W")), # G and S can't return upstream
which(gs.ddl$Psi$stratum %in% c("W") &
gs.ddl$Psi$tostratum %in% c("R")), # W can't return upstream
which(gs.ddl$Psi$stratum %in% c("G","S") &
gs.ddl$Psi$tostratum %in% c("G","S")), #G and S can't switch
which(gs.ddl$Psi$stratum %in% c("R") &
gs.ddl$Psi$tostratum %in% c("B")), # R can't go to B
which(gs.ddl$Psi$stratum %in% c("W") &
gs.ddl$Psi$tostratum %in% c("B")), # W can't go to B
which(gs.ddl$Psi$stratum %in% c("R") &
gs.ddl$Psi$tostratum %in% c("G","S")), # R can't go to G or S
which(gs.ddl$Psi$stratum %in% c("R") &
gs.ddl$Psi$tostratum %in% c("R")), # R can't stay at R
which(gs.ddl$Psi$stratum %in% c("W") &
gs.ddl$Psi$tostratum %in% c("W")), # W can't stay at W
which(gs.ddl$Psi$stratum %in% c("B") &
gs.ddl$Psi$tostratum %in% c("B")) # B can't stay at B
) %>% sort() %>% unique()
psi_fixed_0values = rep(0, length(psi_fixed_0indices))
# Some transitions HAVE to occur: R to W, and G or S to B,
# Set these to 1
psi_fixed_1indices = c(which(gs.ddl$Psi$stratum %in% c("R") &
gs.ddl$Psi$tostratum %in% c("W") &
gs.ddl$Psi$reach == 1),
which(gs.ddl$Psi$stratum %in% c("G","S") &
gs.ddl$Psi$tostratum %in% c("B") &
gs.ddl$Psi$reach == 3)
) %>% sort() %>% unique()
psi_fixed_1values = rep(1, length(psi_fixed_1indices))
psi_comb_indices = c(psi_fixed_0indices,psi_fixed_1indices)
psi_comb_values = c(psi_fixed_0values,psi_fixed_1values)
# Run the Models
model = mark(data = pr.dat,
ddl = gs.ddl,
realvcv = TRUE,
model.parameters = list(S = list(formula = ~reach*stratum,
fixed = list(index = S_stratum_indices,
value = S_stratum_values)),
p = list(formula = ~site*stratum,
fixed = list(index = p_fixed_indices,
value = p_fixed_values)),
Psi = list(formula = ~-1+stratum:tostratum,
fixed = list(index = psi_comb_indices,
value = psi_comb_values)))
)
I wanted to make sure that I had correctly fixed these parameters, logically they make sense to me, but I'm a novice. Additionally, when I use the Transition Matrix, there still seem to be transitions from B to B, R to R, etc. which I was hoping to exclude as these are not possible. I am struggling to figure out how to exclude these in RMark.