A note on P-spline additive models with correlated errors

Summary

We consider additive models with k smooth terms and correlated errors, and use the penalised spline approach of Eilers & Marx (1996) to estimate the smooth functions. We obtain explicit expressions for the hat-matrix of the model and each individual curve. P-splines are represented as mixed models and REML is used to select the smoothing and correlation parameters. The method is applied to the analysis of some time series data.

Appendices

Appendix A

We derive expression (12) for the hat-matrix in the additive model (9) with k terms. First we define ${\boldsymbol{B}}_{ - i}^* = \left( {{\boldsymbol{B}}_1^*: \ldots :{\boldsymbol{B}}_{i - 1}^*:{\boldsymbol{B}}_{i + 1}^*: \ldots :{\boldsymbol{B}}_k^*} \right)$and ${{\boldsymbol{a}}_{ - i}} = \left( {{\boldsymbol{a}}_1^\prime , \ldots ,{\boldsymbol{a}}_{i - 1}^\prime ,{\boldsymbol{a}}_{i + 1}^\prime , \ldots ,{\boldsymbol{a}}_k^\prime } \right)\prime$. We rewrite (9) as

$${\rm{E}}\left( {\boldsymbol{y}} \right) = {\boldsymbol{Ba}}\;{\rm{where}}\;{\boldsymbol{B}}{\rm{ = }}\left( {1:{\boldsymbol{B}}_i^*:{\boldsymbol{B}}_{ - i}^*} \right),\quad {\boldsymbol{a}} = \left( {\alpha ,{\boldsymbol{a}}_i^\prime ,{\boldsymbol{a}}_{ - i}^\prime } \right)\prime ,$$

and obtain ${\boldsymbol{\hat a}}$ by minimising

$$S = \left( {{\boldsymbol{y}} - \alpha 1 - {\boldsymbol{B}}_i^*{{\boldsymbol{a}}_i} - {\boldsymbol{B}}_{ - i}^*{{\boldsymbol{a}}_{ - i}}} \right)\prime \left( {{\boldsymbol{y}} - \alpha 1 - {\boldsymbol{B}}_i^*{{\boldsymbol{a}}_i} - {\boldsymbol{B}}_{ - i}^*{{\boldsymbol{a}}_{ - i}}} \right) + {\boldsymbol{a}}_i^\prime {{\boldsymbol{P}}_i}{{\boldsymbol{a}}_i} + {\boldsymbol{a}}_{ - i}^\prime {{\boldsymbol{P}}_{ - i}}{{\boldsymbol{a}}_{ - i}}.$$

Taking derivatives with respect to the components of a we obtain

$${{\partial S} \over {\partial \alpha }} = 1\prime \left( {{\boldsymbol{y}} - \hat \alpha 1 - {\boldsymbol{B}}_i^*{{{\boldsymbol{\hat a}}}_i} - {\boldsymbol{B}}_{ - i}^*{{{\boldsymbol{\hat a}}}_{ - i}}} \right) = 1\prime {\boldsymbol{y}} - 1\prime 1\hat \alpha = 0$$

(17)

$${{\partial S} \over {\partial {{\boldsymbol{a}}_i}}} = - {\boldsymbol{B}}_i^{*\prime }\left( {{\boldsymbol{y}} - \hat \alpha 1 - {\boldsymbol{B}}_i^*{{{\boldsymbol{\hat a}}}_i} - {\boldsymbol{B}}_{ - i}^*{{{\boldsymbol{\hat a}}}_{ - i}}} \right) + {P_i}{{\boldsymbol{\hat a}}_i} = 0$$

(18)

$${{\partial S} \over {\partial {{\boldsymbol{a}}_{ - i}}}} = - {\boldsymbol{B}}_{ - i}^{*\prime }\left( {{\boldsymbol{y}} - \hat \alpha 1 - {\boldsymbol{B}}_i^*{{{\boldsymbol{\hat a}}}_i} - {\boldsymbol{B}}_{ - i}^*{{{\boldsymbol{\hat a}}}_{ - i}}} \right) + {{\boldsymbol{P}}_{ - i}}{{\boldsymbol{\hat a}}_{ - i}} = 0.$$

