Some Special Spline Models

Example 1. Polynomial Spline. For the polynomial spline of order $m$ on ${\cal T}=[a,b]$, the model space is

$${\cal H}=W_m([a,b])=\{f: ~f, f', \cdots, f^{(m-1)} \mbox{absolutely continuous}, \int_a^b (f^{(m)})^2dt<\infty\}.$$ (6)

Let ${\cal T}=[0,1]$. Define an inner product on $W_m([0,1])$ as

$$(f,g) = \sum_{i=0}^{m-1} \int_0^1 (f^{(i)} dt) \int_0^1 (g^{(i)} dt) + \int_0^1 f^{(m)} g^{(m)} dt.$$

Then we have

$${\cal H}_0 = {\mbox span}\{k_0(t), k_1(t), \cdots, k_{m-1}(t)\},$$

$$R_1(s, t) = k_m(s) k_m(t) + (-1)^{m-1} k_{2m}(s-t) ,$$

$$\vert\vert P_1f\vert\vert^2 = \int_0^1 (f^{(m)})^2 dt,$$ (7)

where $k_r(t)=B_r(t)/r!$ and $B_r$'s are Bernoulli polynomials (Craven and Wahba, 1979).

Let ${\cal T}=[a,b]$. Define inner product on $W_m([a,b])$ as

$$(f,g) = \sum_{i=0}^{m-1} f^{(i)}(a) g^{(i)}(a) + \int_a^b f^{(m)} g^{(m)} dt.$$

Then (Wahba, 1990)

$${\cal H}_0 = {\mbox span}\{1, (t-a), \cdots, (t-a)^{m-1}\},$$

$$R_1(s, t) = \int_a^b\frac{\displaystyle (s-u)_{+}^{m-1}}{\displaystyle (m-1)!} \frac{\displaystyle (t-u)_{+}^{m-1}}{\displaystyle (m-1)!} du .$$ (8)

Table lists statements inside ssr for four simple polynomial splines. We assume variables $k1$, $k2$ and $k3$ have been calculated based on Bernoulli polynomials. All expressions of the reproducing kernels in this table are available in our library. Note that the domain under construction () is restricted to $[0,1]$ while the domain under construction () is an arbitrary interval $[a,b]$. Thus one needs to transform a variable into $[0,1]$ before using rk functions in the third column. The rk functions in the fourth column assume the domain $[0,T]$ for any fixed $T>0$. One can calculate the reproducing kernel on $[a,b]$ by a translation, for example $cubic2(t-a)$.

Table: Statements for fitting simple polynomial splines.

m	splines	under construction ()	under construction ()
1	linear	y~1, rk=linear(t)	y~1, rk=linear2(t)
2	cubic	y~k1, rk=cubic(t)	y~t, rk=cubic2(t)
3	quintic	y~k1+k2, rk=quintic(t)	y~t+t**2, rk=quintic2(t)
4	septic	y~k1+k2+k3, rk=septic(t)	y~t+t**2+t**3, rk=septic2(t)

Example 2. Stein Estimate. A James-Stein shrinkage estimate can be regarded as the solution to () with ${\cal T} =\{1, 2, \cdots, K \}$, ${\cal H}= R^K$ and $L_if=f(t_i)$ (Gu, 2002).

For shrinkage toward a constant,

$${\cal H}_0 = \mbox{span}\{1\},$$

$$R_1(s, t) = I_{[s=t]}-1/K,$$

$$\vert\vert P_1f\vert\vert^2 = \sum_{i=1}^K [f(i)-\sum_{j=1}^K f(j)/K]^2.$$ (9)

For shrinkage toward zero,

$${\cal H}_0 = \mbox{empty set},$$

$$R_1(s, t) = I_{[s=t]},$$

$$\vert\vert P_1f\vert\vert^2 = \sum_{i=1}^K f^2(i).$$ (10)

The corresponding rk expressions for shrinkage toward a constant and shrinkage toward zero are shrink1(t) and shrink0(t) respectively. We use these shrinkage methods to model discrete covariates in smoothing spline ANOVA models as will be discussed later.

Example 3. Periodic Spline. For the $m$-th order periodic spline on ${\cal T}=[0,1]$ (Wahba, 1990),

$$\begin{eqnarray*} {\cal H}&=&W_m(per) \\ &=& \{f: f^{(j)}~ \mbox{abs. cont.}, ~ ... ...s-t),\\ \vert\vert P_1f\vert\vert^2 &=&\int_0^1 (f^{(m)})^2 dt. \end{eqnarray*}$$

The function periodic in our library calculates $R_1$ evaluated at specified points. The order $m$ is specified by the argument order with default order=2. For example, one may fit a cubic periodic spline ($m=2$) by

    ssr(y~1, rk=periodic(t))

Example 4. Thin plate spline (TPS) (Wahba, 1990). For a TPS of order $m$ on ${\cal T}=R^d$ with $2m-d>0$,

$${\cal H} = \{ f:~J^d_m(f) < \infty \},$$

$${\cal H}_0 = \{f:~ J^d_m(f)=0 \},$$

$$\vert\vert P_1f\vert\vert = J^d_m(f),$$ (11)

where

$$J^d_m(f)=\sum_{\alpha_1+\cdots+\alpha_d=m}\frac{m!} {\alpha_... ...\alpha_1}\cdots\partial t_d^{\alpha_d}}\right)^2\prod_j dt_j .$$

