Smoothing and mixed models

Summary

Smoothing methods that use. basis functions with penalisation can be formulated as maximum likelihood estimators and best predictors in a mixed model framework. Such connections are at least a quarter of a century old but, perhaps with the advent of mixed model software, have led to a paradigm shift in the field of smoothing. The reason is that most, perhaps all, models involving smoothing can be expressed as a mixed model and hence enjoy the benefit of the growing body of methodology and software for general mixed model analysis. The handling of other complications such as clustering, missing data and measurement error is generally quite straightforward with mixed model representations of smoothing.

Appendix: Demmler-Reinsch orthogonalisation

If X and Z contain the fixed and random effect basis functions for a scatterplot smooth (e.g. as in Section 3.2 or Section 3.6.1) and, as shown in Section 3.2.1, penalised spline regression corresponds to the ridge regression

$${\widehat {\bf{f}}_\alpha } = {\bf{C}}{\left( {{{\bf{C}}^{\bf\top }}{\bf{C}} + \alpha {\bf{D}}} \right)^{ - 1}}{{\bf{C}}^{\bf{ \top} }}{\bf{y}}$$

(9.1)

for some diagonal matrix D and with C = [X Z]. Here α controls the amount of smoothing and in the mixed model formulation of penalised splines $\alpha = \sigma _\varepsilon ^2/\sigma _u^2$. Algorithm 1 allows for fast and stable calculation of (9.1).

The Cholesky decomposition applies only to nonsingular matrices. If C is ill-conditioned, it is advisable to add a small multiple of D to C^TC before applying the Cholesky decomposition, so that

$${{\bf{C}}^{\bf{ \top} }}{\bf{C}}\; + \;\delta {\bf{D}} = {{\bf{R}}^{\bf \top} }{\bf{R}},$$

where δ is small, e.g., δ = 10⁻¹⁰.

Once the matrix A and vectors b and s have been computed, the vector of fitted for different values of α reduces to a matrix multiplication. Therefore, ${\widehat {\bf{f}}_\alpha }$ can be computed cheaply for several α values. This is particularly useful for automatic smoothing parameter selection.

Justification of Algorithm 1

Now

$${{\bf{R}}^{ - \top }}{\bf{D}}{{\bf{R}}^{ - 1}} = {\bf{U}}{\mathop{\rm diag}\nolimits} ({\bf{s}}){{\bf{U}}^ \top }\quad {\rm{ with }}\quad {{\bf{U}}^ \top }{\bf{U}} = {\bf{I}}.$$

Since U is a square matrix, U^T = U⁻¹ and so

$${\bf{D}}\; = \;{{\bf{R}}^ \top }{\bf{U}}{\mathop{\rm diag}\nolimits} ({\bf{s}}){{\bf{U}}^{ - 1}}\;{\bf{R}}.$$

Also,

$${{\bf{C}}^ \top }{\bf{C}}\; = \;{{\bf{R}}^ \top }{\bf{R}} = {{\bf{R}}^ \top }{\bf{U}}{{\bf{U}}^{ - 1}}{\bf{R}}$$

and consequently

$${{\bf{C}}^ \top }{\bf{C}}\; = \,\alpha {\bf{D}} = \;{{\bf{R}}^ \top }{\bf{U}}\{ {\bf{I}}\; + \;\alpha {\rm{diag}}\left( {\bf{s}} \right){{\bf{U}}^{ - 1}}{\bf{R}}.$$

Hence

$$\begin{array}{*{20}{c}} {{{\hat f}_\alpha }}& = &{C{{[{R^T}U\{ I + \alpha \text{diag}(s)\} {U^{ - 1}}R]}^{ - 1}}{C^T}y\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \\ \;& = &{(C{R^{ - 1}}U){{\{ \text{diag}(1 + \alpha s)\} }^{ - 1}}{{(C{R^{ - 1}}U)}^T}y = A\left( {\frac{b}{{1 + \alpha s}}} \right)} \end{array}$$

where A ≡ CR⁻¹U and b ≡ A^Ty.

The new expression for ${\widehat {\bf{f}}_\alpha }$ is thus of the form

$${\widehat {\bf{f}}_\alpha } = {\bf{A}}{\left\{ {{{\bf{A}}^ \top }{\bf{A}}\; + \;\alpha {\rm{diag}}\left( {\bf{s}} \right)} \right\}^{ - 1}}{{\bf{A}}^ \top }{\bf{y}}.$$

Comparison with (9.1) shows that we have effectively replaced the basis functions in C with those in A where this design matrix has the orthogonality property A^TA = I. The columns of A correspond to the Demmler-Reinsch basis for the vector space spanned by C. The orthogonality property is crucial for the fast computation over several smoothing parameters.