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=penspline` and 3 for `METHOD=pspline` |

`DIFFORDER` = scalar |
Differencing order for P-spline penalty; default 2 |

`X1LOWER` = scalar |
Specifies the lower boundary in the `X1` dimension when `KMETHOD=equal` ; default takes the minimum value in `X1` |

`X1UPPER` = scalar |
Specifies the upper boundary in the `X1` dimension when `KMETHOD=equal` ; default takes the maximum value in `X1` |

`X2LOWER` = scalar |
Specifies the lower boundary in the `X2` dimension when `KMETHOD=equal` ; default takes the minimum value in `X2` |

`X2UPPER` = scalar |
Specifies the upper boundary in the `X2` dimension when `KMETHOD=equal` ; default takes the maximum value in `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 |
---|---|

`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 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 2*d*+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