We present a new method to infer the mass composition of ultra-high energy cosmic rays performing simultaneous modifications of the simulated depth of shower maximu ($X_{\rm max}$) and the ground signal. We introduce a likelihood fit to two-dimensional distributions of "hybrid observables", $X_{\rm max}$ and the ground signal induced by muons and electromagnetic particles at different zenith angles. With this method, the composition fractions of primary particles are simultaneously estimated together with a shift of $X_{\rm max}$ and multiplicative rescaling parameters of the muon and electromagnetic signals in the ground detector. The performance of the method is presented using Monte Carlo simulations showing the importance of a global description of the hybrid data from ground and optical detectors resulting in reduced systematic uncertainties stemming from the modelling of hadronic interactions.