Out-of-Sample Prediction in Multidimensional P-Spline Models

Abstract

The prediction of out-of-sample values is an interesting problem in any regression model. In the context of penalized smoothing using a mixed-model reparameterization, a general framework has been proposed for predicting in additive models but without interaction terms. The aim of this paper is to generalize this work, extending the methodology proposed in the multidimensional case, to models that include interaction terms, i.e., when prediction is carried out in a multidimensional setting. Our method fits the data, predicts new observations at the same time, and uses constraints to ensure a consistent fit or impose further restrictions on predictions. We have also developed this method for the so-called smooth-ANOVA model, which allows us to include interaction terms that can be decomposed into the sum of several smooth functions. We also develop this methodology for the so-called smooth-ANOVA models, which allow us to include interaction terms that can be decomposed as a sum of several smooth functions. To illustrate the method, two real data sets were used, one for predicting the mortality of the U.S. population in a logarithmic scale, and the other for predicting the aboveground biomass of Populus trees as a smooth function of height and diameter. We examine the performance of interaction and the smooth-ANOVA model through simulation studies.

1. Introduction

For more than twenty years, use of penalized splines with B-splines bases (or simply P-splines) [1] has become one of the most popular smoothing techniques in many application fields, see [2] for a detailed review, and also from a theoretical perspective [3,4]. The possibility to reformulate a P-spline as a maximum likelihood estimator and best linear unbiased predictors (BLUPs) in a mixed model framework have led to a wider generalization of the method to longitudinal data, spatial smoothing, correlated errors and Bayesian approaches [5,6]. In a multidimensional setting, P-splines have also been used to model mortality life tables or spatio-temporal data [7,8,9]. This method is not only used to fit data, but also to obtain out-of-sample predictions. For instance, ref. [10] used P-splines to forecast mortality life tables and [11,12] also used them in a three-dimensional spatio-temporal model for predicting the risk of cancer death in the next few years. However, there are some important issues that still need to be addressed, for example, simultaneous prediction in more than one dimension, or in models that include smooth main effects and interactions terms.

One of the main contributions of this paper is an in depth study of the method proposed in [10], presenting the relationship between the regression coefficients to determine the fit and the coefficients to determine the prediction and extending the method to the mixed model framework. For the one-dimensional case, this research has been completed in [13]. In this work, we study the properties of the method proposed by [10] and show that the fit changes when the fit and the prediction are obtained simultaneously. The difference between the fits can sometimes be small, but as we shown in the examples, the impact on the prediction can be very important.

We also show that for specific cases where only one of the two covariates is extended, the extended penalty matrix can be modified to maintain the fit. However, when the two covariates are extended, the penalty matrix cannot be modified because the matrix involved in obtaining the estimated parameters is singular.

As a general solution to maintain the fit, we propose imposing constraints on the P-splines regression coefficients. We achieve this by using the standard Lagrangian formulation of the least squares minimization problem [14]. In this article, we mainly use restrictions to ensure that the fit must be maintained, but our approach can be used for other purposes—for example, combining known information or expert opinions (or upper/lower bounds) about the value of the out-of-sample predicted observations.

In other contexts, such as spatial or spatio-temporal modelling, we might be interested in prediction at new locations, this would involve out-of-sample prediction for two covariates (latitude and longitude), and so, we extend the proposal of [10] to predict when the two covariates are extended and to predict with smooth-ANOVA models. These models have been formulated as mixed models by [8,15] and they allow us to decompose multidimensional smooth functions in terms of additive terms and interactions.

The rest of the paper is organized as follows: in Section 2, we give a brief introduction to multidimensional P-splines, present a general approach for out-of-sample prediction with and without constraints and show the properties satisfied, under certain conditions, by the coefficients that determine the prediction. We apply the methodology proposed to the analysis and forecast of US mortality data. In Section 3, we extend the results to the case of smooth-ANOVA models and, in particular, to the case in which P-splines are reformulated as mixed models. We close the section with the analysis and prediction of aboveground biomass data. Section 4 is dedicated to a simulation study to compare the methods proposed under two different scenarios. In Section 5, we present some conclusions of this work.

2. Multidimensional P-Splines

Consider a multidimensional smoothing problem of a response variable $\mathit{y}$ with two covariates, $\mathit{x}$ and $\mathit{z}$.

$$\mathit{y} = f(\mathit{z},\mathit{x}) + \mathit{ϵ},\mathit{ϵ}\sim \mathcal{N}(0,\mathit{R})$$

(1)

where $\mathit{ϵ}$ is a Gaussian i.i.d error term with zero mean and variance covariance matrix $\mathit{R} = {\sigma }^2\mathit{W}$, where $\mathit{W}$ is a diagonal matrix when data are not correlated. The function $f(·)$ is a $2–$multidimensional smooth function that depends on the two explanatory variables $\mathit{z} = {(z_1,\dots ,z_{n_{\mathit{z}}})}^{\prime }$ and $\mathit{x} = {(x_1,\dots ,x_{n_{\mathit{x}}})}^{\prime }$, and each of them have lengths $n_{\mathit{z}}$ and $n_{\mathit{x}}$, respectively. Although, here, we assume i.i.d. for simplicity, the results can be easily extended to a more general variance–covariance matrix $\mathit{R}$. The model in (1), can be represented in terms of basis functions with

$$f(\mathit{z},\mathit{x}) = \mathit{B}\mathit{\theta },$$

(2)

with $\mathit{B}$, a B-spline regression basis, and $\mathit{\theta }$, the vector of coefficients. If we consider data in a regular array structure, a smooth multidimensional surface is constructed via the Kronecker product of the marginal B-spline basis of each covariate, the basis in model (2) is,

$$\mathit{B} = {\mathit{B}}_{\mathit{x}} \otimes {\mathit{B}}_{\mathit{z}},$$

(3)

where ⊗ is the Kronecker product of two matrices, and ${\mathit{B}}_{\mathit{x}} = \mathit{B}(\mathit{x})$ and ${\mathit{B}}_{\mathit{z}}=\mathit{B}(\mathit{z})$, of dimensions $n_{\mathit{x}} \times c_{\mathit{x}}$ and $n_{\mathit{z}} \times c_{\mathit{z}}$, are the marginal B-spline bases for $\mathit{x}$ and $\mathit{z}$, respectively. Then, the dimension of (3) is $n_{\mathit{x}}{n}_{\mathit{z}} \times c_{\mathit{x}}{c}_{\mathit{z}}$. On the other hand, if we consider scattered data, the basis is constructed by the tensor product of the marginal B-spline basis defined in [7] as the so-called box-product, denoted by symbol □:

$$\mathit{B} = {\mathit{B}}_{\mathit{x}}\square {\mathit{B}}_{\mathit{z}} = ({\mathit{B}}_{\mathit{x}} \otimes {\mathbf{1}}_{c_{\mathit{z}}}^{\prime })\odot ({\mathbf{1}}_{c_{\mathit{x}}}^{\prime } \otimes {\mathit{B}}_{\mathit{z}}),$$

where the operator ⊙ is the element-wise matrix product and ${\mathbf{1}}_{c_{\mathit{z}}}$ and ${\mathbf{1}}_{c_{\mathit{x}}}$ are column vectors of ones of lengths $c_{\mathit{z}}$ and $c_{\mathit{x}}$.

In both cases (array structured or scattered data), the vector of coefficients $\mathit{\theta }$ can be arranged into a $c_{\mathit{z}} \times c_{\mathit{x}}$ matrix $\mathsf{\Theta }$, that is

$$\mathsf{\Theta } = \left[\begin{array}{cccc}{\mathit{\theta }}_{11}& {\mathit{\theta }}_{12}& \cdots & {\mathit{\theta }}_{1c_{\mathit{x}}}\\ {\mathit{\theta }}_{21}& {\mathit{\theta }}_{22}& \cdots & {\mathit{\theta }}_{2c_{\mathit{x}}}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathit{\theta }}_{c_{\mathit{z}}1}& {\mathit{\theta }}_{c_{\mathit{z}}2}& \cdots & {\mathit{\theta }}_{c_{\mathit{z}}{c}_{\mathit{x}}}\end{array}\right],$$

(4)

and so the two-dimensional P-spline model can be written as

$$f(\mathit{z},\mathit{x}) = \left({\mathit{B}}_{\mathit{x}} \otimes {\mathit{B}}_{\mathit{z}}\right)\mathit{\theta } = \mathrm{vec}({\mathit{B}}_{\mathit{z}}\mathsf{\Theta }{\mathit{B}}_{\mathit{x}}^{\prime }),$$

(5)

where $\mathrm{vec}(·)$ denotes the vectorization operator.

In the two-dimensional case, a matrix $\mathit{P}$ penalizes the difference between adjacent coefficients of the rows and columns of the coefficient matrix $\mathsf{\Theta }$ (4). Therefore, the penalty matrix $\mathit{P}$ has the following expression:

$$\mathit{P} = {\mathsf{\lambda }}_{\mathit{z}}\underset{{\mathit{P}}^{\mathit{z}}}{\underbrace{{\mathit{I}}_{c_{\mathit{x}}}\otimes {\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}}}} + {\mathsf{\lambda }}_{\mathit{x}}\underset{{\mathit{P}}^{\mathit{x}}}{\underbrace{{\mathit{D}}_{\mathit{x}}^{\prime }{\mathit{D}}_{\mathit{x}}\otimes {\mathit{I}}_{c_{\mathit{z}}}}},$$

(6)

where ${\mathit{D}}_{\mathit{z}}$ and ${\mathit{D}}_{\mathit{x}}$ are the difference matrices acting on the rows and columns of $\mathsf{\Theta }$, respectively. ${\mathsf{\lambda }}_{\mathit{z}}$ and ${\mathsf{\lambda }}_{\mathit{x}}$ are the smoothing parameters for each dimension of the model. Since ${\mathsf{\lambda }}_{\mathit{z}}$ and ${\mathsf{\lambda }}_{\mathit{x}}$ are not necessary equal, penalty (6) allows for anisotropic smoothing. To estimate the regression coefficients, we minimize the following function for the penalized sum of squares:

$$S(\mathit{\theta }) = {(\mathit{y} - \mathit{B}\mathit{\theta })}^{\prime }(\mathit{y} - \mathit{B}\mathit{\theta }) + {\mathit{\theta }}^{\prime }\mathit{P}\mathit{\theta }.$$

(7)

Therefore, for given values of the smoothing parameters ${\mathsf{\lambda }}_{\mathit{z}}$ and ${\mathsf{\lambda }}_{\mathit{x}}$, the minimizer of the penalized sum of squares (7) is:

$$\widehat{\mathit{\theta }} = {({\mathit{B}}^{\prime }\mathit{B} + \mathit{P})}^{-1}{\mathit{B}}^{\prime }\mathit{y}.$$

(8)

Optimal ${\mathsf{\lambda }}_{\mathit{z}}$ and ${\mathsf{\lambda }}_{\mathit{x}}$ can be estimated using a information criteria (such as Akaike or Bayesian criteria) or via cross-validation (see [1,7] for details).

Once we have given this brief introduction to multidimensional P-splines, we generalize the approach given in [10] to obtain out-of-sample predictions not only in one, $f({\mathit{x}}_{+},\mathit{z})$, but in both independent variables, $f({\mathit{x}}_{+},{\mathit{z}}_{+})$ (the subscript ‘$+$’ denotes the dimension for which we are interested in for out-of-sample prediction). For the particular case in which a single covariate is extended, we will state several properties of the prediction method. Since the natural extension of the penalty matrix provides a change in the fit, in order to overcome this problem, we will suggest using constraints on the model coefficients.

2.1. Out-of-Sample Prediction of a Single Covariate

For the multidimensional P-spline model (1), we will consider first the case where the interest is on predicting outside the range of one of the covariates. This was already discussed in [10] for mortality life tables, where death counts were collected at different ages and years, and the aim was to forecast future mortality across ages.

Given a vector of $n_{\mathit{z}} \times n_{\mathit{x}}$ observations $\mathit{y}$ of the response variable, suppose we want to predict $n_p = n_{\mathit{z}} \times n_{{\mathit{x}}_p}$ new values ${\mathit{y}}_p$ at $\mathit{z}$ and ${\mathit{x}}_p$, where subscript p denotes the ‘predicted‘ range of the ${\mathit{x}}_p = {(x_1,x_2,\cdots ,x_n,x_{n+1},\cdots ,x_{n+p})}^{\prime }$ covariate and the response vector is ${\mathit{y}}_{+} = {(\mathit{y},{\mathit{y}}_p)}^{\prime }$, where $\mathit{y}$ is the vector of observed values and ${\mathit{y}}_p$ is the new unknown observations to be predicted that can take arbitrary values—if we arrange the vector of observations into a matrix $\mathit{Y}$ of dimension $n_{\mathit{z}} \times n_{\mathit{x}}$, such that $\mathrm{vec}(\mathit{Y}) = \mathit{y}$ of length $n_{\mathit{z}}{n}_{\mathit{x}} \times 1$ the new matrix of observations can be written as:

$${\mathit{Y}}_{+} = \left[\begin{array}{c}\mathit{Y}|{\mathit{Y}}_{\mathit{z}{\mathit{x}}_p}\end{array}\right],$$

(9)

