Re: Hand-calculation of predicted values with splines

From: <arthur.gilmour_at_DPI.NSW.GOV.AU>
Date: Wed, 7 Nov 2007 12:30:30 +1100

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}
\special{ pdf:outline 2 << /Title (Cubic Splines) /Dest [@thispage /FitH
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
$$\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) }
$$\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
 returns $\tilde{\bmu}$ (rather than $\tilde{\bmu}_s$) in the solution
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}$ (=

The interpolated spline at $x$ between two knot points
$x_l$ and $x_h$ is
\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 +

where $\bmy= \X\hat{\bmb}+\Z\tilde{\bmu}$ is thus the vector of predicted
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

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

Received on Sun Nov 07 2007 - 12:30:30 EST

This webpage is part of the ASReml-l discussion list archives 2004-2010. More information on ASReml can be found at the VSN website. This discussion list is now deprecated - please use the VSN forum for discussion on ASReml. (These online archives were generated using the hypermail package.)