Geographic information systems (GISs) offer a powerful tool to geographers, foresters, statisticians, public health officials, and other users of spatially referenced regional data sets. However, as useful as they are for data display and trend detection, they typically feature little ability for statistical inference, leaving the user in doubt as to the significance of the various patterns and 'hot spots' identified. Unfortunately, classical statistical methods are often ill suited for this complex inferential task, dealing as it does with data which are multivariate, multilevel, misaligned, and often nonrandomly missing. In this paper we describe a Bayesian approach to this inference problem which simultaneously allows interpolation of missing values, estimation of the effect of relevant covariates, and spatial smoothing of underlying causal patterns. Implemented via Markov-chain Monte Carlo (MCMC) computational methods, the approach automatically produces both point and interval estimates which account for all sources of uncertainty in the data. After describing the approach in the context of a simple, idealized example, we illustrate it with a data set on leukemia rates and potential geographic risk factors in Tompkins County, New York, summarizing our results with numerous maps created by using the popular GIS Arc/INFO.