Re: Hand-calculation of predicted values with splines

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