Combatting edge effects in spatial smoothing

David L. Miller (joint work with Esther Jones and Jason Matthiopoulos)

ISEC 2012, Oslo, Norway

Spatial smoothing with splines

• Just talking about 2-D smoothing in a “box”
• No fancy complex regions
• This should be easy…
• BUT: still have edge effects

What’s going on?

$f(x,y) = {\color{red}\sum_{i=1}^n \delta_i \eta(r_i)} + \sum_{j=1}^3 \alpha_j \phi_j(x,y),$ where $r_i=\sqrt{(x-x_i)^2+(y-y_i)^2}$

What’s going on?

$f(x,y) = \sum_{i=1}^n \delta_i \eta(r_i) + {\color{red}\sum_{j=1}^3 \alpha_j \phi_j(x,y)},$

Can we get around this?

• change to a local basis (e.g. P-splines)
• computationally expensive – tensor products
• treat this as a complex region (e.g. soap)
• wasteful, complex and not well motivated
• remove the offending parts of the spline
• $f$ not guaranteed continuous
• introduce extra penalty for the linear parts
• numerically icky, complicated, not guaranteed to work

Duchon splines

• same basis functions as thin plate splines
• modify penalty
• remove the linear terms
• but keep it mathematically consistent

Mathematical interlude

Thin plate penalty (in 2D, with $2^\text{nd}$ order derivative penalty): $P_{2,2} = \int \int_{\mathbb{R}^2} \left (\frac{\partial^2 f(x,y)}{\partial x^2} \right )^2 + 2\left (\frac{\partial^2 f(x,y)}{\partial x \partial y} \right )^2 + \left (\frac{\partial^2 f(x,y)}{\partial y^2} \right )^2 \text{d}x \text{d}y$ More generally: $P_{m,2} = \int \int_{\mathbb{R}^2} \sum_{\nu_1 + \nu_2=m} \frac{m!}{\nu_1! \nu_2!} \left( \frac{\partial^m f \left (x, y \right )}{\partial x^{\nu_1} \partial y^{\nu_2}} \right)^2 \text{d} x \text{d} y$ to ensure that $f$ remains continuous, $2m>$ # dimensions (here 2).

$m$ also dictates # linearly independent polynomials

Mathematical interlude

Take Fourier transform of derivatives: $P_{m,2} = \int \int_{\mathbb{R}^2} \sum_{\nu_1 + \nu_2=m} \frac{m!}{\nu_1! \nu_2!} \left ( \mathfrak{F} \frac{\partial^m f}{\partial x^{\nu_1} \partial y^{\nu_2}} \left ( \boldsymbol{\tau}\right ) \right )^2 \text{d} \boldsymbol{\tau}.$

Mathematical interlude

What about weighting on the frequencies? $\int \int_{\mathbb{R}^2} w(\boldsymbol{\tau}) \sum_{\nu_1 + \nu_2=m} \frac{m!}{\nu_1! \nu_2!} \left ( \mathfrak{F} \frac{\partial^m f}{\partial x_1^{\nu_1} \partial y^{\nu_2}} \left (\boldsymbol{\tau} \right ) \right )^2 \text{d} \boldsymbol{\tau}$

New penalty

• Letting $w(\boldsymbol{\tau})=\lvert \boldsymbol{\tau} \rvert^{2s}$ $\breve{P}_{m,2} = \int \int_{\mathbb{R}^2} \lvert \boldsymbol{\tau} \rvert^{2s} \sum_{\nu_1 + \nu_2=m} \frac{m!}{\nu_1! \nu_2!}\left ( \mathfrak{F} \frac{\partial^m f}{\partial x^{\nu_1} \partial y^{\nu_2}} \left (\boldsymbol{\tau} \right ) \right )^2 \text{d} \boldsymbol{\tau}$
• We can now set $s$ to “make up” for reducing $m$ : $2m>d-2s$
• $s>d/2-m \quad (s \in \mathbb{Z}/2) \quad \Rightarrow$ $s=1/2, m=1$
• $\breve{P}_{1,2} = \int \int_{\mathbb{R}^2} \lvert \boldsymbol{\tau} \rvert \left\{ \left ( \mathfrak{F} \frac{\partial f}{\partial x} \left(\boldsymbol{\tau} \right) \right )^2 + \left (\mathfrak{F} \frac{\partial f}{\partial y}\left(\boldsymbol{\tau} \right) \right )^2 \right\} \text{d}\boldsymbol{\tau}.$

Software

• already implemented in mgcv
• bs="ds", then set m=c(1,0.5) and you’re good to go
• From the above:
b <- gam(z~s(x, y, k=100), data=dat)
• becomes:
b <- gam(z~s(x, y, bs="ds", m=c(1,0.5), k=100), data=dat)
• where m=c($m$,$s$)

Conclusion

• Edge effects can be nasty
• Duchon splines offer an easy solution
• Already implemented in mgcv
• Simple change to code to see if they help
• Paper coming soon Combatting edge effects in spatial models Miller, Lane and Matthiopoulos.