  Determining spatial distributions of species and communities are key
objectives of ecology and conservation. Joint species distribution models use
multi-species detection-nondetection data to estimate species and community
distributions. The analysis of such data is complicated by residual
correlations between species, imperfect detection, and spatial autocorrelation.
While methods exist to accommodate each of these complexities, there are few
examples in the literature that address and explore all three complexities
simultaneously. Here we developed a spatial factor multi-species occupancy
model to explicitly account for species correlations, imperfect detection, and
spatial autocorrelation. The proposed model uses a spatial factor dimension
reduction approach and Nearest Neighbor Gaussian Processes to ensure
computational efficiency for data sets with both a large number of species
(e.g., > 100) and spatial locations (e.g., 100,000). We compare the proposed
model performance to five candidate models, each addressing a subset of the
three complexities. We implemented the proposed and competing models in the
spOccupancy software, designed to facilitate application via an accessible,
well-documented, and open-source R package. Using simulations, we found
ignoring the three complexities when present leads to inferior model predictive
performance. Using a case study on 98 bird species across the continental US,
the spatial factor multi-species occupancy model had the highest predictive
performance among the candidate models. Our proposed framework, together with
its implementation in spOccupancy, serves as a user-friendly tool to understand
spatial variation in species distributions and biodiversity metrics while
addressing common complexities in multi-species detection-nondetection data.
