Closed Captures-behaviour model

[Owen]

I have a 5-event closed capture expt. (Schnabel-type expt.). I'm new to RMark and am trying to duplicate MARK results with RMark. The model I'd like to fit is a p(time), c(time) model, where the c(time) is additively different to the p(time). I've managed to get the model p(time),c(time), with no difference, using a staement like -> p.time=list(formula=~time,share=TRUE). Now I need to include an additive difference for the p and c. Not sure how that's done.

I've looked in the Cooch Appendix C and can't see much in the way of using RMark and closed captures. Thanks for any input.

[jlaufenb]

Owen

I have used the following statement to model p varying by time with an additive behavioral effect:

ptime.c = list(formula=~time+c,share=TRUE)

Jared

[Owen]

Thanks Jared. That did work. I wouldn't have come up with it.

If I'm now not interested in a behaviour component, but I want to bin times together, can you give adice on that? For example, there was high water in the last 2 events (of 5 events) such that I want to set p1=p2=p3 (and c2=p2=c3=p3) and p4=p5 (and c4=p4=c5=p5). I'm not seeing how I'd invoke the formula/SHARE component to make that work. Thanks very much for any advice. I really appreciate it.

[jlaufenb]

Owen

I'll admit that I am still somewhat of a novice with RMark and have never binned capture periods like you are suggesting, but I'll take a crack at it. Your situation is very similar to the dipper example for design covariates described on page C-27 of Appendix C in the MARK book. Based on the example code, I would suspect something like the following might work. Try it and see.

yourdata.ddl$p$Flood=0
yourdata.ddl$p$Flood[yourdata.ddl$p$time==4 | yourdata.ddl$p$time==5]=1

pFlood.c=list(formula=~Flood+c,share=TRUE)

[jlaake]

Glad to see this list working where others are answering RMark questions. Thanks Jared. Jared's advice is correct except that he specified a model with a behavioral component which you didn't specify. If you define the Flood variable as he suggested and use

pFlood.c=list(formula=~Flood,share=TRUE)

you should get what you want. It will help you to understand RMark more if you look at the design data list (ddl). With the exception of indivdual covariates,this is the data which are used to build the design matrix from the formula. That is the primary job that RMark accomplishes and it isn't a mystery if you look at the design data. To create a particular design, it must work off of fields in the design data. I can understand your frustration with the closed model because I don't have much about it the documentation but had you looked at the design data you would have seen the dichotomous c field which is 0 in the p design data and 1 in the c design data. Also, there is ?edwards.eberhardt example that shows the +c approach and there are posts about this issue in phidot -- many from Jared who has worked is way through it. But there is no documenation section for closed models so I can understand your frustration. However, the question you asked here his described elsewhere under different models in Appendix C and the workshop notes because it a simple model constraint on the parameters that is used in any type of model. It can be done with a 0/1 variable as used in the dipper example if there are only 2 periods or you can bin the times into 2 or more intervals and use the binned time in the formula. See ?make.design.data and ?add.design.data or you can add design data like Jared showed with standard R commands. The ddl is simply a list of dataframes and the fields can be redefined however you want. There is only one limitation -- NEVER sort the dataframes in the ddl because it assumes they are in their original order when it relates the design data to the PIMS.

It will also help if you do a search on model.matrix in Appendix C and the workshop notes and study the examples it gives of how it produces design matrices from the design data and the formula. Remember that if you were using MARK you would be creating those by hand. Thus, to understand RMark, you need to understand how it creates the design matrices for you so you can understand how to construct the models that you want.

One of these days I'll get to creating more complete documentation for each model, but currently the little time I can devote to RMark is used to add models in MARK that are not in RMark at present and that list seemingly never stops growing.

regards --jeff

[Owen]

Excellent. Thanks Jared and Jeff; that code gave results identical to those I got in MARK. I am compeltely new to RMark but am slowly seeing the bigger picture of how the ddl can be customized and then desired design matrix constructed from the formula component. Will check out the Edwards.Ebehardt example and Jared's posts. Thanks.

[Owen]

I should add that I found I needed to do the same manipulations to the yourdata.ddl$c component. That is,
yourdata.ddl$c$Flood=0
yourdata.ddl$c$Flood[yourdata.ddl$c$time==4|yourdta.ddl$c$time==5]=1

Otherwise I get this:
Error: for a shared p&c model, their design data columns must match
. If you add design data to p it must also be added to c.
Columns of p: group time Time c High
Columns of c: group time Time c
Error in make.mark.model(data.proc, title = title, covariates = covariates, :
Function terminated

[jlaake]

That is correct. I forget which is why I created the error message. Glad to see you catching onto this. --jeff

[agoldberg]

I have a similar problem but I am working in a CJS model. I want to bin time into year and I have 2 trapping intervals within each year. I tried the following code but this is different the DM's I would make in MARK.

Code: Select all

multigroup.process=process.data(multigroup,model="CJS",time.intervals=c(.8,10.03,1.008,10.303,1.097),begin.time=0,
                                groups=c("sex","trt"))

multigroup.ddl=make.design.data(multigroup.process)

multigroup.ddl$Phi$year=0
multigroup.ddl$Phi$year[multigroup.ddl$Phi$time==10.83]=1
multigroup.ddl$Phi$year[multigroup.ddl$Phi$time==11.838]=1
multigroup.ddl$Phi$year[multigroup.ddl$Phi$time==22.141]=2

multigroup.ddl$p$year=0
multigroup.ddl$p$year[multigroup.ddl$p$time==10.83]=1
multigroup.ddl$p$year[multigroup.ddl$p$time==11.838]=1
multigroup.ddl$p$year[multigroup.ddl$p$time==22.141]=2
multigroup.ddl$p$year[multigroup.ddl$p$time==23.238]=2


I have read ?add.design.data and it seems like I could do something with:

Code: Select all

ddl=add.design.data(multigroup.process,ddl,parameter="Phi",type="time",
                    bins=c(0,.8,11.838,23.238),name="yearbin")

but I am not quite sure how to make that work or if I am on the correct path.

Thanks,
Amanda

[jlaake]

Amanda-

You are close. Currently you have defined year as a numeric variable which is fine because you only have 2 years so whether you treat it as numeric or factor will not affect the model fit but will make the design matrix (DM) differently. When year is numeric for the formula ~year you'll get a DM that looks like

Code: Select all

Intercept    year
1                  1
1                  2


If you code it as factor then for the formula ~year the DM will be

Code: Select all

Intercept    year2
1                 0
1                 1


If I understand correctly that you only have 2 years then use

Code: Select all

multigroup.ddl$Phi$year=1
multigroup.ddl$Phi$year[multigroup.ddl$Phi$time==22.141]=2
multigroup.ddl$Phi$year=factor(multigroup.ddl$Phi$year)

multigroup.ddl$p$year=1
multigroup.ddl$p$year[multigroup.ddl$p$time%in%c(22.141,23.238)]=2
multigroup.ddl$p$year=factor(multigroup.ddl$p$year)


Run each line above and look to see what it is doing in the ddl, so you understand. You could also use ifelse function like

Code: Select all

multigroup.ddl$p$year=ifelse(multigroup.ddl$p$time%in%c(22.141,23.238),2,1)

Then they will be factor variables and will be more like what you would code in MARK. The add.design.data function would create a factor variable directly with the cut function but you want the cut points to only specify 2 groups and you would need to use Time instead of time because Time is numeric and time is a factor variable. See ?cut in R.

Spend some time reading about numeric and factor variables with the R language and how these differ with regard to formula. In R, see ?factor. Note that I didn't test code so may have typos but hopefully not.

--jeff