From: Butler, David <David.Butler_at_DPI.QLD.GOV.AU>

Date: Wed, 7 Nov 2007 14:48:35 +1000

Date: Wed, 7 Nov 2007 14:48:35 +1000

Kevin,
Further to Arthur's comments, there used to be a helper function that
returned the design matrix etc for a spline term. I have let this
'lapse' since predict() but I can include it in future distributions if
required.
Dave.
________________________________
From: ASReml users discussion group [mailto:ASREML-L_at_AGRIC.NSW.GOV.AU]
On Behalf Of arthur.gilmour_at_DPI.NSW.GOV.AU
Sent: Wednesday, 7 November 2007 11:31 AM
To: ASREML-L_at_AGRIC.NSW.GOV.AU
Subject: 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
><><><><><><><><><><><><><><><><><><><><><><><>
********************************DISCLAIMER****************************
The information contained in the above e-mail message or messages
(which includes any attachments) is confidential and may be legally
privileged. It is intended only for the use of the person or entity
to which it is addressed. If you are not the addressee any form of
disclosure, copying, modification, distribution or any action taken
or omitted in reliance on the information is unauthorised. Opinions
contained in the message(s) do not necessarily reflect the opinions
of the Queensland Government and its authorities. If you received
this communication in error, please notify the sender immediately and
delete it from your computer system network.

Received on Sun Nov 07 2007 - 14:48:35 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.)