Plot of predicted covariate effect on occupancy probability with stPGOcc object?

[Christopher Dennison]

Good afternoon everyone,

I hope this email finds you well. I am working with spOccupancy, and am wondering if I might ask for guidance on plotting.

I have successfully fit a spatial model using stPGOcc, and would like to use said model to create a plot of the predicted effect of the primary occupancy covariate in the model on occupancy of the species in question.

I have done this previously with unmarked using the modavgPred function in AICcmodavg, and assume that I would use the  predict.stPGOcc() function in spOccupancy to create a data frame of predicted values then ggplot?

I have attached a figure of the type of plot I am gunning for.

This humble graduate student thanks you for your time and support. I have thoroughly enjoyed using and learning through spOccupancy thus far.

Warm regards,

- Chris

[Jeffrey Doser]

Hi Chris,

Yes, your intuition is correct on how to make the plot. You will need to predict across a sequence of values for the covariate you're interested in (while holding other covariates constant), and then you can make the plot in ggplot using the prediction output. I've attached a script that does this with stPGOcc and an example bird species in the attached data. Hopefully you can adapt the code in there to make the plot for your data, but let me know if you run into any other problems.

Cheers,

Jeff

--
You received this message because you are subscribed to the Google Groups "spOccupancy and spAbundance users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to spocc-spabund-u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/spocc-spabund-users/e96ef15d-a805-4410-aab4-427162bc4f4an%40googlegroups.com.

--

Jeffrey W. Doser, Ph.D.

Assistant Professor

Department of Forestry and Environmental Resources

North Carolina State University

Statistical Ecology and Forest Science Lab

Pronouns: he/him/his

[Christopher Dennison]

Thank you so much, Jeff! Your code was extremely helpful, and I was able to get it working with my data.

Thank you again for all of your help. I truly appreciate your time and support.

[Christopher Dennison]

Hey again, Jeff. Sorry for the flurry of posts. I did have one question:

I am using a random effect of year in my occupancy covariates (INBU.occ.formula <- ~ scale(vl100m) + (1 | years)), and I am wondering how I should define this in my pred.def object, or if I need to?

I hope that makes some sense! Thank you so much.

[Christopher Dennison]

I believe the headache might be in the following code from your example. Specifically, when I try to run this code:

pred.df <- as.matrix(data.frame(intercept = 1, regen = regen.pred.vals.scale,
                                basalArea = 0, basalArea.2 = 0,
                                percentForest = 0))
X.0 <- array(pred.df, dim = c(nrow(pred.df), 1, ncol(pred.df)))
dimnames(X.0)[[3]] <- c('(Intercept)', 'scale(regen)', 'scale(basalArea)',
                            'I(scale(basalArea)^2)', 'scale(percentForest)',
                            'scale(year)')

I get the error message: Error in dimnames(X.0)[[3]] <- c("(Intercept)", "scale(regen)", "scale(basalArea)",  :
  length of 'dimnames' [3] not equal to array extent.

Perhaps this relates to the random effect?

[Jeffrey Doser]

Hi Chris, 

Oops, that error you encountered in the script is a typo on my end. You can get rid of the 'scale(year)') in the last line of the code there (that was just left over from changing the code from something else). 

You are correct that you will need to do something slightly different to account for the fact that you have a random effect in the model. You will want to do one of two things: 

1. Include "years" as a variable in the prediction "design matrix" and just set it equal to one of the years in your data set. So, if you adapt the code I sent in that script, you could add "years" into the "pred.df" object and just set it equal to one of the years in your data set. 
2. Alternatively, you can set the "ignore.RE" argument in predict() to TRUE. This will generate the predictions only with your covariate of interest. 

Given that your goal is to create a marginal effects plot for your covariate effect, they both should work fine, but the second one will likely be a bit easier. 

Cheers,

Jeff