Fixing Multistrata values in a spatial CJS analysis

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

Fixing Multistrata values in a spatial CJS analysis

Postby tspaulding » Thu Dec 22, 2022 2:03 pm

Hello,
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:
Image

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.
tspaulding
 
Posts: 4
Joined: Wed Nov 04, 2020 6:35 pm

Re: Fixing Multistrata values in a spatial CJS analysis

Postby jlaake » Thu Dec 22, 2022 9:51 pm

So if I understand correctly most of your transitions are 0 and if a critter lasted throughout the experiment the only possible true values are
RWSB or RWGB. One problem that I see is that the data do not contain a G observation. I assume this is because p=0. Regardless if a stratum
does not show up in the data then it must be specified to process.data. So you should do the following:

Code: Select all
pr.dat=process.data(ch,model="Multistrata",strata.labels=c("B","G","R","S","W"))

Now you need to specify which strata transitions you want to specify for subtraction because the Psi must sum to 1 and one in each set is
computed as 1-sum of the others.

I suggest

Code: Select all
gs.ddl = make.design.data(pr.dat,
                          parameters = list(S = list(pim.type = "time"),
                                            p = list(pim.type = "time"),
                                            Psi = list(pim.type = "time",subtract.stratum=c("B","B","W","S","G"))))

You can see what this does by doing

Code: Select all
> with(gs.ddl$Psi,table(stratum,tostratum))
       tostratum
stratum B G R S W
      B 0 3 3 3 3
      G 0 3 3 3 3
      R 3 3 3 3 0
      S 3 3 3 0 3
      W 3 0 3 3 3



You can see which are missing in each row. A row is the stratum and the column is the tostratum (the stratum into which it will transition).

Now since most are 0 and we have set the ones for subtraction to be fixed to 1 as done by subtraction then you should do the following

Code: Select all
gs.ddl$Psi$fix=0
gs.ddl$Psi$fix[gs.ddl$Psi$stratum=="W"&gs.ddl$Psi$tostratum=="S"]=NA

Now you can see that the only ones that are estimated (NA value) are W to S and then W to G is computed by subtraction.
Code: Select all
> with(gs.ddl$Psi,table(stratum,tostratum,fix,useNA="ifany"))
, , fix = 0

       tostratum
stratum B G R S W
      B 0 3 3 3 3
      G 0 3 3 3 3
      R 3 3 3 3 0
      S 3 3 3 0 3
      W 3 0 3 0 3

, , fix = NA

       tostratum
stratum B G R S W
      B 0 0 0 0 0
      G 0 0 0 0 0
      R 0 0 0 0 0
      S 0 0 0 0 0
      W 0 0 0 3 0


Then I ran the default model - model=mark(pr.dat,gs.ddl) and got the TransitionMatrix

Code: Select all
Psilist=get.real(model,"Psi",vcv=TRUE)
Psivalues=Psilist$estimates
> TransitionMatrix(Psivalues[Psivalues$time==1,])
  B         G R         S W
B 1 0.0000000 0 0.0000000 0
G 1 0.0000000 0 0.0000000 0
R 0 0.0000000 0 0.0000000 1
S 0 0.0000000 0 1.0000000 0
W 0 0.7298889 0 0.2701111 0

I can't take it any further than that because I'm not sure if that was your real data or not. If they are then p should be fixed to 0 for G and
that might help but the resultsfrom the default model din't look great with large std errors.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA


Return to RMark

Who is online

Users browsing this forum: Google [Bot] and 7 guests