where ${\mathit{Y}}_{\mathit{z}{\mathit{x}}_p}$ is a matrix of dimensions $n_{\mathit{z}} \times (n_x + n_{{\mathit{x}}_p})$. Hence, matrix ${\mathit{Y}}_{+}$ is, in this case, of dimensions $n_z \times n_{x_p}$. In order to fit and predict ${\mathit{y}}_{+}$ we consider the extended model:

$${\mathit{y}}_{+} = {\mathit{B}}_{+}{\mathit{\theta }}_{+} + {\mathit{ϵ}}_{+},{\mathit{ϵ}}_{+}\sim \mathcal{N}(\mathbf{0},{\sigma }_{\mathit{ϵ}}^2{\mathit{R}}_{+}),$$

(10)

where ${\mathit{R}}_{+}$ is the extended error covariance-matrix defined as ${\mathit{R}}_{+} = {\sigma }^2{\mathit{W}}_{+}$, with ${\mathit{W}}_{+} = {\mathit{W}}_{{\mathit{x}}_{+}} \otimes {\mathit{W}}_{\mathit{z}}$, and ${\mathit{W}}_{{\mathit{x}}_{+}}$ is a diagonal matrix with 1’s for the first $n_x$ entries (when the data are observed) and infinite values for the last $n_{x_p}$ (when the data are missing). For the $\mathit{z}$ covariate, as we are not predicting new observation, ${\mathit{W}}_{\mathit{z}}$ of a diagonal matrix of one’s of dimension $n_z \times n_z$.

Following [10], the extended basis is defined as:

$${\mathit{B}}_{+} = {\mathit{B}}_{{\mathit{x}}_{+}} \otimes {\mathit{B}}_{\mathit{z}} = \left[\begin{array}{cc}{\mathit{B}}_{\mathit{x}}& \mathit{O}\\ {\mathit{B}}_{{\mathit{x}}_{(1)}}& {\mathit{B}}_{{\mathit{x}}_{(2)}}\end{array}\right] \otimes {\mathit{B}}_{\mathit{z}} = \left[\begin{array}{cc}{\mathit{B}}_{\mathit{x}}\otimes {\mathit{B}}_{\mathit{z}}& \mathit{O}\\ {\mathit{B}}_{{\mathit{x}}_{(1)}}\otimes {\mathit{B}}_{\mathit{z}}& {\mathit{B}}_{{\mathit{x}}_{(2)}}\otimes {\mathit{B}}_{\mathit{z}}\end{array}\right],$$

(11)

where ${\mathit{B}}_{{\mathit{x}}_{+}} = \mathit{B}({\mathit{x}}_{+})$ and ${\mathit{B}}_{\mathit{z}} = \mathit{B}(\mathit{z})$ are the regression bases with ${\mathit{x}}_{+}$ and $\mathit{z}$, the two regressors. The new extended marginal B-spline basis for the $\mathit{x}$-covariate, ${\mathit{B}}_{{\mathit{x}}_{+}}$, is constructed from a new set of knots that consists of the original knots covering $x_i$, $i = 1, \dots , n_{\mathit{x}}$, and extended to the new range of the $n_{{\mathit{x}}_p}$ values of $x_{p_j}$, $j = 1, \dots , n_{{\mathit{x}}_p}$. In other words, ${\mathit{B}}_{{\mathit{x}}_{+}}$ is a direct extension of ${\mathit{B}}_{\mathit{x}}$. To estimate the new regression coefficients ${\mathit{\theta }}_{+}$, we minimize the extended penalized sum of squares:

$$\mathcal{L}({\mathit{\theta }}_{+}) = {({\mathit{y}}_{+} - {\mathit{B}}_{+}{\mathit{\theta }}_{+})}^{\prime }{\mathit{W}}_{+}^{-1}({\mathit{y}}_{+} - {\mathit{B}}_{+}{\mathit{\theta }}_{+}) + {\mathit{\theta }}_{+}^{\prime }{\mathit{P}}_{+}{\mathit{\theta }}_{+},$$

(12)

where the extended penalty matrix is given by:

$${\mathit{P}}_{+} = {\mathsf{\lambda }}_{\mathit{z}}{\mathit{P}}_{+}^{\mathit{z}} + {\mathsf{\lambda }}_{\mathit{x}}{\mathit{P}}_{+}^{{\mathit{x}}_{+}}.$$

(13)

since only the covariate $\mathit{x}$ is extended, ${\mathit{P}}_{+}^{\mathit{z}}$ and ${\mathit{P}}_{+}^{{\mathit{x}}_{+}}$ are:

$$\begin{array}{ccc}\hfill {\mathit{P}}_{+}^{\mathit{z}}& =& {\mathit{I}}_{{c_{\mathit{x}}}_{+}} \otimes {\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}} = \left[\begin{array}{cc}{\mathit{I}}_{c_{\mathit{x}}}& \mathit{O}\\ \mathit{O}& {\mathit{I}}_{c_{{\mathit{x}}_p}}\end{array}\right] \otimes {\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}} = \left[\begin{array}{cc}{\mathit{I}}_{c_{\mathit{x}}}\otimes {\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}}& \mathit{O}\\ \mathit{O}& {\mathit{I}}_{c_{{\mathit{x}}_p}}\otimes {\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}}\end{array}\right]\hfill \\ & =& \left[\begin{array}{cc}{\mathit{P}}_{{+}_{11}}^{\mathit{z}}& \mathit{O}\\ \mathit{O}& {\mathit{P}}_{{+}_{22}}^{\mathit{z}}\end{array}\right],\hfill \end{array}$$

(14)

and

$$\begin{array}{ccc}\hfill {\mathit{P}}_{+}^{{\mathit{x}}_{+}}& =& {{\mathit{D}}_{\mathit{x}}}_{+}^{\prime }{{\mathit{D}}_{\mathit{x}}}_{+} \otimes {\mathit{I}}_{c_{\mathit{z}}} = \left[\begin{array}{cc}({\mathit{D}}_{\mathit{x}}^{\prime }{\mathit{D}}_{\mathit{x}} + {{\mathit{D}}_{\mathit{x}}^{\prime }}_{(1)}{{\mathit{D}}_{\mathit{x}}}_{(1)})\otimes {\mathit{I}}_{c_{\mathit{z}}}& {{\mathit{D}}_{\mathit{x}}^{\prime }}_{(1)}{{\mathit{D}}_{\mathit{x}}}_{(2)}\otimes {\mathit{I}}_{c_{\mathit{z}}}\\ {{\mathit{D}}_{\mathit{x}}^{\prime }}_{(2)}{{\mathit{D}}_{\mathit{x}}}_{(1)}\otimes {\mathit{I}}_{c_{\mathit{z}}}& {{\mathit{D}}_{\mathit{x}}^{\prime }}_{(2)}{{\mathit{D}}_{\mathit{x}}}_{(2)}\otimes {\mathit{I}}_{c_{\mathit{z}}}\end{array}\right]\hfill \\ & =& \left[\begin{array}{cc}{\mathit{P}}_{{+}_{11}}^{{\mathit{x}}_{+}}& {\mathit{P}}_{{+}_{12}}^{{\mathit{x}}_{+}}\\ {\mathit{P}}_{{+}_{21}}^{{\mathit{x}}_{+}}& {\mathit{P}}_{{+}_{22}}^{{\mathit{x}}_{+}}\end{array}\right],\hfill \end{array}$$

(15)

where ${\mathit{D}}_{{\mathit{x}}_{+}}$ and ${\mathit{D}}_{\mathit{z}}$ are the difference matrices acting on the columns and rows of the matrix formed by the extended vector of coefficients:

$${\mathsf{\Theta }}_{+} = \left[\begin{array}{cccccccc}{\mathit{\theta }}_{11}& {\mathit{\theta }}_{12}& \cdots & {\mathit{\theta }}_{1c_{\mathit{x}}-1}& {\mathit{\theta }}_{1c_{\mathit{x}}}& {\mathit{\theta }}_{1c_{\mathit{x}}+1}& {\mathit{\theta }}_{1c_{\mathit{x}}+2}& \cdots \\ {\mathit{\theta }}_{21}& {\mathit{\theta }}_{22}& \cdots & {\mathit{\theta }}_{2c_{\mathit{x}}-1}& {\mathit{\theta }}_{2c_{\mathit{x}}}& {\mathit{\theta }}_{2c_{\mathit{x}}+1}& {\mathit{\theta }}_{2c_{\mathit{x}}+2}& \cdots \\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathit{\theta }}_{c_{\mathit{z}}1}& {\mathit{\theta }}_{c_{\mathit{z}}2}& \cdots & {\mathit{\theta }}_{c_{\mathit{z}}{c}_{\mathit{x}}-1}& {\mathit{\theta }}_{c_{\mathit{z}}{c}_{\mathit{x}}}& {\mathit{\theta }}_{c_{\mathit{z}{c}_{\mathit{x}}+1}}& {\mathit{\theta }}_{c_{\mathit{z}}{c}_{\mathit{x}}+2}& \cdots \end{array}\right].$$

(16)

The minimizer of (12) is given by:

$${\mathit{\theta }}_{+} = {({\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{B}}_{+} + {\mathit{P}}_{+})}^{-1}{\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{y}}_{+}.$$

(17)

The fact that we extend just one of the two covariates, and the particular form of the penalties based on differences between adjacent coefficients, yields certain important properties satisfied by the model. These properties are the direct result of the following theorem:

Theorem 1.

The coefficients obtained from minimization of (12) with the extended basis (11), extended error covariance matrix (10) and extended penalty matrix (13) meet the following properties:

I. 

The first c, $c = c_{\mathit{z}} \times c_{\mathit{x}}$, coefficients of ${\widehat{\mathit{\theta }}}_{+}$, are:

$$\begin{array}{c}\hfill {\widehat{\mathit{\theta }}}_{{+}_{1,\dots ,c}} = {\left({\mathit{B}}^{\prime }\mathit{B} + {\mathsf{\lambda }}_{\mathit{x}}{\mathit{P}}_{{+}_{11}}^{{\mathit{x}}_{+}} + {\mathsf{\lambda }}_{\mathit{z}}{\mathit{P}}_{{+}_{11}}^{\mathit{z}} - {\mathsf{\lambda }}_{\mathit{x}}^2{\mathit{P}}_{{+}_{12}}^{{\mathit{x}}_{+}}{({\mathsf{\lambda }}_{\mathit{x}}{\mathit{P}}_{{+}_{22}}^{{\mathit{x}}_{+}} + {\mathsf{\lambda }}_{\mathit{z}}{\mathit{P}}_{{+}_{22}}^{\mathit{z}})}^{-1}{\mathit{P}}_{{+}_{21}}^{{\mathit{x}}_{+}}\right)}^{-1}{\mathit{B}}^{\prime }\mathit{y},\\ & \end{array}$$

(18)

where ${\mathit{P}}_{{+}_{11}}^{\mathit{z}}$ and ${\mathit{P}}_{{+}_{22}}^{\mathit{z}}$ are defined in (14) and ${\mathit{P}}_{{+}_{11}}^{{\mathit{x}}_{+}}$, ${\mathit{P}}_{{+}_{12}}^{{\mathit{x}}_{+}}$, ${\mathit{P}}_{{+}_{21}}^{{\mathit{x}}_{+}}$ and ${\mathit{P}}_{{+}_{22}}^{{\mathit{x}}_{+}}$ defined in (15).

II. 

The coefficients for the $n_p = n_{\mathit{z}} \times n_{{\mathit{x}}_p}$ predicted values are

$${\widehat{\mathit{\theta }}}_p = -{\left(\frac{{\mathsf{\lambda }}_{\mathit{z}}}{{\mathsf{\lambda }}_{\mathit{x}}}{\mathit{P}}_{{+}_{22}}^{\mathit{z}} + {\mathit{P}}_{{+}_{22}}^{{\mathit{x}}_{+}}\right)}^{-1}{\mathit{P}}_{{+}_{21}}^{{\mathit{x}}_{+}}{\widehat{\mathit{\theta }}}_{{+}_{1,\dots ,c}},$$

(19)

where ${\mathit{P}}_{{+}_{22}}^{\mathit{z}}$ is defined in (14) and ${\mathit{P}}_{{+}_{22}}^{{\mathit{x}}_{+}}$, ${\mathit{P}}_{{+}_{21}}^{{\mathit{x}}_{+}}$ is defined in (15).

The proof of Theorem 1 is given in the Supplementary Materials.

Therefore, according to the previous theorem, when we predict in two dimensions extending one covariate, the predicted values, ${\mathit{y}}_p$, obtained by using the new coefficients, ${\mathit{\theta }}_p$; depend on the ratio ${\mathsf{\lambda }}_{\mathit{z}}/{\mathsf{\lambda }}_{\mathit{x}}$, unless ${\mathsf{\lambda }}_{\mathit{x}} = {\mathsf{\lambda }}_{\mathit{z}}$ (i.e., an isotropic smoothing). Therefore, while in the case of one-dimensional smooth models, the smoothing parameter does not play any role in the prediction (see [13] for details), we have found that in two dimensions the smoothing parameters in both directions, ${\mathsf{\lambda }}_{\mathit{x}}$ and ${\mathsf{\lambda }}_{\mathit{z}}$, determine the prediction.

Notice that we have proved that the coefficients that give the fit when the fit and the prediction are obtained at the same time (18) are different from the solution we obtain only fitting the data (8) (this property is verified in the case of one-dimensional smooth models [13]). Although, in the one-dimensional case, the extended penalty is not a direct extension of the penalty used to fit the data, the blocks of the extended penalty are simplified and the fit remains the same. This does not happen in the case of two dimensions, unless the block ${\mathit{P}}_{+22}^{\mathit{z}}$ in (14) is equal to zero, as shown by the following corollary.

