The same original data set as in [5] is used in this analysis. In total, 6,460 water samples from drinking water supplies from 99 Austrian districts were analysed for lithium. The average was 65.3 samples per district, with a range from 1 to 312 samples. The lowest measurable threshold lithium level by inductively coupled plasma optical emission spectrometry was 0.0033 mg/l. For the statistical calculations, lithium levels were averaged per district with a mean lithium level of 0.0113 mg/l (SD ± 0.027). Further, data on suicide rates adjusted as standardized mortality ratios (SMR; [9]), proportion of Roman Catholics, population density, and per capita income were obtained from Statistics Austria. Density of psychiatrists, general practitioners, and psychotherapists were obtained from the Austrian Medical Chamber and the Austrian Institute of Health (ÖBIG). Unemployment rates per district were obtained from the Austrian Public Employment Service (AMS) and were averaged for the available years 2005–2008. Further details were described previously in [5].

First, exploratory spatial data analysis is applied to demask spatial patterns which may contradict fundamental model assumptions of the subsequent regression. To start off, spatial autocorrelation and non-stationarity in the variables are explored using the global Moran’s *I* and the local *G**-statistic. In this context “global” refers to testing for spatial autocorrelation for the entire study area at once and deriving a single value indicating whether spatial autocorrelation exists and what its strength is.

The Moran’s *I*[10] is a widely used measure of global spatial autocorrelation, which tests whether there are some relationships between location and attribute values. A significant positive statistic indicates that nearby locations of similar attribute values are more spatially clustered than randomly distributed. In contrast, a significant negative statistic shows dissimilar values at nearby locations showing a more dispersed pattern.

This global measure is not capable to explore distinctive local features as well as the non-stationarity of a spatial process. Often, a non-stationary response variable is a first clue for spatially-varying relationships in a multivariate regression model. Therefore, Getis and Ord [11] have introduced the *G**-statistic to detect clusters of high or low values by location. High positive values refer to “hot spots” and high negative values to “cold spots”. As such, “hot spots” can be described as areas, where districts with high levels are surrounded by other districts with high levels. In contrast, “cold spots” are clusters with low level districts surrounded by other low level districts. A crucial step in spatial modeling is the choice of an appropriate representation of space. Common choices for such representations are, for instance, contiguity, *k*-nearest neighbors, among others. For a summary see [12].

The classical ordinary least squares (OLS) model is widely used to model the global relationship between a response and one or more explanatory variables. OLS assumes, among other things that residuals are spatially independent. Residual autocorrelation captures unexplained similarities between neighboring districts, which can be a results of omitted variables or a misspecification of the regression model in general [13]. Accounting for spatial effects reduces the magnitude of the prediction error, removes most of the systematic error, and leads to more reliable estimates [8]. In this study the OLS model serves as the reference model.

Global models that account for spatial effects are spatial autoregressive models, introduced by Anselin [6]. They comprise of two special cases, namely the spatial lag and the spatial error model. The spatial lag model extends the standard OLS regression model by including a spatially lagged dependent variable, which can be mostly interpreted as spill-over effects. Ignoring this lagged dependent variable leads the OLS model to be biased and inconsistent [14]. This model will be referred to as the spatial autoregressive (SAR) model in the reminder of this article. The spatial error model addresses the presence of spatial autocorrelation by defining a spatial autoregressive process for the error term and, by doing so, captures unexplained similarities. Not considering this spatial process in the error term yields the OLS estimates to be inefficient, although unbiased [14]. The Lagrange multiplier test statistic (LM) on the estimated OLS residuals helps to decide between these two alternative model specifications.

Global regression models assume a homogeneous behavior of the estimated parameters across space, which has often proven to be unrealistic. A way out is the application of the Geographically Weighted Regression (GWR), a local regression model that explores local spatial variations of the parameters [15]. The GWR models spatial autocorrelation and spatial heterogeneity for subsets of the entire data set. Each subset is established around a regression point with near data points exhibiting a higher influence than more distant data points. The weighting is often based on a bisquare kernel function (e.g. [16]). The setting of an appropriate bandwidth length of the weighting function is crucial. The most common is the adaptive bandwidth, where its length is allowed to vary across space, depending on the density of the data points. In densely populated areas the kernel possesses a shorter bandwidth in contrast to regions with larger inter-point distances, where the bandwidth is longer [15]. A kernel function with an adaptive bandwidth improves the goodness-of-fit, if data points are irregularly distributed across space. Despite the critique that the GWR is more suitable for exploratory analysis (see e.g. [17, 18] and references therein), it is a flexible model type to investigate spatially varying relationships, which is the focus of this research. GWR estimations are modeled for the same level of aggregation or spatial units as the original data and can be interpreted in the same way as OLS regression estimates.