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 |
---|---|

`X` 2 = 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

*t _{j}* = (

*t*

_{j}_{1}…

*t*)ʹ,

_{jr}*j*=1,2

and input variates

*x _{i}* = (

*x*

_{i}_{1}…

*x*)ʹ,

_{in}*i*=1,2

the random basis functions are calculated via a function η (Green & Silverman, 1994), where

η(*z*) = *z*^{2} × log(*z*^{2}) / (16 × π) for *z* > 0

= 0 for *z* = 0.

The *r* random basis functions then take the form

*b _{l}*(

*x*

_{1},

*x*

_{2}) = η( [ (

*x*

_{1}–

*t*

_{1l})

^{2}+ (

*x*

_{2}–

*t*

_{2l})

^{2}]

^{1/2}) for

*l*= 1…

*r*.

These columns can be concatenated into a *n* × *r* matrix *E _{t}*. The corresponding

*r*×

*r*penalty matrix

*K*has entries

*K*[*i*,*j*] = ζ( [ (*t*_{1i}–*t*_{1j})^{2} + (*t*_{2i}–*t*_{2j})^{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 *x*_{1} *x*_{2}]) corresponding to zero eigenvalues, with corresponding transformation of the matrix functions as

*E _{u}* =

*E*

_{t}*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