the topic of spatial autocorrelation (SA) within the context of generalized additive models has already been discussed in several posts within this forum, see e.g.
Why does including latitude and longitude in a GAM account for spatial autocorrelation?
Modelling spatiotemporal data with GAMs
A point which however continues to confuse me is whether GAMs "perform well" for data sets, which display a high degree of SA. This confusion arises because several people have advised me against the use of GAMs for this purpose and indeed in papers such as e.g. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review, you can find lines such as:
Trend surface GAM (Hastie and Tibshirani 1990, Wood2006) does not address the problem of spatial autocorrelation, but merely accounts for trends in the data across larger geographical distances".
The background for my question is that I am currently working on fitting a valuation model based a on data set of spatially autocorrelated house prices (which I have already described in other posts such as Model checking for generalized additive models). My model is based on fitting a GAM to the data set, but the points above make me doubt if this is the best strategy to pursue or whether I should at least think about coupling it with an approach that more directly handles the SA. Below I describe in detail my current model and results and at the end of the post I attempt to summarize some of the things that still puzzle me.
Current model
The model I currently use is obtained by fitting a GAM to training set of 200.000 house sales, and I subsequently evaluate its performance on an independent validation set of 41253 house sales. The response variable and covariates are listed below:
Response variable: sales_price
Covariates:
Geographical region of house: region
Coordinates of house: (coordinate_east, coordinate_north)
Sales date (converted to a numeric): sales_date_numeric
Building year: year_built
Living space: living_space
Lot size: lot_size
Distance to coast: coast_distance
Distance to lake: lake_distance
Distance to train station: train_station_distance
Distance to road: road_distance
Roof type (categorial variable): roof_type
Wall type (categorial variable): wall_type
Heating type (categorical variable): heating_type
The current model structure I use is the following:
model <- mgcv::bam(sales_price ~ s(coordinate_east, coordinate_north, k = 4000) +
s(sales_date_numeric, k = 10, pc = 0, by = region) +
s(year_built, k = 8, pc = 1970, by = region) +
s(living_space, k = 5, pc = 140, by = region) +
s(lot_size, k = 5, pc = 800, by = region) +
s(coast_distance, k = 5, pc = 1500, by = region) +
s(lake_distance, k = 5, pc = 500, by = region) +
s(road_distance, k = 5, pc = 2000, by = region) +
s(train_station_distance, k = 5, pc = 2000, by = region) +
roof_type*region +
wall_type*region +
heating_type*region,
data = training_set,
family = tw(link=log),
method = "fREML")
Modelling the spatial autocorrelation
Does the model structure above capture the SA in my dataset sufficiently well? This is still somewhat an open question to me, but what I have done to adress this so far is to examine the performance of the model above for different number of basis functions in the spatial smooth s(coordinate_east, coordinate_north, k). The result of this is shown in the plot below which plots the validation RMSE as a function of the number of basis functions k:
Evidently the validation error drops substantially with increasing number of basis functions and continues to do so at large numbers, where fitting the model starts to become too computationally intense. In my previous post Model checking for generalized additive models Gavin Simpson commented that my choice of k=4000 might be excessive, but if my results in the plot above are valid, it is indeed necessary to employ such a large number of basis functions.
Comparison with other approaches for modelling spatial autocorrelation
My own interpretation of the plot above is that my GAM-model is able to capture the spatial autocorrelation in my data set but that this comes at the cost of having to fit a very large and computationally demanding model. It is therefore natural to ask whether it would make more sense to employ a more local fitting approach.. I am not an expert of these, but the simplest I could think of is a KNN-model which predicts the house price as:
$$ P_{i}^{KNN} = \frac{1}{N}\sum_{j=1}^{N}P_{j}\tag 1\label 1 $$
where the sum runs over the $N$ nearest neighbors of each property. In itself $\eqref 1$ will be less accurate than the GAM-model, as it does not directly include information on living space, lot size like the GAM-model does. However, denoting the price estimates of my GAM-model by $P_{i}^{GAM}$ I can form as combined KNN-GAM model estimate as:
$$ P_{i}^{\text{Combined}} = P_{j}^{KNN} + (P_{i}^{GAM} - \frac{1}{N}\sum_{j=1}^{N}P_{j}^{GAM}) \tag 2\label 2 $$
The interpretation of $\eqref 2$ is that one first estimates the local price level based on KNN and then corrects for the difference of the property with those in its neighborhood using the GAM-model. I have plotted the results from this combined model in the plot below. Evidently this vastly improves the convergence with respect to the number of basis functions.
Summary and questions
This post has adressed the modelling of spatial autocorrelation with GAMs. Concretely it has focused on modelling a data set of spatially autocorrelated house prices with a gam model and shown that the accuracy of this could be improved by combining it with a local smoothing approach. Below are some of the questions that I am still struggling with and where I could use an expert's point of view:
Q1: Do the results from my convergence analysis with respect to the number of basis function seem reasonable? In another post Model checking for generalized additive models where I posted results from model checking my GAM-model it became clear that my model has some issues with its distributional assumptions, which I still have not managed to fix. My concern is therefore that the extremely larged number of basis functions needed somehow relates to the fact that the model is ill-defined to start with and that a ton of degrees of freedom are needed to correct for this.
Q2: Combining the GAM-model with a KNN model seemed to vastly improve the convergence results. What would be the downside of such an approach? And is there an approach with GAMs that can achieve something similar - i.e. a GAM model that more directly incorporates SA.
