trap dependence in RD model

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

trap dependence in RD model

Postby monicaarso » Wed May 31, 2023 10:43 am

Hi,
I've successfully fitted trap dependence (td) in CJS models in the past by means of a time-varying covariate, as explained in section C16 of MARK book. I am now wanting to fit td in a robust design modelling framework. I have modified the code (see below) to add a sequence of covariates that contain the capture history entry for each secondary occasion, so that when p is estimated it can be done as a function of this covariate td, and account for trap-dependence. The data consists of 21 primary sampling occasions with 17 to 28 secondary occasions within each.

The function is written as described in Appendix C:
Code: Select all
create.td=function(ch,varname="td",begin.time=1)
  #
  # Arguments:
  # ch - capture history vector (0/1 values only)
  # varname - prefix for variable name
  # begin.time - time for first occasion
  #
  # Value:
  # returns a datframe with trap-dependent variables
  # named varnamet+1,...,varnamet+nocc-1
  # where t is begin.time and nocc is the
  # number of occasions
#
{
  # turn vector of capture history strings into a vector of characters
  char.vec=unlist(strsplit(ch,""))
  # test to make sure they only contain 0 or 1
  if(!all(char.vec %in% c(0,1)))
    stop("Function only valid for CJS model without missing values")
  else
  {
    # get number of occasions (nocc) and change it into a matrix of numbers
    nocc=nchar(ch[1])
    tdmat=matrix(as.numeric(char.vec),ncol=nocc,byrow=TRUE)
    # remove the last column which is not used
    tdmat=tdmat[,1:(nocc-1)]
    # turn it into a dataframe and assign the field (column) names
    tdmat=as.data.frame(tdmat)
    names(tdmat)=paste(varname,(begin.time+1):(begin.time+nocc-1),sep="")
    return(tdmat)
  }
}


And then the data are processed, design data added and models specified as follows (for the purpose of this exercise I am only running no movement, no heterogeneity models):

Code: Select all
#Process data specifying primary and secondary capture occasions (trip sec occ)
time.intervals.SAC_trip=c(rep(0,27),1,rep(0,23),1,rep(0,22),1,rep(0,20),1,rep(0,18),1,rep(0,27),1,rep(0,25),1,rep(0,23),1,rep(0,28),1,rep(0,23),1,rep(0,19),1,rep(0,17),1,rep(0,21),1,rep(0,20),1,rep(0,19),1,rep(0,20),1,rep(0,22),1,rep(0,20),1,rep(0,18),1,rep(0,19),1,rep(0,17))

do.SAC.td=function()
{
  # get data and add the td time-varying covariate, process the data
  df=cbind(trip_markdata,create.td(trip_markdata$ch,begin.time=1))
  dp2=process.data(df,begin.time=2001,model="RDFullHet",time.intervals=time.intervals.SAC_trip)
  # and create the design data to constrain gammas
  ad.ddl=make.design.data(dp2)
ad.ddl=add.design.data(dp2,ad.ddl,parameter="GammaPrime",type="time",bins=c(2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2021),right=FALSE, name="constrain",replace=TRUE)  ad.ddl=add.design.data(dp2,ad.ddl,parameter="GammaDoublePrime",type="time",bins=c(2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2021),right=FALSE, name="constrain",replace=TRUE)
 
  # define parameters
  S.dot=list(formula=~1)
  p.td=list(formula=~td) #trap dependance
  pi.null=list(formula=~1,fixed=1)
  GammaDoublePrime.0=list(formula=~1,fixed=0)
  GammaPrime.0=list(formula=~1,fixed=1)
 
  # create model list
  cml=create.model.list("RDFullHet")
  # run and return models
  return(mark.wrapper(cml,data=dp2,ddl=ad.ddl))
}

td.results.nomov=do.SAC.td()
td.results.nomov


When I run this, I can see that the covariate has been added to the data frame of capture histories, named as td2, td3, td4 etc, but I get the following error:

td.results.nomov=do.SAC.td()
Using default formula for c
Using default formula for f0
S.dot.GammaDoublePrime.0.GammaPrime.0.pi.null.p.td
Error in make.mark.model(data.proc, title = title, parameters = model.parameters, :
Error: Variable td used in formula is not defined in data
Error in mark(model.parameters = model.parameters, initial = initial, :
Misspecification of model or internal error in code
No mark models found


I could not find an example using RD. I wonder if the issue is in the naming of the covariate, given that in the CJS example is says they need to match the times for the capture probabilities? Or is it that the covariate should be for each session instead of each occasion within session? (I don't completely see the logic there though). Any guidance would be appreciated.
Many thanks,
Monica
monicaarso
 
Posts: 33
Joined: Wed Feb 22, 2012 2:58 pm

Re: trap dependence in RD model

Postby jlaake » Wed May 31, 2023 9:52 pm

It should be named tdps where p is primary number (session) and s is secondary number. Then you would use td in formula.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: trap dependence in RD model

Postby monicaarso » Mon Jun 05, 2023 10:12 am

Hi Jeff,
Thanks for the quick reply. Apologies but I am still a bit lost here on how to actually code this. Do you mean the covariate should be named e.g. td20011, td20022, etc, which would correspond to session 2001 and secondary times 1 and 2 ...? But then still only call td in the formula for parameter definition?
Code: Select all
p.td=list(formula=~td)

I am not sure how to go about modifying the function code to produce the covariates named like that, i.e. combining primary sessions and secondary occasions, unless I just modify the ddl instead and generate a new column with that new name for each combination of session (s, primary) and time (p, secondary)?
Thanks again,
Monica
monicaarso
 
Posts: 33
Joined: Wed Feb 22, 2012 2:58 pm

Re: trap dependence in RD model

Postby jlaake » Mon Jun 05, 2023 10:27 am

This would be an individual time varying covariate which is in the data file and not design data because it varies by individual animal. I assume that was what you were doing with dependence at session level.

Think carefully about what you are doing and why. Since this is RD it is closed within session and estimates a p for each occasion and not sure how you would define a to covariate for that first time in each session. Use 0?

This idea was originally considered for cjs type data where it is conditioned on release (in first "capture") and only subsequent capture probability are modeled.

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

Re: trap dependence in RD model

Postby monicaarso » Mon Jun 05, 2023 10:37 am

Hi Jeff,
Yes, I agree I need to think this carefully as I am not sure how one can incorporate td in the RD framework given that p is estimated at each occasion. I might have to account for this source of heterogeneity another way but I'll come back when I've figure it out.
Many thanks again, as always.
Monica
monicaarso
 
Posts: 33
Joined: Wed Feb 22, 2012 2:58 pm

Re: trap dependence in RD model

Postby monicaarso » Tue Jun 13, 2023 11:46 am

Hi,
I just thought I'd provide an update on progress now that I've hit a new wall. :wink:
After reading different papers/forum posts, I decided to use Fletcher (1994) / Williams et al (2002) approach which is to incorporate an individual covariate for p (recapture probability) to model capture hetergeneity (trap dependency), by counting the number of secondary-capture occasions during the previous primary capture session that an individual had been detected. This is used, for example, in Clark et al 2010 paper on black bears.

The data now have an extra number of columns (same as total number of primary occ - 1) that are named td1990, td1991, td1992, etc, with td1990 having the number of secondary occasions an animal was captured the session before (i.e. in 1989 in this instance).

Then when I model p, I can do so by calling:

Code: Select all
p.td=list(formula=~td) #trap dependance


But that doesn't work and I get an error that
Code: Select all
#model 1: S(.) p(td)
> nomov1.t=mark(SAC_trip.process,SAC_trip.ddl,model="RDFullHet",threads=1,
+               model.parameters=list(S=S.dot,p=p.td,GammaDoublePrime=GammaDoublePrime.0,GammaPrime=GammaPrime.0,pi=pi.null))
Error in make.mark.model(data.proc, title = title, parameters = model.parameters,  :
 
Error: Variable td used in formula is not defined in data

Error in mark(SAC_trip.process, SAC_trip.ddl, model = "RDFullHet", threads = 1,  :
  Misspecification of model or internal error in code


So I suspect the issue is either on the very first year (1989) where there is no info from the previous year OR the way I've set it up with individual covariate per session is not what is required within the RD framework in terms of naming. As very well explained by Jeff in the RMark workshop notes:

The only trick to using time-varying individual covariates is to name them properly. To link the covariates to the occasion (time), they should be named with a common prefix (e.g.,
cov) and the sufix should be the label given to each occasion (time). The label depends on the value you
assign to begin.time in process.data, so you have to be consistent between labeling and the names of the
covariates in your data. (...) When I create a formula, I'll use ~td and RMark will recognize that td is not an individual covariate and is not in the design data. It then looks for individual covariates with the names td1991, td1992, td1993 if it is used for p and for variables named td1990, td1991, td1992 if used in formula for survival because the values of time in the design data differ.


Any thoughts?

Thanks...

References:
Clark, Joseph D., Rick Eastridge, and Michael J. Hooker. "Effects of exploitation on black bear populations at White River National Wildlife Refuge." The Journal of Wildlife Management 74.7 (2010): 1448-1456.
Fletcher, David J. "A mark-recapture model in which sighting probability depends on the number of sightings on the previous occasion." Statistics in ecology and environmental monitoring (1994): 105-110.
Williams, Byron K., James D. Nichols, and Michael J. Conroy. Analysis and management of animal populations. Academic press, 2002. Page 552
monicaarso
 
Posts: 33
Joined: Wed Feb 22, 2012 2:58 pm

Re: trap dependence in RD model

Postby jlaake » Tue Jun 13, 2023 12:32 pm

Yes, not having a value for 1989 would be a problem. I've only used this approach for CJS type models that don't model first capture/release. You may want to ask those folks how they implemented it. One possibility is to use 0 for 1989. Have you tried that? This is why I'm a big proponent of providing code and data with papers in an online format.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: trap dependence in RD model

Postby monicaarso » Wed Jun 28, 2023 7:43 am

Hi,
I've now sorted this out with help from Jeff and lovely people that shared their data/models/code with me (thanks to Dr Joseph Clark and Dr Mary Conner). I'll post what I've done in case it's useful to someone else.
1) The first problem is that I was using the wrong model. In order to model p based on the number of time captured (individual covariate) one needs to use one of the Huggins models. In my case I chose RDHFHet.
2) As Jeff pointed out, there needs to be a covariate for the first primary session, which I set to zero for all individuals (otherwise the model will not run). The covariates were then named with the prefix "td" (you can name it differently, but it has to match the name of the covariate called in the model), followed by a number. This can be sequential from 1 or it can be, like in my case, startiong at a year (1989). Just make sure your begin.time equals that first number.

And here is a bit of example code on importing the data and setting up a model:
Code: Select all
#Import inp file and define covariates. Data runs from 1989 to 2021 (33 years)

chSAC=convert.inp("data/SAC_hug_2week.inp",use.comments = TRUE,covariates=c("td1989","td1990","td1991","td1992","td1993","td1994","td1995","td1996","td1997","td1998","td1999","td2000","td2001","td2002","td2003","td2004","td2005","td2006","td2007","td2008","td2009","td2010","td2011","td2012","td2013","td2014","td2015","td2016","td2017","td2018","td2019","td2020","td2021"))

#Process data specifying primary and secondary capture occasions (2week blocks here)
time.intervals.2week=c(rep(0,3),1,rep(0,6),1,rep(0,7),1,rep(0,8),1,rep(0,7),1,rep(0,8),1,rep(0,7),1,rep(0,6),1,rep(0,6),1,rep(0,6),1,rep(0,8),1,rep(0,5),1,rep(0,9),1,rep(0,8),1,rep(0,10),1,rep(0,9),1,rep(0,10),1,rep(0,10),1,rep(0,10),1,rep(0,10),1,rep(0,10),1,rep(0,9),1,rep(0,9),1,rep(0,9),1,rep(0,10),1,rep(0,10),1,rep(0,11),1,rep(0,11),1,rep(0,11),1,rep(0,10),1,rep(0,11),1,rep(0,8),1,rep(0,9))

SAC.process=process.data(chSAC,begin.time=1989,model="RDHFHet",time.intervals=time.intervals.2week)

#Create the design data (constrain Gammas)
SAC.ddl=make.design.data(SAC.process)
SAC.ddl=add.design.data(SAC.process,SAC.ddl,parameter="GammaPrime",type="time",bins=c(1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2021),right=FALSE, name="constrain",replace=TRUE)
SAC.ddl=add.design.data(SAC.process,SAC.ddl,parameter="GammaDoublePrime",type="time",bins=c(1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2021),right=FALSE, name="constrain",replace=TRUE)

# example model
S.dot=list(formula=~1)
p.td=list(formula=~td,share=TRUE) #trap dependance
pi.null=list(formula=~1,fixed=1)
GammaDoublePrime.0=list(formula=~1,fixed=0)
GammaPrime.0=list(formula=~1,fixed=1)

#model 1: S(.) p(td)
nomov1=mark(SAC.process,SAC.ddl,model="RDHFHet",threads=1,            model.parameters=list(S=S.dot,p=p.td,GammaDoublePrime=GammaDoublePrime.0,GammaPrime=GammaPrime.0,pi=pi.null),silent=T)



Thanks again and happy mark recapture modelling
Monica
monicaarso
 
Posts: 33
Joined: Wed Feb 22, 2012 2:58 pm


Return to RMark

Who is online

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