Corollary 1

(Theorem 1). If ${\mathit{P}}_{{+}_{22}}^{\mathit{z}} = \mathit{O}$ (or equivalently, if ${\mathit{I}}_{c_{{\mathit{x}}_p}} = \mathit{O}$) in (14), the solution of (12) verifies the following properties:

1. 

When performing out-of-sample predictions, the fit remains invariant.

2. 

Consider the coefficient matrix obtained from the model fit, $\widehat{\mathsf{\Theta }}$, and the c matrix of predicted coefficients, ${\widehat{\mathsf{\Theta }}}_p$, each row $j = 1, \dots , c_{\mathit{z}}$, of the additional matrix of coefficients is a linear combination of the last $q_{\mathit{z}}$ old coefficients of that row ($q_{\mathit{z}}$ is the order of the penalty acting on rows and $c_{\mathit{z}}$ is the number of rows of ${\mathit{B}}_{\mathit{z}}$).

For the particular case of second and third order penalties, i.e., $q_{\mathit{x}} = q_{\mathit{z}} = 2$ and $q_{\mathit{x}} = q_{\mathit{z}} = 3$, the proof is given in the Supplementary Materials.

As we have showed in the previous corollary, setting ${\mathit{I}}_{c_{{\mathit{x}}_p}}$ to zero, the fit is preserved and everything is analogous to the one-dimensional case. In the literature, there are some works in which ${\mathit{I}}_{c_{{\mathit{x}}_p}}$ is considered equal to zero, e.g., [12]. However, in practice, we do not set ${\mathit{I}}_{c_{{\mathit{x}}_p}}$ to zero since we would not be imposing the penalty correctly. In addition, we cannot extend it to situations where we want to make out-of-sample predictions in two dimensions. We cannot set ${\mathit{I}}_{c_{{\mathit{z}}_p}}$ and ${\mathit{I}}_{c_{{\mathit{x}}_p}}$ to zero in (14) and (15), since the matrix ${\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{B}}_{+} + {\mathit{P}}_{+}$ would be singular.

It is important to note that, regardless of the value of ${\mathit{P}}_{{+}_{22}}^{{\mathit{x}}_{+}}$, by (19) the new coefficients are determined by the last $q_z$ columns of the coefficients that give the fit. Therefore, we do not need to set ${\mathit{P}}_{{+}_{22}}^{{\mathit{x}}_{+}}$ equal to zero to know which coefficients determine the prediction.

2.2. Out-of-Sample Prediction in Two Dimensions

In this section, we address the case in which we are interested in making predictions outside the range of the observed values of both covariates $\mathit{x}$ and $\mathit{z}$.

Given a vector of $n_{\mathit{z}}{n}_{\mathit{x}}$ observations $\mathit{y}$ of the response variable, now we want to predict $n_p = n_{\mathit{z}}{n}_{{\mathit{x}}_p} + n_{{\mathit{z}}_p}{n}_{\mathit{x}} + n_{{\mathit{z}}_p}{n}_{{\mathit{x}}_p}$ new values at $(\mathit{z},{\mathit{x}}_p)$, $({\mathit{z}}_p,\mathit{x})$ and $({\mathit{z}}_p,{\mathit{x}}_p)$, i.e., if we arrange the observation vector into a matrix $\mathit{Y}$ of dimension $n_{\mathit{z}} \times n_{\mathit{x}}$, the observed and predicted values can be arranged into a matrix of dimension $n_{{\mathit{z}}_{+}} \times n_{{\mathit{x}}_{+}}$ ($n_{{\mathit{z}}_{+}} = n_{\mathit{z}} + n_{{\mathit{z}}_p}$, $n_{{\mathit{x}}_{+}} = n_{\mathit{x}} + n_{{\mathit{x}}_p}$), as:

$${\mathit{Y}}_{+} = \left[\begin{array}{cc}\mathit{Y}& {\mathit{Y}}_{\mathit{z}{\mathit{x}}_p}\\ {\mathit{Y}}_{{\mathit{z}}_p\mathit{x}}& {\mathit{Y}}_{{\mathit{z}}_p{\mathit{x}}_p}\end{array}\right].$$

(20)

Note that the dimensions of ${\mathit{Y}}_{\mathit{z}{\mathit{x}}_p}$, ${\mathit{Y}}_{{\mathit{z}}_p\mathit{x}}$ and ${\mathit{Y}}_{{\mathit{z}}_p{\mathit{x}}_p}$ are $n_{\mathit{z}} \times n_{{\mathit{x}}_p}$, $n_{{\mathit{z}}_p} \times n_{\mathit{x}}$ and $n_{{\mathit{z}}_p} \times n_{{\mathit{x}}_p}$, respectively.

We again use the model in (10), but now the full extended basis is the Kronecker product of the two extended marginal B-spline basis, ${\mathit{B}}_{{\mathit{x}}_{+}} = \mathit{B}({\mathit{x}}_{+})$ and ${\mathit{B}}_{{\mathit{z}}_{+}} = \mathit{B}({\mathit{z}}_{+})$, of dimensions $n_{{\mathit{x}}_{+}} \times c_{{\mathit{x}}_{+}}$ and $n_{{\mathit{z}}_{+}} \times c_{{\mathit{z}}_{+}}$, respectively:

$${\mathit{B}}_{+} = {\mathit{B}}_{{\mathit{x}}_{+}}\otimes {\mathit{B}}_{{\mathit{z}}_{+}} = \left[\begin{array}{cc}{\mathit{B}}_{\mathit{x}}& \mathit{O}\\ {\mathit{B}}_{{\mathit{x}}_{(1)}}& {\mathit{B}}_{{\mathit{x}}_{(2)}}\end{array}\right]\otimes \left[\begin{array}{cc}{\mathit{B}}_{\mathit{z}}& \mathit{O}\\ {\mathit{B}}_{{\mathit{z}}_{(1)}}& {\mathit{B}}_{{\mathit{z}}_{(2)}}\end{array}\right],$$

(21)

where the extended bases ${\mathit{B}}_{{\mathit{x}}_{+}}$ and ${\mathit{B}}_{{\mathit{z}}_{+}}$ are constructed from a new set of knots that consists of the original knots and extended to cover the full range of ${\mathit{x}}_{+}$ and ${\mathit{z}}_{+}$, respectively. The extended covariance matrix ${\mathit{R}}_{+} = {\sigma }_{\mathit{ϵ}}^2{\mathit{W}}_{+}$ with ${\mathit{W}}_{+} = {\mathit{W}}_{{\mathit{x}}_{+}}\otimes {\mathit{W}}_{{\mathit{z}}_{+}}$, where ${\mathit{W}}_{{\mathit{x}}_{+}}$ and ${\mathit{W}}_{{\mathit{z}}_{+}}$ are diagonal matrices of dimensions $n_{{\mathit{x}}_{+}} \times n_{{\mathit{x}}_{+}}$ ($n_{{\mathit{x}}_{+}} = n_{\mathit{x}} + n_{{\mathit{x}}_p}$) and $n_{{\mathit{z}}_{+}} \times n_{{\mathit{z}}_{+}}$ ($n_{{\mathit{z}}_{+}} = n_{\mathit{z}} + n_{{\mathit{z}}_p}$), respectively, with infinity entry if the data are missing and 1 if the data are observed. The infinity quantities express the uncertainty since we do not have any information about the data ${\mathit{y}}_p$. To estimate the extended coefficients, we minimize the penalized sum of squares in (12) with an extended penalty matrix

$${\mathit{P}}_{+} = {\mathsf{\lambda }}_{\mathit{z}}{\mathit{P}}_{+}^{{\mathit{z}}_{+}} + {\mathsf{\lambda }}_{\mathit{x}}{\mathit{P}}_{+}^{{\mathit{x}}_{+}},$$

(22)

with

$$\begin{array}{ccc}\hfill {\mathit{P}}_{+}^{{\mathit{z}}_{+}}& =& {\mathit{I}}_{{c_{\mathit{x}}}_{+}} \otimes {\mathit{D}}_{{\mathit{z}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{z}}_{+}} = \left[\begin{array}{cc}{\mathit{I}}_{c_{\mathit{x}}}\otimes {\mathit{D}}_{{\mathit{z}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{z}}_{+}}& \mathit{O}\\ \mathit{O}& {\mathit{I}}_{c_{{\mathit{x}}_p}}\otimes {\mathit{D}}_{{\mathit{z}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{z}}_{+}}\end{array}\right]\hfill \\ & =& \left[\begin{array}{cc}{\mathit{P}}_{{+}_{11}}^{{\mathit{z}}_{+}}& \mathit{O}\\ \mathit{O}& {\mathit{P}}_{{+}_{22}}^{{\mathit{z}}_{+}}\end{array}\right],\hfill \end{array}$$

(23)

and ${\mathit{P}}_{+}^{{\mathit{x}}_{+}}$ as in (15). Note that ${\mathit{D}}_{{\mathit{z}}_{+}}$ and ${\mathit{D}}_{{\mathit{x}}_{+}}$ are direct extensions of ${\mathit{D}}_{\mathit{z}}$ and ${\mathit{D}}_{\mathit{x}}$ but ${\mathit{P}}_{+}^{{\mathit{z}}_{+}}$ and ${\mathit{P}}_{+}^{{\mathit{x}}_{+}}$ are not direct extensions of ${\mathit{P}}^{\mathit{z}}$ and ${\mathit{P}}^{\mathit{x}}$. Moreover, if $\mathit{\theta }$ is given in (8) and we are extending the two covariates, ${\mathit{\theta }}_{+} = \mathrm{vec}({\mathsf{\Theta }}_{+})$, in matrix form ${\mathsf{\Theta }}_{+}$ is:

$${\mathsf{\Theta }}_{+} = \left[\begin{array}{cccccccc}{\mathit{\theta }}_{11}& {\mathit{\theta }}_{12}& \cdots & {\mathit{\theta }}_{1c_{\mathit{x}}-1}& {\mathit{\theta }}_{1c_{\mathit{x}}}& {\mathit{\theta }}_{1c_{\mathit{x}}+1}& {\mathit{\theta }}_{1c_{\mathit{x}}+2}& \cdots \\ {\mathit{\theta }}_{21}& {\mathit{\theta }}_{22}& \cdots & {\mathit{\theta }}_{2c_{\mathit{x}}-1}& {\mathit{\theta }}_{2c_{\mathit{x}}}& {\mathit{\theta }}_{2c_{\mathit{x}}+1}& {\mathit{\theta }}_{2c_{\mathit{x}}+2}& \cdots \\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathit{\theta }}_{c_{\mathit{z}}1}& {\mathit{\theta }}_{c_{\mathit{z}}2}& \cdots & {\mathit{\theta }}_{c_{\mathit{z}}{c}_{\mathit{x}}-1}& {\mathit{\theta }}_{c_{\mathit{z}}{c}_{\mathit{x}}}& {\mathit{\theta }}_{c_{\mathit{z}}{c}_{\mathit{x}}+1}& {\mathit{\theta }}_{c_{\mathit{z}}{c}_{\mathit{x}}+2}& \cdots \\ {\mathit{\theta }}_{c_{\mathit{z}}+11}& {\mathit{\theta }}_{c_{\mathit{z}}+12}& \cdots & {\mathit{\theta }}_{c_{\mathit{z}}+1c_{\mathit{x}}-1}& {\mathit{\theta }}_{c_{\mathit{z}}+1c_{\mathit{x}}}& {\mathit{\theta }}_{c_{\mathit{z}}+1c_{\mathit{x}}+1}& {\mathit{\theta }}_{c_{\mathit{z}}+1c_{\mathit{x}}+2}& \cdots \\ {\mathit{\theta }}_{c_{\mathit{z}}+21}& {\mathit{\theta }}_{c_{\mathit{z}}+22}& \cdots & {\mathit{\theta }}_{c_{\mathit{z}}+2c_{\mathit{x}}-1}& {\mathit{\theta }}_{c_{\mathit{z}}+2c_{\mathit{x}}}& {\mathit{\theta }}_{c_{\mathit{z}}+2c_{\mathit{x}}+1}& {\mathit{\theta }}_{c_{\mathit{z}}+2c_{\mathit{x}}+2}& \cdots \\ ⋮& ⋮& \cdots & ⋮& ⋮& ⋮& ⋮& \cdots \end{array}\right].$$

Note that in the case of fitting, the extension of the predicted new values depends on the data structure. If we consider scattered data, we set $n_p$ out-of-sample points $(z_{{+}_i},x_{{+}_i})$ at which we want to predict new $y_{p_i}$ values for $i = 1, \dots , n_p$ and ${\mathit{R}}_{+}$ is a diagonal matrix with the first n values equal to 1 and the last $n_p$ values equal to infinity. Everything else is independent of the data structure.

As in the case of only one-dimensional out-of-sample prediction, when using model (1) or model (10), the fit to the observed data is not the same, as one might expect. In the next section, we recommend using constraints to maintain the fit when fitting and prediction are obtained at the same time.

2.3. Out-of-Sample Prediction with Constraints

As we showed in the previous section, the natural expansion of the penalty matrix provides a change in fit. To overcome this problem, and as a possible way to incorporate known information about the prediction, we recommend the use of constrained P-splines. This method will allow us to impose constant and fixed constraints, as well as constraints that depend on observed data.

