Many critical ecological issues require the analysis of large spatial point data sets - for example, modelling species distributions, abundance and spread from survey data. But modelling spatial relationships, especially in large point data sets, presents major computational challenges. We use a novel Bayesian hierarchical statistical approach, 'spatial predictive process' modelling, to predict the distribution of a major invasive plant species, Celastrus orbiculatus, in the northeastern USA. The model runs orders of magnitude faster than traditional geostatistical models on a large data set of c. 4000 points, and performs better than generalized linear models, generalized additive models and geographically weighted regression in cross-validation. We also use this approach to model simultaneously the distributions of a set of four major invasive species in a spatially explicit multivariate model. This multispecies analysis demonstrates that some pairs of species exhibit negative residual spatial covariation, suggesting potential competitive interaction or divergent responses to unmeasured factors.