  We present a method to simultaneously model the dust far-infrared spectral
energy distribution (SED) and the total infrared $-$ carbon monoxide (CO)
integrated intensity $(S_{\rm IR}-I_{\rm CO})$ relationship. The modelling
employs a hierarchical Bayesian (HB) technique to estimate the dust surface
density, temperature ($T_{\rm eff}$), and spectral index at each pixel from the
observed far-infrared (FIR) maps. Additionally, given the corresponding CO map,
the method simultaneously estimates the slope and intercept between the FIR and
CO intensities, which are global properties of the observed source. The model
accounts for correlated and uncorrelated uncertainties, such as those present
in Herschel observations. Using synthetic datasets, we demonstrate the accuracy
of the HB method, and contrast the results with common non-hierarchical fitting
methods. As an initial application, we model the dust and gas on 100 pc scales
in the Magellanic Clouds from Herschel FIR and NANTEN CO observations. The
slopes of the $\log S_{\rm IR}-\log I_{\rm CO}$ relationship are similar in
both galaxies, falling in the range 1.1$-$1.7. However, in the SMC the
intercept is nearly 3 times higher, which can be explained by its lower
metallicity than the LMC, resulting in a larger $S_{\rm IR}$ per unit $I_{\rm
CO}$. The HB modelling evidences an increase in $T_{\rm eff}$ in regions with
the highest $I_{\rm CO}$ in the LMC. This may be due to enhanced dust heating
in the densest molecular regions from young stars. Such simultaneous dust and
gas modelling may reveal variations in the properties of the ISM and its
association with other galactic characteristics, such as star formation rates
and/or metallicities.
