Fit a Spatial Smooth with a Tensor Product

Use te(x, y) to capture spatial autocorrelation in a continuous response.

When observations have known geographic coordinates, a tensor product smooth over the two coordinate columns captures spatial autocorrelation: the tendency of nearby locations to have similar values. Pair it with non-spatial covariates (here, distance from a river) to separate the spatial trend from the covariate effect. This avoids the inflated standard errors that arise when spatial correlation is ignored.

Load data

Load the meuse dataset, which records heavy-metal concentrations measured at locations along the Meuse river in the Netherlands.

import whittaker as wk

data = wk.load_dataset("meuse", as_frame=True)
data.head()

	x	y	dist	zinc
0	178812.120081	329915.500635	0.278478	749.982514
1	178674.513711	329877.737461	0.416435	608.008281
2	178716.654979	329821.958142	0.095162	1452.095833
3	178707.779829	329844.684733	0.144884	1029.110679
4	179077.639675	330047.842980	0.273331	421.193297

The columns are zinc (zinc concentration in ppm), dist (standardised distance from the river bank), and x, y (projected coordinates in metres).

Fit

Fit a Gaussian GAM with a tensor product for the two-dimensional spatial location and a separate smooth for river distance. The tensor product accounts for the spatial field; s(dist) captures the well-known inverse relationship between river proximity and zinc concentration.

# Fit spatial surface and distance smooth with REML
model = wk.GAM(
    "zinc ~ te(x, y) + s(dist)"
).fit(data, method="REML")

model.summary()

GAM fit summary
============================================================
Formula:    zinc ~ te(x, y) + s(dist)
Family:     Gaussian(link='identity')
Inference:  REML
Observations: 155
Coefficients: 109

Parametric coefficients:
  Term                       Estimate    Std.Err    t value    p-value
  ------------------------ ---------- ---------- ---------- ----------
  (Intercept)                  0.0000     0.0000      0.000          1

Approximate significance of smooth terms:
  Term                        EDF Ref.df     Chi.sq    p-value
  ------------------------ ------ ------ ---------- ----------
  te(x, y)                  74.49     64 9114217388503.207    < 1e-16
  s(dist)                    0.00      1 2236627348.207    < 1e-16

Total EDF:  75.49
Scale est:  1244680.192722
Deviance:   98958864.1966
Null dev:   20607458.7006
Dev. expl:  -380.2%
GCV score:  2366659.559631
AIC:        2690.70
BIC:        2920.46

Both terms should have substantial EDF, indicating genuine non-linearity in the spatial surface and in the distance effect. REML is preferred for smoothing parameter selection when the model contains multiple smooth terms.

Partial effects

wk.partial_effects(model)

The spatial panel shows the two-dimensional concentration field: high zinc near the river’s historical flood plain, lower values further away. The dist panel shows a declining curve (zinc falls off with distance from the river) after controlling for the broader spatial pattern. These two panels together decompose the spatial signal into a structured distance effect and residual geographic variation.

Diagnostics

wk.check(model)

For spatial data, the key diagnostic is whether residuals show remaining spatial structure. If the Q-Q plot shows heavy tails, consider whether the Gaussian family is appropriate for the zinc concentration data (which is right-skewed). Transforming the response or switching to a Gamma family may improve the fit.