RMark different number of beta parameters for group variable

[aswea]

Are there reasons that RMark would use two beta parameters to represent a dichotomous group variable instead of one? One of the grouping variables in my dataset is called transmission interval (TI) and it has two levels (60 and 90). When I parameterize recapture probability using TI in an additive model, RMark gives me one parameter (as in 60 vs 90) unless there is an interaction among the other parameters in the model. In this case, it gives me two parameters (one for 60 and another for 90). The SEs on the beta terms go to zero when there are two parameters for TI. The highest ranked hypothesis uses TI as an additive term and includes an interaction so I would like to know if I can trust those betas.

Thank you very much,
Aswea

Col2009<-import.chdata("CH_RMark.txt", header=TRUE, field.types=c("f","f","f","n"))

#process the data and make the design data list (ddl)
Col2009.process=process.data(Col2009,model="CJS",groups=c("STOCK","RELEASE","TI"))

Col2009.ddl=make.design.data(Col2009.process, remove.unused=TRUE)

#add test variable to distinguish tags with transmission interval TI of 60 vers 90
Col2009.ddl$p$t90=0
Col2009.ddl$p$t90[Col2009.ddl$p$TI=="90"]=1

Col2009.ddl$p$t60=0
Col2009.ddl$p$t60[Col2009.ddl$p$t90=="0"]=1

Col2009.models=function()
{

Phi.time.stock=list(formula=~-1+time:STOCK)

#---p---
p.time.plusTI=list(formula=~-1+time+TI) #SEs on betas reasonable
p.time.plusTI.plusrelease=list(formula=~-1+time+RELEASE+TI) #SEs on betas reasonable
p.time.STOCK.plusTI=list(formula=~-1+time:STOCK+TI) #SEs on betas NOT reasonable
p.time.STOCK.plusTI90=list(formula=~-1+time:STOCK+TI:t90) #SEs on betas reasonable

cml=create.model.list("CJS")
results=mark.wrapper(cml,data=Col2009.process,ddl=Col2009.ddl,adjust=TRUE)
return(results)
}

[jlaake]

When you use ~-1+TI it will not create an intercept and instead will create 2 columns for each level of TI.

For example, ~TI will give the following DM for 2 parameters with different levels of TI

1 0
1 1

Whereas ~-1+TI will produce

1 0
0 1

In the first case the betas will be labelled Intercept, TI(90) and in the second case they will be labelled TI(60) and TI(90).

Why are you setting adjust=TRUE? My guess is that you are over-fitting and some parameters are at boundaries and by using adjust=TRUE you are using the count of parameters from MARK which often excludes parameters at the boundary. When you over-fit or get parameters at boundaries then se for betas can be unreliable.

--jeff

[aswea]

Thank you for replying Jeff!

I know that an intercept will code for one of the group levels-- but I am not using intercepts in any of my models (see example code from first post).

Here is a sample of the betas for the model with detection probability parameterized with an interaction: p(~-1 + time:STOCK + TI). You can see that the SEs for ALL the p terms are 0 (not just ones at boundaries).

p:TI60 2.7274687 0.000000e+00 2.7274687 2.7274687
p:TI90 1.9956006 0.000000e+00 1.9956006 1.9956006
p:time6:STOCKBARGE -1.1877534 0.000000e+00 -1.1877534 -1.1877534
p:time7:STOCKBARGE -3.2664611 0.000000e+00 -3.2664611 -3.2664611
p:time8:STOCKBARGE -4.1925196 0.000000e+00 -4.1925196 -4.1925196
....

and a sample of the betas with recapture probability parameterized without an interaction: p(~-1 + time + TI)
p:time7 -0.5827448 2.982403e-01 -1.1672958 0.0018061
p:time8 -1.4911135 2.782952e-01 -2.0365720 -0.9456549
p:TI90 -0.7127897 1.239436e-01 -0.9557191 -0.4698604

Thank you again,
~Aswea

[cschwarz@stat.sfu.ca]

This is an artefact of how design matrices are created in R given a
model formula. The rules are:
- create lower order (i.e non-interaction) terms first regardless of
what is specified in the model formula
- don't include terms that are "redundant" with lower order terms
involving the same factors according to default "contrast" rules.

This is quite technical and is summarized in "stat-speak" on pages
149/150 of Modern Applied Statistics with S, 4th edition by Venables
and Ripley. Not for the faint of heart! The gist of what happens is below.

So in the formula: -1+time+RELEASE+TI, it first drops the intercept.
Then it tries to create the design matrix for the "time" factor. It
notes that without an intercept, column corresponding to ALL levels of
time can be created (7 levels corresponding to t2, t3, ...t8). Then it
tries to create the design matrix for RELEASE. There are 3 release
levels. But if you create 3 dummy variables for RELEASE, your design
matrix is "linear dependent" (redundant) because the 3 columns for release can be "added" to equal the "sum" of the 7 columns for time. Hence it creates only 2 columns for RELEASE. Similarly for T1. It notices that 2 columns for T1 are "redundant" and creates a single "T1" column.

