Though in the last decade many works have appeared in the literature dealing with model-based extensions of the classical (univariate) geostatistical mapping methodology based on linear Kriging, very few authors have concentrated, mainly for the inferential problems they pose, on model-based extensions of classical multivariate geostatistical techniques like the linear model of coregionalization, or the related ‘factorial kriging analysis’. Nevertheless, in presence of multivariate spatial non-Gaussian data, in particular count data, as in many environmental applications, the use of these classical techniques can lead to incorrect predictions about the underling factors. To overcome this problem, here we discuss a hierarchical geostatistical factor model that extends, following a model-based geostatistical approach, the classical geostatistical proportional covariance model. For this model we investigate a likelihood-based inferential procedure using the Monte Carlo EM algorithm. In particular, we discuss some of its theoretical properties and show, through some thorough simulation studies, its sampling performances.