MultiNMix: A Package for Multispecies N-Mixture Models

MultinMix is an R package designed for fitting Multispecies N-Mixture (MNM) Models (Mimnagh, Niamh, et al. (2022)), a powerful tool for estimating abundance and occurrence of multiple species in a hierarchical Bayesian framework.

Features

• Bayesian Modeling: Fit hierarchical Bayesian MNM models using Nimble.
• Customisable Priors: Define prior distributions easily for each parameter.
• Comprehensive Outputs: Includes posterior summaries, convergence diagnostics, and model fit statistics (log-likelihood, AIC, BIC).
• User-Friendly API: Simple interface to specify data, initial values, and model parameters.
• Visualisation: Built-in methods for producing density plots and traceplots, for model diagnostics.

Installation

To install the development version of MultiNMix, use the following commands in R:

 if (!requireNamespace("devtools", quietly = TRUE)) {
 install.packages("devtools")
 
devtools::install_github("niamhmimnagh/MultiNMix")

Getting Started

Here is a quick example to get you started with MultiNMix:

 library(MultinMix)
 
 # Example data
x <- simulateData(model = "MNM")
 R<-x$R
 T<-x$T
 S<-x$S
 K<-x$K
 
 Xp <- array(rnorm(R * S * 2), dim = c(R, S, 2)) # creating 2 detection probability covariates
 Xn <- array(rnorm(R * S *3), dim = c(R, S, 3)) # creating 3 abundance covariates
 
 # Fit 
fit <- MNM_fit(
 Y = species_counts,
 Xp = Xp,
 Xn = Xn,
 Hurdle=FALSE,
 AR = FALSE,
 iterations = 5000, # Number of iterations
 burnin = 1000, # Burn-in period
 thin = 10, # Thinning interval
 prior_detection_probability="dnorm(0,0.01)" # user-defined normal prior distribution
)
 
 # Summarize results
fit@summary
 
 # Plot diagnostic results by specifying the model and the parameter
 tracePlot(fit, param="N[8,1]")
 density(fit, param="N[8,1]")
 
 # A list of all available diagnostic plots can be found:
 View(y@plot)

Functions

Main Function

• MNM_fit(): Fits a Multispecies N-Mixture Model using specified data and parameters.

Utility Functions

• tracePlot(): Generates traceplots of monitored parameters.
• density(): Generates density plots of monitored parameters.
• logLik(): Extracts the log-likelihood of the model.
• AIC(), BIC(): Computes AIC and BIC values for model comparison.
• check_convergence(): Assesses model convergence using Gelman-Rubin diagnostics.

Documentation

Detailed documentation and vignettes are available in the package. After installation, access them using:

??MultiNMix

Datasets

There are two datasets available in the package birds and the zero-inflated birds_ZI. Both are a subset of the North American Breeding Bird Survey dataset (https://www.pwrc.usgs.gov/BBS/). birds is a dataframe with 2,880 observations and 13 columns (R=24, T=10, S=20, K=6) while birds_ZI is a dataframe with 600 observations and 13 columns (R=15, T=10, S=10, K=4).

In this vignette, we will show the birds dataset, the processing steps required and a worked example of it.

The birds Dataset

 data(birds)
 head(birds)

The birds dataset is currently a data frame of dimension (600, 10). It needs to be reformatted into an array of dimension (R=15, T=10, S=10, K=4) before it can be used with the MultiNMix functions.

 # Data must first be reformatted to an array of dimension (R,T,S,K)
 R <- 15
 T <- 10
 S <- 10
 K <- 4
 
 # Ensure data is ordered consistently
 birds <- birds[order(birds$Route, birds$Year, birds$English_Common_Name), ]
 
 # Create a 4D array with proper dimension
 Y <- array(NA, dim = c(R, T, S, K))
 
 # Map route, species, and year to indices
 route_idx <- as.numeric(factor(birds$Route))
 species_idx <- as.numeric(factor(birds$English_Common_Name))
 year_idx <- as.numeric(factor(birds$Year))
 
 # Populate the array
 stop_data <- as.matrix(birds[, grep("^Stop", colnames(birds))])
 
 for (i in seq_len(nrow(birds))) {
 Y[route_idx[i], , species_idx[i], year_idx[i]] <- stop_data[i, ]
 }
 
 # Assign dimnames
 dimnames(Y) <- list(
 Route = sort(unique(birds$Route)),
 Stop = paste0("Stop", 1:T),
 Species = sort(unique(birds$English_Common_Name)),
 Year = sort(unique(birds$Year)))

The function MNM_fit in the MultiNMix package allows for easy implementation of a multi-species N-mixture model using data of this format.

model<-MNM_fit((Y=Y, AR=FALSE, Hurdle=FALSE))

We can then access elements of the model as follows:

model@summary # outputs the mean estimate, standard deviation, standard error, 95% credible interval, effective sample size and gelman rubin statistic for each monitored variable
 
model@estimates$N # outputs the estimated mean N
 
 logLik(model) # estimates the log likelihood of the model
 
 AIC(model)/BIC(model) # outputs the AIC or BIC values
 
 tracePlot(model, param="N[1,1]") # outputs the traceplot of the N[1,1] parameter
 
 density(model, param="N[1,1]") #outputs the density plot for the N[1,1] parameter

Contributions

Contributions are welcome! If you encounter any issues or have suggestions for improvement, please submit a report or a pull request.

References

Mimnagh, Niamh, et al. "Bayesian multi-species N-mixture models for unmarked animal communities." Environmental and Ecological Statistics 29.4 (2022): 755-778.

Acknowledgements

MultiNMix was developed as part of research into multispecies abundance modeling. Special thanks to the creators of Nimble (r-nimble.org) for their invaluable tools in Bayesian modeling.