Our proposal is to derive the solution of the extended model (10) subject to a set of l linear constraints given by $\mathit{C}{\mathit{\theta }}_{+} = \mathit{r}$, where $\mathit{C}$ is a known matrix of dimension $l \times c_{+}$ acting on all coefficients, and $\mathit{r}$ is the restrictions vector of dimension $l \times 1$. Therefore, we have the following extended regression model with constraints:

$${\mathit{y}}_{+} = {\mathit{B}}_{+}{\mathit{\theta }}_{+}^{*} + {\mathit{ϵ}}_{+},{\mathit{ϵ}}_{+}\sim \mathcal{N}(\mathbf{0},{\mathit{R}}_{+})$$

(24)

subject to

$$\mathit{C}{\mathit{\theta }}_{+}^{*} = \mathit{r},$$

(25)

Depending on whether we are predicting out-of-sample in one or two dimensions, we define ${\mathit{y}}_{+}$, ${\mathit{B}}_{+}$ and ${\mathit{R}}_{+}$ as shown in Section 2.1 and Section 2.2. The superscript ${(}^{*})$ refers to the use of constraints.

We extend the results of [14] to the case of penalized least squares and obtain the Lagrange formulation of the penalized least squares minimization problem:

$$\mathcal{L}({\mathit{\theta }}_{+}^{*},\mathit{\omega }) = {({\mathit{y}}_{+} - {\mathit{B}}_{+}{\mathit{\theta }}_{+}^{*})}^{\prime }{\mathit{W}}_{+}^{-1}({\mathit{y}}_{+} - {\mathit{B}}_{+}{\mathit{\theta }}_{+}^{*}) + {\mathit{\theta }}_{+}^{{*}^{\prime }}{\mathit{P}}_{+}{\mathit{\theta }}_{+}^{*} + 2{\mathit{\omega }}^{\prime }(\mathit{C}{\mathit{\theta }}_{+}^{*} - \mathit{r}),$$

(26)

where ${\mathit{W}}_{+}^{-1}$ is defined depending on whether we extend one or two covariates, ${\mathit{P}}_{+}$ is the extended penalty matrix (13) or (22), ${\mathit{\theta }}_{+}^{*}$ denotes the restricted penalized least squares (RPLS) estimator and $\mathit{\omega }$ is a $l \times 1$ vector of Lagrange multipliers.

Differentiating the expression in (26), we obtain:

$$\begin{array}{ccc}\hfill \frac{\partial \mathcal{L}}{\partial {\mathit{\theta }}_{+}^{*}}& =& -2{\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}({\mathit{y}}_{+} - {\mathit{B}}_{+}{\mathit{\theta }}_{+}^{*}) + 2{\mathit{P}}_{+}{\mathit{\theta }}_{+}^{*} + 2{\mathit{C}}^{\prime }\mathit{\omega },\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill \frac{\partial \mathcal{L}}{\partial \mathit{\omega }}& =& \mathit{C}{\mathit{\theta }}_{+}^{*} - \mathit{r}.\hfill \end{array}$$

(28)

The system of equations can be written as:

$$\left[\begin{array}{cc}{\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{B}}_{+}+{\mathit{P}}_{+}& {\mathit{C}}^{\prime }\\ \mathit{C}& \mathit{O}\end{array}\right]\left[\begin{array}{c}{\widehat{\mathit{\theta }}}_{+}^{*}\\ \widehat{\mathit{\omega }}\end{array}\right] = \left[\begin{array}{c}{\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{y}}_{+}\\ \mathit{r}\end{array}\right].$$

yielding the following solutions:

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\theta }}}_{+}^{*}& =& {\widehat{\mathit{\theta }}}_{+} - {({\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{B}}_{+}+{\mathit{P}}_{+})}^{-1}{\mathit{C}}^{\prime }\widehat{\mathit{\omega }},\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill \widehat{\mathit{\omega }}& =& {[\mathit{C}{({\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{B}}_{+}+{\mathit{P}}_{+})}^{-1}{\mathit{C}}^{\prime }]}^{-1}(\mathit{C}{\mathit{\theta }}_{+} - \mathit{r}).\hfill \end{array}$$

(30)

where ${\widehat{\mathit{\theta }}}_{+}$ is the unconstrained penalized least squares estimator as in (17).

Therefore, the constrained vector of coefficients ${\widehat{\mathit{\theta }}}_{+}^{*}$, are obtained by computing the vector of Lagrange multipliers (30) and substituting in (29), i.e., ${\widehat{\mathit{\theta }}}_{+}^{*}$ is the unconstrained solution, ${\widehat{\mathit{\theta }}}_{+}$, plus a multiple of the discrepancy vector.

Using these results, we can obtain the constrained fitted and predicted values as

$${\widehat{\mathit{y}}}_{+}^{*} = {\mathit{B}}_{+}{\widehat{\mathit{\theta }}}_{+}^{*},$$

defining the matrices

$$\begin{array}{ccc}\hfill {\mathit{A}}_{\mathbf{1}}& =& {({\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{B}}_{+}+{\mathit{P}}_{+})}^{-1}{\mathit{B}}_{+}^{\prime }\hfill \\ \hfill {\mathit{A}}_{\mathbf{2}}& =& {({\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{B}}_{+}+{\mathit{P}}_{+})}^{-1}{\mathit{C}}^{\prime }{[\mathit{C}{({\mathit{B}}_{+}^{\prime }{\mathit{W}}_{+}^{-1}{\mathit{B}}_{+}+{\mathit{P}}_{+})}^{-1}{\mathit{C}}^{\prime }]}^{-1},\hfill \end{array}$$

the response vector of new values ${\widehat{\mathit{y}}}_{+}^{*}$ can be written as:

$${\widehat{\mathit{y}}}_{+}^{*} = {\mathit{B}}_{+}({\mathit{A}}_{\mathbf{1}}{\mathit{W}}_{+}^{-1}{\mathit{y}}_{+} - {\mathit{A}}_{\mathbf{2}}\mathit{C}{\mathit{A}}_{\mathbf{1}}{\mathit{W}}_{+}^{-1}{\mathit{y}}_{+} + {\mathit{A}}_{\mathbf{2}}\mathit{r}).$$

The variance of ${\mathit{y}}_{+}^{*}$ depends on the following conditions:

(a)

If the constraint $\mathit{r}$ is a constant, i.e., it does not depend on the data, the variance is:

$$\mathrm{Var}\left[{\widehat{\mathit{y}}}_{+}^{*}\right] = {\sigma }_{\mathit{ϵ}}^2{\mathit{B}}_{+}{\mathit{A}}_{\mathbf{3}}{\mathit{W}}_{+}^{-1}{\mathit{A}}_{\mathbf{3}}^{\prime }{\mathit{B}}_{+}^{\prime },$$

with ${\mathit{A}}_{\mathbf{3}} = {\mathit{A}}_{\mathbf{1}} - {\mathit{A}}_{\mathbf{2}}\mathit{C}{\mathit{A}}_{\mathbf{1}}$.

(b)

If the constraint depends on the data, we must consider the variability of $\mathit{r}$. For instance, if the constraint implies that the fit has to be maintained, i.e., $\mathit{r} = \widehat{\mathit{\theta }}$, the variance is:

$$\mathrm{Var}\left[{\widehat{\mathit{y}}}_{+}^{*}\right] = {\sigma }_{\mathit{ϵ}}^2{\mathit{B}}_{+}{\mathit{A}}_{\mathbf{4}}{\mathit{W}}_{+}^{-1}{\mathit{A}}_{\mathbf{4}}^{\prime }{\mathit{B}}_{+}^{\prime },$$

with

$${\mathit{A}}_{\mathbf{4}}={\mathit{A}}_{\mathbf{1}} - {\mathit{A}}_{\mathbf{2}}\mathit{C}{\mathit{A}}_{\mathbf{1}} + {\mathit{A}}_{\mathbf{2}}\left[\begin{array}{cc}{({\mathit{B}}^{\prime }\mathit{B}+\mathit{P})}^{-1}& \\ & \mathit{O}\end{array}\right],$$

where $\mathit{O}$ is a null matrix of dimension, $c_p \times n_p$, $c_p = c_{\mathit{z}}{c}_{{\mathit{x}}_p} + c_{{\mathit{z}}_p}{c}_{\mathit{x}} + c_{{\mathit{z}}_p}{c}_{{\mathit{x}}_p}$ is the number of new coefficients, and $n_p$ is the number of new observations.

To help readers understand the process, we now explain how to impose constraints on the fit in practice. Suppose we just make out-of-sample predictions in one of the two covariates and, for example, the estimated coefficients matrix obtained when only fitting the observed data has dimensions of $4 \times 3$, and the one obtained when fitting and forecasting simultaneously has dimension $4 \times 5$, then

$$\widehat{\mathsf{\Theta }} = \left[\begin{array}{ccc}{\widehat{\theta }}_1& {\widehat{\theta }}_5& {\widehat{\theta }}_9\\ {\widehat{\theta }}_2& {\widehat{\theta }}_6& {\widehat{\theta }}_{10}\\ {\widehat{\theta }}_3& {\widehat{\theta }}_7& {\widehat{\theta }}_{11}\\ {\widehat{\theta }}_4& {\widehat{\theta }}_8& {\widehat{\theta }}_{12}\end{array}\right],{\mathsf{\Theta }}_{+}^{*} = \left[\begin{array}{ccccc}{\stackrel{•}{\theta }}_1& {\stackrel{•}{\theta }}_5& {\stackrel{•}{\theta }}_9& {\stackrel{\circ }{\theta }}_{13}& {\stackrel{\circ }{\theta }}_{17}\\ {\stackrel{•}{\theta }}_2& {\stackrel{•}{\theta }}_6& {\stackrel{•}{\theta }}_{10}& {\stackrel{\circ }{\theta }}_{14}& {\stackrel{\circ }{\theta }}_{18}\\ {\stackrel{•}{\theta }}_3& {\stackrel{•}{\theta }}_7& {\stackrel{•}{\theta }}_{11}& {\stackrel{\circ }{\theta }}_{15}& {\stackrel{\circ }{\theta }}_{19}\\ {\stackrel{•}{\theta }}_4& {\stackrel{•}{\theta }}_8& {\stackrel{•}{\theta }}_{12}& {\stackrel{\circ }{\theta }}_{16}& {\stackrel{\circ }{\theta }}_{20}\end{array}\right],$$

where in ${\mathsf{\Theta }}_{+}^{*}$ the coefficients that determine the fit have superscript ${}^{•}$ and the coefficients that determine the forecast have superscript symbol ${}^{\circ }$. If we impose the constraint that the fit has to be maintained, we use the constraint equation in (25) with ${\mathit{\theta }}_{+}^{*} = \mathrm{Vec}({\mathsf{\Theta }}_{+}^{*})$, $\mathit{C} = \left[{\mathit{I}}_{12\times 12}\right|{\mathit{O}}_{12\times 8}]$ (${\mathit{I}}_{12\times 12}$ is an identity matrix of dimension 12 and ${\mathit{O}}_{12\times 8}$ a zero matrix of dimension $12 \times 8$) and $\mathit{r} = \widehat{\mathit{\theta }} = \mathrm{vec}(\widehat{\mathsf{\Theta }})$.

On the other hand, if we extend the two covariates, the coefficient matrices obtained when only fitting and when fit and prediction is carried out simultaneously are:

$$\widehat{\mathsf{\Theta }} = \left[\begin{array}{ccc}{\widehat{\theta }}_1& {\widehat{\theta }}_5& {\widehat{\theta }}_9\\ {\widehat{\theta }}_2& {\widehat{\theta }}_6& {\widehat{\theta }}_{10}\\ {\widehat{\theta }}_3& {\widehat{\theta }}_7& {\widehat{\theta }}_{11}\\ {\widehat{\theta }}_4& {\widehat{\theta }}_8& {\widehat{\theta }}_{12}\end{array}\right],{\mathsf{\Theta }}_{+}^{*} = \left[\begin{array}{ccccc}{\stackrel{•}{\theta }}_1& {\stackrel{•}{\theta }}_7& {\stackrel{•}{\theta }}_{13}& {\stackrel{\circ }{\theta }}_{19}& {\stackrel{\circ }{\theta }}_{25}\\ {\stackrel{•}{\theta }}_2& {\stackrel{•}{\theta }}_8& {\stackrel{•}{\theta }}_{14}& {\stackrel{\circ }{\theta }}_{20}& {\stackrel{\circ }{\theta }}_{26}\\ {\stackrel{•}{\theta }}_3& {\stackrel{•}{\theta }}_9& {\stackrel{•}{\theta }}_{15}& {\stackrel{\circ }{\theta }}_{21}& {\stackrel{\circ }{\theta }}_{27}\\ {\stackrel{•}{\theta }}_4& {\stackrel{•}{\theta }}_{10}& {\stackrel{•}{\theta }}_{16}& {\stackrel{\circ }{\theta }}_{22}& {\stackrel{\circ }{\theta }}_{28}\\ {\stackrel{\circ }{\theta }}_5& {\stackrel{\circ }{\theta }}_{11}& {\stackrel{\circ }{\theta }}_{17}& {\stackrel{\circ }{\theta }}_{23}& {\stackrel{\circ }{\theta }}_{29}\\ {\stackrel{\circ }{\theta }}_6& {\stackrel{\circ }{\theta }}_{12}& {\stackrel{\circ }{\theta }}_{18}& {\stackrel{\circ }{\theta }}_{24}& {\stackrel{\circ }{\theta }}_{30}\end{array}\right],$$

