Estimating Range Parameter ($\rho$) for GAMs in 'mgcv' R Package

Question

A while back, I asked a question regarding the fitting of Gaussian Process (GP) smooths within a GAM framework that garnered some interest:

Gaussian Process smooths in mgcv: choosing between spherical and exponential covariance functions

A very recent article (published yesterday on BioRxiv) by @GavinSimpson expands on general ideas mentioned in the response to the above post through application to paleoecological and paleolimnological climate data.

In fitting GPs within the R package 'mgcv', it is best to specify the estimated range parameter ($\rho$), above which observations are no longer correlated with one another. 'mgcv' uses a default for $\rho$ equal to the maximum distance between pairs of observations, which may be adequate for some purposes (i.e., smooth, non-stochastic trends = non-time series).

The most straightforward method to go about selecting the optimal $\rho$ is to compare against some model selection criterion (AIC/UBRE/ML/REML) to see where a global minimum occurs.

My question is: how should one go about selecting a suitable interval of values to test for the range parameter? @GavinSimpson tested $\rho \in$ [15, 500] in his recent work, whereas Simon Wood uses 1:10*10 (10, 20, 30, ...) in his book (2nd ed. p. 362).

Any thoughts are greatly appreciated and welcomed.

Answer 1

$h$ is a distance in the units of the variable defining distance.

In my example (the paper you linked too) I was using data where the time ordering variable I was using was in units of years, so I search for separation distances of 15 years through 500 years (given that the example time series was 3000 years long, 500 years upper limit seemed reasonable)**.

In Simon's example, the spatial distance between samples is in KM (if I'm following the entire example properly) hence 10, 20, etc km separation relative to distances on the whole of Portugal see reasonable values to search over.

In summary, you need to use domain knowledge and the limits of the data to decide what values to search over. You need to have observations closer than 10 km for example if you are going to fit a model with basis functions based on a correlation function with an effective range parameter $\phi$ set to 10 km.

** In hindsight, that lower limit was likely too low for a chunk of the series and is perhaps what is causing the odd behaviour in the trace of the REML score at low values of $h$ — many of the observations earlier in the series are separated by more than 15 years.