GAM on values per cell in a grid cell

Question

I am running a GAM model on grid-cell data at a 1Km spatial resolution ~89K cells. My response variable is species richness, and I use a tw() distribution as it works better than Poisson and nb. The predictor variables are 2 Dimensions of a PCA analysis I conducted on bird abundance data. I am modelling spatial autocorrelation by including a smooth of lat long. Still, the plots show problems in the model and significance in the number of k values used for the spatial smooth.

This is the model I am running:

mod2 <- gam( epm_ON.x ~ s(x, y,  bs = "ts", k = 100) + s(Dim.1, Dim.2), data = data_model_ON_reduced, family = tw(), method = "REML")

I am thinking of running the spatial autocorrelation with the bs = 'mrf", but I am unsure what else I might be missing.

Does anyone have a recommendation regarding using the PCA axes?

Answer 1

I wouldn't use Tweedie for this given what is shown. Species richness should be an integer value and the output is showing that you are effectively fitting a Poisson GAM (the $p = 1.01). I've rarely had good diagnostic plots when modelling species richness data with a Poisson model. Instead try the nb() family for the negative binomial; that should help some with the extra-Poisson dispersion evident in the QQ plot.

I assume the response is not species richness of those birds that you computed the PCA on?

In general, I'm not a fan of using PC axes like this in models, mainly because you're going to run into problems explaining what the smooth of the two axes means in terms of the original data (bird abundances). I would also suggest that this smooth should be a tensor product, not a 2d thin plate spline, as the latter assumes the same wiggliness in all directions, which need not be the case with PCs.

As for the $k$ issue; the heuristic basis of the test starts to run into problems when your data are in the form of a time series or spatial layout. Even small amounts of autocorrelation left over in the residuals will cause the test to reject the null hypothesis (in my experience). Simon Wood, in his GAM book (2nd ed.) suggests an alternative:

1. take the deviance residuals from your model,

2. add them as the response in a GAM with the same covariate terms, the following family

   family = quasi(link = "identity", variance = "constant")
   

   and double the k on the smooth terms,

3. look at the model summary and see if the EDF for the smooth terms in greater than 1 (the minimum),

4. Go back to your original model and double k for any term that had a greater than 1 EDF in the model for the deviance residuals,

5. rinse and repeat if necessary

I probably wouldn't use the the shrinkage smooth for the spatial effect as my starting point, get the main features figured out first, then switch on extra penalties if you need them.

The "mrf" basis will also work for this kind of spatial grid layout, but you'll either run into computational problems (if you run the full MRF you'll need a 89,000 x 89,000 penalty matrix) or you'll need to set k lower to get the low-rank version. But give it a go if you want to try it. As the MRF has a different definition of wiggliness, it might work better in the spatial setting, but it will only do so by being very wiggly. Whether a smoother fit or a wigglier fit is desirable isn't clear from your question.