i.e., $c_{\mathit{z}} = 4$, $c_{{\mathit{z}}_p} = 2$, $c_{\mathit{x}} = 3$ and $c_{{\mathit{x}}_p} = 2$. To impose the restriction the fit has to be maintained, we set ${\mathit{\theta }}_{+}^{*} = \mathrm{Vec}({\mathsf{\Theta }}_{+}^{*})$, $\mathit{C} = \mathrm{blockdiag}({\mathit{I}}_{4\times 4},\left[{\mathit{O}}_{2\times 2}\right|{\mathit{I}}_{4\times 4}],\left[{\mathit{O}}_{2\times 2}\right|{\mathit{I}}_{4\times 4}])$ (${\mathit{I}}_{4\times 4}$ is an identity matrix of dimension 4 and ${\mathit{O}}_{2\times 2}$ is a zero matrix of dimension 2) and $\mathit{r} = \widehat{\mathit{\theta }} = \mathrm{vec}(\widehat{\mathsf{\Theta }})$ in Equation (25).

In general, if the constraint implies that the fit has to be maintained and $c_{{\mathit{x}}_p} \ge c_{{\mathit{z}}_p}$, $\mathit{C}$ is a block diagonal matrix with the first block an identity matrix of dimensions $c_{\mathit{z}} \times c_{\mathit{z}}$ and $c_{{\mathit{x}}_p}$ blocks equal to $[{\mathit{O}}_{c_{{\mathit{z}}_p}\times c_{{\mathit{z}}_p}}|{\mathit{I}}_{c_{\mathit{z}}\times c_{\mathit{z}}}]$, i.e.,

$$\mathit{C} = \left[\begin{array}{ccccc}{\mathit{I}}_{c_{\mathit{z}}\times c_{\mathit{z}}}& & & & \\ & [{\mathit{O}}_{c_{{\mathit{z}}_p}\times c_{{\mathit{z}}_p}}|{\mathit{I}}_{c_{\mathit{z}}\times c_{\mathit{z}}}]& & & \\ & & [{\mathit{O}}_{c_{{\mathit{z}}_p}\times c_{{\mathit{z}}_p}}|{\mathit{I}}_{c_{\mathit{z}}\times c_{\mathit{z}}}]& & \\ & & & [{\mathit{O}}_{c_{{\mathit{z}}_p}\times c_{{\mathit{z}}_p}}|{\mathit{I}}_{c_{\mathit{z}}\times c_{\mathit{z}}}]& \\ & & & & \ddots \end{array}\right],$$

and $\mathit{r} = \mathrm{vec}(\widehat{\mathsf{\Theta }})$.

2.4. Example 1: Prediction of Mortality Data

We illustrate our constrained out-of-sample prediction method with the analysis of male mortality data in the US. Although mortality data are usually analyzed by Poisson distribution, for the simple purpose of illustrating the proposed method, we will assume that the logarithm of the mortality rates are normally distributed. The data come from the Human Mortality Database (Link: https://www.mortality.org/ (last accessed date: 1 July 2021). We observed the log of mortality rates from ages 0 to $110+$ over the period 1960–2014, and forecasted up to 2050, i.e., we obtained out-of-sample predictions in one of the two covariates, in this case, year.

The Lee–Carter method [16] is one of the most commonly used methods for estimating and forecasting mortality data; however, this method may produce an unwanted cross-over of forecasted mortality (higher mortality rates for younger ages than for older ages). The Lee-Carter model is defined as:

$$log(m_{x,y}) = {\alpha }_x + {\beta }_x{\kappa }_y + {ϵ}_y,$$

where $m_{x,y}$ is the central rate of mortality at age x and year y, ${\alpha }_x$, ${\beta }_x$ and ${\kappa }_y$ are parameters to be estimated, and ${ϵ}_y$ is the error term with mean zero and variance ${\sigma }^2$. This model is fitted to historical data and the resulting estimated ${\kappa }_y$s are then modelled and projected as a stochastic time series using standard Box–Jenkins methods. Ref. [17] proposed a modification to the Lee–Carter model based on the estimation of smooth versions of ${\beta }_x$ to avoid the cross-over of projected mortality rates. In order to compare our method with the solution given by [17], we obtained fits and predictions through four different models:

(a)

Model 1: an unconstrained model:

$${\mathit{y}}_{+} = f(\mathit{age},{\mathit{year}}_{+}) + {\mathit{ϵ}}_{+},{\mathit{ϵ}}_{+}\sim \mathcal{N}(\mathbf{0},{\mathit{R}}_{+})$$

where ${\mathit{y}}_{+}$ is the log mortality rate and ${\mathit{R}}_{+}$ is defined as in (10).

(b)

Model 2: The model defined in (a) subject to the constraint that the fit is maintained when fitting and forecasting at the same time.

(c)

Model 3: The model defined in (a) is subject to two constraints:

-

The fit is maintained;

-

The inter-age structure is preserved. We impose this constraint to avoid age crossing, that is, to prevent the death rate of younger ages from being higher than that of older ages. To achieve this, we use the pattern of the estimated coefficients in the past few years and forecast it. In particular, we ensure that the difference between the coefficients of every two consecutive projections is constant and equal to the difference between the corresponding last coefficient in the fitting.

Let us use an example to explain how to impose these two constraints at the same time. Suppose that the dimension of the coefficient matrix from the fitting is $4 \times 3$, and the dimensions of the coefficient matrix given the fitting and prediction is $4 \times 5$, that is,

$$\widehat{\mathsf{\Theta }} = \left[\begin{array}{ccc}{\widehat{\theta }}_1& {\widehat{\theta }}_5& {\widehat{\theta }}_9\\ {\widehat{\theta }}_2& {\widehat{\theta }}_6& {\widehat{\theta }}_{10}\\ {\widehat{\theta }}_3& {\widehat{\theta }}_7& {\widehat{\theta }}_{11}\\ {\widehat{\theta }}_4& {\widehat{\theta }}_8& {\widehat{\theta }}_{12}\end{array}\right],{\mathsf{\Theta }}_{+}^{*} = \left[\begin{array}{ccccc}{\stackrel{•}{\theta }}_1& {\stackrel{•}{\theta }}_5& {\stackrel{•}{\theta }}_9& {\stackrel{\circ }{\theta }}_{13}& {\stackrel{\circ }{\theta }}_{17}\\ {\stackrel{•}{\theta }}_2& {\stackrel{•}{\theta }}_6& {\stackrel{•}{\theta }}_{10}& {\stackrel{\circ }{\theta }}_{14}& {\stackrel{\circ }{\theta }}_{18}\\ {\stackrel{•}{\theta }}_3& {\stackrel{•}{\theta }}_7& {\stackrel{•}{\theta }}_{11}& {\stackrel{\circ }{\theta }}_{15}& {\stackrel{\circ }{\theta }}_{19}\\ {\stackrel{•}{\theta }}_4& {\stackrel{•}{\theta }}_8& {\stackrel{•}{\theta }}_{12}& {\stackrel{\circ }{\theta }}_{16}& {\stackrel{\circ }{\theta }}_{20}\end{array}\right],$$

in ${\mathsf{\Theta }}_{+}^{*}$, the coefficients that determine the fit have the superscript • and the coefficients that determine the forecast have the superscript ∘. Now, in the restriction Equation (25), ${\mathit{\theta }}_{+}^{*} = \mathrm{Vec}({\mathsf{\Theta }}_{+}^{*})$, and $\mathit{C} = \mathrm{blockdiag}({\mathit{I}}_{12\times 12},\mathit{U},\mathit{U})$ with

$$\mathit{U} = \left[\begin{array}{cccc}1& -1& 0& 0\\ 0& 1& -1& 0\\ 0& 0& 1& -1\end{array}\right]\mathrm{and}\mathit{r} = \left[\begin{array}{c}\widehat{\mathit{\theta }}\\ {\widehat{\theta }}_9-{\widehat{\theta }}_{10}\\ {\widehat{\theta }}_{10}-{\widehat{\theta }}_{11}\\ {\widehat{\theta }}_{11}-{\widehat{\theta }}_{12}\\ {\widehat{\theta }}_9-{\widehat{\theta }}_{10}\\ {\widehat{\theta }}_{10}-{\widehat{\theta }}_{11}\\ {\widehat{\theta }}_{11}-{\widehat{\theta }}_{12}\end{array}\right]$$

and $\widehat{\mathit{\theta }} = \mathrm{Vec}(\widehat{\mathsf{\Theta }})$, i.e., we are imposing that the coefficients that determine the fit (when fitting and forecasting simultaneously) have to be the ones we obtain when only fitting the data, and that the difference between the coefficients that determine the forecast of two consecutive rows has to be equal to the difference between the last coefficients from the fit of those rows.

(d)

Model 4: is the Lee–Carter model proposed in [16], but modified to get regular projected life tables [17].

Figure 1 shows the fit and the forecast obtained with all four models. In order to highlight the difference between the results, we have selected four ages: 20, 40, 60, and 80. Figure 2 shows the fit and forecast for these ages obtained by model 1 (red line), model 2 (green line), model 3 (blue line) and model 4 (orange line). The first three models (models 1, 2 and 3) provide almost the same fit. However, the fit given by model 4 is quite different and worse than the fits of other models, at 20 and 40 years. This is to be expected because the Lee–Carter model is not flexible enough to capture the increase in mortality in the early 1980s.

The predictions with models 1 and 2 are almost identical. Although it gave very similar results in the fitting, model 3 provided completely different results in the predictions, as expected. For 60 years of age, model 1 and model 2 predict that the logarithmic death rate will increase during 2020–2050, because they predict an upward trend in death rates from 2010 to 2014. This is inconsistent with the expected result. In the case of model 4, forecasts are close to model 3 for ages 40 and 80, but in the case of age 20 it seems to clearly overestimate the mortality. In model 3, the increasing trend is corrected after the initial years to maintain the structure of the adjacent ages, thus avoiding irregular predictions (as happens at age 60).

As we have seen, when we impose the two constraints, we obtain the most realistic results: maintain the fit and the cross-age structure, that is, those given by model 3. If we do not maintain the coefficient pattern, age crossover may occur. We further illustrate this fact in Figure 3, we plot the projections obtained with model 2 (green line) and model 3 (blue line) for ages 66 and 67. We can see that the fit of the two models is the same, and within the known data range, the log mortality of age 66 is lower than that of age 67. However, in the forecast, model 2 does not avoid the cross-over for ages 66 and 67 and so, in 2050, the log mortality rate is larger for age 66 than for age 47, this cross-over is not biologically acceptable since it is expected that mortality rates of older people are higher. This does not occur with model 3 where the restrictions imposed preserve the cross-age structure among ages.

3. Component-Wise Prediction with P-Spline Smooth-ANOVA Models

In some cases, the multidimensional model in (1) is not flexible enough, and it may result in unnecessary complex models (that is, use too many degrees of freedom to fit the data). In order to add more flexibility and avoid unnecessary and irrelevant terms, the authors of [8] propose the use of P-spline smooth-ANOVA models (or S-ANOVA). This method decomposes the sum of smooth functions in a similar way to the analysis of variance decomposition, i.e., main and interaction effects. In this section, we will show how to obtain out-of-sample predictions in this setting. First, we briefly review how to fit the S-ANOVA model and show how the P-spline mixed model representation can avoid the identifiability problem that occurs when multiple smoothing terms are used in the same model.

Suppose we have data laid out as a rectangular array, and $\mathit{y}$ is the response vector of length $n \times 1$, where $n = n_{\mathit{z}}{n}_{\mathit{x}}$, and $\mathit{z} = {(z_1,z_2,\dots z_{n_{\mathit{z}}})}^{\prime }$ and $\mathit{x} = {(x_1,x_2,\dots x_{n_{\mathit{x}}})}^{\prime }$ are the covariates. Let us consider the following S-ANOVA model:

$$\mathit{y} = \mathsf{\gamma } + f_1(\mathit{z}) + f_2(\mathit{x}) + f_{1,2}(\mathit{z},\mathit{x}) + \mathit{ϵ} = \mathit{B}\mathit{\theta } + \mathit{ϵ},\mathit{ϵ}\sim \mathcal{N}(\mathbf{0},\mathit{R})$$

(31)

where $\mathit{R} = {\sigma }_{\mathit{ϵ}}^2\mathit{I}$, i.e., the errors are independent and identically distributed and the B-spline basis $\mathit{B}$ is defined as:

$$\mathit{B} = \left[{\mathbf{1}}_n\right|{\mathbf{1}}_{n_{\mathit{x}}}\otimes {\mathit{B}}_{\mathit{z}}|{\mathit{B}}_{\mathit{x}}\otimes {\mathbf{1}}_{n_{\mathit{z}}}|{\mathit{B}}_{\mathit{x}}\otimes {\mathit{B}}_{\mathit{z}}],$$

(32)

of dimension $n \times (1 + c_{\mathit{z}} + c_{\mathit{x}} + c_{\mathit{z}}{c}_{\mathit{x}})$, and where ${\mathbf{1}}_{n_{\mathit{z}}}$ and ${\mathbf{1}}_{n_{\mathit{x}}}$ are column vectors of ones of length $n_{\mathit{z}}$ and $n_{\mathit{x}}$ respectively, and the vector of regression coefficients are $\mathit{\theta } = {(\mathsf{\gamma },{\mathit{\theta }}_{\mathit{z}}^{\prime },{\mathit{\theta }}_{\mathit{x}}^{\prime },{\mathit{\theta }}_{\mathit{s}}^{\prime })}^{\prime }$, where ${\mathit{\theta }}_{\mathit{z}}$ and ${\mathit{\theta }}_{\mathit{x}}$ are the vectors of coefficients for the main effects, of dimension $c_{\mathit{z}} \times 1$ and $c_{\mathit{x}} \times 1$, respectively, and ${\mathit{\theta }}_s$ is the vector of coefficients for the interaction, of dimension $c_{\mathit{z}}{c}_{\mathit{x}} \times 1$. Therefore, model (31) is written as:

$$\mathit{y} = \mathsf{\gamma }{\mathbf{1}}_n + ({\mathbf{1}}_{n_{\mathit{x}}}\otimes {\mathit{B}}_{\mathit{z}}){\mathit{\theta }}_{\mathit{z}} + ({\mathit{B}}_{\mathit{x}}\otimes {\mathbf{1}}_{n_{\mathit{z}}}){\mathit{\theta }}_{\mathit{x}} + ({\mathit{B}}_{\mathit{x}}\otimes {\mathit{B}}_{\mathit{z}}){\mathit{\theta }}_{\mathit{s}},$$

(33)

with penalty equal to the following block-diagonal matrix:

$$\mathit{P} = \left[\begin{array}{cccc}0& & & \\ & {\mathsf{\lambda }}_{\mathit{z}}{\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}}& & \\ & & {\mathsf{\lambda }}_{\mathit{x}}{\mathit{D}}_{\mathit{x}}^{\prime }{\mathit{D}}_{\mathit{x}}& \\ & & & {\tau }_{\mathit{z}}{\mathit{I}}_{c_{\mathit{x}}}\otimes {\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}}+{\tau }_{\mathit{x}}{\mathit{D}}_{\mathit{x}}^{\prime }{\mathit{D}}_{\mathit{x}}\otimes {\mathit{I}}_{c_{\mathit{z}}}\end{array}\right],$$

