Language Reference

TPSPLINE Call

CALL TPSPLINE (fitted, coeff, adiag, gcv, x, y <*>, lambda );

The TSPLINE subroutine fits a thin-plate smoothing spline (TPSS) to data. The generalized cross validation (GCV) function is used to select the smoothing parameter.

The TPSPLINE subroutine returns the following values:

fitted: is an $n\times 1$ vector of fitted values of the TPSS fit evaluated at the design points x. The n is the number of observations. The final TPSS fit depends on the optional lambda.
coeff: is a vector of spline coefficients. The vector contains the coefficients for basis functions in the null space and the representer of evaluation functions at unique design points. (see Wahba (1990) for more detail on reproducing kernel Hilbert space and representer of evaluation functions.) The length of coeff vector depends on the number of unique design points and the number of variables in the spline model. In general, let nuobs and k be the number of unique rows and the number of columns of x respectively. The length of coeff equals to $k+\mbox{nuobs}+1$ . The coeff vector can be used as an input to the TPSPLNEV subroutine to evaluate the resulting TPSS fit at new data points.
adiag: is an $n \times 1$ vector of diagonal elements of the "hat" matrix. See the "Details" section.
gcv: If lambda is not specified, then gcv is the minimum value of the GCV function. If lambda is specified, then gcv is a vector (or scalar if lambda is a scalar) of GCV values evaluated at the lambda points. It provides you with both the ability to study the GCV curves by plotting gcv against lambda and the chance to identify a possible local minimum.

The input arguments to the TPSPLINE subroutine are as follows:

x: is an $n \times k$ matrix of design points on which the TPSS is to be fit. The k is the number of variables in the spline model. The columns of x need to be linearly independent and contain no constant column.
y: is the $n \times 1$ vector of observations.
lambda: is a optional $q \times 1$ vector that contains $\lambda$ values in $log_{10}(n\lambda )$ scale. If lambda is not specified (or lambda is specified and $q > 1$ ) the GCV function is used to choose the "best" $\lambda$ and the returning fitted values are based on the $\lambda$ that minimizes the GCV function. If lambda is specified and $q=1$ , no minimization of the GCV function is involved and the fitted, coeff and adiag values are all based on the TPSS fit that uses this particular lambda.

Aside from the values returned, the TPSPLINE subroutine also prints other useful information such as the number of unique observations, the dimensions of the null space, the number of parameters in the model, a GCV estimate of $\lambda$ , the smoothing penalty, the residual sum of square, the trace of $(I-A(\lambda ))$ , an estimate of $\sigma ^2$ , and the sum of squares for replication.

No missing values are accepted within the input arguments. Also, you should use caution if you want to specify small lambda values. Since the true $\lambda =(10^{\log _{10} \mi{lambda}})/n$ , a very small value for lambda can cause $\lambda$ to be smaller than the magnitude of machine error and usually the returned gcv values from such a $\lambda$ cannot be trusted. Finally, when using TPSPLINE be aware that TPSS is a computationally intensive method. Therefore a large data set (that is, a large number of unique design points) will take a lot of computer memory and time.

For convenience, the TPSS method is illustrated with a two-dimensional independent variable $\bX =(x^1,x^2)$ . More details can be found in Wahba (1990), or in Bates et al. (1987).

Assume that the data are from the model

$y_ i = f({x}_ i) + \epsilon _ i,$

where $({x}_ i, y_ i), i=1,\ldots ,n$ are the observations. The function f is unknown and you assume that it is reasonably smooth. The error terms $\epsilon _ i, i=1,\ldots ,n$ are independent zero-mean random variables.

You measure the smoothness of f by the integral over the entire plane of the square of the partial derivatives of f of total order 2, that is

$J_2(f) = \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } \left[ \frac{\partial ^2 f}{\partial {x_1}^2} \right]^2 + 2 \left[ \frac{\partial ^2 f}{\partial {x_1} \partial {x_2}} \right]^2 + \left[ \frac{\partial ^2 f}{\partial {x_2}^2} \right]^2 ~ dx_1 dx_2$

Using this as a smoothness penalty, the thin-plate smoothing spline estimate $f_{\lambda }$ of f is the minimizer of

$S_\lambda (f) = \frac{1}{n} \sum ^ n_{i=1} (y_ i-f(x_ i))^2 + \lambda J_2(f).$

Duchon (1976) derived that the minimizer $f_\lambda$ can be represented as

$f_\lambda (x) = \sum _{i=1}^{3} \beta _ i \phi _ i(x) + \sum _{i=1}^ n \delta _ i E_2(x-x_ i),$

where $(\phi _1(x),\phi _2(x),\phi _3(x)) = (1,x^1,x^2)$ and $E_2(s)=\frac{1}{2^3 \pi } \| s\| ^2ln(\| s\| )$ .

Let matrix $\mb{K}$ have entries $(\mb{K})_{ij} = E_2(x_ i-x_ j)$ and matrix $\mb{T}$ have entries $(\mb{T})_{ij} = \phi _ j(x_ i)$ . Then the minimization problem can be rewritten as finding coefficients $\beta$ and $\delta$ to minimize

$S_\lambda (\beta , \delta ) = \frac{1}{n} \| y - \bT \beta - \bK \delta \| ^2 + \lambda \delta ^ T \bK \delta$

The final TPSS fits can be viewed as a type of generalized ridge regression estimator. The $\lambda$ is called the smoothing parameter, which controls the balance between the goodness of fit and the smoothness of the final estimate. The smoothing parameter can be chosen by minimizing the generalized cross validation function (GCV). If you write

$\hat{y} = \bA (\lambda ) y$

and call the $\bA (\lambda )$ as the " $hat$ " matrix, the GCV function $V(\lambda )$ is defined as

$V(\lambda ) = \frac{(1/n) \| (\bI - \bA (\lambda ) y \| ^2}{[(1/n) \mbox{tr}(\bI - \bA (\lambda ))]^2}$

The returned values from this function call provide the $\hat{y}$ as fitted, the $(\beta ,\delta )$ as coeff, and $diag(A(\lambda ))$ as adiag.

To evaluate the TPSS fit $f_\lambda (x)$ at new data points, you can use the TPSPLNEV call.

Suppose $\bX ^{\mbox{new}}$ , a $m \times k$ matrix, contains the m new data points at which you want to evaluate $f_\lambda$ . Let $(\mb{T}^{\mbox{new}}_{ij}) = \phi _ j(x^{\mbox{new}}_ i)$ and $(\mb{K}^{\mbox{new}}_{ij})=E_2(x^{\mbox{new}}_ i-x_ j)$ be the (i, j) elements of $\mb{T}^{\mbox{new}}$ and $\mb{K}^{\mbox{new}}$ respectively. The prediction at new data points $\bX ^{\mbox{new}}$ is

$y_{\mbox{pred}}=T^{\mbox{new}} \beta + K^{\mbox{new}} \delta$

Therefore, the $y_{\mbox{pred}}$ can be easily evaluated by using the coefficient $(\beta ,\delta )$ obtained from the TPSPLINE call.

An example is given in the documentation for the TPSPLNEV call .