A pseudo reproducing kernel (conditional positive definite rather than positive definite) for ${\cal H}_1$ is $R_1(s, t) = E(\vert s-t\vert)$, where $\vert s-t\vert$ is the Euclidean distance, and

$$E(u)=\left\{\begin{array}{clcr} &\frac{\displaystyle (-1)^{... ...vert u\vert^{2m-d}, &d \ \ {\mbox odd}. \end{array} \right.$$

The function tp.psuedo in our library calculates this pseudo kernel. The argument order of this function specifies $m$ with default order=2. For example, for $d=2$ and $m=2$, one may fit the TSP model by

    ssr(y~t1+t2, rk=tp.pseudo(list(t1,t2)))

The true kernel discussed in Gu and Wahba (1993a) is calculated by the function tp. It takes longer to compute the true kernel and is only necessary for calculating posterior variances.

Example 5. Spline on the sphere is an extension of both the periodic spline defined on the unit circle and the TPS on $R^2$. Let the domain be ${\cal T}={\cal S}$, where ${\cal S}$ is the unit sphere. Any point $t$ on ${\cal S}$ can be represented as $t=(\theta,\phi)$, where $\theta$ ( $0\le \theta \le 2\pi$) is the longitude and $\phi$ ( $-\pi/2 \le \phi \le \pi/2$) is the latitude. Define

$$J(f)=\int_{t\in {\cal S}} (\Delta^{m/2} f)^2 dt,$$

where $\Delta f$ is the surface Laplacian on the unit sphere

$$\Delta f = \frac{1}{\cos^2 \phi} f_{\theta\theta} + \frac{1}{\cos \phi} (\cos \phi f_{\phi})_{\phi} .$$

The model space ${\cal H}={\cal H}_0 \oplus {\cal H}_1$, where ${\cal H}_0=\mbox{span} \{ 1 \}$ and

$${\cal H}_1 = \{ f \in {\cal L}_2({\cal S}):~~ J(f)<\infty \} .$$

${\cal H}$ is a RKHS with RK $R(s,t)=1+R_1(s,t)$, where

$$R_1(s,t) = \sum_{i=1}^\infty \frac{2i+1}{4\pi} \frac{1}{[i(i+1)]^m} L_i(\cos \gamma (s,t))$$

is the rk of ${\cal H}_1$, $\gamma(s,t)$ is the angle between $s$ and $t$, and $L_i$'s are the Legendre polynomials. $\vert\vert P_1f\vert\vert=J(f)$. $R_1$ is in the form of an infinite series which is inconvenient to compute. Closed form expressions are only available for $m=2$ and $m=3$. Wahba (1981) proposed replacing $J$ by a topologically equivalent norm $Q$ under which closed form of rk's can be derived. See Wahba (1981), Wahba (1982), Wahba (1990) and Wahba and Luo (1996) for more details. The function sphere in our library calculates $R_1$ under the norm $Q$ for $2 \le m \le 6$. The argument order of this function specifies $m$ with default order=2. For example, for $m=3$, one may fit a spline on the sphere by

    ssr(y~t1+t2, rk=sphere(cbind(t1,t2),order=3))

where $t1$ and $t2$ are longitude and latitude respectively.

Example 6. $L$-spline. The penalty term, $\vert\vert P_1f\vert\vert$, is usually used to penalize the roughness of the function $f$. However, sometimes it is advantageous to use other forms of penalty. For example, prior information may be incorporated or even estimated by a penalty to the departure of the non-parametric function $f$ from a specific parametric model (Wahba, 1990; Heckman and Ramsay, 2000). Let $L$ be a linear differential operator $L=D^m+\sum_{j=0}^{m-1}\omega_j D^j$, where $D^j$ denotes the $j$th derivative operator and the $\omega$'s are continuous real-valued weight functions. The spline estimate with the penalty $\vert\vert P_1f\vert\vert^2=\int_a^b (Lf(t))^2 dt$ is called an $L$-spline. See Heckman and Ramsay (2000) and Gu (2002) for more details about the $L$-spline. The lspline function in our library calculates reproducing kernels for the following four $L$-spline models.

(a) Suppose that ${\cal T}=[0,1]$. If prior knowledge suggests that $f$ is close to a linear combination of $\sin 2\pi t$ and $\cos 2\pi t$, one may use ${\cal H}=W_2(per) \ominus \{ 1 \}$, and $L=D^2+(2\pi)^2$. Then

$$\begin{eqnarray*} {\cal H}_0 &=& \mbox{span} \{ \sin 2\pi t,~ \cos 2\pi t \}, \\... ...2}^{\infty} \frac{2}{(2\pi)^4(1-\nu^2)^2} \cos 2\pi \nu (s-t). \end{eqnarray*}$$

The statement for fitting such a model is

    ssr(y~sin(2*pi*t)+cos(2*pi*t)-1, rk=lspline(t,type="sine0"))

If we want to include the constant in the model space, then ${\cal H}=W_3(per)$, $L=D[D^2+(2\pi)^2]$,

$$\begin{eqnarray*} {\cal H}_0 &=& \mbox{span} \{ 1,~ \sin 2\pi t,~ \cos 2\pi t \}... ...fty} \frac{2}{(2\pi)^6 \nu^2 (1-\nu^2)^2} \cos 2\pi \nu (s-t). \end{eqnarray*}$$

The statement for fitting such a model is

    ssr(y~sin(2*pi*t)+cos(2*pi*t), rk=lspline(t,type="sine1"))

(b) Suppose that ${\cal T}=[0,T]$. If prior knowledge suggests that $f$ is close to a linear combination of $1$ and $\exp (-t)$, one may use ${\cal H}=W_2([0,T])$, and $L=D^2+D$. Then

$${\cal H}_0 = \mbox{span} \{ 1,~ \exp (-t) \},$$

$$R_1(s, t) = \min(s, t)+e^{-t}+e^{-s}-e^{\min(s,t)-s} -e^{\min(s, t)-t}- \frac{e^{-(s+t)}}{2} +\frac{e^{2\min(s, t)-s-t}}{2} .$$ (12)

The statement for fitting such a model is

    ssr(y~exp(-t), rk=lspline(t,type="exp"))

(c) Suppose that ${\cal T}=[0,T]$. If prior knowledge suggests that $f$ is close to the logistic function $\exp(t)/(1+\exp(t))$, one may use ${\cal H}=W_1([0,T])$, and $L=D-\frac{1}{1+e^t}I$. Then

$$\begin{eqnarray*} {\cal H}_0 &=& \mbox{span} \{ \exp (t) / (1+\exp (t)) \}, \\ ... ...)- 2e^{-\min(s, t)}-\frac{1}{2}e^{-2\min(s, t)} + \frac{5}{2}] . \end{eqnarray*}$$

The statement for fitting such a model is

    ssr(y~I(exp(t)/(1+exp(t)))-1, rk=lspline(t,type="logit"))

(d) Suppose that ${\cal T}=[0,T]$. If prior knowledge suggests that $f$ is close to a linear combination of $1$, $t$, $\sin 2\pi t$ and $\cos 2\pi t$, one may use ${\cal H}=W_4([0,T])$, and $L=D^4+D^2$. Then

$$\begin{eqnarray*} {\cal H}_0 &=& \mbox{span} \{ 1, ~ t, ~ \sin (t),~ \cos (t) \}... ...}\sin (s-t)-\frac{1}{4}\sin (s+t), &s\le t . \end{array}\right. \end{eqnarray*}$$

The statement is

    ssr(y~t+cos(t)+sin(t), rk=lspline(t,type="linSinCos"))