Constant survival for one group, time interaction for others

[jlaake]

This is a query I received via email. I'm posting here and then a reply. I believe something like this has been posted before in possibly a different form.

Dear Jeff Laake,

I found your email address on the R-Mark documentation, I hope you don’t mind me contacting you.

I am a beginner to R and R-Mark. I need to use MARK Pradel temporal symmetry models, and the R-MARK package is very appealing so as to avoid the design matrix headaches!

I have read most of the documentation available, including the appendix from the Mark manual, yet I can’t seem to find the answer to my question and I suspect that despite investing quite a bit of time in trying to figure this out, I may be missing the point! Initially I need to run a CJS analysis to determine the appropriate structure for survival to then apply in the Pradel models. My first stumbling block is that I have 3 groups, Male Female and Unsexed birds. I need to have the unsexed birds constant and males and females time dependent because in previous analyses we found that the unsexed birds were creating problems due to low numbers in this group, however I need to keep them ‘in there’ for the Pradel analysis.

So if I was doing this in MARK my model structure would be:
survival (Males time, females time, unsexed constant) (males and females different but time dependent, unsexed constant)
survival (Males and females time, unsexed constant) (males and females bunched together but time dependent, unsexed constant)
survival (males constant, females constant, unsexed constant) (all 3 groups different but constant)
survival (males and females constant, unsexed constant). (males and females the same but constant, unsexed constant)

Sorry if this sounds really stupid, but I can’t figure out how to apply this when I get to writing the model functions in R-MARK; from the examples and documentation I have seen I don’t know how to ‘drill down’ to give individual groups different temporal structure in this way. I know how to write it in GEMACO but M/E-surge doesn’t have the nice parameterisations of the Pradel model available, hence the need to switch over to R-MARK.

Hope you can help and apologies if I am ‘missing the point’ any nudges in the right direction would be very gratefully received as RMARK looks pretty cool if only I could get past this block!

Thanks!
Lianne

[jlaake]

Lianne-

I don't mind being contacted individually but would prefer that you use phidot listserver because the question and answer are archived for others to see, I get an email for all phidot listings so it is equivalent to emailing me directly. I will post this query and solution to phidot. If you have any further queries about it I ask that you use phidot so it is archived.

As with any linear model, there are many ways to specify it. The design data are the key to understanding how to make RMark work. Any time you want to do something in a model you want to think about how to create or adapt the design data to accomplish it. In this case, what you are describing is a separate group that is treated differently than the other groups with an interaction. The simplest way is to create a 0/1 variable that limits the interaction to those that have a value of 1. So once you have created your design data list with make.design.data you add another field I'll call KnownSex which is 1 for Male and Female group and 0 for unsexed group. I'll assume it is called ddl and one group is called "U" and the grouping field is called sex

Code: Select all

ddl$Phi$KnownSex=ifelse(ddl$Phi$sex=="U",0,1)

Then the formula can be

Code: Select all

Phi=list(formula=~KnownSex:time:sex)

to get constant value (intercept) for unknown sex, and time variation for M and F.

To get additive difference for sex with time variation you need to create an additional variable say male that is also 0/1 as follows:

Code: Select all

ddl$Phi$Male=ifelse(ddl$Phi$sex=="M",1,0)
Phi=list(formula=~KnownSex:time+Male)


Below is an example with the dipper data.

Code: Select all

data(dipper)
# modify sex field to include unknowns
dipper$sex=substr(as.character(dipper$sex),1,1)
dipper$sex[seq(10,200,10)]="U"
dipper$sex=factor(dipper$sex)
# process data/make design data
dp=process.data(dipper,model="CJS",groups="sex")
ddl=make.design.data(dp)
ddl$Phi$KnownSex=ifelse(ddl$Phi$sex=="U",0,1)
#look at design data
summary(ddl$Phi)
# get summary of model
summary(mark(dp,ddl,model.parameters=list(Phi=list(formula=~KnownSex:sex:time))))
# add male variable
ddl$Phi$Male=ifelse(ddl$Phi$sex=="M",1,0)
# get summary of model
summary(mark(dp,ddl,model.parameters=list(Phi=list(formula=~KnownSex:time+Male))))


It helps if you look at the design matrix it is creating:

Code: Select all

dm=model.matrix(~KnownSex:sex:time,ddl$Phi)

If you look at dm it has columns for unknown sex-time interaction but they contain all 0's and those columns are not used. To see the columns that are used

Code: Select all

dm[,colSums(dm)>0]


Second model:

Code: Select all

dm=model.matrix(~Male+KnownSex:time,ddl$Phi)
dm[,colSums(dm)>0]

To see how this matches up with design data use:

Code: Select all

cbind(ddl$Phi,dm)


There are other ways to specify this model in which the unknown sex is the mean of the female/male using sum contrasts. I believe this is described in Cooch and White. Something like that can be done in RMark by specifying columns in the design data that are essentially columns of the design matrix.

regards --jeff

[Lianne Concannon]

Thank you for the help with this, Jeff. I will give this a try.
Many thanks
Lianne