www.phidot.org

Hi,
I am running a CRDMS analysis in RMark using 740ish capture histories with 18 primary occassions*5 secondary occassions each. I would like to modify the structure of the design data to merge some estimates of "p" and "c".
When trying to modify the design structure of primary session estimates of "p" ("session" in the design data) , I am finding problems with using some functions related to the fact that the functions expect a field of type "time" instead of type "session" (SEE BELOW). However, when modifying the design structure of secondary sessions in "c", I have no problems because there is a field called "time" (which actually stands for secondary session or visit for that parameter). I'm sure that I am missing something, but I have hanged around through postings and manuals but haven't found anything related.

Preliminar code:

Code: Select all

grazing=convert.inp("Grazing_Robust_RMARK",group.df=data.frame(sex=c(rep("Male",4),rep("Female",4)),grid=rep(c("A2","A4","E2","E4"),2)))
summary(grazing)                   
#set time intervals with 18 primary periods each with 5 secondary occasions
t.int<-c(rep(c(0,0,0,0,1),17),c(0,0,0,0))
grazing.processed=process.data(grazing,model="CRDMS",time.interval=t.int,groups=c("sex","grid"))
grazing.ddl=make.design.data(grazing.processed)

Problems start when I try to merge some primary sessions by using bins through the function "add.design.data"

Code: Select all

grazing.ddl=add.design.data(grazing.processed, grazing.ddl, parameter="p", type="session",
bins=c(0,1,2,3,4,5,6,13,14,15,16,17,18),name="newsession",right="FALSE",replace=TRUE)


I get the following message

Code: Select all

#Error in add.design.data(grazing.processed, grazing.ddl, parameter = "p",  :
#invalid type


If I try it in a different way, by using "merge_design.covariates" and a new categorical covariate

Code: Select all

Monthdf=data.frame(session=1:18,Month=c("May06","June06","July06","Aug06","Sept06","Oct06","Winter06","Winter06","Winter06","Winter06","Winter06","Winter06","May07","June07","July07","Aug07","Sept07",
"Oct07"))
grazing.ddl$p= merge_design.covariates(grazing.ddl$p,Monthdf)


I get another error message (below)

Code: Select all

#Error in merge_design.covariates(grazing.ddl$p, Monthdf) :
#Monthdf does not contain field named time


Now, when defining a covariate "night" in parameter "c" that accounts for different capture rates morning vs. night shifts, everything works fine

Code: Select all

grazing.ddl$c$night=0  #dummy field "night" in c with all 0s
grazing.ddl$c$night[grazing.ddl$c$time==2|grazing.ddl$c$time==4]=1   #dummy night (=1) affects captures 2 and 4 of every session


This problem also arises when using formulas with "p"

Code: Select all

p.session.stratum.sex.grid.coffset=list(formula=session*stratum*sex*grid+c:night, share=TRUE)
#Error: object 'session' not found


Can anyone help me to see what am I missing?

Cheers,
N

p and c are session-specific parameters. session is the primary time value. time is secondary time value and it is the same for p and c so I'm not sure why you thought they differed. What does differ between p and c is that p starts at time 1 and c starts at time 2 because it is recapture probability. merge_occasion.data was designed for non-robust designs. In retrospect, I should have kept session as time and made the current time secondary.time. But I didn't. You can make it work with some modification of the design data. The design data created by make.design.data are the default data. You can change it many different ways but there are some restrictions. Namely, never do anything that sorts the design data in any way other than the way it was created. Also, parameters that are shared (eg p and c) have to have the exact same column definitions in the design data because they are appended (rbind) to each other when the parameters are shared. You are not restricted to using add.design.data and merge_occasion.data. Those are simply convenience functions. You can use any R you want to modify the design data as long as you meet the restrictions above.

For example, you can rename time to say secondary.time and session to time and then you should find that merge_occasion.data. If you can create examples with one of the example data sets in RMark it helps me help you. For example, here is an example with the crdms data.

Code: Select all

