One of the tenets of geostatistical modelling is that close things in space are more similar than distant things, a principle also known as “the first law of geography”. However, this may be questionable when unmeasured covariates affect, not only the mean of the underlying process, but also its covariance structure. In this paper we go beyond the assumption of stationarity and propose a novel modelling approach which we justify in the context of disease mapping. More specifically, our goal is to incorporate spatially referenced risk factors into the covariance function in order to model non-stationary patterns in the health outcome under investigation. Through a simulation study, we show that ignoring such non-stationary effects can lead to invalid inferences, yielding prediction intervals whose coverage is well below the nominal confidence level. We then illustrate two applications of the developed methodology for modelling anaemia in Ethiopia and Loa loa risk in West Africa. Our results indicate that the non-stationary models give a better fit than standard geostatistical models yielding a lower value for the Akaike information criterion. In the last section, we conclude by discussing further extensions of the new methods.