www.phidot.org

Hi, I am relatively new to RMark but have been running my data using MARK for a bit longer. I'm modeling 4 primary occasions with various secondary lengths. We are trapping prairie dogs essentially for two primary periods a year. We have finished two field seasons (4 primary sessions) and I have one more to go). I want to move it into RMark to make it easier to add new data to models and run heterogeneity models as my design matrices were getting very complicated to produce. I am interested in looking at survival rates but I also need abundance estimates to report to the parks I am trapping on. I've been attempting to write my code by fudging around a friends code from a previous project but I seem to be running into a couple of problems:
1) For now I want to just model no movement (when I feel more comfortable I'll probably add in random movement which I had back in MARK but I'm trying to keep things as simple as I can for now), but with the code I am currently running, I am still getting Gamma' to come up .5 while Gamma'' is 0. What am I doing wrong?

2) Am I modeling p=c correctly? I think I read from others coding on here this is how I should be doing it but I just wanted to double check? We only catch about 80-120 animals each primary occasion so our p's tend to be low and as a result I need to keep p=c if I'm going to model time for sake of having too many parameters.

I am attaching the code that I am working with and any help would be greatly appreciated. Thanks so much,
Amanda

#RMark analysis
#Load RMark package
library(RMark);
rm(list=ls())
cleanup(ask=FALSE)

#inport .inp
bentsadults=convert.inp("C:/R_robust/inpfiles/bentsadultsrobust.inp",
group.df=data.frame(sex=c("Male", "Female")))

bentsadults.robust.run=function()
{

#Process the data so models can be analyzed
#Time intervals for 4 primary occasions (300+ collectors)
#use the format from MARK where you put in the x/60 for time intervals
time=c(0,0,0,0,0,0,.733,0,0,0,0,4.23,0,0,0,0,0,0,0,0,0,0,0,0,1.47,0,0,0,0,0,0,0,0,0,0,0,0,0)
bentsadults.proc=process.data(bentsadults,model="RDFullHet",time.intervals=time,
groups="sex")
summary(bentsadults.proc)

#Create the design data
bentsadults.ddl=make.design.data(bentsadults.proc)

#Create parameter specifications for formulas
#S.dot=list(formula=~1)
#S.group=list(formula=~group)
S.time=list(formula=~time)
#S.grpbytime=list(formula=~group*time)
#S.grpplustime=list(formula=~group+time)
#GammaDoublePrime.0=list(formula=~1,fixed=0)
#GammaPrime.0=list(formula=~1,fixed=0)
GammaDoublePrime.dot=list(formula=~1,fixed=0,share=TRUE)
#p.dot=list(formula=~1)
#p.time=list(formula=~time)
#p.group=list(formula=~group)
#p.time.group=list(formula=~time*group)
#p.timeplusgroup=list(formula=~time+group)
#p.dotsession=list(formula=~session)
p.dotc=list(formula=~session,share=TRUE)
p.dotgroupc=list(formula=~session:group,share=TRUE)
p.dottimec=list(formula=~session:time,share=TRUE)
p.timegroupc=list(formula=~session:time*group,share=TRUE)
p.timeplusgroupc=list(formula=~session:time+group,share=TRUE)
#c.dot=list(formula=~1)
#c.time=list(formula=~time)
#c.group=list(formula=~group)
#c.time.group=list(formula=~time*group)
#c.timeplusgroup=list(formula=~time+group)
#c.dotsession=list(formula=~session)
pi.dot=list(formula=~1)
pi.group=list(formula=~group)
pi.session=list(formula=~session)
pi.none=list(formula=~1, fixed=0)
#N.dot=list(formula=~1)
#N.group=list(formula=~group)
N.session=list(formula=~session)
N.grpbysession=list(formula=~group*session)

# create model list
cml=create.model.list("RDFullHet")
#row.names(cml)=c(1:length(cml$p))

# run and return models
results=mark.wrapper(cml,data=bentsadults.proc,ddl=bentsadults.ddl,adjust=TRUE)
return(results)
}

robust.results=bentsadults.robust.run()
robust.results

results.table=robust.results$model.table

With regard to #1, see example on page C-99 in Appendix C. When you use create.model.list it uses each pairing of models and when you don't list one it uses the default. What you want to do is to combine the model specifications for both Gamma parameters. The example shows how to do that for the parameters p1 and p2 in that model example.

With regard to #2, it looks fine to me. When you use share=T for p, it sets up p=c and you can also add +c to the formula to create an additive difference between c and p.

--jeff

Jeff,
Thanks that was very helpful. Sorry but I have one more question for you. I am running my data now and looking more closely at it and it appears that my pi's are not factoring in? I've really only been running full heterogeneity models in closed capture not robust but when I did that I would get double the number of p's and c's as I would without heterogeneity. But with Rmark I am not getting any difference? Instead I am just getting that my pi's are .5. Something must be up as I know it should be much different than this (I expect that pi's should be around .7 or . and I'm sure it's just something that I'm coding wrong or forgetting to put in there. I tried reading the help and running the same set-up as the help file in R (?robust) but that seems to give me the same answer.

Thanks,
Amanda
--

See ?edwards.eberhardt for an example. You want to use the word mixture in your formula for p to get it vary it for each mixture. ~time + mixture creates a model with time differences and an additive difference for mixtures; whereas, ~time*mixture has all separate time-varying p's for each mixture.

--jeff

Jeff,
Thanks so much, that was really helpful.

Amanda