At higher energies the uncertainty in the estimated cosmic ray mass composition, extracted from the observed distributions of the depth of shower maximum $X_{\mathrm{max}}$, is dominated by uncertainties in the hadronic interaction models. Thus, the estimated composition depends strongly on the particular model used for its interpretation. To reduce this model dependency in the interpretation of the mass composition, we have developed a novel approach which involves fitting some of the model parameters to the data. Using this approach we extracted cosmic ray mass composition information from the published Pierre Auger $X_{\mathrm{max}}$ distributions. We have obtained a consistent mass composition interpretation for EPOS-LHC and QGSJetII-04. The Sibyll2.3 model provides two possible composition interpretations, one of them in agreement with EPOS-LHC and QGSJetII-04. Our fit suggests a composition dominated by protons below $10^{18.5}$ eV, while between $10^{18.5}$ and $10^{19.4}$ eV it becomes dominated by helium, and in the two energy bins above $10^{19.4}$ eV a growing nitrogen component is hinted. The fit also suggests that the normalisation level for $\langle X_{\mathrm {max}}\rangle$ is consistent with the QGSJetII-04 model prediction, but the normalization level for $\sigma(X_{\mathrm max})$ is consistent with the EPOS-LHC prediction. These results remain sensitive to the other model parameters that we keep fixed, such as the <Xmax> elongation rate and the $\langle X_{\mathrm {max}}\rangle$ separation between p and Fe.