||How to orthogonalize the random basis (
||Scaling of the
||Coordinates in the first dimension for which spline values are required|
||Coordinates in the second dimension for which spline values are required|
||Saves the design matrix to define the fixed terms (excluding the constant) for fitting the radial spline|
||Saves the design matrix to define the random terms for fitting the radial spline|
||Specifies the coordinates in the first dimension of the internal knots used to form the basis for the spline|
||Specifies the coordinates in the second dimension of the internal knots used to form the basis for the spline|
||Specifies the coordinates in the first dimension at which to predict|
||Specifies the coordinates in the second dimension at which to predict|
||Saves the design matrix for the fixed terms (excluding the constant) for the radial spline at the prediction points|
||Saves the design matrix for the random terms for the radial spline at the prediction points|
RADIALSPLINE generates the fixed and random terms required to fit a radial-spline surface as a linear mixed model, using
REML estimation of the smoothing parameter. The coordinates at which the spline is to be calculated are specified in two variates using
X2 parameters. The coordinates to be used as knots must be specified (in variates) using the
ORTHOGONALIZATION option specifies whether the components of the spline to be fitted as random terms should be made orthogonal to the components to be fitted as fixed. The default action (
ORTHOGONALIZATION=fixed) is to perform the orthogonalization, and this means that all of the polynomial trend associated with the fixed terms will be captured in the fixed part of the model. When
ORTHOGONALIZATION=none, some of this trend may be contained within the random terms.
The fixed and random components of the radial-spline terms are saved separately. The terms required to be fitted as fixed terms can be saved (in a matrix) using
XFIXED parameter. This matrix does not include the constant term as this is added by default as part of a mixed model. The terms to be fitted as random can be saved (in a matrix) using the
The random terms can be scaled so that, for a random spline matrix
TRACE(Z *+ T(Z)) = NROWS(Z)
This ensures that the average contribution of each component to the variance of an observation is equal to one. This improves interpretability of the spline variance components.
The radial-spline terms required for prediction can be saved using the
PXRANDOM parameters. The
PX2 parameters provide the coordinates at which predictions are to be made.
This procedure calculates a low-rank thin-plate spline in two dimensions, following an approach equivalent to that of Section 13.5 in Ruppert, Wand & Carroll (2003). The fixed terms comprise the input variates,
X2. For r knots at co-ordinates
tj = (tj1 … tjr)ʹ, j=1,2
and input variates
xi = (xi1 … xin)ʹ, i=1,2
the random basis functions are calculated via a function η (Green & Silverman, 1994), where
η(z) = z2 × log(z2) / (16 × π) for z > 0
= 0 for z = 0.
The r random basis functions then take the form
bl(x1,x2) = η( [ (x1–t1l)2 + (x2–t2l)2 ]1/2 ) for l = 1…r.
These columns can be concatenated into a n × r matrix Et. The corresponding r × r penalty matrix K has entries
K[i,j] = ζ( [ (t1i–t1j)2 + (t2i–t2j)2 ]1/2 ) for i,j = 1…r.
The matrix K can be transformed to full rank as
H = Cʹ K C
where matrix C contains the eigenvectors of X Xʹ (X = [1 x1 x2]) corresponding to zero eigenvalues, with corresponding transformation of the matrix functions as
Eu = Et C
This is translated to a set of independent random effects via post-multiplication by H-1/2.
The design matrices for use in prediction are calculated by evaluating the same set of basis functions at the predict points.
Input structures must not be restricted.
Green, P.J. & Silvername, B.W. (1994). Nonparametric Regression and Generalized Linear Models. Chapman & Hall, London.
Ruppert, D., Wand, M.P. & Carroll, R.J. (2003). Semiparametric Regression. Cambridge University Press, Cambridge.
CAPTION 'RADIALSPLINE example'; STYLE=meta " read the data " SPLOAD '%gendir%/data/Slatehall.gsh' " compute row and column coordinates " CALCULATE nrow,ncol = NLEVELS(fieldrow,fieldcolumn) & x1,x2 = fieldrow,fieldcolumn " predict points - on much finer grid " CALCULATE rx1 = !(10,11...100)/10 & rx2 = !(10,11...150)/10 & nr1,nr2 = NVALUES(rx1,rx2) VARIATE [VALUES=#nr2(#rx1)] px1 & [VALUES=(#rx2)#nr1] px2 GROUPS px1,px2; FACTOR=prow,pcol & px1; FACTOR=cprow; LIMITS=!(1.5...9.5) & px2; FACTOR=cpcol; LIMITS=!(1.5...14.5) " radial spline - default settings " RADIALSPLINE X1=x1; X2=x2; XFIXED=XYFa; XRANDOM=XYRa; PX1=px1; PX2=px2;\ X1KNOTS=!(8(1,2.5,4,5.5,7,8.5,10));\ X2KNOTS=!((1,3,5,7,9,11,13,15)7);\ PFIXED=PFa; PRANDOM=PRa " fit spline as spatial trend as part of mixed model " VCOMPONENTS [FIXED=XYFa+variety] RANDOM=XYRa; CONSTRAINT=positive REML [PRINT= model,components] yield " predict surface at observed points and interpolated points " VPREDICT [PRINT=description; PREDICTIONS=Tpa1] XYFa,XYRa;\ LEVELS=XYFa,XYRa; PARALLEL=XYFa & [PREDICTIONS=Tpa2] XYFa,XYRa; LEVELS=PFa,PRa; PARALLEL=XYFa " plot fitted spatial surfaces " CALCULATE np1,np2 = nval(unique(px1,px2)) MATRIX [ROWS=SORT(UNIQUE(x1)); COLUMNS=SORT(UNIQUE(x2))]\ Predicta1; VALUES=Tpa1 MATRIX [ROWS=SORT(UNIQUE(px1)); COLUMNS=SORT(UNIQUE(px2))]\ Predicta2; VALUES=Tpa2 FRAME [RESET=yes] 1...4 XAXIS [RESET=yes] 1...4 YAXIS [RESET=yes] 1...4 FFRAME [ROWS=2; COLUMNS=2] FRAME 1...4; XMLOWER=0.075; XMUPPER=0.075; YMLOWER=0.075; YMUPPER=0.075 DSTART [TITLE=\ 'Fitted surface at low & high resolution for radial spline with 20 knots'] DSHADE [WINDOW=1; KEY=3; GRID=*] Predicta1; PEN=!(4,7); LIMITS=!(9...22) DSHADE [WINDOW=2; KEY=4; SCREEN=keep; GRID=*] Predicta2; PEN=!(4,7);\ LIMITS=!(9...22) DFINISH