RMark - covariate.predictions

[vin]

I encountered following problem when plotting survival against an individual covariate using the function ‘covariate.predictions’:

I have as a group variable: population; coded as “pop1” and “pop2” and as individual covariate: “speed”. We expect survival to be higher when an individuals speed is higher and that the effect differs between the two populations. I build following model:

do.model=function(){
Phi.speed=list(formula=~population+speed)
p.time=list(formula=~time)
cml=create.model.list("CJS")
results=mark.wrapper(cml,data=hag.process,ddl=hag.ddl,adjust=FALSE)
return(results)}
modelcov=do.model()

When I look at the output file, the PIM’s for Phi are :

INPUT --- group=1 Phi rows=4 cols=4 Triang ;
INPUT --- 1 1 1 1 ;
INPUT --- 1 1 1 ;
INPUT --- 1 1 ;
INPUT --- 1 ;

INPUT --- group=2 Phi rows=4 cols=4 Triang ;
INPUT --- 2 2 2 2 ;
INPUT --- 2 2 2 ;
INPUT --- 2 2 ;
INPUT --- 2 ;

Now when I use

Phispeed=covariate.predictions(modelcov,data=data.frame(speed=speed.values),indices=c(1,2))

I get EXACTLY the same survival estimates in both populations (index 1 and 2 if I understand correctly). While they should differ, based on the real beta estimates. I even tried to recode my population as 0 and 1, but with no improvement. Can anyone suggest what I do wrong here?


Many thanks in advance,
Vincent.

[dhewitt]

Vincent,

Jeff can correct me if I screw up here, but I think what you've run into is a very important difference between the simplified PIMS (parameter indices) and the all-different PIMS as they relate to RMark internal functions like covariate.predictions(). I'm not sure this difference is as well known as it ought to be among RMark users, so I'm glad you brought it up.

RMark simplifies the PIMS for a given model whenever possible before sending to MARK. So when you look at the output file (generated by MARK), it only shows you what it knows -- the simplified indices. **Within RMark, you need the all-different indices for the function covariate.predictions().**

Without seeing your actual RMark-created objects, I'm not sure if this is your issue. But, you can get the all-different indices with the PIMS() function (see ?PIMS). Try something like this:

PIMS(modelcov[[1]], "Phi", simplified=FALSE)

Dropping the simplified argument goes with the default of TRUE and you'll get what was sent to MARK (and what is shown in the output file). If these two sets of PIMS are different, that will probably put you on the right track.

Please re-post with a resolution, or let the list know if I am off track.

- Dave Hewitt

[vin]

Thanks for the quick response Dave. Your comment is wright on track and very helpfull. Indeed when I use the PIM function with the “simplified=FALSE” option I got a different set of PIMS.

INPUT --- group=1 Phi rows=4 cols=4 Triang ;
INPUT --- 1 2 3 4 ;
INPUT --- 5 6 7 ;
INPUT --- 8 9 ;
INPUT --- 10 ;

INPUT --- group=2 Phi rows=4 cols=4 Triang ;
INPUT --- 11 12 13 14 ;
INPUT --- 15 16 17 ;
INPUT --- 18 19 ;
INPUT --- 20 ;

Meaning I have to use

Phispeed=covariate.predictions(modelcov,data=data.frame(speed=speed.values),indices=c(1,11)) to get the survival estimates for the two different populations.

Vincent

[jlaake]

Thanks to Dave for answering this post. One clarification is that the PIMS function may not work with every type of model. It is fine for CJS and similar models. PIMS was created initially as a teaching tool as a demostration for those used to working wth pims in MARK. The reason that covariate.predictions uses the all-different indices is because it will do model averaging and the only way to specify parameters uniquely is with the all-different formulation; otherwise for example parameter 2 (using the simplified indices) may be for the second time for group 1 in one model and then in a different model it may be the constant value for group 2 like in the proposed model. But the all-different indices are constant across all models.

Another more efficient way to get the all-different indices is to use the numeric equivalents of the row.names of the design data. This allows you to select the indices based on the design data without looking them up. For example,

as.numeric(row.names(ddl$Phi[ddl$Phi$group=="M",]))

would pick out all the indices for Phi for group=="M" if that value was used as a grouping variable.

This can be a bit confusing and I have plans to write some functions to help simplify this. regards --jeff