Fixing the relationship among a set of experimental quantities is a fundamental issue in many scientific disciplines. In the two-dimensional case, the classical approach is to compute the linear correlation coefficient from a scatterplot. This method, however, implicitly assumes a linear relationship between the variables. Such assumption is not always correct. With the use of the partial correlation coefficients, an extension to the multi-dimensional case is possible. However, the problem of the assumed mutual linear relationship among the variables still remains. A relatively recent approach which permits to avoid this problem is modeling the joint probability density function (PDF) of the data with copulas. These are functions which contain all the information on the relationship between two random variables. Although in principle this approach can work also with multi-dimensional data, theoretical as well computational difficulties often limit its use to the two-dimensional case. In this paper, we consider an approach, based on so-called vine copulas, which overcomes this limitation and at the same time is amenable to a theoretical treatment and feasible from the computational point of view. We apply this method to published data on the near-IR and far-IR luminosities and atomic and molecular masses of the Herschel Reference Sample. We determine the relationship among the luminosities and gas masses and show that the far-IR luminosity can be considered as the key parameter which relates all the other three galaxy properties. Once removed from the 4D relation, the residual relation among the other three is negligible. This may be interpreted as that the correlation between the gas masses and near-IR luminosity is driven by the far-IR luminosity, likely by the star-formation activity of the galaxy.