(34)

of dimension $(1 + c_{\mathit{z}} + c_{\mathit{x}} + c_{\mathit{z}}{c}_{\mathit{x}}) \times (1 + c_{\mathit{z}} + c_{\mathit{x}} + c_{\mathit{z}}{c}_{\mathit{x}})$, where each block corresponds to the penalty over each set of coefficients in the model. Ref. [15] pointed out that model (33) and penalty (34) should be modified in order to preserve identifiability. Their proposal is to construct identifiable bases and penalties by reparameterizing the model as a mixed model instead of imposing numerical constraints as suggested by [18]. Ref. [15] defines the following transformation matrix $\mathsf{\Omega } = \left[{\mathsf{\Omega }}_f\right|{\mathsf{\Omega }}_r]$ with dimension $(1 + c_{\mathit{z}} + c_{\mathit{x}} + c_{\mathit{z}}{c}_{\mathit{x}})$, where:

$${\mathsf{\Omega }}_f = \left[\begin{array}{cccc}1& & & \\ & 1\otimes {\mathit{u}}_{\mathit{z}f}^{(2)}& & \\ & & {\mathit{u}}_{\mathit{x}f}^{(2)}\otimes 1& \\ & & & {\mathit{u}}_{\mathit{x}f}^{(2)}\otimes {\mathit{u}}_{\mathit{z}f}^{(2)}\end{array}\right],$$

(35)

$${\mathsf{\Omega }}_r = \left[\begin{array}{cccc}1& & & \\ & 1\otimes {\mathit{U}}_{\mathit{z}r}& & \\ & & {\mathit{U}}_{\mathit{x}r}\otimes 1& \\ & & & {\mathit{u}}_{\mathit{x}f}^{(2)}\otimes {\mathit{U}}_{\mathit{z}r}|{\mathit{U}}_{\mathit{x}r}\otimes {\mathit{u}}_{\mathit{z}f}^{(2)}|{\mathit{U}}_{\mathit{x}r}\otimes {\mathit{U}}_{\mathit{z}r}\end{array}\right],$$

(36)

and ${\mathit{U}}_{\mathit{z}f}$ and ${\mathit{U}}_{\mathit{x}f}$ are the eigenvectors corresponding to the positive values of the SVD of the marginal penalties, ${\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}}$ and ${\mathit{D}}_{\mathit{x}}^{\prime }{\mathit{D}}_{\mathit{x}}$, and ${\mathit{u}}_{\mathit{z}f}^{(2)}$ and ${\mathit{u}}_{\mathit{x}f}^{(2)}$ are the second columns of ${\mathit{U}}_{\mathit{z}r}$ and ${\mathit{U}}_{\mathit{x}r}$, respectively.

Given the previous transformation matrix, ref. [15] showed that the fixed and random effects matrices $\mathit{X}$ and $\mathit{Z}$ are:

$$\begin{array}{ccc}\hfill \mathit{X}& =& \left[{\mathbf{1}}_n\right|{\mathbf{1}}_{n_{\mathit{x}}}\otimes \tilde{\mathit{z}}|\tilde{\mathit{x}}\otimes {\mathbf{1}}_{n_{\mathit{z}}}|\tilde{\mathit{x}}\otimes \tilde{\mathit{z}}]\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill \mathit{Z}& =& \left[{\mathbf{1}}_n\right|{\mathbf{1}}_{n_{\mathit{x}}}\otimes {\mathit{Z}}_{\mathit{z}}|{\mathit{Z}}_{\mathit{x}}\otimes {\mathbf{1}}_{n_{\mathit{z}}}|\tilde{\mathit{x}}\otimes {\mathit{Z}}_{\mathit{z}}|{\mathit{Z}}_{\mathit{x}}\otimes \tilde{\mathit{z}}|{\mathit{Z}}_{\mathit{x}}\otimes {\mathit{Z}}_{\mathit{z}}]\hfill \end{array}$$

(38)

where $\tilde{\mathit{z}} = {\mathit{B}}_{\mathit{z}}{\mathit{u}}_{\mathit{z}f}^{(2)}$, $\tilde{\mathit{x}} = {\mathit{B}}_{\mathit{x}}{\mathit{u}}_{\mathit{x}f}^{(2)}$, ${\mathit{Z}}_{\mathit{z}} = {\mathit{B}}_{\mathit{z}}{\mathit{U}}_{\mathit{z}r}$ and ${\mathit{Z}}_{\mathit{x}} = {\mathit{B}}_{\mathit{x}}{\mathit{U}}_{\mathit{x}r}$. Moreover, the mixed model penalty (or equivalently, the precision matrix of random effects) is:

$$\mathit{F} = \mathrm{blockdiag}({\mathit{F}}_{(1)},{\mathit{F}}_{(2)},{\mathit{F}}_{(1,2)}),$$

(39)

where for a second order penalty, it has size $(c_{\mathit{z}}{c}_{\mathit{x}} - 4) \times (c_{\mathit{z}}{c}_{\mathit{x}} - 4)$, and with:

$$\begin{array}{ccc}\hfill {\mathit{F}}_{(1)}& =& {\mathsf{\lambda }}_{\mathit{z}}{\tilde{\Sigma }}_{\mathit{z}},\hfill \\ \hfill {\mathit{F}}_{(2)}& =& {\mathsf{\lambda }}_{\mathit{x}}{\tilde{\Sigma }}_{\mathit{x}},\hfill \\ \hfill {\mathit{F}}_{(1,2)}& =& \mathrm{blockdiag}({\tau }_{\mathit{z}}{\tilde{\Sigma }}_{\mathit{z}},{\tau }_{\mathit{x}}{\tilde{\Sigma }}_{\mathit{x}},{\tau }_{\mathit{z}}{\mathit{I}}_{c_{\mathit{x}}-2}\otimes {\tilde{\Sigma }}_{\mathit{z}} + {\tau }_{\mathit{x}}{\tilde{\Sigma }}_{\mathit{x}}\otimes {\mathit{I}}_{c_{\mathit{z}}-2})\hfill \end{array}$$

where ${\tilde{\Sigma }}_{\mathit{z}}$ and ${\tilde{\Sigma }}_{\mathit{x}}$ are the nonzero eigenvalues of the SVD of ${\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}}$ and ${\mathit{D}}_{\mathit{x}}^{\prime }{\mathit{D}}_{\mathit{x}}$, respectively. Then, model (31) becomes:

$$\mathit{y} = \mathit{X}\mathit{\beta } + \mathit{Z}\mathit{\alpha } + \mathit{ϵ}\mathit{\alpha }\sim N\left(\mathbf{0},{\mathit{F}}^{-1}\right).$$

With the previous representation of model (31), ref. [15] avoid the identifiability problem removing the column vector of $\mathbf{1}$s in the random effects matrix (38). For more details, see [15].

Although the out-of-sample prediction will be performed in the context of a mixed model, we introduced the method in the original P-spline formulation because the reparameterization required for out-of-sample prediction will be based on this definition. In particular, we need to know the extended penalty matrix in order to calculate the precision matrix of random effects.

In model (31), given a vector of $n_{\mathit{z}}{n}_{\mathit{x}}$ observations $\mathit{y}$ of the response variable, suppose that we want to predict $n_p = n_{\mathit{z}}{n}_{{\mathit{x}}_p} + n_{{\mathit{z}}_p}{n}_{\mathit{x}} + n_{{\mathit{z}}_p}{n}_{{\mathit{x}}_p}$ new values at $(\mathit{z},{\mathit{x}}_p)$, $({\mathit{z}}_p,\mathit{x})$ and $({\mathit{z}}_p,{\mathit{x}}_p)$ of observed and predicted values can be arranged as in (20). Then, we consider the following extended smooth-ANOVA model:

$${\mathit{y}}_{+} = \mathsf{\gamma } + f_1({\mathit{z}}_{+}) + f_2({\mathit{x}}_{+}) + f_{1,2}({\mathit{z}}_{+},{\mathit{x}}_{+}) + {\mathit{ϵ}}_{+},{\mathit{ϵ}}_{+}\sim \mathcal{N}(\mathbf{0},{\sigma }_{\mathit{ϵ}}^2{\mathit{R}}_{+})$$

(40)

where ${\mathit{R}}_{+}$ is defined as in (10) with an extended B-spline basis ${\mathit{B}}_{+}$ defined as:

$${\mathit{B}}_{+} = [{\mathbf{1}}_{n_{+}}|{\mathbf{1}}_{n_{{\mathit{x}}_{+}}}\otimes {\mathit{B}}_{{\mathit{z}}_{+}}|{\mathit{B}}_{{\mathit{x}}_{+}}\otimes {\mathbf{1}}_{n_{{\mathit{z}}_{+}}}|{\mathit{B}}_{{\mathit{x}}_{+}}\otimes {\mathit{B}}_{{\mathit{z}}_{+}}],$$

(41)

of dimension $n \times (1 + c_{{\mathit{z}}_{+}} + c_{{\mathit{x}}_{+}} + c_{{\mathit{z}}_{+}}{c}_{{\mathit{x}}_{+}})$, and where ${\mathbf{1}}_{n_{{\mathit{z}}_{+}}}$ and ${\mathbf{1}}_{n_{{\mathit{x}}_{+}}}$ are column vectors of ones of length $n_{{\mathit{z}}_{+}}$ and $n_{{\mathit{x}}_{+}}$, respectively.

Therefore, model (40) is written as:

$${\mathit{y}}_{+} = \mathsf{\gamma }{\mathbf{1}}_{n_{+}} + ({\mathbf{1}}_{n_{{\mathit{x}}_{+}}}\otimes {\mathit{B}}_{{\mathit{z}}_{+}}){\mathit{\theta }}_{{\mathit{z}}_{+}} + ({\mathit{B}}_{{\mathit{x}}_{+}}\otimes {\mathbf{1}}_{n_{{\mathit{z}}_{+}}}){\mathit{\theta }}_{{\mathit{x}}_{+}} + ({\mathit{B}}_{{\mathit{x}}_{+}}\otimes {\mathit{B}}_{{\mathit{z}}_{+}}){\mathit{\theta }}_{{\mathit{s}}_{+}},$$

(42)

with extended penalty the following block-diagonal matrix:

$${\mathit{P}}_{+} = \left[\begin{array}{cccc}0& & & \\ & {\mathsf{\lambda }}_{\mathit{z}}{\mathit{D}}_{{\mathit{z}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{z}}_{+}}& & \\ & & {\mathsf{\lambda }}_{\mathit{x}}{\mathit{D}}_{{\mathit{x}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{x}}_{+}}& \\ & & & {\tau }_{\mathit{z}}{\mathit{I}}_{c_{{\mathit{x}}_{+}}}\otimes {\mathit{D}}_{{\mathit{z}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{z}}_{+}}+{\tau }_{\mathit{x}}{\mathit{D}}_{{\mathit{x}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{x}}_{+}}\otimes {\mathit{I}}_{c_{{\mathit{z}}_{+}}}\end{array}\right],$$

(43)

