# Fit a Spatial Smooth with a Tensor Product

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.


``` python
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.


``` python
# 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


``` python
wk.partial_effects(model)
```


<style>
  #altair-viz-f318a50f604c47e7bd8cff8ec8485c50.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-f318a50f604c47e7bd8cff8ec8485c50.vega-embed details,
  #altair-viz-f318a50f604c47e7bd8cff8ec8485c50.vega-embed details summary {
    position: relative;
  }
</style>


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


``` python
wk.check(model)
```


<style>
  #altair-viz-7e155de7ec8848bba73e7b0b61777318.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-7e155de7ec8848bba73e7b0b61777318.vega-embed details,
  #altair-viz-7e155de7ec8848bba73e7b0b61777318.vega-embed details summary {
    position: relative;
  }
</style>


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.