www.phidot.org

Greetings forum,

I need some help with figuring out how to extract nest survival estimate results when using model averaging, which may involve using the covariate.predictions function.

I've run model selection in RMark and my covariates include Year, Habitat, AgeDay1, Incub (a time-dependent nest stage variable as explained here http://www.phidot.org/forum/viewtopic.php?f=21&t=2613&p=8286&hilit=nest+stage#p8286), AgeFound, NestAge, Time, and Edge)

I performed model averaging on these results:

Code: Select all

## Model averaging
my.avg=model.average(myrexs.results, vcv=TRUE)$estimates
names(my.avg)
my.avg$estimate  # (returns list of >1000 numbered estimates)
my.avg$Hab$results$real  #  (returns "NULL" output)

I am confused about how to view specific model averaged results to obtain my nest survival estimates. What are the commands to see estimates based on factor variables like Habitat? Before I used model averaging, the last line of code would work for this.

I saw that Jeff recommended using covariate.predictions to obtain the nest stage results, but I haven't been able to execute it. The examples in the RMark guide and workshop notes involve estimating individual covariates at specific values or their means, but I don't have numeric individual covariates and I've not been able to figure out how to apply this to get an estimate for the incubation (<17 days) and nestling stages.

Any guidance here would be awesome
Thanks,
Deb

I realize there is only one parameter type for this model but as the help explains for "model.average.marklist", if parameter=NULL (the default), the design data are not copied over with the estimates. You can append (cbind) the estimates you have now to the design data but if you use code below it should work as well

Code: Select all

## Model averaging
my.avg=model.average(myrexs.results, parameter="S",vcv=TRUE)$estimates
names(my.avg)
my.avg$estimate  # (returns list of >1000 numbered estimates)


The following only works with a model list and not with the result of model.average which averages the estimates across all models in the list. When you did names(my.avg) it should have shown only 2 and Hab wouldn't have been one of the names so the code below had no chance to work. Also, you don't want to use the code below in practice if you scale for overdispersion. Instead use summary(model)$reals which will scale the std errors for over-dispersion.

Code: Select all

my.avg$Hab$results$real  #  (returns "NULL" output)


You are correct that you only need covariate.predictions if you have individual numeric covariates which is what that person had in my response. I think it was % bare ground. See if the above gets you to where you want and if not post back.

Hi Jeff,

Thank you so much for the quick reply! I've corrected the mistake of leaving out the parameter specification. Thanks for pointing that out. Unfortunately, I still cannot get the output to subset usefully. I hope this isn't too basic of a question! For example, I want to see the nest survival estimates for each Habitat type. I'm also not clear on how I should proceed for more complicated (time-dependent) variables like nest stage. I know how to plot DSR by nest age values (without model averaging) to view the results, so do I also need to do something similar for nest stage? I named the stage variable "Incub" but I don't see it under the names:

Code: Select all

>my.avg=model.average(myrexs.results, parameter="S", vcv=TRUE)$estimate
my.avg> names(my.avg)
 [1] "par.index"   "estimate"    "se"          "lcl"         "ucl"       
 [6] "fixed"       "note"        "par.index"   "model.index" "group"     
[11] "age"         "time"        "Age"         "Time"        "Habitat"   
[16] "Year

My attempt at subsetting gives me every possible estimate, but I haven't figured out how to get them to group, for example by Habitat:

Code: Select all

> head(subset(my.avg,select=c("Habitat","estimate","se","lcl","ucl")))
                   Habitat  estimate         se       lcl       ucl
S gPenin2010 a0 t1   Penin 0.9009308 0.04570073 0.7692454 0.9612522
S gPenin2010 a1 t2   Penin 0.9008257 0.04539183 0.7704053 0.9609194
S gPenin2010 a2 t3   Penin 0.9007200 0.04507956 0.7715716 0.9605811
S gPenin2010 a3 t4   Penin 0.9006135 0.04476386 0.7727441 0.9602372
S gPenin2010 a4 t5   Penin 0.9005063 0.04444472 0.7739227 0.9598876
S gPenin2010 a5 t6   Penin 0.9003984 0.04412212 0.7751074 0.9595322

I'd greatly appreciate any further help,
Deb

For example, I want to see the nest survival estimates for each Habitat type.

I don't know what you mean with the above. Are you saying the estimate of surviving the entire interval? From what you showed, those are the survival estimates by time.

If you have a variable that you added to the design data that is not a grouping variable then you should be able to extract that column from the design data you used and append (cbind) it to my.avg. Or you can get the estimates like you did originally and then append them to your design data.

Maybe someone else on the list that has used nest survival analysis with RMark can help you further.

regards --jeff