of dimension $(1 + c_{{\mathit{z}}_{+}} + c_{{\mathit{x}}_{+}} + c_{{\mathit{z}}_{+}}{c}_{{\mathit{x}}_{+}}) \times (1 + c_{{\mathit{z}}_{+}} + c_{{\mathit{x}}_{+}} + c_{{\mathit{z}}_{+}}{c}_{{\mathit{x}}_{+}})$, where each block corresponds to the penalty over each of the coefficients of the model.

In order to reformulate the extended S-ANOVA model (40) as a mixed model, we need to define an extended transformation matrix. The obvious choice would be to use the natural extended transformation matrix based on the SVD of the extended difference matrices; however, this method cannot ensure the invariance of the fit when performing out-of-sample predictions. In the following section, we propose the use of an extended transformation matrix, ${\mathsf{\Omega }}_{+}^{*}$, that ensures coherent predictions.

3.1. Coherent Prediction with S-ANOVA Model

In order to use the S-ANOVA model for prediction under the constraint that the fit must be maintained, an extended transformation matrix that preserves the model matrix must be used. When the two covariates are extended, we define the following extended transformation matrix ${\mathsf{\Omega }}_{+}^{*} = [{\mathsf{\Omega }}_{+f}^{*}|{\mathsf{\Omega }}_{+r}^{*}]$ with dimension $(1 + c_{{\mathit{z}}_{+}} + c_{{\mathit{x}}_{+}} + c_{{\mathit{z}}_{+}}{c}_{{\mathit{x}}_{+}})$, and:

$${\mathsf{\Omega }}_{+f}^{*} = \left[\begin{array}{cccc}1& & & \\ & 1\otimes {\mathit{u}}_{{\mathit{z}}_{+}f}^{*(2)}& & \\ & & {\mathit{u}}_{{\mathit{x}}_{+}f}^{*(2)}\otimes 1& \\ & & & {\mathit{u}}_{{\mathit{x}}_{+}f}^{*(2)}\otimes {\mathit{u}}_{{\mathit{z}}_{+}f}^{*(2)}\end{array}\right],$$

(44)

$${\mathsf{\Omega }}_{+r}^{*} = \left[\begin{array}{cccc}1& & & \\ & 1\otimes {\mathit{U}}_{{\mathit{z}}_{+}r}^{*}& & \\ & & {\mathit{U}}_{{\mathit{x}}_{+}r}^{*}\otimes 1& \\ & & & {\mathit{u}}_{{\mathit{x}}_{+}f}^{*(2)}\otimes {\mathit{U}}_{{\mathit{z}}_{+}r}^{*}|{\mathit{U}}_{{\mathit{x}}_{+}r}^{*}\otimes {\mathit{u}}_{{\mathit{z}}_{+}f}^{*(2)}|{\mathit{U}}_{{\mathit{x}}_{+}r}^{*}\otimes {\mathit{U}}_{{\mathit{z}}_{+}r}^{*}\end{array}\right].$$

(45)

where ${\mathit{u}}_{{\mathit{z}}_{+}f}^{*(2)}$ and ${\mathit{u}}_{{\mathit{x}}_{+}f}^{*(2)}$ are the second columns of ${\mathit{U}}_{{\mathit{z}}_{+}f}^{*}$ and ${\mathit{U}}_{{\mathit{x}}_{+}f}^{*}$, respectively,

$${\mathit{U}}_{{\mathit{z}}_{+}f}^{*} = \left[\begin{array}{c}{\mathit{U}}_{\mathit{z}f}\\ -{\mathit{D}}_{\mathit{z}(2)}^{-1}{\mathit{D}}_{\mathit{z}(1)}{\mathit{U}}_{\mathit{z}f}\end{array}\right],\mathrm{and}$$

(46)

${\mathit{U}}_{\mathit{z}f}$ are the eigenvectors corresponding to the positive values of the SVD of ${\mathit{D}}_{\mathit{z}}^{\prime }{\mathit{D}}_{\mathit{z}}$ and ${\mathit{D}}_{\mathit{z}(2)}$ and ${\mathit{D}}_{\mathit{z}(1)}$ are blocks of (14). ${\mathit{U}}_{{\mathit{z}}_{+}f}^{*}$ ${\mathit{U}}_{{\mathit{x}}_{+}f}^{*}$, ${\mathit{U}}_{\mathit{x}f}$, ${\mathit{D}}_{\mathit{x}(2)}$ and ${\mathit{D}}_{\mathit{x}(1)}$ have analogous definitions. Then, the fixed effects matrix is:

$${\mathit{X}}_{+}^{*} = {\mathit{B}}_{+}{\mathsf{\Omega }}_{+f}^{*} = [{\mathbf{1}}_{n_{+}}|{\mathbf{1}}_{n_{{\mathit{x}}_{+}}}\otimes {\tilde{\mathit{z}}}_{+}|{\tilde{\mathit{x}}}_{+}\otimes {\mathbf{1}}_{n_{{\mathit{z}}_{+}}}|{\tilde{\mathit{x}}}_{+}\otimes {\tilde{\mathit{z}}}_{+}],$$

(47)

where ${\tilde{\mathit{z}}}_{+}$ and ${\tilde{\mathit{x}}}_{+}$ are ${\mathit{B}}_{{\mathit{z}}_{+}}{\mathit{u}}_{\mathit{z}f}^{*(2)}$ and ${\mathit{B}}_{{\mathit{x}}_{+}}{\mathit{u}}_{{\mathit{x}}_{+}f}^{*(2)}$, respectively, and the random effects matrix is:

$${\mathit{Z}}_{+}^{*} = {\mathit{B}}_{+}{\mathsf{\Omega }}_{+r}^{*} = [{\mathbf{1}}_{n_{+}}|{\mathbf{1}}_{n_{{\mathit{x}}_{+}}}\otimes {\mathit{Z}}_{{\mathit{z}}_{+}}^{*}|{\mathit{Z}}_{{\mathit{x}}_{+}}^{*}\otimes {\mathbf{1}}_{n_{{\mathit{z}}_{+}}}|{\tilde{\mathit{x}}}_{+}\otimes {\mathit{Z}}_{{\mathit{z}}_{+}}^{*}|{\mathit{Z}}_{{\mathit{x}}_{+}}^{*}\otimes {\tilde{\mathit{z}}}_{+}|{\mathit{Z}}_{{\mathit{x}}_{+}}^{*}\otimes {\mathit{Z}}_{{\mathit{z}}_{+}}^{*}],$$

(48)

where ${\mathit{Z}}_{{\mathit{z}}_{+}}^{*}$ and ${\mathit{Z}}_{{\mathit{x}}_{+}}^{*}$ are ${\mathit{B}}_{{\mathit{z}}_{+}}{\mathit{U}}_{{\mathit{z}}_{+}r}^{*}$ and ${\mathit{B}}_{{\mathit{x}}_{+}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*}$, respectively. The covariance matrix of random effects is ${\mathit{G}}_{+}^{*} = {\sigma }_{\mathit{ϵ}}^2{\mathit{F}}_{+}^{*-1}$, with ${\mathit{F}}_{+}^{*}$ given by the following theorem.

Theorem 2.

The extended precision matrix of the random effects for S-ANOVA model in (40) with extended transformation matrix given in (45), and extended penalty given in (43), is a block-diagonal matrix given by:

$${\mathit{F}}_{+}^{*} = \mathrm{blockdiag}({\mathit{F}}^{(1)},{\mathit{F}}_{+}^{(2)},{\mathit{F}}_{+}^{(1,2)}),$$

(49)

where:

$$\begin{array}{ccc}\hfill {\mathit{F}}^{(1)}& =& {\mathsf{\lambda }}_{\mathit{z}}{\mathit{U}}_{{\mathit{z}}_{+}r}^{*\prime }{\mathit{D}}_{{\mathit{z}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{z}}_{+}}{\mathit{U}}_{{\mathit{z}}_{+}r}^{*},\hfill \\ \hfill {\mathit{F}}_{+}^{(2)}& =& {\mathsf{\lambda }}_{\mathit{x}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*\prime }{\mathit{D}}_{{\mathit{x}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{x}}_{+}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*},\hfill \\ \hfill {\mathit{F}}_{+}^{(1,2)}& =& \left[\begin{array}{ccc}{\mathit{F}}_{+11}^{(1,2)}& \mathit{O}& {\mathit{F}}_{+13}^{(1,2)}\\ \mathit{O}& {\mathit{F}}_{+22}^{(1,2)}& {\mathit{F}}_{+23}^{(1,2)}\\ {\mathit{F}}_{+13}^{(1,2)\prime }& {\mathit{F}}_{+23}^{(1,2)\prime }& {\mathit{F}}_{+33}^{(1,2)}\end{array}\right],\hfill \end{array}$$

with

$$\begin{array}{ccc}\hfill {\mathit{F}}_{+11}^{(1,2)}& =& {\tau }_{\mathit{z}}{\mathit{u}}_{{\mathit{x}}_{+}f}^{*(2)\prime }{\mathit{u}}_{{\mathit{x}}_{+}f}^{*(2)}\otimes {\mathit{U}}_{{\mathit{z}}_{+}r}^{*\prime }{\mathit{D}}_{{\mathit{z}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{z}}_{+}}{\mathit{U}}_{{\mathit{z}}_{+}r}^{*}\hfill \\ \hfill {\mathit{F}}_{+13}^{(1,2)}& =& {\tau }_{\mathit{z}}{\mathit{u}}_{{\mathit{x}}_{+}f}^{*(2)\prime }{\mathit{U}}_{{\mathit{x}}_{+}r}^{*}\otimes {\mathit{U}}_{{\mathit{z}}_{+}r}^{*\prime }{\mathit{D}}_{{\mathit{z}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{z}}_{+}}{\mathit{U}}_{{\mathit{z}}_{+}r}^{*},\hfill \\ \hfill {\mathit{F}}_{+22}^{(1,2)}& =& {\tau }_{\mathit{x}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*\prime }{\mathit{D}}_{{\mathit{x}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{x}}_{+}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*}\otimes {\mathit{u}}_{{\mathit{z}}_{+}r}^{*(2)\prime }{\mathit{u}}_{{\mathit{z}}_{+}r}^{*(2)},\hfill \\ \hfill {\mathit{F}}_{+23}^{(1,2)}& =& {\tau }_{\mathit{x}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*\prime }{\mathit{D}}_{{\mathit{x}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{x}}_{+}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*}\otimes {\mathit{u}}_{{\mathit{z}}_{+}f}^{*(2)\prime }{\mathit{U}}_{{\mathit{z}}_{+}r}^{*}.\hfill \\ \hfill {\mathit{F}}_{+33}^{(1,2)}& =& {\tau }_{\mathit{z}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*\prime }{\mathit{U}}_{{\mathit{x}}_{+}r}^{*}\otimes {\mathit{U}}_{{\mathit{z}}_{+}r}^{*\prime }{\mathit{D}}_{{\mathit{z}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{z}}_{+}}{\mathit{U}}_{{\mathit{z}}_{+}r}^{*} + {\tau }_{\mathit{x}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*\prime }{\mathit{D}}_{{\mathit{x}}_{+}}^{\prime }{\mathit{D}}_{{\mathit{x}}_{+}}{\mathit{U}}_{{\mathit{x}}_{+}r}^{*}\otimes {\mathit{U}}_{{\mathit{z}}_{+}r}^{*\prime }{\mathit{U}}_{{\mathit{z}}_{+}r}^{*}.\hfill \end{array}$$

The proof is given in the Supplementary Materials.

3.2. Example 2: Prediction of Aboveground Biomass

In this section, we apply the proposed 2D interaction P-spline and S-ANOVA models to another real data sets. The data come from an agricultural experiment conducted in Spain [19] to evaluate the economic viability of Populus trees prior to harvesting. In this context, it is important to estimate the aboveground biomass and obtain accurate predictions using only a minimal set of readily available information, i.e., diameter and height (denoted by $\mathit{z}$ and $\mathit{x}$). Ref. [20] proposed the use of smooth additive mixed models for predicting aboveground biomass. Here, we analyzed the data of one of the nine clones included in the experiment. The observation data consist of 315 observations for the diameter value measured at a breast height of 1.30 m. The data are shown in Figure 4, where it can be seen that the weight is a nonlinear function of diameter and height. We propose predicting out-of-sample aboveground biomass as a smooth function of height and diameter using the following two models:

• 2D interaction P-spline model:

  $${\mathit{y}}_{+} = f({\mathit{z}}_{+},{\mathit{x}}_{+}) + {\mathit{ϵ}}_{+},{\mathit{ϵ}}_{+}\sim \mathcal{N}(\mathbf{0},{\mathit{R}}_{+})$$
• Smooth-ANOVA model:

  $${\mathit{y}}_{+} = f({\mathit{z}}_{+}) + f({\mathit{x}}_{+}) + f({\mathit{z}}_{+},{\mathit{x}}_{+}) + {\mathit{ϵ}}_{+},{\mathit{ϵ}}_{+}\sim \mathcal{N}(\mathbf{0},{\mathit{R}}_{+})$$

We have predicted weights for 10 new out-of-sample values for diameter and height. In Figure 5, we show the estimated smooth main effects for height (left panel) and diameter (right panel) obtained after fitting and predicting with the S-ANOVA model imposing that the fit is maintained. As shown, the effect of height is weaker than diameter, but both smooth terms are significantly different from zero. Figure 6 shows the fitted and predicted interaction function ($f({\mathit{z}}_{+},{\mathit{x}}_{+})$) for the restricted S-ANOVA model.

