Direct imaging searches for exoplanets around stars detect many spurious candidates that are in fact background field stars. To help distinguish these from genuine companions, multi-epoch astrometry can be used to identify a common proper motion with the host star. Although this is frequently done, many approaches lack an appropriate model for the motions of the background population, or do not use a statistical framework to properly quantify the results. For this study we used Gaia astrometry combined with 2MASS photometry to model the parallax and proper motion distributions of field stars around exoplanet host stars as a function of candidate magnitude. We developed a likelihood-based method that compares the positions of a candidate at multiple epochs with the positions expected under both this field star model and a co-moving companion model. Our method propagates the covariances in the Gaia astrometry and the candidate positions. True companions are assumed to have long periods compared to the observational baseline, so we currently neglect orbital motion. We applied our method to a sample of 23 host stars with 263 candidates identified in the B-Star Exoplanet Abundance Study (BEAST) survey on VLT/SPHERE. We identified seven candidates in which the odds ratio favours the co-moving companion model by a factor of 100 or more. Most of these detections are based on only two or three epochs separated by less than three years, so further epochs should be obtained to reassess the companion probabilities. Our method is publicly available as an open-source python package from 1em