2. Space-time GAM contruction and considertations
2.1 Overview
This section works through an example of constructing a space-time GAM to illustrate a workflow that allows for spatio-temporally varying terms, using the Chaco dataset. It fits an initial complex model and then, through investigation of the results, refines the how the model is specified. In so doing a number of important considerations are illustrated relating the bases that are specified, their dimensions, the GAM fitting method, the model family and the number of knots used in the smooths.
In the GAMs created in this vignette and this package, the Intercept is treated as an explicitly addressable term.
By way of example, consider the model m created
below:
m <- gam(ndvi~s(X,Y) + s(month) + pr, data = chaco)
summary(m)
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> ndvi ~ s(X, Y) + s(month) + pr
#>
#> Parametric coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 5.595e-01 3.391e-03 165.01 <2e-16 ***
#> pr 5.366e-04 1.948e-05 27.54 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df F p-value
#> s(X,Y) 21.648 26.055 5.444 <2e-16 ***
#> s(month) 8.925 8.998 13.355 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.335 Deviance explained = 34.5%
#> GCV = 0.011302 Scale est. = 0.011117 n = 2000This includes the terms s(X,Y) + s(month). These model
the spatially varying and temporally varying Intercept plus error. Note
that if this was a spatial model (without time), then just
s(X,Y) would be included for the Intercept and similarly
just s(month) for a purely temporal model.
The model can be reformulated as follows with an addressable Intercept term:
m <- gam(ndvi ~
0 + Intercept + s(X,Y,by=Intercept) + s(month, by=Intercept) +
pr,
data = chaco |> mutate(Intercept = 1))
summary(m)
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> ndvi ~ 0 + Intercept + s(X, Y, by = Intercept) + s(month, by = Intercept) +
#> pr
#>
#> Parametric coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> Intercept 5.595e-01 3.391e-03 165.01 <2e-16 ***
#> pr 5.366e-04 1.948e-05 27.54 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df F p-value
#> s(X,Y):Intercept 21.648 26.055 5.444 <2e-16 ***
#> s(month):Intercept 8.925 8.998 13.355 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.335 Deviance explained = 97.3%
#> GCV = 0.011302 Scale est. = 0.011117 n = 2000The model summary now includes Intercept as a parametric
term (rather than (Intercept)) with
s(X,Y):Intercept and s(month):Intercept as
smooth terms, rather than s(X,Y)and s(month),
but the values in both model summaries are all the same.
The rationale for treating the Intercept as addressable term is because, in varying coefficient modelling, specifying a varying intercept over space and or time detects situations when the response variable (the \(y\)-value) may show notably low or high values in a space-time locality that were not accounted for by variation in the explanatory values (Harris 2019). This could be caused, for example, if variables not recorded in the study were causing otherwise unexplained space-time patterns in the \(y\) variable. As such, a varying intercept commonly tends to capture autocorrelated errors (Harris 2019) and background variation that has the potential to affect other varying coefficient estimates. It also simplifies the interpretation of other varying coefficient estimates, particularly if the predictor variables are not mean-centred.
The models in this vignette include the Intercept as an addressable term and the code below creates this variable in the input data:
2.2 Space-time GAMs
An initial space-time GAM is fitted using the code below. Notice that
for each predictor variable and the explicitly specified intercept, the
model specifies Tensor Product (TP) smooths to capture the variations in
NDVI (ndvi) in both space and time simultaneously. These
are tensor product smooths, specified using the te()
function, and are described below. Notice also that global
parametric parameters are also specified for tmax,
pr and Intercept. This is to provide a
multiplicative constant as well as a function around the coefficient
term, and potentially an offset term.
# initial model
m1 <- gam(
ndvi ~ 0 +
Intercept + te(X,Y,month, by=Intercept,bs = c("tp", "cr"), d = c(2,1)) +
tmax + te(X,Y,month, by = tmax, bs = c("tp", "cr"), d = c(2,1)) +
pr + te(X,Y,month, by = pr, bs = c("tp", "cr"), d = c(2,1)),
data = chaco,
method = "REML",
family = gaussian()
)Each of the space-time smooths is specified with Tensor Product (TP)
smooth (te()) to model interactions between space and time
margins. These construct smooth surfaces over space and time without
random effects (Wood 2025). They enable
variables that are measured over different scales (i.e. that are
anisotropic), such as location coordinates and time, to be combined. The
basis dimensions of each component are specified through the
d parameter, here specifying a 2-dimensional smooth basis
for location (X and Y) and a 1-dimensional one
for time (month). These are combined via the TPs. Thin
plate regression splines (tp) were specified to model
location and low-rank cubic regression splines (cr) were
specified for time. Thin plate regression splines are ideal for 2-D
smooths because they are isotropic (they are rotation invariant
with have no directional bias) and they perform well even if the
observations are irregularly located over space. Cubic regression
splines are efficient in 1-D and known to be good for seasonal or
long-term effects. The low-rank basis avoids overfitting and improves
computational efficiency. Finally, the model is specified with a
Gaussian family and the model is fitted using a REML (Relaxed Maximum
Likelihood) method. Model estimation in GAMs can be conducted in two
ways: (a) predictive (i.e., out-of-sample minimization) via Generalized
Cross-Validation (GCV) and (b) Bayesian (i.e. attach priors on basis
coefficients) via some maximum likelihood approach. These tend to
provide more stable estimates of the smoothing parameters (i.e., reduce
over-fitting) and tend to provide more robust estimates of the variance.
GCV is more computationally efficient but can over-fit, especially when
errors are correlated. Note that GCV is set as the default in
mgcv, but REML based approaches are preferred because they
provide less biased estimates of variance and smoothing parameters, and
are more robust during the optimization process.
The model can be examined using the appraise function in
the gratia package (Simpson 2024;
Simpson 2025) which generates diagnostic plots of the model
residuals (the difference between the observed value of \(y\) (in this case ndvi) and
the predicted value of \(y\), also
referred to as \(\hat{y}\) :
Diagnostic plots of m1.
The diagnostic plots show that the top and bottom and of the QQ plot
are close to the theoretical fit line and that the modelled quantiles
are close to the theoretical one. This indicates that there are a low
number of residuals - observed quantiles that are not on the theoretical
fit line. If there were a high number - i.e. the QQ plot had an S-shape
- then this would suggest residual heterogeneity, potentially due to
un-modelled interactions or due to an inappropriate assumption of a
Gaussian model. In such cases other response should be investigated for
example through histograms and specified accordingly via the
family parameter that is passed to the GAM - see below.
2.2 Space-time GAM Refinements I: Model family
The assumptions behind a Gaussian model are that the target variable is normally distributed, meaning that it is continuous, has constant variance and is unbounded (i.e. it can theoretically take any value from \(-\propto\) to \(+\propto\)).
If the target variable has a non-normal distribution (a long tail, a skewed distribution, etc) then a different model family should be specified. For example, a target variable such as house price which has values that cannot be negative and may have some very high values, would suggest specifying a a Gamma family with a log link. This is used when the response is continuous, strictly positive, right-skewed (i.e., with a long tail to the right) and when the variance increases with the mean (i.e. is heteroscedastic). What this effectively does to provide a log-scale regression for positive, skewed data.
The point is that there are a number options for specifying the model
family in a mgcv GAM, depending on the
distribution of the response variable, as summarised in the table below.
Further information can also be found in the help for the
family function (enter ?family).
| Family | Response | Link | Example |
|---|---|---|---|
| Gaussian | continuous, symmetric | identity | residuals roughly normal |
| Gamma | continuous, positive, skewed | log | costs, durations, prices |
| Inverse Gaussian | continuous, positive, heavy-tailed | log | reaction times, dispersions |
| Poisson | count data | log | integer counts |
| Binomial | 0–1 or proportion data | logit | success/failure outcomes |
This vignette and the current version of the stgam
package considers only a Gaussian, normally distributed response. Future
versions will address other responses.
2.3 Space-time GAM Refinement II: Effect size
Having determined that correct model family has been specified for the response variable, the model summary can be examined:
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> ndvi ~ 0 + Intercept + te(X, Y, month, by = Intercept, bs = c("tp",
#> "cr"), d = c(2, 1)) + tmax + te(X, Y, month, by = tmax, bs = c("tp",
#> "cr"), d = c(2, 1)) + pr + te(X, Y, month, by = pr, bs = c("tp",
#> "cr"), d = c(2, 1))
#>
#> Parametric coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> Intercept 1.562e+00 4.069e-02 38.385 < 2e-16 ***
#> tmax -1.821e-02 5.604e-03 -3.249 0.00118 **
#> pr 5.020e-04 4.887e-05 10.273 < 2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df F p-value
#> te(X,Y,month):Intercept 5.01 5.014 6.435 6.35e-06 ***
#> te(X,Y,month):tmax 34.83 42.522 4.347 < 2e-16 ***
#> te(X,Y,month):pr 10.48 12.178 1.686 0.0591 .
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Rank: 375/377
#> R-sq.(adj) = 0.478 Deviance explained = 49.1%
#> -REML = -1758.3 Scale est. = 0.0087294 n = 2000This indicates that pr specified in a space-time smooth
is not significant but that all of the other model components are.
To unpick this further and to provide a further basis for model
evaluation, it is also instructive to examine the size of the effect of
each model component. The code below generates numeric summaries of the
effects of each smooth by extracting the effect of each term from the
model using the predict function. The typical size of each
smooth’s contribution to predictions are summarised from their standard
deviations (sd) and the total amplitude of the smooth
effects (range), where higher values indicate a larger
effect on the target variable. This is undertaken below and the code
snippets have been wrapped into the effect_size function in
stgam.
# extract the terms
f1 <- predict(m1, type = "terms")
# generate numeric summaries
apply(f1, 2, function(v) c(sd = sd(v), range = diff(range(v)))) |> round(3)#> Intercept tmax pr te(X,Y,month):Intercept te(X,Y,month):tmax te(X,Y,month):pr
#> sd 0 0.033 0.062 0.239 0.234 0.011
#> range 0 0.242 0.282 1.639 1.518 0.121This shows the size of the effect of each model component.
Considering the ranges, we can see the zero effect of the parametric
Intercept and pr term and the low partial
effect of the TP smooth for the pr term, suggesting a lack
of interaction over space-time for this this predictor variable. A
possible refinement in this case could be to remove the fixed
Intercept and pr terms but we would like to
retain them as constants and offsets. The range value of the combined
space-time smooth for pr is relatively low
(te(X,Y,month):pr) and it may be useful to examine
pr over space and time separately. In the code below the
individual space and time smooths are specified as thin plate regression
smooths, in s() function in the mgcv
package.
# second model with separate space-time effects for pr
m2 <- gam(
ndvi ~ 0 +
Intercept + te(X,Y,month, by=Intercept, bs = c("tp", "cr"), d = c(2,1)) +
tmax + te(X,Y,month, by = tmax, bs = c("tp", "cr"), d = c(2,1)) +
pr + s(X,Y, by = pr) + s(month, by = pr),
data = chaco,
method = "REML",
family = gaussian()
)Then effect size can be examined using the effect_size
function:
#> Intercept tmax pr te(X,Y,month):Intercept te(X,Y,month):tmax s(X,Y):pr s(month):pr
#> sd 0 0.037 0 0.219 0.214 0.046 0.107
#> range 0 0.271 0 1.437 1.356 0.209 0.437Now the size of the effect of pr is split and is larger
over both space (tmax s(X,Y):pr) and time
(s(month):pr), compared to the m1 GAM model.
This demonstrate how evaluating and comparing models with different
space-time smooths can help to determine whether to model them
interactively or not, and what the changes in the magnitude of the
modelled response surface might be.
2.4 Space-time GAM Refinement III: k the number of
knots
In the mgcv GAM implementation the number of knots in
each smooth are automatically determined. These control how ‘wiggly’
(non-linear) the the smooth can be. The number of knots in each smooth
can be examined using the gam.check function. This provides
some diagnostic information summarising the model fitting procedures can
also be generated and also prints out some diagnostic plots (although
not as informative as the ones created by the appraise
function in the gratia R package).
This checks whether the number of knots used in the smooths are too
low. This is indicated when the effective degrees of freedom (EDF in the
edf column) are close to the number of knots
(k) and some of the k-index value are somewhat
less than 1. In this case, the number of knots used in the temporal
smooth for pr may be too low.
The number of knots in a used in a smooth are automatically
determined but can be specified. In mgcv GAM smooths,
k controls the basis dimension of the smooth term —
i.e. the maximum possible complexity of the smooth, not the actual
complexity. The effective degrees of freedom (edf) after
fitting should be much smaller. If k is too small, the
model fails to capture the true pattern of the relationships in the data
(i.e. it fails to capture non-linearities) and can result in biased fits
and residual structure. There are some approximate guides for specifying
k that are summarised in the table below.
| Smooth type | Heuristic | Notes |
|---|---|---|
| 1D temporal smooth (e.g., days, years) |
k ≈ min(n/10, 20–40)
|
Typically start around 10 basis functions; if long time series (> 1000 points), may go to 50–100 |
| 2D spatial smooth (e.g., te(X, Y)) |
k ≈ 50–200 total basis functions
|
Typically start around 50 basis functions; choose higher if data cover large/complex region |
| 3D space–time smooth (e.g., te(X, Y, time)) |
k = product of marginal bases,
e.g. k = c(50, 20)
|
mgcv handles tensor product decomposition; but do not let any single dimension be too small |
| Factor-by smooths (e.g., by = region) |
k ≈ same as main smooth per level
|
You can use fewer if levels are many and data per level small |
As a general rule of thumb, the k value should be
doubled if the check reveals it to be too low. This is done for the
temporal smooth specified for pr in the code below to
create m3 and note that as k increases so does
the time to fit the model:
# third model with increased k for the pr temporal smooth
m3 <- gam(
ndvi ~ 0 +
Intercept + te(X,Y,month, by=Intercept, bs = c("tp", "cr"), d = c(2,1)) +
tmax + te(X,Y,month, by = tmax, bs = c("tp", "cr"), d = c(2,1)) +
pr + s(X,Y, by = pr) + s(month, k=20, by = pr),
data = chaco,
method = "REML",
family = gaussian()
)The gam.check function calls the k.check
function, provides some interpretation notes and diagnostic plots.
However, the k.check function can be called directly:
#> k' edf k-index p-value
#> te(X,Y,month):Intercept 124 11.297948 0.8969713 0
#> te(X,Y,month):tmax 125 34.676559 0.8971304 0
#> s(X,Y):pr 30 7.268772 0.7575468 0
#> s(month):pr 20 14.121735 0.5536068 0This shows some improvement, with the EDF now increased but still
close to k for the smooth and the k-index
still much lower than 1. The model below tries to improve the fit by
increasing k once more:
# fourth model with increased k for the pr temporal smooth
m4 <- gam(
ndvi ~ 0 +
Intercept + te(X,Y,month, by=Intercept, bs = c("tp", "cr"), d = c(2,1)) +
tmax + te(X,Y,month, by = tmax, bs = c("tp", "cr"), d = c(2,1)) +
pr + s(X,Y, by = pr) + s(month, k=40, by = pr),
data = chaco,
method = "REML",
family = gaussian()
)
# check the model smooths
k.check(m4)#> k' edf k-index p-value
#> te(X,Y,month):Intercept 124 9.442720 0.9193043 0
#> te(X,Y,month):tmax 125 36.175682 0.9107176 0
#> s(X,Y):pr 30 7.941069 0.7444885 0
#> s(month):pr 40 36.853976 0.6355069 0Again some improvement is shown but this can be nudged further:
# fifth model with increased k for the pr temporal smooth
m5 <- gam(
ndvi ~ 0 +
Intercept + te(X,Y,month, by=Intercept, bs = c("tp", "cr"), d = c(2,1)) +
tmax + te(X,Y,month, by = tmax, bs = c("tp", "cr"), d = c(2,1)) +
pr + s(X,Y, by = pr) + s(month, k=80, by = pr),
data = chaco,
method = "REML",
family = gaussian()
)
# check the model smooths
k.check(m5)#> k' edf k-index p-value
#> te(X,Y,month):Intercept 124 15.71355 0.9177279 0
#> te(X,Y,month):tmax 125 37.46644 0.9162767 0
#> s(X,Y):pr 30 7.65725 0.7242695 0
#> s(month):pr 80 68.65270 0.8199559 0Some further improvement is evident but this process of increasing
k should continue. However, the number of knots cannot be
increased infinitely due the data observations, the number of terms in
the model and how they combine to determine the effective degrees of
freedom (edf). The code below increases the number of knots
by 40 as in this case there are only 120 unique values for months.
# sixth model with increased k for the pr temporal smooth
m6 <- gam(
ndvi ~ 0 +
Intercept + te(X,Y,month, by=Intercept, bs = c("tp", "cr"), d = c(2,1)) +
tmax + te(X,Y,month, by = tmax, bs = c("tp", "cr"), d = c(2,1)) +
pr + s(X,Y, by = pr) + s(month, k=120, by = pr),
data = chaco,
method = "REML",
family = gaussian()
)
# check the model smooths
k.check(m6)#> k' edf k-index p-value
#> te(X,Y,month):Intercept 124 15.047285 0.9195204 0
#> te(X,Y,month):tmax 125 37.547213 0.9281455 0
#> s(X,Y):pr 30 7.533317 0.7161696 0
#> s(month):pr 120 80.443605 0.8443115 0In this case the check shows that there is sufficiently a large
difference between k (120) and edf (~80) and
the k-index is closer to 1 (0.84). The final model can be
examined the usual way:
The effects of the terms can be examined as before:
#> Intercept tmax pr te(X,Y,month):Intercept te(X,Y,month):tmax s(X,Y):pr s(month):pr
#> sd 0 0.023 0 0.253 0.249 0.060 0.134
#> range 0 0.169 0 1.326 1.381 0.305 0.663The refined model above, m6, shows improvement, with the
effect of pr over time relatively high
(s(month):pr). Its effect over over space
(s(X,Y):pr) is relatively low and potentially the spatial
smooth for pr could be dropped. Some confirmation for this
decision is provided when the module summary is examined:
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> ndvi ~ 0 + Intercept + te(X, Y, month, by = Intercept, bs = c("tp",
#> "cr"), d = c(2, 1)) + tmax + te(X, Y, month, by = tmax, bs = c("tp",
#> "cr"), d = c(2, 1)) + pr + s(X, Y, by = pr) + s(month, k = 120,
#> by = pr)
#>
#> Parametric coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> Intercept 1.238686 0.059797 20.715 <2e-16 ***
#> tmax -0.012671 0.005584 -2.269 0.0234 *
#> pr 0.000000 0.000000 NaN NaN
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df F p-value
#> te(X,Y,month):Intercept 15.047 20.21 2.287 0.000885 ***
#> te(X,Y,month):tmax 37.547 43.99 4.380 < 2e-16 ***
#> s(X,Y):pr 7.533 10.47 1.230 0.288856
#> s(month):pr 80.444 86.84 12.649 < 2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Rank: 399/402
#> R-sq.(adj) = 0.67 Deviance explained = 69.3%
#> -REML = -1994.6 Scale est. = 0.0055213 n = 2000However in this case, the m6 model is retained as the
final model as the effect sizes are relatively strong for both
components. It specifies a space-time TP smooth for the Intercept and
for tmax, and separate spatial and temporal smooth for
pr. All of the parametric terms are retained:
The size of the effects can be examined and compared with the
original model m1:
#> Intercept tmax pr te(X,Y,month):Intercept te(X,Y,month):tmax te(X,Y,month):pr
#> sd 0 0.033 0.062 0.239 0.234 0.011
#> range 0 0.242 0.282 1.639 1.518 0.121#> Intercept tmax pr te(X,Y,month):Intercept te(X,Y,month):tmax s(X,Y):pr s(month):pr
#> sd 0 0.023 0 0.253 0.249 0.060 0.134
#> range 0 0.169 0 1.326 1.381 0.305 0.663The ranges of the final model (m6) still indicate the
strong space-time interaction of NDVI with the tmax in the
TP smooth (te(X,Y,month):tmax) and the interactions of
pr in separate space and time smooths
(s(X,Y):pr and s(month):pr)).
The AICs of the 2 models can be compared and the final model shows a dramatic improvement (i.e. reduction in value) from the original:
#> [1] -3751.768
#> [1] -4583.2882.5. Extracting, summarising and plotting the coefficients
The spatially and temporally varying coefficients estimates generated
by the model can be extracted in a number of different ways using the
calulate_vcs function in the stgam package.
There are basically 2 options for generating the varying coefficient
estimates: i) over the original data points or ii) over spatial surfaces
for specific time slices.
For the first option, the original data, the GAM and a vector of the
model terms (i.e. predictor variables names) are simply passed to the
calculate_vcs function. The second option requires slices
of the space-time data to be passed to the calculate_vcs
function and a bit more work to process the results. In both cases the
results can be summarised and mapped.
Considering first the original observations:
A summary over space and time of the coefficient estimates shows that
the Intercept is large and varying and that tmax and
pr also have varying relationships with the target variable
over space and time, with negative and positive values in some places
and / or at some times. Their values reflect the association of degrees
Centigrade and millimetres of rainfall. Considering the mean values,
each additional 1 degree of tmax is associated with a
0.0222 decrease in NDVI (biomass), and each additional millimetre of
rainfall is associated with a 0.0013 increase in NDVI.
vcs |>
# drop the geometry if the object is spatial
st_drop_geometry() |>
select(starts_with("b_")) |>
apply(2, summary) |>
round(4)
#> b_Intercept b_tmax b_pr
#> Min. 0.4448 -0.0436 -0.0077
#> 1st Qu. 1.0610 -0.0284 0.0004
#> Median 1.2679 -0.0229 0.0007
#> Mean 1.2387 -0.0222 0.0013
#> 3rd Qu. 1.4349 -0.0165 0.0020
#> Max. 1.7712 0.0023 0.0082These can be explicitly examined over time and notice the increasingly positive relationship with NDVI over time with the Intercept and the increasingly negative relationship of Maximum temperature. There is much more variation in the temporal trend of the relationship between NDVI and Precipitation.
vcs |>
st_drop_geometry() |>
select(date, starts_with("b_")) |>
rename(`Intercept` = b_Intercept,
`Max Temperature` = b_tmax,
`Precipitation` = b_pr) |>
pivot_longer(-date) |>
group_by(date, name) |>
summarise(
lower = quantile(value, 0.25),
median = median(value),
upper = quantile(value, 0.75)
) |>
ggplot(aes(x = date, y = median)) +
geom_point(col = "blue", alpha = 0.2) +
geom_smooth() +
facet_wrap(~name, scale = "free_y") +
theme_bw() + xlab("") + ylab("") +
theme(strip.background = element_rect(fill="white"))
#> `summarise()` has grouped output by 'date'. You can override using the `.groups` argument.
#> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Standard ggplot approaches can be used to map the
coefficient estimates as in the figure below, which summarises the
coefficient estimates for pr (b_pr) by year.
The code uses the year variable to facets the maps of the
spatially varying coefficients. Notice the increasingly negative
relationship of Maximum temperature with NDVI over time, with a distinct
South West (low) to North East (high) gradient.
# title
tit <-expression(paste(""*beta[`tmax`]*""))
# plot
ggplot() +
geom_sf(data = vcs, aes(col = b_tmax)) +
scale_colour_continuous_c4a_div(palette="brewer.rd_yl_bu",name = tit) +
facet_wrap(~year) +
theme_bw() +
theme(
legend.position = "inside",
legend.position.inside = c(0.9, 0.17),
strip.background = element_rect(fill="white", colour = "white"),
strip.text = element_text(size = 8, margin = margin(b=4)),
axis.title=element_blank(),
axis.text=element_blank(),
axis.ticks=element_blank()) 
tmax (Maximum
temperature) coefficient estimates over space and time.A second approach is to estimate the varying coefficients over a
grid. The code below creates the chaco_grid spatial object,
in the same way as was done in the introductory vignette, but this time
also adding dummy variables:
# create a grid object from the Chaco data
chaco_grid <-
st_make_grid(chaco, square=FALSE,n=20) |>
st_sf()
# rename the geometry, sort row names
st_geometry(chaco_grid) = "geometry"
rownames(chaco_grid) <- 1:nrow(chaco_grid)
# create and add coordinates X and Y / 1000
coords <- chaco_grid |> st_centroid() |> st_coordinates()
chaco_grid <- chaco_grid |> bind_cols(coords/1000)
# add dummy variables to the grid
chaco_grid <-
chaco_grid |>
mutate(Intercept = NA_real_,
tmax = NA_real_,
pr = NA_real_)You may wish to inspect this object and what it contains :
In order to generate annual gridded estimates of the coefficients
over grid space, the model passed to the calculate_vcs
function requires the space and time inputs that were used in the model:
X, Y and month. The
chaco_grid object has locational attributes but the
temporal attributes need to be created. The for loop below
does this, extracts the coefficient estimates for the specified month -
here the sixth month in each year - and then combines them.
# create time slices
years <- sort(unique(vcs$year))
# calculate over the grid for each time slice
res_out <- matrix(nrow = nrow(chaco_grid), ncol = 0)
for (i in 1:length(years)){
# convert years to months - here getting month 6
month.val <- ((years[i]-min(years)) * 12) + 6
res.i <- calculate_vcs(input_data = chaco_grid |> mutate(month = month.val),
mgcv_model = m6,
terms = c("Intercept", "tmax", "pr"))
# select all the coefficient estimates
res.i <-
res.i |>
st_drop_geometry() |>
select(starts_with("b_"),
starts_with("se_"),
starts_with("t_"))
# rename coefficients estimates
names(res.i) <- paste0(names(res.i), "_", years[i])
# bind to the result
res_out <- cbind(res_out, res.i)
cat(years[i], "\t")
}The result can be joined to the grid and mapped. The code below does
this for tmax (Maximum temperature) as this has space-time
TP smooth. The maps describe the variations over time and space of the
relationship between the tmax predictor variable and
ndvi. The maps indicate the increasingly negative
relationship over time and the South West to North East gradient as
before.
# join to the grid
chaco_grid |>
cbind(res_out) |>
# select the variables and pivot longer
select(starts_with("b_tmax")) |>
# rename and select with transmute
transmute(`2012` = b_tmax_2012, `2013` = b_tmax_2013,
`2014` = b_tmax_2014, `2015` = b_tmax_2015,
`2016` = b_tmax_2016, `2017` = b_tmax_2017,
`2018` = b_tmax_2018, `2019` = b_tmax_2019,
`2020` = b_tmax_2020, `2021` = b_tmax_2021,
`2022` = b_tmax_2022
) |>
pivot_longer(-geometry) |>
# make the new time object a factor (to enforce plotting order)
mutate(name = factor(name, levels = years)) |>
# 4. and plot
ggplot() +
geom_sf(aes(fill = value), col = NA) +
# adjust default shading
scale_fill_continuous_c4a_div("brewer.rd_yl_bu", name = tit) +
# facet
facet_wrap(~name, ncol = 4) +
# apply and modify plot theme
theme_bw() +
theme(
legend.position = "inside",
legend.position.inside = c(0.9, 0.17),
strip.background = element_rect(fill="white", colour = "white"),
strip.text = element_text(size = 8, margin = margin(b=4)),
axis.title=element_blank(),
axis.text=element_blank(),
axis.ticks=element_blank()) 
tmax (Maximum
temperature) coefficient estimates over space and time and the grid
surface.2.6 Summary
The above subsections describes a sequence of model construction, investigation and refinement, which, after specifying an initial Gaussian model that incorporated space-time smooths for each predictor variable, required a number of GAM-related considerations to be expanded upon.
- The initial model used Tensor Product smooths, under an assumption
of space-time interactions, rather than separate space and time effects.
These required the basis dimensions of each component to be specified
(through the
dparameter), here specifying a 2-dimensional basis for location and a 1-dimensional one for time). - Additionally thin plate regression smooths (
tp) were used to model location because they are isotropic with no directional bias) and can handle irregularly located observations. Cubic regression smooths (cr) were used to model time effects because there known to be good at handling seasonal or long-term effects. - The model was fitted using REML as this tends to provide more stable estimates of the smoothing parameters and robust estimates of the variance. However, model diagnostic investigations suggested residual heterogeneity, potentially due to the incorrect assumption of a Gaussian response.
- A sequence of refinements were then undertaken and a number of considerations to be discussed: different model families and when they should be used, quantifying the magnitude of the effect size of the model terms, and checking and adjusting the knots used in model fitting procedure in order to determine the final model form.
In this case, the final model specified was Gaussian with a combined
space-time TP smooth for tmax and a spatial smooth for
pr with the number of knots set to 120.
So although statistical practice tends towards fitting the more complex model, penalizing the estimated effects and then evaluating the magnitude of those estimated effects, as was done here, it also suggests that there may be some benefit in undertaking some kind of automatic model selection (see Section 3).
Next they can be mapped over space and time using the original
observation locations and the 