data(crdms)
#set time intervals
#4 primary periods each with 3 secondary occasions
t.int<-c(rep(c(0,0,1),3),c(0,0))
#process data for RMark
crdms.data<-process.data(crdms,model="CRDMS",time.interval=t.int,
strata.labels=c("1","U"))
#change Psi parameters that are obtained by subtraction
crdms.ddl<-make.design.data(crdms.data,
parameters=list(Psi=list(subtract.stratum=c("1","1"))))
# change names of design data in p and c
pnames=names(crdms.ddl$p)
pnames[pnames=="time"]="secondary.time"
pnames[pnames=="session"]="time"
names(crdms.ddl$p)=pnames
names(crdms.ddl$c)=pnames
# change definition of Time trend variable to be for primary periods
crdms.ddl$p$Time=as.numeric(as.character(crdms.ddl$p$time))-1
crdms.ddl$c$Time=as.numeric(as.character(crdms.ddl$c$time))-1


Now you have to recognize that with the definitions being different the models in the example for crdms will either not work or they will be different because I've changed the names of the variables. Again, if you don't like the design data I created change it.

--jeff

Thanks Jeff, very instructive.

I have fiddled around with the crdms data you suggested. As you suggested, renaming works fine, but I have found other problems

Code: Select all

data(crdms)
    #set time intervals
    #4 primary periods each with 3 secondary occasions
    t.int<-c(rep(c(0,0,1),3),c(0,0))
    #process data for RMark
    crdms.data<-process.data(crdms,model="CRDMS",time.interval=t.int,
    strata.labels=c("1","U"))
    #change Psi parameters that are obtained by subtraction
    crdms.ddl<-make.design.data(crdms.data,
    parameters=list(Psi=list(subtract.stratum=c("1","1"))))
    #the resulting design data
    summary(crdms.ddl)
    summary(crdms.ddl$S)
    summary(crdms.ddl$p)
    summary(crdms.ddl$c)
    summary(crdms.ddl$Psi)
    summary(crdms.ddl$N)
    #say we want to constrain 2 bins into primary sessions e.g. [0,3) and [3,4]
    crdms.ddl=add.design.data(crdms.data, crdms.ddl, parameter="p", type="session",
    bins=c(0,1,3,4),name="newsession",right="FALSE",replace=TRUE)
    #Error in add.design.data(crdms.data, crdms.ddl, parameter = "p", type = "session",  :
    #invalid type
    #We need to change names of design data in p and c
    pnames=names(crdms.ddl$p)
    pnames[pnames=="time"]="secondary.time"
    pnames[pnames=="session"]="time"
    names(crdms.ddl$p)=pnames
    names(crdms.ddl$c)=pnames
    #Now we can add the bins correctly to both p and c
    crdms.ddl=add.design.data(crdms.data, crdms.ddl, parameter="p", type="time",
    bins=c(0,3,4),name="newsession",right="FALSE",replace=TRUE)
    crdms.ddl=add.design.data(crdms.data, crdms.ddl, parameter="c", type="time",
    bins=c(0,3,4),name="newsession",right="FALSE",replace=TRUE)
    summary(crdms.ddl$p)
    summary(crdms.ddl$c)


We could also modify secondary occassions ("secondary.time") to make occasion 2 different from others e.g. "night"

Code: Select all

crdms.ddl$p$night=0  #dummy field "night" in p with all 0s
    crdms.ddl$p$night[crdms.ddl$p$secondary.time==2]=1
    crdms.ddl$c$night=0  #dummy field "night" in p with all 0s
    crdms.ddl$c$night[crdms.ddl$c$secondary.time==2]=1
    summary(crdms.ddl$p)
    summary(crdms.ddl$c)


However say we also wanted to add bins to "secondary.time" to both p and c. The following shows error, because again type="secondary.time" does not work with "add.design.data".

Code: Select all

crdms.ddl=add.design.data(crdms.data, crdms.ddl, parameter="p", type="secondary.time",
    bins=c(0,2,3),name="newsecond",right="FALSE",replace=TRUE)
    #Error in add.design.data(crdms.data, crdms.ddl, parameter = "p", type = "secondary.time",  :
    #invalid type


We could try to change "secondary.time" back to "time", and "time" to "session" (or whatever name you wanted) for it to work. However, this messes around the fields "newsession" and "night" we've just defined because we loose their consistency with "time" and "secondary.time"

Code: Select all

pnames[pnames=="time"]="session"
    pnames[pnames=="secondary.time"]="time"
    names(crdms.ddl$p)=pnames
    names(crdms.ddl$c)=pnames
    summary(crdms.ddl$p)
    summary(crdms.ddl$c)


Any ideas? Cheers,

N

