plot density as a function of covariates

[HL Y]

Dear All, 

 

I am trying to plot density as a function of each of the covariates. The density was calculated using function get.real() in the oSCR, based on model D~SESSION + covariate 1 + covariate 2, where covariate 1 is continuous and covariate 2 is categorical. So, now I have densities (per pixel per year for each sex), and lower and upper values of the densities. I think the density to the categorical covariate can be easily plotted with a boxplot? But how to plot a density curve as a function of the continuous covariate (like Fig.2 in doi: 10.1016/j.biocon.2022.109875)? If anyone knows this, could you share the code with me please? Thank you very much in advance.

 

Best wishes,

 

Hongli

Dr Hongli Yu

Nottingham Trent University

hongli...@my.ntu.ac.uk; hongli...@gmail.com

[Daniel Linden]

Hi Hongli,

What you want is a marginal effects plot, where you illustrate the change in expected density with different values of a given covariate (holding other covariates constant).  You need to create a new data frame where the covariate of interest varies from the min to the max observed values (using the seq function), and other covariates are given some average value (e.g., 0 for standardized variables).  Below is an example using the red-backed salamander models in oSCR.  It is illustrating the marginal effect of a detection covariate, but the concept is the same for density covariates.

library(oSCR); library(ggplot2)
data(rbs_ecography_mods)

newdf <- data.frame(jday=seq(-1,7,by=1),b=0,jday2=0)
# use "m1" from the rbs models

realdf <- get.real(m1,type="det",newdata=newdf)

realdf %>%
  ggplot(aes(x = jday, y = estimate, ymax = upr, ymin = lwr)) +
  geom_ribbon(alpha = 0.25) + geom_line(lwd=1)

NOTE!  This plot does not make much sense because Julian day was fit as a quadratic relationship in this model.  In this case, the proper way to view the marginal effects is to include the square of the values as well:

newdf <- data.frame(jday=seq(-1,7,by=1),b=0,jday2=seq(-1,7,by=1)^2)
realdf <- get.real(m1,type="det",newdata=newdf)

realdf %>%
  ggplot(aes(x = jday, y = estimate, ymax = upr, ymin = lwr)) +
  geom_ribbon(alpha = 0.25) + geom_line(lwd=1)

[HL Y]

Dear Daniel,

It finally worked! Thank you so much for your help!

Best regards,

Hongli

--
You received this message because you are subscribed to the Google Groups "oSCR" group.
To unsubscribe from this group and stop receiving emails from it, send an email to oscr_package...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/oscr_package/41eacf0d-7f4b-4571-8ff0-9cc98abdb185n%40googlegroups.com.