We developed an open source, extensible Python-based framework, that we call the Versatile Modeling of Deformation (VMOD), for forward and inverse modeling of crustal deformation sources. VMOD abstracts from specific source model implementations, data types and inversion methods. We implement the most common geodetic source models which can be combined to model and analyze multi-source deformation. VMOD supports Global Navigation Satellite System (GNSS), InSAR, electronic distance measurement, Leveling and tilt data. To infer source characteristics from observations, VMOD implements non-linear least squares and Markov Chain Monte-Carlo Bayesian inversions, including joint inversions using different sources of data. VMOD's structure allows for easy integration of new geodetic models, data types, and inversion strategies. We benchmark the forward models against other published results and the inversion approaches against other implementations. We apply VMOD to analyze deformation at Unimak Island, Alaska, observed with continuous and campaign GNSS, and ascending and descending InSAR time series generated from Sentinel-1 satellite radar acquisitions. These data show an inflation pattern at Westdahl volcano and subsidence at Fisher Caldera. We use VMOD to test a range of source models by jointly inverting the GNSS and InSAR data sets. Our final model simultaneously constrains the parameters of two sources. Our results reveal a depressurizing spheroid under Fisher Caldera ∼4–6 km deep, contracting at a rate of ∼2–3 Mm³/yr, and a pressurizing spherical source underneath Westdahl volcano ∼6–8 km deep, inflating at ∼5 Mm³/yr. This and past applications of VMOD to volcanic unrest benefit from an extensible framework which supports jointly inversions of data sets for parameters of easily composable multi-source models.