Calculates design matrices to fit a tensor-spline surface as a linear mixed model (S.J. Welham & P.H.C. Eilers).
Options
METHOD = string token |
Type of spline to use to construct the basis (pspline, penalizedspline); default pspl |
|---|---|
PENALTYMETHOD = string token |
Which tensor-spline penalty to use (isotropic, semiconstrained, unconstrained); default unco |
NX1SEGMENTS = scalar |
Specifies the number of segments between boundaries in the X1 dimension; default * obtains a value automatically |
NX2SEGMENTS = scalar |
Specifies the number of segments between boundaries in the X2 dimension; default * obtains a value automatically |
DEGREE = scalar |
Degree of polynomial used to form the underlying spline basis functions; default 1 for METHOD=pena and 3 for METHOD=pspl |
DIFFORDER = scalar |
Differencing order for P-spline penalty; default 2 |
X1LOWER = scalar |
Specifies the lower boundary in the X1 dimension; default takes the minimum value of X1 |
X1UPPER = scalar |
Specifies the upper boundary in the X1 dimension; default takes the minimum value of X1 |
X2LOWER = scalar |
Specifies the lower boundary in the X2 dimension; default takes the minimum value of X2 |
X2UPPER = scalar |
Specifies the upper boundary in the X2 dimension; default takes the minimum value of X2 |
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 tensor spline |
XRANDOM = pointers |
Saves the design matrices to define the random terms for fitting the tensor spline |
X1KNOTS = variates |
Saves the coordinates in the first dimension of the internal knots used to form the basis for the spline |
X2KNOTS = variates |
Saves 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 tensor spline at the prediction points |
PRANDOM = pointers |
Saves the design matrices for the random terms for the tensor spline at the prediction points |
Description
TENSORSPLINE generates the fixed and random terms required to fit a tensor-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 full range of the spline in the X1 dimension can be defined using the X1LOWER and X1UPPER options; by default the lower limit is equal to the minimum values of X1 and the upper limit is equal to the maximum value. The range in the X2 dimension is defined similarly by the X2LOWER and X2UPPER options. In each dimension, the region between these bounds is divided into a number of equal segments, specified by the NX1SEGMENTS and NX2SEGMENTS options. The boundaries of these segments provide the set of knots used to form the spline terms in each dimension, and can be saved in variates by the X1KNOTS and X2KNOTS parameters. If the number of segments is unset in either dimension, the number is determined automatically as
min([p/4], 35) + 1
(Ruppert 2002) where p is the number of unique values of the variate (X1 or X2), and [r] denotes the integer part of the number r.
The METHOD option specifies whether to use P-splines (the default) or penalized splines to construct the basis. The degree of polynomial used to form the underlying spline basis functions is specified by the DEGREE option. This has a default of 3 for P-spline models ,and 1 for penalized spline models.
The DIFFORDER option specifies the differencing order to be used with P-spline models. This determines the strength of the penalty (for a given smoothness parameter). The default is to use second-order differencing. For a P-spline model, the underlying fixed polynomial in each dimension has degree d equal to DIFFORDER-1. For a penalized spline model, the underlying fixed polynomial in each dimension has degree d equal to the value specified by the DEGREE option.
The tensor-spline basis is constructed via interactions of the one-dimensional spline bases, as detailed in the methods section.
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. Experience suggests that ORTHOGONALIZATION=none can lead to computational instability when the model is fitted using REML, especially for more complex models (d>1).
The fixed and random components of the tensor-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. For P-spline models with DIFFORDER=1, this is a null term and no matrix is saved.
The terms to be fitted as random can be saved (in a matrix) using the XRANDOM parameter. The terms to be fitted as random can be saved using the XRANDOM parameter. This is a saves a pointer, containing a number of matrices that depends on the setting of the PENALTYMETHOD option. The components of the random terms consist of interactions between:
● the underlying fixed polynomials in the first dimension with the random basis functions in the second dimension (d+1 terms)
● the random basis functions in the first dimension with the underlying fixed polynomials in the second dimension (d+1 terms); and
● the interaction of the two sets of random basis functions from each dimension (1 term).
The default is an unconstrained penalty, where a separate smoothing parameter is allowed for each term. In this case, the XRANDOM pointer has 2d+3 matrices, one for each term. For a semi-constrained penalty (PENALTYMETHOD=semiconstrained), the same smoothing parameter is imposed across the interaction of polynomials in the first dimension with random terms in the second, and for the interaction of random terms in the first dimension with polynomials in the second dimension. This is implemented by combining terms, so the XRANDOM pointer then has 3 matrices. For an isotropic penalty, which uses a single common penalty (PENALTY=isotropic), the terms are combined into a single matrix, so the pointer XRANDOM has a single element.
The random terms can be scaled so that, for each component matrix Z in XRANDOM,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.
For PENALTYMETHOD=unconstrained, the fitted model is invariant to the scale of X1 and X2. This is not the case for the semi-constrained and isotropic penalties. Full (automatic) scaling is imposed by default, but this can be avoided by setting option SCALING=none. An intermediate option is available (SCALING=standardize) where the polynomial components of the random terms are standardized before being added into the random terms.
The tensor-spline terms required for prediction can be saved using the PXFIXED (for P-spline models, provided d is greater than 1) and PXRANDOM parameters. The PX1 and PX2 parameters provide the coordinates at which predictions are to be made.
Options: METHOD, PENALTYMETHOD, NX1SEGMENTS, NX2SEGMENTS, DEGREE, DIFFORDER, X1LOWER, X1UPPER, X2LOWER, X2UPPER, ORTHOGONALIZATION, SCALING.
Parameters: X1, X2, XFIXED, XRANDOM, X1KNOTS, X2KNOTS, PX1, PX2, PFIXED, PRANDOM.
Method
For each dimension, the appropriate one-dimensional P-spline or penalized spline basis functions are calculated (see procedures PSPLINE and PENSPLINE for details). Where the degree of the underlying fixed polynomials is equal to d (as explained under Description), there are d+1 polynomial terms for each dimension. To explain the construction of the tensor-splines, we represent these polynomials as vectors
X1[0...d] = X1**(0...d)
and
X2[0...d] = X2**(0...d)
The fixed part of the tensor-spline comprises all (d+1)2 products of the form
X1X2[i][j] = X1[i]*X2[j]
for i=0…d, j=0…d. These vectors are copied into the columns of matrix XFIXED, with the exception of X1X2[0][0] which is the constant term, and added into the model automatically by the VCOMPONENTS directive. We represent the random parts of the two one-dimensional spline bases as matrices Z1 and Z2, with n1 and n2 columns, respectively. The random parts of the tensorspline are then created as interactions between the polynomial and spline terms, calculated as follows (i,j=0…d):
| Term | Matrix calculation |
X1[i].Z2 |
KRONECKER(X1[i]; ROW1(n2)) * Z2 |
Z1.X2[j] |
Z1 * KRONECKER(ROW1(n1); X2[j]) |
Z1.Z2 |
KRONECKER(Z1; ROW1(n2)) * KRONECKER(ROW1(n1); Z2) |
When SCALING=automatic or SCALING=standardize, the terms X1[] and X2[] are standardized before calculation of the random components of the model. When SCALING=none, no transformation is used.
For PENALTYMETHOD=unconstrained, the 2d+3 terms generated are kept separate, each with a separate smoothing parameter, which is estimated via the variance component associated with each random term. Because each of the terms has a separate variance component, this model is invariant to rescaling of the X1 and X2 variates. For PENALTY=semiconstrained, the matrices X1[i].Z2 (i=0…d) are concatenated into a single matrix and fitted as a single random term, with a common smoothing parameter. Similarly the matrices Z1.X2[j] (j=0…d) form a single term, and Z1.Z2 forms the third random term in the model. When SCALING=automatic or SCALING=standardized, this model is invariant to rescaling of X1 and X2, as all of the polynomial terms are standardized before formation of the matrices – this corresponds to a particular (and arbitrary) choice of scaling. When SCALING=none, the fitted model depends on the scale of the input variates. For PENALTYMETHOD=isotropic, all the random terms are concatenated into a single matrix and fitted with a common smmothing parameter. Again, with SCALING=none, the fit will depend on the scale of the input variates.
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.
Reference
Ruppert, D. (2002). Selecting the number of knots for penalised splines. Computational & Graphical Statistics, 11, 735-757.
See also
Directives: VCOMPONENTS, REML.
Procedures: NCSPLINE, PENSPLINE, PSPLINE, RADIALSPLINE, SPLINE.
Function: SSPLINE.
Commands for: Calculations and manipulation, Regression analysis, REML analysis of linear mixed models.
Example
CAPTION 'TENSORSPLINE 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
" tensortspline based on Pspline - default settings "
TENSORSPLINE X1=x1; X2=x2; XFIXED=XYFa; XRANDOM=XYRa; PX1=px1; PX2=px2;\
PFIXED=PFa; PRANDOM=PRa; X1KNOTS=xk; X2KNOTS=yk
" 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 P-splines with k=3, d=2']
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
" tensor spline based on linear penalized spline - default knots "
TENSORSPLINE [METHOD=penalized] X1=x1; X2=x2; XFIXED=XYFb; XRANDOM=XYRb;\
PX1=px1; PX2=px2; PFIXED=PFb; PRANDOM=PRb
" fit spline as spatial trend as part of mixed model "
VCOMPONENTS [FIXED=XYFb+variety] RANDOM=XYRb[]; CONSTRAINT=positive
REML [PRINT=model,components] yield
" predict surface at observed points and interpolated points "
VPREDICT [PRINT=description; PREDICTIONS=Tpb1] XYFb,XYRb[];\
LEVELS=XYFb,XYRb[]; PARALLEL=XYFb
& [PREDICTIONS=Tpb2] XYFb,XYRb[]; LEVELS=PFb,PRb[]; PARALLEL=XYFb
" plot fitted spatial surfaces "
MATRIX [ROWS=SORT(UNIQUE(x1)); COLUMNS=SORT(UNIQUE(x2))]\
Predictb1; VALUES=Tpb1
MATRIX [ROWS=SORT(UNIQUE(px1)); COLUMNS=SORT(UNIQUE(px2))]\
Predictb2; VALUES=Tpb2
DSTART [TITLE=\
'Fitted surface at low & high resolution for linear penalized splines']
DSHADE [WINDOW=1; KEY=3; GRID=*] Predictb1; PEN=!(4,7); LIMITS=!(9...22)
DSHADE [WINDOW=2; KEY=4; SCREEN=keep; GRID=*] Predictb2; PEN=!(4,7);\
LIMITS=!(9...22)
DFINISH
" tensor spline based on Pspline - knots at data "
TENSORSPLINE [NX1SEGMENTS=9; NX2SEGMENTS=14] X1=x1; X2=x2; XFIXED=XYFc;\
XRANDOM=XYRc; PX1=px1; PX2=px2; PFIXED=PFc; PRANDOM=PRc;\
X1KNOTS=xk; X2KNOTS=yk
" Fit spline as spatial trend as part of mixed model "
VCOMPONENTS [FIXED=XYFc+variety] RANDOM=XYRc[]; CONSTRAINT=positive
REML [PRINT=model,components] yield
" predict surface at observed points and interpolated points "
VPREDICT [PRINT=description; PREDICTIONS=Tpc1] XYFc,XYRc[];\
LEVELS=XYFc,XYRc[]; PARALLEL=XYFc
& [PREDICTIONS=Tpc2] XYFc,XYRc[]; LEVELS=PFc,PRc[]; PARALLEL=XYFc
" plot fitted spatial surfaces "
MATRIX [ROWS=SORT(UNIQUE(x1)); COLUMNS=SORT(UNIQUE(x2))]\
Predictc1; VALUES=Tpc1
MATRIX [ROWS=SORT(UNIQUE(px1)); COLUMNS=SORT(UNIQUE(px2))]\
Predictc2; VALUES=Tpc2
DSTART [TITLE=\
'Fitted surface at low/high res for P-splines with k=3, d=2 & knots at data']
DSHADE [WINDOW=1; KEY=3; GRID=*] Predictc1; PEN=!(4,7); LIMITS=!(9...22)
DSHADE [WINDOW=2; KEY=4; SCREEN=keep; GRID=*] Predictc2; PEN=!(4,7);\
LIMITS=!(9...22)
DFINISH