www.phidot.org

I'm new to RMark (but know both R an MARK well) and I'm working with a student analyzing a data set where most birds are marked as chicks. We want a model with a different Phi for the first year of life, but the same Phi thereafter (i.e., age has two levels, 'juv' and 'adult'). There should be plenty of information in the data to fit a Phi(time*age) model except that Phi for 'adult' in the first year is not estimable (there is only one 'adult' individual released in the first year). We can fix the real parameter for this Phi, but we still get one redundant beta parameter in the model matrix.

If I remember correctly MARK does not object to having a row of all zeros in the model matrix, so the way I would do this in MARK is to construct a model matrix like this

1:1997:juv: 1 0 0 0 0
2:1998:juv: 1 1 0 0 0
3:1999:juv: 1 0 1 0 0
4:1997:adult: 0 0 0 0 0
5:1998:adult: 1 1 0 1 0
6:1999:adult: 1 0 1 1 1

...and then fix real parameter 4 (I guess one can put whatever values in row 4, but the point is that I have 5 columns whereas the equivalent model matrix in RMark has 6 columns).

What is the best way to do the equivalent in RMark?

If it helps to answer, here is some code to base it on:

GLOB.proc = process.data(GLOB, model="CJS", group="age.class", begin.time=1997, age.var=1, initial.age=c(1,1,2,4,4))

GLOB.ddl = make.design.data(GLOB.proc, parameters=list(Phi=list(age.bins=c(1,2,18)), p=list(age.bins=c(2,3,4,19))), right=FALSE)
levels(GLOB.ddl$Phi$age)=c("1K","2Kpluss")
levels(GLOB.ddl$p$age)=c("2K","3K","4Kpluss")

phis.to.fix = GLOB.ddl$Phi$par.index[GLOB.ddl$Phi$time == 1997 & GLOB.ddl$Phi$age=="2Kpluss"]

time.age.int = list(formula= ~time*age, fixed=list(index=phis.to.fix, value=0), link="loglog")
td.time.age=list(formula=~td + time + age)

fit = mark(GLOB.proc, GLOB.ddl, model.parameters = list(Phi=time.age.int, p=td.time.age))

Best regards,

Torbjørn Ergon

This is discussed in section C.15 of the appendix in GIM for the case when only young are marked. A simple solution for your case is to delete the one adult capture history. You are already assigning its survival to 0 so not much difference. Then use ~-1+time:age which fits a separate beta for each valid time and age combination. I've not tried it but that may work with the fixed parameter as well.

You'll want to review C.15 because it describes how to restrict time effects to say young with a constant adult survival. The : operator is a good way to restrict effects and interactions to specific subsets. For example, you could create different time categories (eg bin first and second year) and then apply that to adults by creating a 0/1 adult field and using adult:alt.time and young:time.

The trick in moving from MARK to RMark is thinking about how to create and use design data to build the model you want. You can then use model.matrix (except with indiv covariates) to see that it is creating the DM that you envisioned.

Setting an all 0 row in MARK for Phi is different than fixing the real Phi to 0. With a logit link I think MARK will compute Phi=0.5.

Let me know if you have any further questions. --jeff

One further point. If you build the DM with model.matrix, you may seem columns that are removed by RMark. It removes any column that is all 0. That is discussed in the appendix. Also, note the use of remove.intercept for the cases where you want to use something like ~-1+time:age+sex. That will not remove the intercept because it simply creates a male and female effect. The solution is to use ~time:age+sex and add remove.intercept=TRUE in the model specification. Don't use -1 and remove.intercept.

--jeff

Many thanks Jeff! I think I understand how it works now.

The approach of using ~-1+time:age works when I delete the rows corresponding to survival for the one adult individual released in the first year in ddl$Phi – and when I do that it doesn’t matter whether I remove this individual from the data or not (ddl$Phi still contains the parameters (rows) for adults in the first year even though there are no such individuals in the data).

Before I read your reply I tried obtaining the DM I wanted in another less elegant way:

> time.age.int = list(formula= ~
I(time=="1998") +
I(time=="1999") +
I(time=="2000") +
I(time=="2001") +
I(time=="2002") +
I(time=="2003") +
I(time=="2004") +
I(time=="2005") +
I(time=="2006") +
I(time=="2007") +
I(time=="2008") +
I(time=="2009") +
I(time=="2010") +
I(time=="2011") +
I(age=="2Kpluss") +
I(age=="2Kpluss" & time=="1999") +
I(age=="2Kpluss" & time=="2000") +
I(age=="2Kpluss" & time=="2001") +
I(age=="2Kpluss" & time=="2002") +
I(age=="2Kpluss" & time=="2003") +
I(age=="2Kpluss" & time=="2004") +
I(age=="2Kpluss" & time=="2005") +
I(age=="2Kpluss" & time=="2006") +
I(age=="2Kpluss" & time=="2007") +
I(age=="2Kpluss" & time=="2008") +
I(age=="2Kpluss" & time=="2009") +
I(age=="2Kpluss" & time=="2010") +
I(age=="2Kpluss" & time=="2011"),
link="loglog")

This construct will produce the DM I want in model.matrix(), and it works in other R functions like lm().

However, this did not work in mark() and I got the following output:

STOP NORMAL EXIT
Error in (npar + 1):nbeta : NA/NaN argument
In addition: Warning message:
In Ops.factor(npar, 1) : + not meaningful for factors

I then tried a simpler example:

> age.bf = list(formula= ~ I(age=="2Kpluss"), link="loglog")
> fit.age.bf = mark(GLOB.proc, GLOB.ddl, model.parameters = list(Phi=age.bf, p=td.time.age))

This works, but the heading in the output looks funny:

Note: only 1 parameters counted of 20 specified parameters
AICc and parameter count have been adjusted upward

Output summary for CJS model
Name : Phi(~I(age == "2Kpluss"))p(~td + time + age)

Npar : 20 (unadjusted== "2Kpluss"))p(~td + time + age) } = 20 )
-2lnL: 1
AICc : NA (unadjusted== "2Kpluss"))p(~td + time + age) } = 21706.817 )

...and I get the following warning messages:

1: In extract.mark.output(out, model, adjust, realvcv, vcvfile) :
Not all parameters were estimated but not able to find non-estimable parameters

2: In Ops.factor(lnl, 2 * nbeta) : + not meaningful for factors

However, the parameter estimates and the vcov are *exactly* the same as I get with

> age = list(formula= ~ age, link="loglog")
> fit.age = mark(GLOB.proc, GLOB.ddl, model.parameters = list(Phi=age, p=td.time.age))


This is not a very elegant way of constructing the DM, but it is an easy and quick way if one are used to thinking about DM’s in the MARK editor. And it would be useful to be able to use the I() construct so that one didn’t have to add variables to the ddl object. Is it a quick fix?

Cheers,

Torbjørn

Note that when you delete the design data for Phi, you are setting Phi=1 (the default value) for that parameter which should affect the LnL value you get if you leave the adult in the data.

I have to fiddle with the result of model.matrix before passing it to MARK to handle things like individual covariates and it is likely causing problems with I(). I'll add it to my list of tasks to look at for RMark but I think it is more straightforward and clear to create variables in the design data rather than long formula using I(). Although, I don't do it or recommend it, you can essentially build the dm by adding variables to the design data with 0 and 1 values and specifying that variable in the formula. I do use columns like that to restrict interactions or to restrict the use of a covariate to specific times or groups, etc.

If you want to explore the problem, type debug(make.mark.model) and step through the process to see where the issue with I() is occurring.

regards --jeff