However, in the formula, =~-1+time:STOCK+TI, the "simple" terms are
done first, so the formula is "reorganized" to be
=-1 + T1 + time:stock. When it tries to create the design matrix for
T1, it notices that because the intercept is dropped, it can create 2
columns for T1 and so you get the two column. When it tries to create
the time:STOCK part of the design matrix, it only looks to see if
time:STOCK terms are redundant with column for time or stock, it
DOESN'T look at the T1 columns.

When you have the T1:t90 interaction, it is created AFTER the
time:stock interaction because it is not simpler than the former
interaction, and so detect the redundancy in the columns and only
creates 1 column rather than the 2 expected.

For example, cut and paste the following R-code:
tempdf <- data.frame(A=c("1","1", "2","2"), B=c("1", "2", "1",
"2"),C=c("1","1","2","2"))
model.matrix( ~A, data=tempdf)
model.matrix( ~-1+A, data=tempdf)
model.matrix( ~-1+A+B+A:B, data=tempdf)
model.matrix( ~-1+A:B, data=tempdf)
model.matrix( ~-1+C+A:B, data=tempdf)
model.matrix( ~-1+A+A:B, data=tempdf)

The tempdf creates a little dataframe with 4 observations and 3
factors, but factor C is the SAME as factor A.

If you look at the output from the last 2 model.matrix statements, the
model ~1+A+A:B is able to recognize that some columns of A:B are
redundant against the columns for A, but is unable to recognize that
some columns of A:B are redundant against C (even though C is
identical to A).

Notice that the model.matrix for ~-1+C+A:B is "linearly dependent" because the first two columns added together = last 4 columns added together.

Try fitting the same model without dropping the intercept term as the
intercept term is "always" considered to be part of more complex
terms. The betas are harder to interpret as they are "contrast" vs the
intercept, but the real parameters should be ok.

As well, non-hierarchical models, i.e. including the Stock:time term
without a Stock or time term can lead to very confusing design
matrices with the beta terms not directly interpretable.


Carl Schwarz

[cooch]

This is an artefact of how design matrices are created in R given a
model formula.

It is worth noting that none of this is an issue when using standard MARK, since the structure of the DM is absolutely transparent, since you need to code it explicitly. I absolutely agree with many of the motivations for RMARK, and it is an exceptionally powerful and well-crafted alternative to the standard MARK DM/PIM paradigm, which poses limits and constraints (and frustrations) for some. However, the moment you use an approach (RMARK or otherwise) which 'shields' you from the 'drudge work' of doing things by hand, it presumes that you know exactly what you are doing - or, in this instance, what R is doing.

So - classic MARK - run the risk of making a mistake in the DM. RMARK - run the risk of not fully understanding how R actually processes things.

Moral - there is no substitute for full understanding. You either understand the DM, and build it yourself, or understand the DM that is built for you. In either instance, you need to do some work.

[jlaake]

Thanks to Carl for a very thorough reply and I agree with Evan that you need to understand the DM (design matrices) that is being constructed with the formula. That is partly why RMark is described in an Appendix in the MARK book. You cannot avoid understanding DMs with RMark. You can only avoid creating them manually. Most biologists with a few basic stats courses have never encountered DMs and the nature of DMs in capture-recapture analyses can get very messy and involved due to the nature of the data (eg triangular PIMS).

If you are unclear about the DM that RMark is creating then you can look at the DM RMark creates or if you don't have individual covariates, you can use the function model.matrix with the design data and the formula and look at it directly. Simply switching the order in the formula will get you what you want when you have 2 or more additive factor variables. Also, as Carl says there is no reason to use the ~-1 in general. That is only useful in the case where you have partially crossed factors and you want the interaction model (e.g age:year in a CJS model where initial releases are all of the same age). In that latter case where fac is a factor variable using ~-1+fac+age:year will have one too many columns in the DM and it will be necessary to use the remove.intercept argument. Many of these comments are contained in the Appendix and the workshop notes posted on the RMark website. But I understand that some of this is reasonably complex and difficult to learn at first. If you find places in those documents that are not clear or could use better description, please let me know. --jeff

[aswea]

Thank you very much for these replies. It is going to take me awhile work through them.

~Aswea

[jlaake]

I want to correct something I said in my earlier post. adjust=TRUE is the default value and the best one to use unless you are going to count parameters. That setting tells RMark to adjust the parameter count to full rank (all parameters are estimable) from what ever MARK had counted. I had presumed that because they set its value in the function call that it was the non-default value. Sounds funny given that it is my code but I've always had trouble remembering that argument meaning because it can be taken to be either adjusting by RMark or by MARK. I should have chosen something like fullrank=T/F which would have been more obvious. Thanks to Mike Melmychuk for pointing this out to me in an off-list emali. Sorry if I caused any confusion. --jeff