Model checking for generalized additive models

Question

I am currently constructing a GAM model to describe house prices. The dataset is a collection of roughly 200K house sales in a particular geographic region. Based on previous suggestions in this forum, the mean house price $p_{i}$ is modelled as:

$$ \ln (p_{i}) = \sum_{j}s_{j}(x_{i,j}) \ \ , \ \ P_{i} \sim Tweedie(p_{i}) $$

, where $s_{j}$ denotes the jth smooth for my set of covariates, which are geograpgic location, house size, lot size etc. etc. Due to the strict positivity of prices, I have found that the combination of a log-link and a Tweedie-distribution provides best results measured by out-of-sample error (however, a Gamma-distribution with log-link also does quite well - see my previous post Simulating from fitted t-distribution mgcv for details on this).

In search for the best model I have tried a range of different covariates and arrived at a model, that I believe to be the optimal in terms of balance between predictive capability and parsimony. I am however concerned about the results I get when applying the model checking routine offered by mgcv (or appraise() in the gratia library) as shown below:

Evidently the results based on the QQ-plot suggest, that my model are not able to capture the fat tails in my distribution. Also the residual vs linear predictor plot is not a flat band as one would hope but has some systematic linear component in the first part of the negative residuals - I suspect that this has something to do with the fact that the model can only output positive values.

As I do have limited experience with model checking of GAMs, I therefore want to ask the following question to those that have more experience in the field:

Do the above model diagnostics indicate a serious problem in my model, or are they within the tolerance of what one generally would accept? And if so, what are the general rules of thumb that people use when assesing model quality based on QQ-plots etc.

Addition to original question based on answer from Gavin Simpson:

As Gavin correctly points out the above question is not sufficiently detailed for one to be able to provide concrete suggestions. While I am not allowed to share the data I am working with I below add details on the exact model used to obtain diagnostics described in the original question.

The model I used is based on fitting a data set of 222299 house sales with the following information:

1. Response variable: House price per square meter, P

2. Covariates:

Coordinates, (x,y)
Sales date, t (converted to an integer)
Building size, A
Lot size, L
Construction year, Y
Distance to coast, D1
Distance to lake, D2
Distance to railway station, D3
Distance to highway, D4
Roof type (categorial variable), R
Wall type (categorial variable), W
Heating source (categorical variable), H

The distance variables are special. These are either in the interval [0, 1000] or 10000 (a large number to indicate they are far away from the entity in question).

The model I specify in mgcv has the following form:

model <- mgcv::bam(P ~ s(x, y, k = 4000) +
                       s(t, k = 10, pc = 0) +
                       s(Y, k = 8, pc = 1970) +
                       s(A, k = 5, pc = 140) +                                   
                       s(L, k = 5, pc = 800) +
                       s(D1, k = 5, pc = 10000) +                   
                       s(D2, k = 5, pc = 10000) +
                       s(D3, k = 5, pc = 10000) +
                       s(D4, k = 5, pc = 10000) +
                       R +
                       W +
                       H,
                       data = data,
                       family = tw(link=log),
                       method = "fREML")

Answer 1

The problem with large data sets is that heterogeneities often become visible. With a location-only model, you are fitting a model that also includes constant "effects" for $p$ and $\varphi$, the power and scale parameters of the Tweedie distribution. As a result, if you are going to fit a location-only model, the model can only fit the data with a constant mean-variance relationship. This is often inflexible to sufficiently model the heterogeneities in large data sets.

There are clear problems with this model; the data are more dispersed than expected under a Tweedie with $\hat{p} = 1.242$ and $\hat{\varphi}$ constant.

This could be modelled with extra terms in the location linear predictor; as you don't actually specify what the model that you fitted it is (why don't you ever post this useful information?), it's impossible to give concrete advice. If the data are spatially distributed, adding a spatial smooth might help. If the data are clustered in some way to neighbourhoods, then a random effect or Markov random field smooth could be used to account for effects related to desirable vs less-desirable neighbourhoods, the latter accounting for the spatial relationship too.

Beyond that, if the diagnostics remain bad, you'll need to move to a more complex model and start modelling the $p$ and $\varphi$ parameters with linear predictors. For this you'll need to fit with gam() and use the twlss() family (or you could also try the gammals() family if you prefer to fit a Gamma location-scale model).

As there is less information in the data about the power and scale parameters (and thence variance and other higher moments) than there is for the mean (location), you should start the model off with simple effects in the linear predictors for $p$ and $\varphi$ (if fitting with twlss()) or $\varphi$ (if fitting with gammals()). For example, include any random effects or Markov random field smooths to capture potentially gross differences due to (spatial)clustering in the data - i.e. the neighbourhood effects.

The mgcViz package has additional diagnostic plots that can be used to help diagnose the need for a location scale shape (aka "distributional model").