AbstractThe investigation of spatial variation in disease rates is a standard epidemiological practice used to describe the geographic clustering of diseases which is helpful for making hypotheses about the possible ‘factors’ responsible for differences in risk. Up to the most recent statistical and computational developments, studies have almost entirely focused on the spatial modeling of univariate distributions of cases, that is, on the spatial modeling of single diseases. However, many diseases show similar patterns of geographical variation which may suggest the existence of common underlying risk factors, whether these are related to the environment, to particular local food habits, or to the clustering of a particular population (genetic origin). In this work, for multivariate categorical data pointwise geo‐referenced in a ‘geostatistical’ fashion, we propose a model for the study of the joint spatial variation of more diseases. Our approach is based on a hierarchical (generalized linear mixed) multivariate model where the underlying latent structure is given by a Gaussian geostatistical spatial factor model. The methodology proposed can be seen as an extension of the geostatistical linear model of coregionalization, and of the related ‘factorial kriging analysis’, to the case of geo‐referenced, in general multi‐way, contingency tables. An application of the proposed methodology is shown on an epidemiological data set coming from an extensive survey on diabetes mellitus patients which involved the majority of the family practitioners of the region of Umbria in central Italy in 1990. Attention is centered on the study of nephropathy and retinopathy, two of the chronic diabetic complications affecting life quality and expectancy. Copyright © 2004 John Wiley & Sons, Ltd.