(19)

Equation (17) yields $\hat \alpha = {\left( {1\prime 1} \right)^{ - 1}}1\prime {\boldsymbol{y}} = {\boldsymbol{\bar y}}$. From (18) we obtain

$$\left( {{\boldsymbol{B}}_i^{*\prime }{\boldsymbol{B}}_i^* + {{\boldsymbol{P}}_i}} \right){{\boldsymbol{\hat a}}_i} = {\boldsymbol{B}}_i^{*\prime }\left( {{\boldsymbol{y}} - {\boldsymbol{B}}_{ - i}^*{{{\boldsymbol{\hat a}}}_{ - i}}} \right)$$

and so

$${\boldsymbol{B}}_i^*{{\boldsymbol{\hat a}}_i} = {\boldsymbol{S}}_i^*\left( {{\boldsymbol{y}} - {\boldsymbol{B}}_{ - i}^*{{{\boldsymbol{\hat a}}}_{ - i}}} \right)$$

(20)

where

$${\boldsymbol{S}}_i^* = {\boldsymbol{B}}_i^*{\left( {{\boldsymbol{B}}_i^{*\prime }{\boldsymbol{B}}_i^* + {{\boldsymbol{P}}_i}} \right)^ - }{\boldsymbol{B}}_i^{*\prime }$$

(21)