Figure 7 shows the predictions when the fit and the forecast are obtained with the S-ANOVA model (top left panel), with the S-ANOVA model subject to the constraint that the fit is maintained (top right panel), with the 2D interaction P-spline model (bottom left panel) and with the 2D interaction P-spline models imposing that the fit is maintained (bottom right panel). As shown in the figure, the results obtained by the S-ANOVA model and the restricted S-ANOVA model are almost equal. However, the solution obtained from the 2D interaction P-spline model is different depending on whether the constraint of maintaining the fit is imposed. However, the solutions obtained from the 2D interaction P-spline models are different depending on whether the restriction the fit is maintained is imposed or not. In particular, if no constraints are imposed, the fit will change significantly.

In the prediction, the most significant difference is that weight values predicted by the 2D interaction P-spline model at the largest diameter and height values are lower than the ones obtained by the S-ANOVA model. It is important to note that, in the fit, the degrees of freedom are 32 for the 2D interaction P-spline model and 12 for the S-ANOVA model, i.e., with the 2D interaction P-spline model, we are adding unnecessary complexity.

Furthermore, comparing the S-ANOVA model and the 2D interaction P-spline model, we conclude that the S-ANOVA model gives the most coherent solution because it provides the largest weight value for the highest values of height and diameter.

4. Simulation Study

We propose two simulation studies to examine the performance of the interaction 2D model and smooth-ANOVA model compared to their constrained versions (imposing the invariance of the fit when performing out-of-sample predictions). For this propose we have simulated the data in two different scenarios:

(a)

Scenario 1. From an interaction model:

$${\mathit{S}}^1 = f_{1,2}(\mathit{z},\mathit{x}).$$

(b)

Scenario 2. From a two main effects with interaction model:

$${\mathit{S}}^2 = f_1(\mathit{z}) + f_2(\mathit{x}) + f_{1,2}(\mathit{z},\mathit{x}).$$

where $f_1(\mathit{z}) = \sin (2\pi \mathit{z})$, $f_2(\mathit{x}) = \cos (3\pi \mathit{x})$ and $f_{1,2}(\mathit{z},\mathit{x}) = 3\sin (2\pi \mathit{z})(2\mathit{x} - 1)$. To simulate the data, we have generated a grid of 4900 values. Both covariates, $\mathit{z}$ and $\mathit{x}$, take 70 equidistant values in the interval $[0,1]$, and the errors are i.i.d., with mean 0 and variance ${\sigma }_{\mathit{ϵ}}^2 = 0.25$. Figure 8 shows $f_1(\mathit{z})$, $f_2(\mathit{x})$ and the surfaces proposed in scenarios 1 and 2.

The models compared are:

• 2D interaction P-spline model, i.e.,

  $$\begin{array}{ccc}\hfill {\mathit{y}}_{+}& =& f({\mathit{z}}_{+},{\mathit{x}}_{+}) + {\mathit{ϵ}}_{+}{\mathit{ϵ}}_{+}\sim \mathcal{N}(\mathbf{0},{\sigma }_{\mathit{ϵ}}^2{\mathit{R}}_{+}).\hfill \end{array}$$
• Smooth-ANOVA model, i.e.,

  $$\begin{array}{ccc}\hfill {\mathit{y}}_{+}& =& f_1({\mathit{z}}_{+}) + f_2({\mathit{x}}_{+}) + f_{1,2}({\mathit{z}}_{+},{\mathit{x}}_{+}) + {\mathit{ϵ}}_{+}{\mathit{ϵ}}_{+}\sim \mathcal{N}(\mathbf{0},{\sigma }_{\mathit{ϵ}}^2)\hfill \end{array}$$
• Constrained 2D interaction P-spline model;
• Constrained smooth-ANOVA model.

The constraint imposed is that the coefficients of the fit remain invariant when out-of-sample-prediction is performed. In all cases, smoothing parameters are estimated using REML.

For each scenario and each model, we repeated the following steps 100 times:

(1)

Start with a data set of observations arranged into a matrix ${\mathit{Y}}_{+}$ of dimension $n_{{\mathit{z}}_{+}} \times n_{{\mathit{x}}_{+}}$:

$${\mathit{Y}}_{+} = \left[\begin{array}{cc}\mathit{Y}& {\mathit{Y}}_{\mathit{z}{\mathit{x}}_p}\\ {\mathit{Y}}_{{\mathit{z}}_p\mathit{x}}& {\mathit{Y}}_{{\mathit{z}}_p{\mathit{x}}_p}\end{array}\right],$$

where $\mathit{Y}$, ${\mathit{Y}}_{\mathit{z}{\mathit{x}}_p}$, ${\mathit{Y}}_{{\mathit{z}}_p\mathit{x}}$ and ${\mathit{Y}}_{{\mathit{z}}_p{\mathit{x}}_p}$ have dimensions $n_{\mathit{z}} \times n_{\mathit{x}}$, $n_{\mathit{z}} \times n_{{\mathit{x}}_p}$, $n_{{\mathit{z}}_p} \times n_{\mathit{x}}$ and $n_{{\mathit{z}}_p} \times n_{{\mathit{x}}_p}$, respectively;

(2)

Split our data into two groups: the training data $\mathit{Y}$ and the test data ${\mathit{Y}}_{\mathit{z}{\mathit{x}}_p}$, ${\mathit{Y}}_{{\mathit{z}}_p\mathit{x}}$ and ${\mathit{Y}}_{{\mathit{z}}_p{\mathit{x}}_p}$;

(3)

Use the training data to predict $4900 - n$ new observations ($n = n_{\mathit{z}} \times n_{\mathit{x}}$);

(4)

Check the accuracy of each model.

The marginal B-spline bases ${\mathit{B}}_{\mathit{z}}$ and ${\mathit{B}}_{\mathit{x}}$ are constructed with 15 knots and cubic splines and then extended to cover the whole range of ${\mathit{z}}_{+}$ and ${\mathit{x}}_{+}$, the penalty orders are two.

As is known, a model that fits the data well does not necessarily yields good predictions, so we compared the fit performance, the forecast performance, and the overall method performance. In order to check the accuracy of the methods proposed, we followed [21] and used the difference between the function values from which we simulated the data and the fit and forecast produced using only the data in the training set:

$${\mathit{E}}_{+} = {\mathit{S}}^k - {\widehat{\mathit{Y}}}_{+},k = 1,2$$

i.e, we have the errors matrix ${\mathit{E}}_{+} = \left[\begin{array}{cc}\mathit{E}& {\mathit{E}}_{\mathit{z}{\mathit{x}}_p}\\ {\mathit{E}}_{{\mathit{z}}_p\mathit{x}}& {\mathit{E}}_{{\mathit{z}}_p{\mathit{x}}_p}\end{array}\right],$ and therefore the vectors containing the errors in the fit, in the forecast and in the overall performance:

• Vector with the fit errors: ${\mathit{e}}^{(f)} = \mathrm{vec}(\mathit{E})$;
• Vector with the prediction errors: ${\mathit{e}}^{(p)} = (\mathrm{vec}({\mathit{E}}_{\mathit{z}{\mathit{x}}_p}),\mathrm{vec}({\mathit{E}}_{{\mathit{z}}_p\mathit{x}}),\mathrm{vec}({\mathit{E}}_{{\mathit{z}}_p{\mathit{x}}_p}))$;
• Vector with the total errors: ${\mathit{e}}^{(t)} = \mathrm{vec}({\mathit{E}}_{+})$.

The errors measure that we use is the mean absolute error because, as it is concluded in [21], it is less sensitive to outliers than the root mean square error.

Mean absolute error: MAE = $\frac{{\sum }_{i=1}^N|e_i^{(l)}|}{N},N = \mathrm{length}({\mathit{e}}^{(l)})$ for $l = f,p,t$, i.e., for the errors in the fit, in the forecast or in total.

4.1. Simulations Results for Scenario 1

Here, we present results from Scenario 1. We use different ranges for the out of sample prediction, $n_{{\mathit{x}}_p}/n_{{\mathit{z}}_p} = 0,5,10,15,20$ and fit and predict simultaneously with the four models described above. Note that for the different values of $n_{{\mathit{z}}_p}$ and $n_{{\mathit{x}}_p}$ we consider as the training data an observations matrix of dimension $(70 - n_{{\mathit{z}}_p}) \times (70 - n_{{\mathit{x}}_p})$ and we predict $(70 - n_{{\mathit{z}}_p}) \times n_{{\mathit{x}}_p} + n_{{\mathit{z}}_p} \times (70 - n_{{\mathit{x}}_p}) + n_{{\mathit{z}}_p} \times n_{{\mathit{x}}_p}$ new values.

In Figure 9, we report MAE values from the fit, the forecast, and the overall performance. The x-axis label of boxplots refers to the model: unrestricted(2D interaction P-spline model), restricted (2D interaction P-spline model with constraints), ANOVA (S-ANOVA Model) and restricted S-ANOVA (S-ANOVA model with constraints).

Plots for the different values of $n_{{\mathit{z}}_p}$ and $n_{{\mathit{x}}_p}$ are available in the Supplementary Materials, and below we show the particular case of $n_{{\mathit{z}}_p} = 0$ and $n_{{\mathit{x}}_p} = 5$.

In most cases, the fit remains relatively invariant in the case of the 2D P-spline models, with and without constraints, and in the restricted S-ANOVA model, while the fit in unrestricted S-ANOVA model is modified in most scenarios.

In the prediction, the three models behave similarly; however, the performance of the restricted S-ANOVA model is a bit worse and its MAE error values increase with an increasing prediction horizon (see Figures S1 and S2 of the Supplementary Materials).

Based on this simulated scenario, we conclude that simultaneous fit and prediction with S-ANOVA models produces significant changes in the fit and that in the overall performance (including fit and prediction), the 2D P-spline models (with and without restriction) and the restricted S-ANOVA models are the most accurate.

4.2. Simulations Results for Scenario 2

In the Supplementary Materials (Figure S12), we show the results obtained for Scenario 2, i.e., in the case in which the true surface is constructed from a model with two main effects and an interaction. Figure 10 shows the boxplots for the particular case of $n_{{\mathit{z}}_p} = 0$ and $n_{{\mathit{x}}_p} = 5$.

In this case, the conclusion is that there is a significant difference between the results obtained with the interaction models and the S-ANOVA models: two-dimensional P-spline models (restricted and unrestricted) are always less accurate than the S-ANOVA models. This is due to the fact that 2D interaction models try to fit the data without considering the main effects.

When fitting the observed data, the performance of 2D P-spline interaction models, with and without constraints, is quite similar. However, the restricted interaction model is better than the unrestricted model in the out-of-sample prediction, the difference between the two models is apparent for most scenarios.

In the case of the S-ANOVA models, the restricted model is always better than the unrestricted one, in the fit, in the prediction and, therefore, in the overall performance. Although, in the particular case in which we only predict in one dimension, the prediction errors for both models are almost equal (see Supplementary Materials, scenario 2). The conclusion for this scenario is that the restricted S-ANOVA model is clearly the most accurate.

Therefore, considering all test scenarios, our recommendation is that, as a general rule, restricted S-ANOVA models should be used to fit and predict models that require interaction terms.

5. Conclusions

In this article, we propose a general framework for out-of-sample prediction in a smooth additive model with interaction terms. We build our proposal from the method proposed in [10], and then we extend their method to situations where out-of-sample predictions are required in both directions of the interaction term. The method proposed in [10] considers out-of-sample predictions as missing values with a weight of 0, and performs fitting and prediction at the same time. We have shown that this method produces different fitted values, depending on whether we only fit or fit and predict simultaneously.

In order to solve this inconsistency, we recommend maximizing the penalized likelihood under linear constraints to ensure that the coefficients obtained in the data fitting remain unchanged when the fitting and out-of-sample prediction are performed at the same time. For this purpose, we use Lagrange multipliers. This general method can be used to impose other related constraints, such as mortality forecasting when the cross-age structure needs to be preserved.

The methodology proposed has been extended to the case of smooth mixed models since it will allow us to predict out-of-sample in a wide range of models. We have shown that care should be taken when formulating the models as mixed models, since the matrices of fixed and random effects for out-of-sample prediction need to be direct extensions of the matrices used in the fit. In this context, we have developed a method for out-ofsample prediction for the S-ANOVA model proposed in [8]. This is an important step forward since, until now, there was no general method available for prediction in models that include main effects and interactions terms. A simulation study has been carried out to compare restricted and unrestricted out-of-sample predictions when using 2D interaction models and S-ANOVA models. From the results of the simulation study, we conclude that, in most cases, when the prediction is made on the two dimensions of the interaction term, the restricted S-ANOVA model performs better in fitting and out-of-sample prediction. The results for only one of the covariates depends on the simulation scenario, although in most cases, the S-ANOVA model is better than the full interaction (restricted or unrestricted) model.

The multidimensional prediction P-splines method with constraints proposed has also been used in two real data examples, one in which mortality rates are forecasted over the years, and the importance of imposing constraints on the coefficients to ensure coherent forecast is shown. The other example predicts tree biomass as a function of the tree height and diameter using the S-ANOVA model; this example shows the ability of this model to produce separate out-of-sample predictions for main effects and interaction terms in the model. Although here, due to the nature of the data used, we have only shown examples of forward out-of-sample prediction, our methodology allows for backwards out-of-sample predictions.

Currently, the work presented here is being extended to non-Gaussian responses in the context of generalized additive and mixed models.