Model averaging with covariate predictions

[tasha]

Hello,
I am modeling survival estimates with groups and individual covariates. I have read and tried to follow the examples given in the Rmark chapter of the Mark book as well as the Rmark help files on covariate.predictions but I still get an error when I try to do model averaging with individual covariates.
I am using R 2.10.1 and Rmark version 1.9.7. My data file consists of individual mark recapture histories, trapping site (a plot nested within a colony), a treatment group (half the plots were food supplemented and half controls), and for each individual we have the sex and a body condition index (BCI). The BCI is the individual covariate. I am using a Robust model with four primary trapping sessions over two years.

Here is an example of the Rmark code I have used:

Survival_Trial_Adult_BC3=import.chdata("C:/Workspace/RMark/RMark All GPs/Ch_Groups_Adults_Final_BCIS3.txt",use.comments=TRUE,header=TRUE,
field.types=c("f","f","f","f","f","n","f"))

Survival_Trial_Adult_BC3.process=process.data(Survival_Trial_Adult_BC3,groups=c("Treatment","Sex","Plot","Colony"),
model="Robust",time.intervals=c(0,0,0,0,0.333,0,0,0,0,1.333,0,0,0,0,0.333,0,0,0,0))

Survival_Trial_Adult_BC3.ddl=make.design.data(Survival_Trial_Adult_BC3.process)

Survival_Trial_Adult_GammasPCTimeSexTrPlotColonyBC3.models=function()
{
S.Global=list(formula=~-1+time:Treatment:Sex:Plot:Colony:BCIS3)
S.dot=list(formula=~1)
S.time=list(formula=~time)
S.time.sex=list(formula=~time+Sex+BCIS3)
S.time.treatment=list(formula=~time+Treatment+BCIS3)
#S.etc (for all survival models)#
GammaDoublePrime.dot=list(formula=~1,share=TRUE)
p.timesessionplot.c=list(formula=~session:time+Plot+c,share=TRUE,remove.intercept=TRUE)
N.group.session=list(formula=~-1+group:session)
cml=create.model.list("Robust")
cml
results=mark.wrapper(cml,data=Survival_Trial_Adult_BC3.process,ddl=Survival_Trial_Adult_BC3.ddl)
return(results)
}
Survival_Trial_Adult_BC3.results=Survival_Trial_Adult_GammasPCTimeSexTrPlotColonyBC3.models()
Survival_Trial_Adult_BC3.results

minBCIS3=min(Survival_Trial_Adult_BC3$BCIS3)
maxBCIS3=max(Survival_Trial_Adult_BC3$BCIS3)
bci.values=minBCIS3+(0:242)*(maxBCIS3-minBCIS3)/242
S.bci3=covariate.predictions(Survival_Trial_Adult_BC3.results,data=data.frame(BCIS3=bci.values),indices="S",drop=FALSE)


And here is the error message:
Error in fixedvalues[fixedparms] = model$fixed$value[match(indices[fixedparms], : replacement has length zero


When I run this without the individual covariate the model.average command works fine and there are no fixed values for survival.

Any help or advice would be greatly appreciated. Or if you need any more information or clarification please ask.

Thanks again,
natasha

[jlaake]

In the call to covariate.predictions, you have specified indices="S" and it needs to be a vector of numeric indices for the all different PIM values because it allows you to get values of parameters of all types (eg S,p etc). See examples with ?covariate.predictions. Also, I think this is covered in the workshop notes that are available on the RMark page on phidot.

--jeff

[tasha]

Thank-you for your response Jeff,

I have tried again following the commands in the covariate.predictions help file, however I still receive an error message when trying to retrieve the all different PIM indices numbers. Here is my code and error messages:

> PIMS(Survival_Trial_Adult_GammasPCTimeSexTrPlotColonyBC3.models,"S",simplified=FALSE)
Error in model$model : object of type 'closure' is not subsettable
> summary.mark(Survival_Trial_Adult_GammasPCTimeSexTrPlotColonyBC3.models,se=TRUE)$real
Error in model$model : object of type 'closure' is not subsettable
>

I apologize if I am missing something simple here. Could this be due to the Robust model type at all?

Thanks again for your time and any more advice,
natasha

[jlaake]

You are using the name of your function as the value of the calling argument to PIMS and summary.mark. You need to specify a model which would be something like

summary.mark(Survival_Trial_Adult_BC3.results[[1]],models,se=TRUE)$real

or

PIMS(Survival_Trial_Adult_BC3.results[[1]],"S","FALSE")

In addition, there is a function extract.indices that you may find useful. I added it recently and have not added it into Appendix C or in workshop notes. If you type ?extract.indices you can see the help file which contains an example at the bottom.

--jeff