Using TransitionMatrix with logit link for Psi in RMark

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

Using TransitionMatrix with logit link for Psi in RMark

Postby jbauder » Wed Jan 15, 2020 5:10 pm

Hello,

I am using a multi-state model to estimate cause-specific mortality rates for nuisance black bears. I have constrained my survival and transition probabilities such that the transition probabilities from the "Alive" state to either a "Harvest Mortality" or "Non-Harvest Mortality" are the cause-specific mortality probabilities. Under this parameterization, Survival is 1 minus the sum of all mortality probabilities. I would like to calculate SEs and CIs for Survival and saw that the TransitionMatrix function in RMark can do just that. However, on page C-86 of the RMark chapter in the Gentle Introduction, it states that TransitionMatrix "...only works if you use the mlogit link for..." Psi. I have used the logit link for Psi in some models to enforce constraints among mortality probabilities (e.g., allow Harvest Mortality to vary by sex and have Non-Harvest Mortality constant). I found that TransitionMatrix does work with my logit link models (I am using RMark v. 2.2.7) and the results look reasonable but I wanted to check to see TransitionMatrix has indeed been updated to accommodate logit link models (knowing that just because something runs doesn't necessarily mean the results are correct!).

Thanks,
Javan
jbauder
 
Posts: 56
Joined: Wed May 25, 2011 12:01 pm

Re: Using TransitionMatrix with logit link for Psi in RMark

Postby jlaake » Thu Jan 16, 2020 3:47 pm

I have not made any changes. Have you checked values? Mlogit and logit are the same with 2 categories. I didn't follow your logic as to why you needed to use logit for some models.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: Using TransitionMatrix with logit link for Psi in RMark

Postby jbauder » Tue Jan 21, 2020 5:17 pm

Hi Jeff,

Sorry for my late reply, I'm just now getting back online. I was using the logit link for models where I wanted one of the other transition probabilities to vary independently of the other. For example, we had reason to suspect that harvest mortality might vary temporally but non-harvest mortality might remain constant (or remain constant across sexes or age classes). It was this type of model, fit with the logit link for transition probability, that I tested with the TransitionMatrix() function. The resulting survival estimates were indeed 1 - the sum of all harvest mortalities so the function was correctly calculating the point estimate for survival. I just wasn't sure if there was any reason I should be skeptical of the standard errors/CIs.

Thanks again,
Javan
jbauder
 
Posts: 56
Joined: Wed May 25, 2011 12:01 pm

Re: Using TransitionMatrix with logit link for Psi in RMark

Postby jlaake » Tue Jan 21, 2020 7:07 pm

Hmm! Was along time ago when I wrote that in Appendix C. Not sure why I wrote that or if I changed the code. That restriction is not in the help file for TransitionMatrix. I looked back through NEWS for the package and couldn't find any documented changes when you first asked about this. I ran the example below from the help with the default mlogit and then logit links and I got nearly identical results.

The code in TransitionMatrix is fairly simple to look at. If you look at it you will see that it computes independent confidence intervals for the cell probabilities using the inverse logit transformation. It takes the real estimates (which are probabilities) and their std errors, converts them to a logit of the probability and its std error, computes a normal 95% conf interval on the logit and then converts end points back to probabilities using the inverse logit transformation. Thus, this approach is not dependent on the mlogit transformation. However because the confidence intervals are computed independently for each cell that means that if you sum the lcl and ucl values across cells(states) they will not necessarily nor likely sum to 1.

Code: Select all
# This example is excluded from testing to reduce package check time
data(mstrata)

# fit the sequence of multistrata models as shown for ?mstrata
run.mstrata=function()
{
  #
  # Process data
  #
  mstrata.processed=process.data(mstrata,model="Multistrata")
  #
  # Create default design data
  #
  mstrata.ddl=make.design.data(mstrata.processed)
  #
  #  Define range of models for S; note that the betas will differ from the output
  #  in MARK for the ~stratum = S(s) because the design matrix is defined using
  #  treatment contrasts for factors so the intercept is stratum A and the other
  #  two estimates represent the amount that survival for B abd C differ from A.
  #  You can use force the approach used in MARK with the formula ~-1+stratum which
  #  creates 3 separate Betas - one for A,B and C.
  #
  S.stratum=list(formula=~stratum)
  S.stratumxtime=list(formula=~stratum*time)
  #
  #  Define range of models for p
  #
  p.stratum=list(formula=~stratum)
  #
  #  Define range of models for Psi; what is denoted as s for Psi
  #  in the Mark example for Psi is accomplished by -1+stratum:tostratum which
  #  nests tostratum within stratum.  Likewise, to get s*t as noted in MARK you
  #  want ~-1+stratum:tostratum:time with time nested in tostratum nested in
  #  stratum.
  #
  Psi.s=list(formula=~-1+stratum:tostratum)
  #
  # Create model list and run assortment of models
  #
  model.list=create.model.list("Multistrata")
  #
  # Add on a specific model that is paired with fixed p's to remove confounding
  #
  p.stratumxtime=list(formula=~stratum*time)
  p.stratumxtime.fixed=list(formula=~stratum*time,fixed=list(time=4,value=1))
  model.list=rbind(model.list,c(S="S.stratumxtime",p="p.stratumxtime.fixed",
                                Psi="Psi.s"))
  #
  # Run the list of models
  #
  mstrata.results=mark.wrapper(model.list,data=mstrata.processed,ddl=mstrata.ddl)
  #
  # Return model table and list of models
  #
  return(mstrata.results)
}
mstrata.results=run.mstrata()
mstrata.results
# for the best model, get.real to get a list containing all Psi estimates
#  and the v-c matrix
Psilist=get.real(mstrata.results[[1]],"Psi",vcv=TRUE)
Psivalues=Psilist$estimates
# call Transition matrix using values from time==1; the call to the function
# must only contain one record for each possible transition. An error message is
# given if not the case
tt=TransitionMatrix(Psivalues[Psivalues$time==1,])
# call it again but specify the vc matrix to get se and conf interval
tt=TransitionMatrix(Psivalues[Psivalues$time==1,],vcv.real=Psilist$vcv.real)

# fit the sequence of multistrata models as shown for ?mstrata
run.mstrata=function()
{
  #
  # Process data
  #
  mstrata.processed=process.data(mstrata,model="Multistrata")
  #
  # Create default design data
  #
  mstrata.ddl=make.design.data(mstrata.processed)
  #
  #  Define range of models for S; note that the betas will differ from the output
  #  in MARK for the ~stratum = S(s) because the design matrix is defined using
  #  treatment contrasts for factors so the intercept is stratum A and the other
  #  two estimates represent the amount that survival for B abd C differ from A.
  #  You can use force the approach used in MARK with the formula ~-1+stratum which
  #  creates 3 separate Betas - one for A,B and C.
  #
  S.stratum=list(formula=~stratum)
  S.stratumxtime=list(formula=~stratum*time)
  #
  #  Define range of models for p
  #
  p.stratum=list(formula=~stratum)
  #
  #  Define range of models for Psi; what is denoted as s for Psi
  #  in the Mark example for Psi is accomplished by -1+stratum:tostratum which
  #  nests tostratum within stratum.  Likewise, to get s*t as noted in MARK you
  #  want ~-1+stratum:tostratum:time with time nested in tostratum nested in
  #  stratum.
  #
  Psi.s=list(formula=~-1+stratum:tostratum,link="logit")
  #
  # Create model list and run assortment of models
  #
  model.list=create.model.list("Multistrata")
  #
  # Add on a specific model that is paired with fixed p's to remove confounding
  #
  p.stratumxtime=list(formula=~stratum*time)
  p.stratumxtime.fixed=list(formula=~stratum*time,fixed=list(time=4,value=1))
  model.list=rbind(model.list,c(S="S.stratumxtime",p="p.stratumxtime.fixed",
                                Psi="Psi.s"))
  #
  # Run the list of models
  #
  mstrata.results=mark.wrapper(model.list,data=mstrata.processed,ddl=mstrata.ddl)
  #
  # Return model table and list of models
  #
  return(mstrata.results)
}
mstrata.results=run.mstrata()
mstrata.results
# for the best model, get.real to get a list containing all Psi estimates
#  and the v-c matrix
Psilist=get.real(mstrata.results[[1]],"Psi",vcv=TRUE)
Psivalues=Psilist$estimates
# call Transition matrix using values from time==1; the call to the function
# must only contain one record for each possible transition. An error message is
# given if not the case
TransitionMatrix(Psivalues[Psivalues$time==1,])
# call it again but specify the vc matrix to get se and conf interval
TransitionMatrix(Psivalues[Psivalues$time==1,],vcv.real=Psilist$vcv.real)

tt
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA


Return to RMark

Who is online

Users browsing this forum: No registered users and 8 guests

cron