is the centred smoother matrix from the model with the single smooth term ${\boldsymbol{B}}_i^*$. Here A⁻ represents a generalised inverse of A (although the generalised inverse is not unique, ${\boldsymbol{S}}_i^*$ is invariant to the choice of the generalised inverse; (see for example Harville 1999, chap. 9). Substitution of (20) in (19) yields

$$\left( {{\boldsymbol{B}}_{ - i}^{*\prime }\left( {{\boldsymbol{I}} - {\boldsymbol{S}}_i^*} \right){\boldsymbol{B}}_{ - i}^* + {{\boldsymbol{P}}_{ - i}}} \right){{\boldsymbol{\hat a}}_{ - i}} = {\boldsymbol{B}}_{ - i}^{*\prime }\left( {{\boldsymbol{I}} - {\boldsymbol{S}}_i^*} \right){\boldsymbol{y}}$$

from which we find

$${\boldsymbol{B}}_{ - i}^*{{\boldsymbol{\hat a}}_{ - i}} = {\boldsymbol{H}}_{ - i}^*{\boldsymbol{y}}$$

where

$${\boldsymbol{H}}_{ - i}^* = {\boldsymbol{B}}_{ - i}^*{\left( {{\boldsymbol{B}}_{ - i}^{*\prime }\left( {{\boldsymbol{I}} - {\boldsymbol{S}}_i^*} \right){\boldsymbol{B}}_{ - i}^* + {{\boldsymbol{P}}_{ - i}}} \right)^ - }{\boldsymbol{B}}_{ - i}^{*\prime }\left( {{\boldsymbol{I}} - {\boldsymbol{S}}_i^*} \right)$$

(22)

and ${\boldsymbol{H}}_{ - i}^*$ is the centred hat-matrix of a weighted additive model (with weights $\left. {\left( {{\boldsymbol{I}} - {\boldsymbol{S}}_i^*} \right)} \right)$. Thus, ${\boldsymbol{\hat y}} = 1\hat \alpha + {\boldsymbol{B}}_i^*{{\boldsymbol{\hat a}}_i} + {\boldsymbol{B}}_{ - i}^*{{\boldsymbol{\hat a}}_{ - i}} = {{\boldsymbol{H}}_k}{\boldsymbol{y}}$ where

$${{\boldsymbol{H}}_k} = 11\prime /n + {\boldsymbol{S}}_i^* + \left( {{\boldsymbol{I}} - {\boldsymbol{S}}_i^*} \right){\boldsymbol{H}}_{ - i}^*$$

(23)

is the hat-matrix for an additive model with k terms, as required.

When ${\boldsymbol{\epsilon }} \sim {\cal N}\left( {0,{\sigma ^2}{\boldsymbol{V}}} \right)$, it is straightforward to show that (13) becomes

$${\hat f_i}\left( {{{\boldsymbol{x}}_i}} \right) = {\boldsymbol{S}}_{iV}^*\left( {{\boldsymbol{I}} - {\boldsymbol{H}}_{ - iV}^*} \right){\boldsymbol{y}}$$

where

$$\begin{array}{*{20}c}{\quad {\boldsymbol{S}}_{iV}^*} \; = \; {{\boldsymbol{B}}_i^*{{\left( {{\boldsymbol{B}}_i^{*\prime }{{\boldsymbol{V}}^{ - 1}}{\boldsymbol{B}}_i^* + {{\boldsymbol{P}}_i}} \right)}^ - }{\boldsymbol{B}}_i^{*\prime }{{\boldsymbol{V}}^{ - 1}}\quad {\rm{as}}\;{\rm{in}}\left( {21} \right)}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \\ {{\boldsymbol{H}}_{ - iV}^*} \; = \; {{\boldsymbol{B}}_{ - i}^*{{\left( {{\boldsymbol{B}}_{ - i}^{*\prime }{{\boldsymbol{V}}^{ - 1}}\left( {{\boldsymbol{I}} - {\boldsymbol{S}}_{iV}^*} \right){\boldsymbol{B}}_{ - i}^* + {{\boldsymbol{P}}_{ - i}}} \right)}^ - }{\boldsymbol{B}}_{ - i}^{*\prime }{{\boldsymbol{V}}^{ - 1}}\left( {{\boldsymbol{I}} - {\boldsymbol{S}}_{iV}^*} \right).} \\ \end{array} $$

Appendix B

The additive P-spline setup chooses a for given λ to minimise

$$S\left( {\boldsymbol{a}} \right) = \left( {{\boldsymbol{y}} - {\boldsymbol{Ba}}} \right)\prime {{\boldsymbol{V}}^{ - 1}}\left( {{\boldsymbol{y}} - {\boldsymbol{Ba}}} \right) + {\boldsymbol{a}}\prime {\boldsymbol{Pa}}$$

(24)

with ${\boldsymbol{B}} = \left( {1:{\boldsymbol{B}}_1^*: \ldots :{\boldsymbol{B}}_k^*} \right)$ and P = blockdiag(0, P₁,…, P_k) where ${{\boldsymbol{P}}_j} = {\lambda _j}{\boldsymbol{D}}_j^\prime {{\boldsymbol{D}}_j}$. The value of a that minimises (24) satisfies

$$\left( {{\boldsymbol{B}}\prime {{\boldsymbol{V}}^{ - 1}}{\boldsymbol{B}} + {\boldsymbol{P}}} \right){\boldsymbol{a}} = {\boldsymbol{B}}\prime {{\boldsymbol{V}}^{ - 1}}{\boldsymbol{y}}.$$

(25)

We show that (25) is the result of estimation in a mixed model. We write ${\boldsymbol{Ba}} = {\boldsymbol{B\tilde G\beta }} + {{\boldsymbol{B}}^*}{\boldsymbol{Zu}}$ with ${\boldsymbol{\tilde G}} = {\rm{blockdiag}}\left( {1,{\boldsymbol{G}}} \right),\;{{\boldsymbol{B}}^*} = \left( {{\boldsymbol{B}}_1^*: \ldots :{\boldsymbol{B}}_k^*} \right)$ and G and Z defined as follows: G = blockdiag(1, G₁,…, G_k) where G_i = $\left( {{{\boldsymbol{g}}_i},{\boldsymbol{g}}_i^2, \ldots ,{\boldsymbol{g}}_i^{{q_i} - 1}} \right)$, q_i is the order of penalty for the ith regressor x_i and ${\boldsymbol{g}}_i^\prime = \left( {1,2, \ldots ,{p_i}} \right)$ where p_i is the rank of B_i; Z = blockdiag(Z₁,…, Z_k) with ${{\boldsymbol{Z}}_i} = {\boldsymbol{D}}_i^\prime {\left( {{{\boldsymbol{D}}_i}{\boldsymbol{D}}_i^\prime } \right)^{ - 1}}$. Substituting ${\boldsymbol{Ba}} = {\boldsymbol{B\tilde G\beta }} + {{\boldsymbol{B}}^{\boldsymbol{*}}}{\boldsymbol{Zu}}$, in (25) and using D_iG_i = 0 we find

$${{\boldsymbol{B}}^\prime }{{\boldsymbol{V}}^{ - 1}}{\boldsymbol{B}}\widetilde{\boldsymbol{G}}{\boldsymbol{\beta }} + {{\boldsymbol{B}}^\prime }{{\boldsymbol{V}}^{ - 1}}{{\boldsymbol{B}}^*}{\boldsymbol{Zu}} + {{\boldsymbol{P}}^*}{\boldsymbol{Zu}} = {{\boldsymbol{B}}^\prime }{{\boldsymbol{V}}^{ - 1}}{\boldsymbol{y}}$$

(26)

with P* =blockdiag(P₁,…,P_k). Multiplying (26) by ${\boldsymbol{\tilde G}}\prime$ gives

$${{\boldsymbol{\tilde G}}^\prime }{{\boldsymbol{B}}^\prime }{{\boldsymbol{V}}^{ - 1}}{\boldsymbol{B\tilde G\beta }} + {{\boldsymbol{\tilde G}}^\prime }{{\boldsymbol{B}}^\prime }{{\boldsymbol{V}}^{ - 1}}{{\boldsymbol{B}}^*}{\boldsymbol{Zu}} = {{\boldsymbol{\tilde G}}^\prime }{{\boldsymbol{B}}^\prime }{{\boldsymbol{V}}^{ - 1}}{\boldsymbol{y}}$$

(27)

(again using D_iG_i = 0) while multiplying (26) by blockdiag(0, Z′) gives

$${Z^\prime }{B^{*\prime }}{V^{ - 1}}B\tilde G\beta + {Z^\prime }{B^{*\prime }}{V^{ - 1}}{B^*}Zu + {Z^\prime }{P^*}Zu = {Z^\prime }{B^{*\prime }}{V^{ - 1}}y.$$

(28)

Let r_i = ndx_i + bdeg_i − pord_i be the number of columns of D_i, the difference matrix for the ith regressor x_i. Then Z′ P* Z = blockdiag(λ₁I_r1,…, λ_kIr_k) = Λ. Then (27) and (28) can be written

$$\left[ {\begin{array}{*{20}c}{{{\tilde G}^\prime }{B^\prime }{V^{ - 1}}B\tilde G} \;\;\;\;\;\;\;\;\;\; {{{\tilde G}^\prime }{B^\prime }{V^{ - 1}}{B^*}Z} \\ {{Z^\prime }{B^{*\prime }}{V^{ - 1}}B\tilde G} \;\;\;\;\;\;\;\; {{Z^\prime }{B^{*\prime }}{V^{ - 1}}{B^*}Z + {\rm{\Lambda }}} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c}{\hat \beta } \\ {\hat u} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c}{{{\tilde G}^\prime }{B^\prime }} \\ {{Z^\prime }{B^{*\prime }}} \\ \end{array} } \right]{V^{ - 1}}y.$$

Thus, $\hat \beta$ and $\hat u$ are estimates that arise from the mixed model

$$y = B\tilde Gb + {B^*}Zu + \epsilon $$

where $u\sim{\cal N}\left( {0,\sigma _u^2} \right),\epsilon \sim{\cal N}\left( {0,{\sigma ^2}V} \right){\rm{ and }}\sigma _u^2 = {\sigma ^2}{{\rm{\Lambda }}^{ - 1}}$