Re: Hand-calculation of predicted values with splines

Dear Kevin,

It will be very hard to do the prediction without using the spl()
coefficients used by ASReml.
These are returned in the .res file with standalone ASReml though I can't
see them in the returned asreml() object..

The coding is based on Verbyla et al Applied Stats 1999 48:269-311 and
described in the following tex code (from the Old ASReml reference manaul)

\section{Cubic\index{Cubic, see splines}\label{DET:SPLINE}
 splines\index{splines}}
\special{ pdf:outline 2 << /Title (Cubic Splines) /Dest [@thispage /FitH
@ypos]>>}
A recent advance in \AAR\ is the implementation of
cubic smoothing splines as described by Verbyla {\it et al.}\footnote{
Verbyla, A.P., Cullis. B.R., Kenward, M.G. and Welham, S.J. (1999).
Smoothing splines in the analysis of designed experiments
and longitudinal data. {\it Applied Statistics} {\bf 48:} 269-311. }
A cubic spline is a smooth curve over an interval $[a,b]$
formed by linking segments of cubic polynomials at {\sl knot}
points such that the whole curve and it's first and second differentials
are continuous over the interval. See Green \& Silverman\footnote{
Green, P.J. and Silverman, B.W. (1994), Nonparametric Regression and
Generalized Linear models. Chapman \& Hall, 182pp}
for a general treatment of splines.
Verbyla {\em et al.} show that a cubic spline can be fitted as a
linear model of the form
$$\bmy=\X\bmb + \Z_s\bmu_s + \bmeta$$
where $\X\bmb$ fits a linear regression
(intercept and slope), $\Z_s=\Delta(\Delta^T\Delta)^{-1},$ both
$\bmu\;\sim$ N$(0, \sigma^2_s\G)$ and $\bmeta\;\sim$ N$(0, \sigma^2\I)$
are vectors of random effects, and
$\G$ and $\Delta$ are functions of the first differences between the
$N$ ordered
knot points. \AAR\ transforms this model to use $\Z=k^{-1.5}\Z_s\G^{1/2}$
so that the fitted effects are
$\tilde{\bmu}=\G^{-1/2}\tilde{\bmu}_sk^{1.5}$ where
$k$ is the averge distance between knots (range/$(N-1)$)
and $k$ standardizes $\Z$ so the variance component is not poorly scaled.

In \AAR, $\G^{1/2}$ is formed as a cholesky factor (unlike the
Genstat and Splus implementations which use an eigen factorization).

For $h_i=x_{i+1}-x_i$ being the differences between five ordered knot
points,
$$\Delta^T =
      {\left( \begin{array}{ccccc}
           1/h_1 & -( 1/h_1 + 1/h_2) & 1/h_2 & 0 & 0 \\
      0 & 1/h_2 & -( 1/h_2 + 1/h_3) & 1/h_3 & 0 \\
  0 & 0 & 1/h_3 & -( 1/h_3 + 1/h_4) & 1/h_4
  \end{array}\right) }
$$
and
$$\G =
      {\left( \begin{array}{ccc}
       {h_1 +h_2\over 3} & {h_2\over 6} &0 \\*[0.15cm]
      {h_2\over 6} & {h_2 +h_3\over 3} & {h_3\over 6} \\*[0.15cm]
  0 & {h_3\over 6} & {h_3 +h_4\over 3}
  \end{array}\right). }
$$
\AAR\ performs the scaling by actually using $h_i=(x_{i+1}-x_i)/k$ in
these expressions
where $k$ is $(x_N-x_1)/(N-1)$.
Note that
\AAR\
 returns $\tilde{\bmu}$ (rather than $\tilde{\bmu}_s$) in the solution
vector.
It returns the knot points ($\bmx$) and the
fitted spline curvature ($\Z\tilde{\bmu}=\Z_s\tilde{\bmu}_s$) values
above mini plots of the curvature against the knot points
% and $\bmx$ (correpsonding to $\tilde{\bmu}$ and $\Z$) is presented
%above the spline plots
in the {\tt .asr} file. $\Z$ and $\bmx$ are
also returned in the {\tt .res} file. A {\tt .pin} file can be
used to predict the full spline $\X\hat{\bmb}+\Z\tilde{\bmu}$ (=
$\X\hat{\bmb}+\Z_s\tilde{\bmu}_s).$

The interpolated spline at $x$ between two knot points
$x_l$ and $x_h$ is
\begin{eqnarray*}
\frac{ (x-x_l)y_h + (x_h-x)y_l}{(x_h-x_l)}
- \hfill\hbox{\hspace*{6cm}}\\\frac{(x-x_l)(x_h-x)}{ 6}\left[
\left(1+\frac{x-x_l}{x_h-x_l}\right)d_h +
\left(1+\frac{x_h-x}{x_h-x_l}\right)d_l\right]
\end{eqnarray*}

where $\bmy= \X\hat{\bmb}+\Z\tilde{\bmu}$ is thus the vector of predicted
points,
and for
 $$\bmdel = \G^{-1}\tilde{\bmu}_s
= \G^{-1}\Delta^T\Delta(\Delta^T\Delta)^{-1}\tilde{\bmu}_s
= \G^{-1}\Delta^T\Z_s\tilde{\bmu}_s,$$
$\bmd= (0\ \ \bmdel'\ \ 0)'$ is vector of second derivatives at the knot
points.

Extract from .res file showing spl coeffs at knot pouints, and then
interpolated at predict points

 spl(x,7) Z matrix for specified KNOT points: Scale factor 0.0061
    118.00 1.3072 1.4651 0.9457 0.4209 0.1926
    484.00 -0.9158 -0.3924 -0.0087 0.0500 0.0352
    664.00 -0.8386 -1.3059 -0.4781 -0.1323 -0.0421
   1004.00 -0.3283 -0.9049 -1.3647 -0.4768 -0.1883
   1231.00 0.0124 -0.1593 -0.5303 -0.7068 -0.2859
   1372.00 0.2240 0.3038 0.1852 -0.1233 -0.3466
   1582.00 0.5391 0.9936 1.2509 0.9684 0.6351
spl(x,7) has 5 levels. Z matrix standardized for scale.
       118.00000 1.30719 1.46505 0.94568 0.42087 0.19262
       150.00000 1.04362 1.31451 0.85826 0.38953 0.17864
       200.00000 0.63940 1.07798 0.72211 0.34043 0.15682
etc

Arthur
><><><><><><><><><><><><><><><><><><><><><><><>
Received on Sun Nov 07 2007 - 12:30:30 EST