The problem you are having is that you are trying to use add.design.data which is not a very smart function. It only works with time,age and cohort. I probably need to make it smarter or better yet just get rid of it because it is making you think in a very limited way. You can use any R statements to add fields to the design data. Look at the code in add.design.data. Type RMark:::add.design.data to see the code. Actually typing add.design.data will work because it is exported but if you include the packagename::: in front of a function you can look at any function (unless it is internal). Anyhow if you do that you'll see it is simply using the cut function which takes any numeric field and breaks it up into bins. Here is the statement from the function for time

Code: Select all

new.data = cut(design.data$time, bins, include.lowest = TRUE,
            right = right)


It is actually shorter to use R than add.design.data

Code: Select all

numeric.session=as.numeric(as.character(crdms.ddl$p$session))
crdms.ddl$p$newsession=cut(numeric.session,breaks=c(1,3,4),include.lowest=TRUE,right=FALSE)
crdms.ddl$p
> crdms.ddl$p
   par.index model.index group time stratum session Time 1 U c newsession
1          1          25     1    1       1       1    0 1 0 0      [1,3)
2          2          26     1    2       1       1    1 1 0 0      [1,3)
3          3          27     1    3       1       1    2 1 0 0      [1,3)
4          4          28     1    1       1       2    0 1 0 0      [1,3)
5          5          29     1    2       1       2    1 1 0 0      [1,3)
6          6          30     1    3       1       2    2 1 0 0      [1,3)
7          7          31     1    1       1       3    0 1 0 0      [3,4]
8          8          32     1    2       1       3    1 1 0 0      [3,4]
...


You can also make it more readable with things like

Code: Select all

> levels(crdms.ddl$p$newsession)=c("sessions1&2","sessions3&4")
> crdms.ddl$p
   par.index model.index group time stratum session Time 1 U c  newsession
1          1          25     1    1       1       1    0 1 0 0 sessions1&2
2          2          26     1    2       1       1    1 1 0 0 sessions1&2
3          3          27     1    3       1       1    2 1 0 0 sessions1&2
4          4          28     1    1       1       2    0 1 0 0 sessions1&2
5          5          29     1    2       1       2    1 1 0 0 sessions1&2
6          6          30     1    3       1       2    2 1 0 0 sessions1&2
7          7          31     1    1       1       3    0 1 0 0 sessions3&4
8          8          32     1    2       1       3    1 1 0 0 sessions3&4


where the labels are anything you want. You'll benefit by spending some time learning more R.

Here are some links you may find useful.
https://github.com/NMML/R-User-Meetings
http://pairach.com/2012/02/26/r-tutorials-from-universities-around-the-world/

Plus all the material you'll find on CRAN. Let me know if you still have further questions. Glad to help.

regards--jeff

Thanks Jeff,
I managed to solve the bin issue by means of using some "smarter" R commands which gave me more flexibility with the design data, so thanks for the tip. I have moved on the analysis and have a quite simple doubt but I haven't actually found a simple solution through the postings.

Recalling, I am running a CRDMS analysis in RMark using 740ish capture histories with 18 primary occassions*5 secondary occassions each (actually I have moved on to HFHetRDMS, but the issue is common to serveral model "families"). I have an issue because I have not data for some of the primary occassions. On Year1, I captured in May, June, July, Aug, Oct and same on Year2. Estimates of monthly transitions and survival rates are calculated fine for May-June, June-July and July-Aug both years, but not for Aug-Sept or Sept-Oct or Oct-Nov Year1, etc...

I have tried several solutions. First, I have fixed p=0 for missing primary occasions, which actually produces more acurate estimates of May-June, June-July and July-Aug estimates, but no effect on the problematic intervals. I have also tried defining new bins trying to reduce the number of unique estimates with missing primary occasions, e.g. simultaneously or separately defining bin1=(Aug-Sept year1 + Sept-Oct year1), bin2=(Aug-Sept year2 + Sept-Oct year2), bin3=(all winter Year 1 + spring Year2 monthly intervals), so that at each end of each bin there is a non-missed capture occasion. I understand that for bin3 it's difficult to obtain estimates because there is plenty of missing data, but I don't see why it is so hard to obtain estimates of the other two bins.

Ideally I would like at least to obtain reasonable estimates for all monthly intervals between May and October, and not throw to the bin the October primary capture data...

Again I must be missing something. Help much apreciated. Cheers,

N

Hopefully someone with more experience with CRDMS than me will be answer your question. I've not used it that much. I would make sure that you are setting all the p's to zero that you need and not setting incorrect ones to 0. Other than that I don't have much to offer. Sorry --jeff