Occupancy probability at average value of 1 covariate (single season single species occupancy model)

[Mickayla Johnston]

Hi everyone,

First, I apologize if this is a silly question-- I'm no statistician by any means. 

I've built single season, single species occupancy models for red-winged blackbirds during their breeding season following the eBird best practices methods (https://cornelllabofornithology.github.io/ebird-best-practices/occupancy.html). These methods use land cover types as occupancy and detection covariates. I've made a model for each ecoregion in North America where red-wings appear to be present (after going through a model selection process using AICs to determine the best model). So, I have 5 total models (1 for each ecoregion), all with the same occupancy / detection covariates. I've run goodness-of-fit tests on them, and they seem to be ok (except for 1 ecoregion). Now, for each model, I'm trying to determine the occupancy probability given the average value of 1 land cover type (and then repeat this for each land cover type). For example, I want to be able to say something like, "the occupancy probability of red-winged blackbirds is X, given 50% wetland land cover". 

Here's an example of code for one of my models-

                                                                                                                                       #detection covs

occ_model_north_for_mount_2019_new <- occu(~ time_observations_started + 

                                         duration_minutes + 

                                         effort_distance_km + 

                                         number_observers + 

                                         protocol_type +

                                         pland_08_woody_savanna +

                                         pland_09_savanna +

                                         pland_10_grassland +

                                         pland_11_wetland +

                                         pland_12_cropland +

                                         pland_13_urban +

                                         pland_14_mosiac

                                       # occupancy covs

                                       ~ pland_08_woody_savanna +

                                         pland_09_savanna +

                                         pland_10_grassland +

                                         pland_11_wetland +

                                         pland_12_cropland +

                                         pland_13_urban +

                                         pland_14_mosiac,

                                       data = occ_um_north_for_mount_new)

I've seen people use backTransform() with linearComb() to get the occupancy probability when all variables are held at their average, so I tried that--

backTransform(linearComb(occ_model_north_for_mount_2019_new,  coefficients = c(1,0,0,0,0,0,0,0), type = 'state'))

I get this as a result--

But this doesn't give me exactly what I'm looking for. 

Or, is there a way to do this using the predict() function instead?

Thank you in advance!

[Mickayla Johnston]

Does anyone have any suggestions for this? Thanks!

[Jeffrey Royle]

hi Mickayla,

 There are quite a few examples in Chapter 10 of the AHM book which show applications of the predict() function for occupancy models, which I think is what you want. You can provide the specific values of each covariate you want a prediction for and the matrix of predictions, SEs, and confidence intervals are returned.  e.g., see p. 572 for one example.

 (AHM book = "Applied hierarchical model in ecology.... volume 1" by Marc Kery and myself).

 Check that out and if you don't have it email me directly and I'll send you some PDF chunks.

regards

andy

[Mickayla Johnston]

Hi Andy,

Thanks so much for getting back to me. I just emailed you!