Plate boundary deformation zones represent a challenge in terms of understanding their underlying geodynamic drivers. Active deformation is well constrained by GNSS observations in the SW Balkans, Greece and W Turkey, and is characterized by variable extension and strike slip in an overall context of slow convergence of the Nubia plate relative to stable Eurasia. Diverse, and all potentially viable, forces have been proposed as the cause of the observed surface deformation, e.g., asthenospheric flow, horizontal gravitational stresses (HGSs) from lateral variations in gravitational potential energy, and rollback of the Hellenic slab. We use Bayesian inference to constrain the relative contribution of the proposed driving and resistive regional forces. Our models are spherical 2D finite element models representing vertical lithospheric averages. In addition to regional plate boundaries, the models include well-constrained fault zones like north and south branches of the North Anatolian Fault, Gulf of Corinth and faults bounding the Menderes Massif. Boundary conditions represent geodynamic processes: (1) far-field relative plate motions (2) resistive fault tractions (3) HGSs mainly from lateral variations in topography and Moho topology (4) slab pull and trench suction at subduction zones. The magnitude of each of these is a parameter in a Bayesian analysis of the models in the context of horizontal GNSS velocities. The search yields a probability distribution over all parameters, allowing us to determine mean/median parameter values, robustly estimate parameter uncertainties, and identify tradeoffs. Significant trench suction forces from the Hellenic slab act on the overriding Aegean Sea, including along the Pliny-Strabo STEP Fault. Resistive tractions on most plate boundaries and faults are low. The best-fitting models compare well with paleomagnetic rotations and fault slip rates from previous studies.