Calculates design matrices to fit a radial-spline surface as a linear mixed model (S.J. Welham & D.B. Baird).
Options
ORTHOGONALIZATION = string token |
How to orthogonalize the random basis (fixed, none); default fixe |
|---|---|
SCALING = scalar |
Scaling of the XRANDOM terms (automatic, none); default auto |
Parameters
X1 = variates or factors |
Coordinates in the first dimension for which spline values are required |
|---|---|
X2 = variates or factors |
Coordinates in the second dimension for which spline values are required |
XFIXED = matrices |
Saves the design matrix to define the fixed terms (excluding the constant) for fitting the radial spline |
XRANDOM = matrices |
Saves the design matrix to define the random terms for fitting the radial spline |
X1KNOTS = variates |
Specifies the coordinates in the first dimension of the internal knots used to form the basis for the spline |
X2KNOTS = variates |
Specifies the coordinates in the second dimension of the internal knots used to form the basis for the spline |
PX1 = variates |
Specifies the coordinates in the first dimension at which to predict |
PX2 = variates |
Specifies the coordinates in the second dimension at which to predict |
PFIXED = matrices |
Saves the design matrix for the fixed terms (excluding the constant) for the radial spline at the prediction points |
PRANDOM = matrices |
Saves the design matrix for the random terms for the radial spline at the prediction points |
Description
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 X1 and X2 parameters. The coordinates to be used as knots must be specified (in variates) using the X1KNOTS and X2KNOTS parameters.
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 XRANDOM parameter.
The random terms can be scaled so that, for a random spline matrix Z,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 PXFIXED and PXRANDOM parameters. The PX1 and PX2 parameters provide the coordinates at which predictions are to be made.
Options: ORTHOGONALIZATION, SCALING.
Parameters: X1, X2, XFIXED, XRANDOM, X1KNOTS, X2KNOTS, PX1, PX2, PFIXED, PRANDOM.
Method
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, X1 and 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.
Action with RESTRICT
Input structures must not be restricted.
References
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.
See also
Directives: VCOMPONENTS, REML.
Procedures: NCSPLINE, PENSPLINE, PSPLINE, SPLINE, TENSORSPLINE.
Function: SSPLINE.
Commands for: Calculations and manipulation, Regression analysis, REML analysis of linear mixed models.
Example
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