SpAbundance - msNMix Prediction Plot Question

[Harry]

Dear Jeffrey,

Hope you are well! I am currently working on a multi species N-Mixture model with two abundance covariates, one abundance factor and three detection covariates, with my MSc supervisor Todd Lewis who posted on here recently.

I am trying to plot the predicted abundance of the two abundance covariates separately, but I am having some trouble doing so. My code seems to combine the covariates together and I can't find a way to explicitly choose which covariate to plot. I think the issue may occur due to both plots using the same mu.med values for the y-axis and having the same quantile values, however I am unsure how to change this.

I have attached the code below, any help would be greatly appreciated!

Thank you! All the best,

Harry

 ############ Pred Modelling Abundance  ############ 

# Define the first covariate values

cov.pred.vals.A <- seq(min(dataNMixSim$abund.covs[, 'abund.cov.1']),

                       max(dataNMixSim$abund.covs[, 'abund.cov.1']),

                       length.out = 50)

 

# Scale the first covariate values

cov.pred.vals.As <- (cov.pred.vals.A - mean(dataNMixSim$abund.covs[, 'abund.cov.1'])) /

  sd(dataNMixSim$abund.covs[, 'abund.cov.1'])

 

# Define the second covariate values

cov.pred.vals.B <- seq(min(dataNMixSim$abund.covs[, 'abund.cov.2']),

                       max(dataNMixSim$abund.covs[, 'abund.cov.2']),

                       length.out = 50)

 

# Scale the second covariate values

cov.pred.vals.Bs <- (cov.pred.vals.B - mean(dataNMixSim$abund.covs[, 'abund.cov.2'])) /

  sd(dataNMixSim$abund.covs[, 'abund.cov.2'])

 

#Put into Matrix

X.A <- as.matrix(data.frame(Intercept = 1, abund.cov.1 = cov.pred.vals.As, abund.cov.2 = cov.pred.vals.Bs))

 

# Check the resulting matrix

head(X.A)

dim(X.A)

 

# Run predictive model for abundance estimation

abund.pred <- predict(mod.NB, X.0 = X.A, ignore.RE = TRUE, type = 'abundance')

 

#Apply confidence interval range

mu.0.quants.A <- apply(abund.pred$mu.0.samples, c(2, 3), quantile,

                     prob = c(0.025, 0.5, 0.975))

 

#specify number of species

n.sp <- 16

 

#Assemble as dataframe for plotting ggplot2

mu.0.plot.A <- data.frame(mu.med = c(mu.0.quants.A[2, , ]),

                           mu.low = c(mu.0.quants.A[1, , ]),

                           mu.high = c(mu.0.quants.A[3, , ]),

                           sp = rep(1:n.sp, nrow(X.A)),

                           abund.cov.1 = rep(cov.pred.vals.A, each = n.sp))

 

# Example species names vector

sp_names <- c("AgCl","BaRu","CrFz","CrNb","CrPs","CrPy","DnPh","DiDa","InCn","LpMl","LpSv","OpPu","RnVa","RhMr","ScEl","SmBd")

 

# Assuming `sp` in your data frame is a numeric vector from 1 to 16

mu.0.plot.A$sp_name <- sp_names[mu.0.plot.A$sp]

 

#Plot prediction with covariate influence

P1 <- ggplot(mu.0.plot.A, aes(x = abund.cov.1, y = mu.med)) +

  geom_ribbon(aes(ymin = mu.low, ymax = mu.high), fill = 'grey70') +

  geom_line() +

  theme_bw() +

  facet_wrap(vars(sp_name), scales = 'free_y') +

  labs(x = 'Covariate', y = 'Expected abundance')

P1

########################################

[Jeffrey Doser]

Hi Harry,

Thanks for the note. Your code is definitely on the right track. What you need to do to generate a conditional probability plot for each covariate effect is do the prediction (i.e., generate mu.0.quants.A) for one covariate at a time while holding the other covariate constant at some value. So, if you want to first generate a covariate for abund.cov.1, then you should use the same approach for generating "cov.pred.vals.As", but when you generate "cov.pred.vals.Bs", you should set that to some constant value (i.e., the most common choice is to set it at the mean, which if you standardized everything is simply 0). Then do the prediction to generate the set of "mu.0.quants" for making the plot for the first covariate. To get a plot for the second covariate, you then need to do the whole process over again, but this time when you form "cov.pred.vals.As" and "cov.pred.vals.Bs", "cov.pred.vals.As" should be constant and then "cov.pred.vals.Bs" should be generated like how you did in your code that you sent. Then the resulting "mu.0.quants" that come from that prediction can be used to make the plot for cov.pred.vals.Bs.

I hope that makes sense! Let me know if I can clarify anything.

Kind regards,

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/f719b5d6-570c-42b1-a8bd-e279185f40edn%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