A general framework for prediction in penalized regression

Abstract

There are two main approaches to carrying out prediction in the context of penalized regression: with low-rank basis and penalties or through the smooth mixed models. In this article, we give further insight in the case of P-splines showing the influence of the penalty on the prediction. In the context of mixed models, we can connect the new predicted values to the observed values through a joint normal distribution, which allows us to compute prediction intervals. In this work, we propose an alternative approach, called the extended mixed model approach, that allows us to fit and predict data simultaneously. The methodology is illustrated with two real datasets, one of them on aboveground biomass and the other on monthly sulphur dioxide ( ${SO}_{2}$ ) levels in a selection of monitoring sites in Europe.

Keywords

mixed models Penalized regression prediction P-splines smooth models

1 Introduction

There are many situations in which prediction of new observations in the context of regression is needed; in particular, when ‘out-of-range’ prediction is required, that is, beyond the range of observed covariates. This problem extends in the framework of smoothing, that is, for models where the regression function is a smooth but otherwise unspecified function. Examples are numerous; for instance, hourly temperatures at a weather station or the yearly number of deaths, where the latter has a major impact in areas such as demography (mortality tables). A scatterplot of data often exhibits patterns, such as an upward or downward trend, or a pattern that repeats, seasonal variation, both of which might be used to predict new values.

The generalities of the problem and the particular characteristics of datasets are one of the main reasons that encourages us to work in the prediction field and base our work on the forecasting method proposed in Currie et al. (2004). They have shown how the method of penalized splines (P-splines), introduced by Eilers and Marx (1996) and extensively discussed in Ruppert et al. (2003), can be extended to smooth and predict two-dimensional mortality tables. In particular, the authors show how to construct the appropriate regression bases and penalty matrices for forecasting.

Most of the existing literature in the area is related to the prediction of new observations in a temporal context, that is, forecast of new observations. Although we have a more general approach, we start by providing a brief review of the main literature related to forecasting in smoothing models by commenting the main approaches of Ba et al. (2012) and Sacks et al. (1989). We also refer to Hyndman et al. (2008) who gave an overview of exponential smoothing methods.

Exponential smoothing or weighted smoothing, respectively, refers to a class of forecasting methods, each of them having the property that forecasts are weighted combinations of past observations, with recent observations given relatively more weight than older observations. Hyndman et al. (2008) provide extensive information about exponential smoothing methods. A summary of the exponential smoothing history shows the great importance of this method. It certainly is the most popular forecasting method used in business and industry since the 1950s.

In this article, we focus on computationally driven methods. A recent strategy is to use penalized splines to fit and forecast time series data (Ba et al., 2012). In this case, one minimizes the penalized least squares criteria:

S = (y - B θ)^{'} M (y - B θ) + λ θ^{'} P θ,

where $y$ is the vector of observed responses, $θ$ is basis coefficient, $B$ is a spline basis that covers the hole range of the explanatory variable and $M$ is a weight matrix that puts exponentially decreasing weights on the samples, according to the order of their arrival. Matrix $P$ is a penalty matrix controlling the smoothness of the fitted function and $λ$ is the smoothing parameter. Ba et al. (2012) proposed a data-driven choice for the exponential weight matrix.

A similar prediction problem is tackled in the framework of global optimization where the interest is to evaluate an unknown function at point $x$ , say. The question is now where to place future values of $x$ to evaluate the function, such that relevant (preferably most) information about the function is achieved. Exemplary, we refer to Sacks et al. (1989) and Jones et al. (1998), who fit a stochastic process to data and predict the process at a new point, given the already observed data. They treat the observations as if they were generated by a constant and an error component which is modelled as a stochastic process. Their approach is called Bayesian global optimization, and the concept is the same as the idea behind the well-known technique in spatial statistics called kriging (Cressie, 1993). Prediction can be addressed in a Bayesian setting using Bayesian P-splines by exploiting the properties of the random walk prior (Besag et al., 1995; Rue and Held, 2005). Computation of the predicted pattern and its credible intervals can be performed via Markov chain Monte Carlo (MCMC).

Beyond the aforementioned articles, we have not found methods to predict in penalized regression. In this article, we propose a general framework for prediction in penalized regression and also delve deeper into our knowledge of the forecasting method proposed in Currie et al. (2004). Although we will use B-spline basis and penalties based on differences, the methodology proposed here can be extended to any basis and quadratic penalty. We have organized the remaining of the article as follows. Section 2 is dedicated to introducing the fundamentals of the method proposed by Currie et al. (2004) and to show some properties that relate the order of the penalty with the shape of the predictions. In Section 3, we describe how forecasting can be carried out in the context of smooth mixed models. It can be done through two-stages procedures, following Gilmour et al. (2004) or based on the conditional distribution of the new values given the observed data. We show how the last procedure allows us to compute prediction intervals (see Section 3.1). In addition, we propose an alternative method for prediction with smooth mixed models that can be done through a one-stage procedure, by extending the proposal of Currie et al. (2004) to the mixed model framework (see Section 3.2). The equivalence of two-stage and one-stage procedures (in terms of predicted values and choice of smoothing parameter) is shown in the particular case of penalties based on the differences between adjacent coefficients (Eilers and Marx, 1996). The proposed methodology is illustrated in Section 4 with the analysis of two real datasets. Finally, concluding remarks are made in Section 5.

2 Prediction with smooth models and quadratic penalties

We start our presentation by giving a brief revision of penalized regression. Consider the case of univariate data with response variable $y$ and regressor $x$ . The smooth model is of the form

y = f (x) + ε, ε \sim N (0, σ_{ε}^{2} I),

(2.1)

where $f (\cdot)$ is an unknown smooth function, subsequently also called trend component. We assume to have the data points $(x_{i}, y_{i})$ drawn from the model, for $i = 1, . . ., n$ , and $ε_{i}$ being independent and identically distributed errors with variance $σ_{ε}^{2}$ . In order to estimate $f (x)$ , we replace the function by a basis representation, i.e. we assume $f (x) = \sum_{k = 1}^{c} B_{k} (x) θ_{k}$ , where $B_{k} (\cdot)$ is appropriately chosen basis function, for example, B-splines (De Boor, 1972). We can then rewrite model (2.1) in matrix form:

y = B θ + ε, ε \sim N (0, σ_{ε}^{2} I),

(2.2)

where $B$ , of size $n \times c$ , is the regression basis constructed row-wise from $B (x_{i})$ for $i = 1, . . ., n$ and $θ$ is the vector of regression coefficients with dimension $c \times 1$ . Rather than estimating the coefficients $θ$ in (2.2) by simple maximum likelihood (ML) methods, we penalize the coefficients through a quadratic penalty, that is, the fit is

\hat{μ} = B \hat{θ} = B (B^{'} B + λ P)^{- 1} B^{'} y,

where $P$ is any quadratic penalty of dimension $c \times c$ that forces the coefficients to vary smoothly, and consequently to obtain a smoothed curve, and $λ$ is the smoothing parameter—it can be estimated using a information criteria or a cross-validation criteria method. There are several alternatives for the choice of the regression basis $B$ and the penalty matrix $P$ . We choose a B-spline basis and set the penalty to $P = D_{q}^{'} D_{q}$ , with $D_{q}$ a difference matrix of order $q$ and dimension $(c - q) \times c$ (P-splines), although the results in Section 3 are valid for any basis and quadratic penalty.

Currie et al. (2004) proposed a method to fit and predict simultaneously in penalized regression models. We call their proposal the missing value approach, subsequently, and give a brief summary of their methodology.

In the framework of model (2.1), given a vector $y$ of $n$ observations of the response variable, suppose that we want to predict $n_{p}$ new values $y_{p}$ associated to the values of the independent variable $x_{p}$ , where $x_{p}$ may be within or, more interestingly, outside of the range of observed values $x_{i}$ . In the following, we focus on the case when $x_{p}$ is not in the convex hull of $x_{i}$ , $i = 1, . . ., n$ and it is to the right of $x$ (forward prediction), but $x_{p}$ could also be to the left of $x$ (backward prediction). We define the new vector of observations as

y_{+} = (y^{'}, y_{p}^{'})^{'},

(2.3)

which contains the observed response $y$ and the unknown values $y_{p}$ to be predicted. A new extended B-spline basis, $B_{+}$ , is built from a new set of knots that consists of the original knots covering $x_{i}$ , $i = 1, . . ., n$ , and extended to the range of the $n_{p}$ values of $x_{p_{j}}$ , $j = 1, . . ., n_{p}$ . This leads to the basis $B_{+} = [\begin{matrix} B & O \\ B_{1} & B_{2} \end{matrix}]$ of size $n_{+} \times c_{+}$ , where $B$ is the $n \times c$ basis used for fitting the trend component, $B_{1}$ and $B_{2}$ are auxiliary bases for prediction up to $n_{+} = n + n_{p}$ values, which are of sizes $n_{p} \times c$ and $n_{p} \times c_{p}$ , respectively, and $c_{+} = c + c_{p}$ . Figure 1 represents an extended splines basis. We show the original basis $B$ in black, the $B_{1}$ component in grey and the $B_{2}$ part with dashed line.

Figure 1:

Example of an extended basis to the right of the data (forward)

Associated to the new basis $B_{+}$ , a new vector of coefficients, $θ_{+} = (θ^{'}, θ_{p}^{'})^{'}$ , is defined, with length $c_{+} \times 1$ . A new quadratic penalty associated with the new set of coefficients needs to be introduced, lets say $P_{+}$ . Similar to $B_{+}$ , we can also decompose $P_{+}$ to

P_{+} = [\begin{matrix} P_{1} & P_{2} \\ P_{2}^{'} & P_{3} \end{matrix}] .

(2.4)

In the case that $P_{+}$ is built from the $q$ th order difference matrices:

P_{+} = {D^{'}}_{+} D_{+} = [\begin{matrix} D^{'} D + D_{1}^{'} D_{1} & D_{1}^{'} D_{2} \\ D_{2}^{'} D_{1} & D_{2}^{'} D_{2} \end{matrix}], with D_{+} = [\begin{matrix} D & O \\ D_{1} & D_{2} \end{matrix}],

(2.5)

that is, $P_{1} = D^{'} D + D_{1}^{'} D_{1}$ , $P_{2} = D_{1}^{'} D_{2}$ and $P_{3} = D_{2}^{'} D_{2}$ . Notice that $D_{+}$ has size $(c_{+} - q) \times c_{+}$ and that $D$ is the difference matrix used to build the penalty matrix for the observed data. Moreover, for a second-order penalty, $D_{+}$ is a banded matrix with three non-zero elements per row. Here, the subscripts do not indicate the order of the penalty, but the blocks of the extended differences matrix.

The model can be fitted and predicted simultaneously by minimizing the following penalized least squares criterion for $θ_{+}$ :

\begin{matrix} S & = & (y_{+} - B_{+} θ_{+})^{'} W (y_{+} - B_{+} θ_{+}) + λ θ_{+}^{'} P_{+} θ_{+}, \end{matrix}

(2.6)

where the unknown $y_{p}$ values of $y_{+}$ are arbitrary and $W$ is a diagonal matrix of dimension $n_{+} \times n_{+}$ with $0$ entries if the data is missing, that is for $y_{p}$ , and $1$ if the data is observed, that is for $y$ . Differentiating with respect to $θ_{+}$ leads to

{\hat{θ}}_{+} = (B_{+}^{'} W B_{+} + λ P_{+})^{- 1} B_{+}^{'} W y_{+},

(2.7)

and ${\hat{μ}}_{+} = H_{+} y_{+}$ with $H_{+} = B_{+} (B_{+}^{'} W B_{+} + λ P_{+})^{- 1} B_{+}^{'} W$ . Note that $λ > 0$ is required so that the matrix inversion in (2.7) exists.

Writing the fit and the forecast as a function of the extended penalty matrix (2.4) and applying Theorem 9.6.1 given in Harville (2000), we get the fit and the prediction

{\hat{μ}}_{+} = B_{+} [\begin{matrix} I \\ - P_{3}^{-} P_{2}^{'} \end{matrix}] (B^{'} B + λ P_{1} - λ P_{2} P_{3}^{-} P_{2}^{'})^{- 1} B^{'} y,

(2.8)

where the superscript ( $^{-}$ ) denotes the generalized inverse. Note that ${\hat{μ}}_{+} = ({\hat{μ}}^{'}, {\hat{μ}}_{p}^{'})^{'}$ is a fitted mean value, that is, ${\hat{μ}}_{+} = \hat{[y_{+}]}$ . Hence, in particular the new values $y_{p}$ are predicted by their fitted mean, where the fit is based on the observed values $y$ . Taking formula (2.8), we can derive the expectation of ${\hat{μ}}_{+}$ which for the new values result to $[{\hat{μ}}_{p}] = (B_{1} - B_{2} P_{3}^{-} P_{2}^{'}) (B^{'} B + λ P_{1} - λ P_{2} P_{3}^{-} P_{2})^{- 1} B^{'} B θ .$ This does not simplify to $[B_{1} | B_{2}] θ_{+}$ . Consequently, a bias is induced in the forecast. This bias cannot be determined in size since $θ_{p}$ is unknown. Hence, we cannot say anything about the quality of the estimate ${\hat{μ}}_{p}$ , nor can we derive any confidence interval statements. This is an unsatisfactory result. But we can change the perspective if we link penalized estimation to mixed models. This is done in Section 3.

When the aim of prediction is to forecast future values, the independent error assumption in (2.2) might not always be appropriate. If this is the case, the estimation of the smoothing parameter will be affected by the correlation that is not being accounted for. The methods presented here are easily extended to the case of non i.i.d. errors where smoothing and correlation parameters can be estimated simultaneously (Durbán and Currie, 2003). See also Krivobokoa and Kauermann (2007) for related results in mixed models.

2.1 Properties of the predictions in the case of P-splines with penalties based on differences

When penalties are based on differences between adjacent coefficients, for interpolating, there is no need for new coefficients and the B-spline coefficients form a polynomial sequence of degree $2 q - 1$ (Eilers and Marx, 2010); for instance when $q = 2$ , we get cubic interpolation. If we extrapolate instead of interpolate, the method above satisfies certain important properties. The subsequent results are based on a basis constructed from equally spaced knot. However, the results also extend to the non-equal spaced knots case, if we define the appropriately scaled penalty matrices. The central results are the following:

The fit remains the same regardless of the forecast horizon (the coefficients that yield the fit do not change).

The shape of the forecast is determined by the order of the penalty.

These properties are an immediate consequence of the following theorems.

Theorem 1. The coefficients from minimizing (2.6) with extended penalty matrix (2.5) satisfy the following properties:

The first $c$ coefficients of ${\hat{θ}}_{+}$ , are those obtained from the fit of $y$ : ${\hat{θ}}_{+_{1, . . ., c}} = \hat{θ} .$

The coefficients for the $n_{p}$ predicted values are ${\hat{θ}}_{p} = - D_{2}^{- 1} D_{1} \hat{θ}$ .

Proof. Substituting the blocks of $P_{+}$ by their specific values in (2.8), we have:

{\hat{μ}}_{+} = B_{+} [\begin{matrix} I \\ - D_{2}^{- 1} D_{1} \end{matrix}] (B^{'} B + λ P)^{- 1} B^{'} y,

that is the first $c$ coefficients of ${\hat{θ}}_{+}$ are $(B^{'} B + λ D^{'} D)^{- 1} B^{'} y,$ the same as the ones that give the fit, and the additional coefficients are

{\hat{θ}}_{p} = - D_{2}^{- 1} D_{1} (B^{'} B + λ D^{'} D)^{- 1} B^{'} y = - D_{2}^{- 1} D_{1} \hat{θ} .

(2.9)

If the knots are not equally spaced, the expression above would be modified, since the penalty would have to account for the difference between the knots (Eilers and Marx, 2010).

Corollary 1. Given penalties of order $q$ , the new coefficients are combinations of order $q - 1$ of the last $q$ fitted coefficients.

As the most popular penalties are of second or third order, the proof of the previous corollary for such cases and for penalties of order 1 is given in the complementary material.

In many situations, the fit is not greatly affected by the order of the penalty. However, there is an immediate connection between the penalty (or prior distribution) on the coefficients on the shape of the out-of-sample prediction shown in the above corollary. This is known in the framework of Bayesian P-splines, where the difference penalty corresponds to assuming a random walk prior on the coefficients (see Lang and Brezger, 2004); however, it is not common knowledge in a non-Bayesian context. We believe this is an important result that needs to be addressed when using this methodology.

3 Prediction with mixed effects smooth models

The connection between penalized smoothing and mixed models was established more than 30 years ago in (Green, 1987; see also Currie and Durbán, 2002; Wand, 2003). The key point of this equivalence is the fact that the smoothing parameter becomes a ratio of variances, $λ = \frac{σ_{ε}^{2}}{σ_{α}^{2}}$ , and both variance components can be estimated through restricted maximum likelihood (REML) procedure (see Patterson and Thompson, 1971; Kauermann, 2005). The idea is extensively discussed in Ruppert et al. (2003). We here exploit the link to mixed models to extend the results of the previous section and derive properties of the predicted value in (2.8).

In order to reparameterize a penalized smooth model, it is necessary to find a new basis that allows the representation of model (2.1) as a mixed model. Like before, we replace the smooth function by a basis representation which is now written as $y = X β + Z α + ε$ . Coefficient $α$ is penalized to achieve smoothness, but the penalty is rephrased as a normal prior on $α$ . This leads to a mixed model of the form

y = X β + Z α + ε, with α \sim N (0, σ_{α}^{2} G) and ε \sim N (0, σ_{ε}^{2} I),

(3.1)

where $X$ and $Z$ are the model matrices and $β$ and $α$ are the fixed and random effects coefficients, respectively. The random effects have the covariance matrix $σ_{α}^{2} G$ , which depends on the variance of the random effects $σ_{α}^{2}$ , as stated above $λ = \frac{σ_{ε}^{2}}{σ_{α}^{2}}$ .

There are different alternatives for the reparameterization of the original smooth model described before which was based on quadratic penalties as a mixed model in (3.1). The idea is to find a transformation $Ω$ such that:

B Ω = [X | Z] and Ω^{'} θ = {[β^{'} | α^{'}]}^{'} to have B θ = X β + Z α,

where $Ω$ is an orthogonal matrix. We split the matrix $Ω$ into two submatrices (for the fixed and the random components respectively), that is, $Ω = [Ω_{f} | Ω_{r}]$ , and achieve that $X = B Ω_{f}$ and $Z = B Ω_{r}$ . Since the fixed effects are unpenalized, the matrix $X$ , may be replaced by any submatrix such that the composed matrix $[X | Z]$ has full rank (this also implies that both $X$ and $Z$ have full column rank). For the submatrix $Ω_{r}$ , there are different alternatives. Following the approach of Currie and Durbán (2002), we use the eigenvalue decomposition of the penalty matrix. We here decompose $P = U \tilde{Σ} U^{'}$ , where $\tilde{Σ}$ is a diagonal matrix that contains the eigenvalues of $P$ , and $U$ is the corresponding matrix of eigenvectors. This allows to define $Ω_{r} = U_{r} Σ^{- 1 / 2},$ where $Σ$ contains the positive eigenvalues and $U_{r}$ contains the span of the decomposition. With this reparameterization, it is straightforward to obtain the relationship between the inverse of the covariance matrix $σ_{α}^{2} G$ of the random effects and the penalty $P$ : $\frac{1}{σ_{α}^{2}} G^{- 1} = \frac{1}{σ_{α}^{2}} Ω_{r}^{'} P Ω_{r} .$

Prediction in the context of mixed models has always been done as a two-stage procedure: first fit and then predict. We will first show how to apply the existing results to the context of smooth mixed models, and then, we will propose an alternative one-stage approach. As we will see, the variance–covariance matrix for the random effects in the two-stage approach is a direct extension of the variance–covariance matrix of the random effects in the fit. This implies that the extended mixed model transformation used has to provide a direct extension of the variance-covariance matrix of the random effects in the fit. One of the advantages of the one-stage approach is that we can use any transformation.

3.1 Two-stage approaches

3.1.1 Standard methodology for prediction

Once we have established the connection between mixed models and P-splines, we can use the results given in Gilmour et al. (2004) to predict new observations. In this case, the prediction is a linear function of the best linear unbiased predictor (BLUP) of random effects and the best linear unbiased estimator (BLUE) of the fixed effects in the model. The results are based on the following augmented mixed model,

y_{+} = X_{+} β + Z_{+} α_{+} + ε_{+},

(3.2)

that is,

[\begin{matrix} y \\ y_{p} \end{matrix}] = [\begin{matrix} X \\ X_{p} \end{matrix}] β + [\begin{matrix} Z & O \\ Z_{1} & Z_{2} \end{matrix}] [\begin{matrix} α \\ α_{p} \end{matrix}] + [\begin{matrix} ε \\ ε_{p} \end{matrix}],

with $(ε^{'}, ε_{p}^{'})^{'} \sim N (0, σ_{ε}^{2} I_{+})$ , $β$ as the fixed effects (the same as the ones that give the fit, the fixed part is linear and therefore no new parameters for the linear part are needed), and $α_{+} = (α^{'}, α_{p}^{'})^{'}$ as the augmented random effects with covariance matrix

Var [α_{+}] = σ_{α}^{2} G_{+} = σ_{α}^{2} [\begin{matrix} G & G_{op} \\ G_{po} & G_{pp} \end{matrix}],

where $G$ is the covariance matrix of the random effects in the model for the observed data, $G_{op}$ is the covariance matrix between the random effects for the observed data and for the unobserved data, and $G_{pp}$ is the covariance matrix of the random effects for the unobserved data. The variance components, $σ_{ε}^{2}$ and $σ_{α}^{2}$ , are the ones estimated in the fit through REML procedure (Patterson and Thompson, 1971).

Now we need to formulate the extended P-spline model into the extended mixed model (3.2). For that, we define the extended transformation matrix $Ω_{r_{+}} = [\begin{matrix} Ω_{r} & O \\ O & Ω_{p_{r}} \end{matrix}],$ where $Ω_{r}$ is the transformation matrix used for the observed data, and $Ω_{p_{r}}$ is the transformation matrix used for the predicted values. Hence, we have $Z_{1} = B_{1} Ω_{r}$ and $Z_{2} = B_{2} Ω_{p_{r}}$ . There are many ways in which $Ω_{r_{+}}$ may be chosen. In the context of penalties based on differences, for simplicity, we chose $Ω_{r} = U_{r} Σ^{- 1 / 2}$ , based on the eigenvalue decomposition of $D^{'} D = U \tilde{Σ} U^{'}$ , and ${Ω_{p}}_{r} = D_{2}^{- 1}$ , with $D$ and $D_{2}$ blocks of the extended difference matrix $D_{+}$ , see Equation (2.5). We choose this extended transformation matrix to obtain an extended variance–covariance matrix of random effects that is a direct extension of $G$ , the variance–covariance matrix of the random effects in the fit.

Then, the new predicted values are

{\hat{μ}}_{p} = X_{p} \hat{β} + Z_{(p)} \hat{α},

(3.3)

with $Z_{(p)} = Z_{1} + Z_{2} G_{po} G^{- 1}$ and $\hat{β}$ and $\hat{α}$ being the BLUE and BLUP, respectively, estimated from the observed data, that is,

[\begin{matrix} \hat{β} \\ \hat{α} \end{matrix}] = Q^{- 1} [\begin{matrix} \frac{1}{σ_{ε}^{2}} X^{'} \\ \frac{1}{σ_{ε}^{2}} Z^{'} \end{matrix}] y,

(3.4)

where $Q = [\begin{matrix} \frac{1}{σ_{ε}^{2}} X^{'} X & \frac{1}{σ_{ε}^{2}} X^{'} Z \\ \frac{1}{σ_{ε}^{2}} Z^{'} X & \frac{1}{σ_{ε}^{2}} Z^{'} Z + \frac{1}{σ_{α}^{2}} G^{- 1} \end{matrix}] .$

It follows that the predicted random effects vector for $α_{p}$ is ${\hat{α}}_{p} = G_{po} G^{- 1} \hat{α}$ .

Therefore, the $95 %$ confidence interval is

{\hat{μ}}_{p} \pm 1.96 \sqrt{diag ([\begin{matrix} X_{p} | Z_{(p)} \end{matrix}] Q^{- 1} [\begin{matrix} X_{p}^{'} \\ Z_{(p)}^{'} \end{matrix}])},

(3.5)

with $Q$ defined above and the variance components estimated through REML procedure; (see Patterson and Thompsonm, 1971). We denote this method as mixed model approach.

3.1.2 Prediction based on the conditional distribution of

y_{p} | y

The previous method can be seen from the point of view of conditional distributions, which allows us to compute prediction intervals. We rewrite (3.2) as

y_{+} = {[y^{'} | y_{p}^{'}]}^{'} \sim N (X_{+} β, V_{+}),

(3.6)

with $Var [y_{+}] = V_{+} = σ_{α}^{2} Z_{+} G_{+} Z_{+}^{'} + σ_{ε}^{2} I_{+}$ . The mixed model formulation connects the new values $y_{p}$ to the observed vector $y$ through a joint normal distribution. It, therefore, appears natural to predict $y_{p}$ given $y$ based on the conditional model resulting from (3.6), that is,

y_{p} | y \sim N (X_{p} β + V_{po} V_{oo}^{- 1} (y - X β), V_{pp} - V_{po} V_{oo}^{- 1} V_{op}),

where $V_{oo}$ , $V_{op}$ and $V_{pp}$ are the submatrices of matrix $V_{+}$ matching to $y$ and $y_{p}$ . The mean value results through

[y_{p} | y] = X_{p} \hat{β} + V_{po} V_{oo}^{- 1} (y - X \hat{β}) = X_{p} \hat{β} + Z_{1} \hat{α} + Z_{2} G_{po} G^{- 1} \hat{α},

which equals (3.3). Note that the first term in the above equation is the result of plugging $X_{p}$ into the regression equation, and it represents the adjustment to this prediction based on the covariance between $y$ and $y_{p} .$ The conditional variance of $y_{p} | y$ gives the prediction error. This follows since

_{y, y_{p}} [({\hat{μ}}_{p} - y_{p})^{2}] =_{y} [_{y_{p}} [({\hat{μ}}_{p} - y_{p})^{2} | y]] =_{y} [Var [y_{p} | y]] = Var [y_{p} | y],

where the latter equality holds, since in the Normal model the conditional variance does not depend on the value of the variable we condition on. With these results, we can construct $95 %$ ‘prediction interval’:

{\hat{μ}}_{p} \pm 1.96 \sqrt{Var [y_{p} | y]},

where observed values are considered as fixed. Note that, as we have mentioned before, this approach allows us to work out the posterior predictive distribution $y_{p} | y$ as a Gaussian distribution and therefore compute the prediction intervals. While we cannot compute prediction intervals with the standard methodology described in 3.1.1 unless we link it with a Gaussian process.

3.2 One-stage approach

The previous method was a two-stage procedure. First, the model is fitted to the available data and, second, based on the fitted model we predict the new observations. As mentioned previously, this approach imposes constraints on the reparameterization used to obtain the smooth mixed model, since the variance–covariance matrix of the extended model (the model that includes out-of-sample prediction) needs to be an extension of the variance–covariance matrix of the original model. Now we propose an alternative approach which can be used with any reparameterization and yields the same results as the two-step approach for appropriate transformations. This approach relates the above results to the method presented in Section 2. We will call it ‘extended mixed model approach’, since we include $y_{p}$ in the model, but with infinite variance (zero weight). In this case, we consider the model

y_{+} = X_{+} β + Z_{+} α_{+} + ε_{+}, α_{+} \sim N (0, σ_{α}^{2} G_{+}), ε_{+} \sim N (0, σ_{ε}^{2} R_{+}),

(3.7)

where $R_{+}$ is a diagonal weight matrix of dimension $n_{+} \times n_{+}$ , with $1$ entry if the data is observed, that is, for $y$ , and infinity if the data is considered to be forecasted, that is, for $y_{p}$ . The quantity infinity expresses that we do not have any information about the data $y_{p}$ . Its estimation is done using the extended mixed model equations of Henderson (1975):

\begin{matrix} \hat{\hat{β}} & = & (X_{+}^{'} {\tilde{V}}_{+}^{-} X_{+})^{- 1} X_{+}^{'} {\tilde{V}}_{+}^{-} y_{+}, \\ {\hat{\hat{α}}}_{+} & = & σ_{α}^{2} G_{+} Z_{+}^{'} {\tilde{V}}_{+}^{-} (y_{+} - X_{+} \hat{\hat{β}}), \end{matrix}

(3.8)

where $y_{+} = (y^{'}, y_{p}^{'})^{'}$ as in (2.3), ${\tilde{V}}_{+} = σ_{α}^{2} Z_{+} G_{+} Z_{+}^{'} + σ_{ε}^{2} R_{+}$ and by Theorem 18.2.8 given in Harville (2000) ${\tilde{V}}_{+}^{-} = \frac{1}{σ_{ε}^{2}} R_{+}^{- 1} - \frac{1}{σ_{α}^{4}} R_{+}^{- 1} Z_{+} (\frac{1}{σ_{α}^{2}} G_{+}^{- 1} + Z_{+}^{'} \frac{1}{σ_{ε}^{2}} R_{+}^{- 1} Z_{+})^{- 1} Z_{+}^{'} R_{+}^{- 1}$ . Finally, $Z_{+} = B_{+} Ω_{+_{r}}$ , with $Ω_{+_{r}}$ as any orthogonal transformation such that

B_{+} [Ω_{+_{f}} | Ω_{+_{r}}] = [X_{+} | Z_{+}] and {[Ω_{+_{f}} | Ω_{+_{r}}]}^{'} θ_{+} = {[β^{'} | α_{+}^{'}]}^{'} .

For the particular case of penalties based on differences, we choose $Ω_{+_{r}}$ based on the eigenvalue decomposition of $D_{+}^{'} D_{+}$ . This yields $B_{+} θ_{+} = X_{+} β + Z_{+} α_{+}$ .

We use the double hat symbol ( $\hat{\hat{}}$ ) here to remark that in this case the estimation is based on the extended version assuming that the unknown values $y_{p}$ have infinity variance. Note that to compute the fixed and random effects, we do not need $R_{+}$ , we need its inverse, $R_{+}^{- 1}$ , with $1$ entry if the data is observed and $0$ entry if the data is considered to be predicted.

Like the above, we have ${\hat{\hat{α}}}_{+} = \hat{\hat{[α_{+} | y]}}$ with the double hat on the expectation denoting that $β$ has been estimated. It is easy to see that taking expectation over $y$ (or $y_{+}$ ), we get $[{\hat{\hat{α}}}_{+}] = 0$ so that the estimate is unbiased in the mixed model.

As we mentioned earlier, the above extension of the missing value approach to the mixed model framework fits and predicts simultaneously, while the approach of Gilmour et al. (2004) is a two-stage method. In order to know the relationship between the two methods, we need to know the relationship between the covariance matrix of the random effects that gives the fit and the extended covariance matrix. This is shown in the following theorem.

Theorem Given model (2.1) penalizing differences between adjacent coefficients, the fit and the prediction of new observations given by extended mixed model approach and mixed model approach are the same if the transformation matrix used to obtain the model components in (3.8) is the direct extension of the original transformation,

Ω_{+_{r}} = [\begin{matrix} Ω_{r} & O \\ O & {Ω_{p}}_{r} \end{matrix}] .

(3.9)

Under the previous hypothesis, the variance components ( $σ_{α}^{2}$ , $σ_{ε}^{2}$ ) that maximize the REML criteria, $l$ , and the REML criteria corresponding to the extended weighted model, $l_{+}$ , are equal,

l (σ_{ε}^{2}, σ_{α}^{2}) = - \frac{1}{2} \log | V | - \frac{1}{2} \log | X^{'} V^{- 1} X | - \frac{1}{2} {(y - X β)}^{'} V^{- 1} (y - X β),

(3.10)

l_{+} (σ_{ε}^{2}, σ_{α}^{2}) = - \frac{1}{2} \log | {\tilde{V}}_{+} | - \frac{1}{2} \log | X_{+}^{'} {\tilde{V}}_{+}^{-} X_{+} | - \frac{1}{2} {(y_{+} - X_{+} β)}^{'} {\tilde{V}}_{+}^{-} (y_{+} - X_{+} β) .

(3.11)

Equivalent results would be obtained using ML.

The proof can be found in the complementary material. Since the fitted and predicted values in both approaches do not depend on the transformations used, the previous theorem states a stronger result: The approaches always give the same solution, regardless of the used transformations. We have stated the theorem for a particular transformation to establish a relationship between both the methods. The last statement of the previous theorem means that the variance parameters used to predict are the same as the ones used for estimating the original fit. In other words, prediction in and out of sample cannot only be done simultaneously but also the optimal smoothing parameter is the same.

4 Applications

In this section, we apply the proposed methods to two real datasets. The first one, allows us to show an example of predicting within the framework of additive models and in the case when out-of-sample prediction is needed to the left and right of the interval where the covariate is observed. The second dataset illustrates a classic example where forecasting is needed, when data are collected over time. Although, as we have shown, all the methods give us the same result, in order to obtain the prediction intervals we use the two-stage approach.

4.1 Prediction of aboveground biomass

All the results presented in the previous sections are obtained in the case of smooth models with a single covariate. However, it is immediate to extend these results to the case of semiparametric or additive models. In this section, we apply the proposed methodology to such dataset. The dataset corresponds to an agricultural trial carried out in Spain (Rivas-Martínez et al., 2002) with the aim of evaluating the economic viability of Populus trees prior to harvesting. In this context, it is essential to estimate aboveground biomass and obtain accurate predictions using only a minimum set of easily obtainable information, that is diameter and height. Sánchez-González et al. (2016) proposed the use of smooth additive mixed models for predicting aboveground biomass. Here, we analyse data for a single clone of the nine included in the trials. The aim is to estimate the productions (measured as aboveground dry biomass) as a function of diameter and height, and give out-of-sample predictions. The observed data consists of $315$ observations for diameter values measured at $1.30$ m breast height. This data is illustrated in Figure 2.

Figure 2:

Plot of weight versus diameter (left panel) and plot of weight versus height (right panel)

From the plot, it is immediate to see that the weight is a non-linear function of diameter and height. Therefore, we fit the following model:

[y_{i} | x_{i}, y_{i}] = f (x_{i}) + f (z_{i}),

(4.1)

where $x$ is the diameter and $z$ is the height, that is, $f (x)$ and $f (z)$ are the functions that represent the main effects of the diameter and of the height, respectively. The regression matrix is then defined by blocks as $B = [B_{x} | B_{z}]$ with marginal B-spline bases of degree three of the covariates diameter and height, $B_{x}$ and $B_{z}$ , respectively. The penalty matrix associated with model (4.1) has a block diagonal form: $P = blockdiag (λ_{x} P_{x}, λ_{z} P_{z}),$ where $P_{x}$ and $P_{z}$ are the marginal second-order difference penalties for diameter and height, respectively.

We predict weight for $6$ new out-of-sample values for diameter and $15$ new values for height where $5$ are to the left and $10$ to the right of the range of the observed height values. Applying the previous methodology, we extend the basis and the penalty: $B_{+} = [B_{x_{+}} | B_{z_{+}}], P_{+} = blockdiag (λ_{x} P_{x_{+}}, λ_{z} P_{z_{+}}),$ with $B_{x_{+}} = [\begin{matrix} B_{x} & O \\ B_{x_{(1)}} & B_{x_{(2)}} \end{matrix}]$ and $B_{z_{+}} = [\begin{matrix} B_{z_{(0)}} & O & O \\ O & B_{z} & O \\ O & B_{z_{(1)}} & B_{z_{(2)}} \end{matrix}]$ .

Once we have extended the basis and the penalty, it is straightforward to obtain the fit and the prediction applying Equation (2.7). But in order to obtain confidence and prediction intervals and to avoid identifiability problems, since the column of $1$ s is contained in the space spanned by the columns of $B_{x_{+}}$ and $B_{z_{+}}$ , we reparameterize the model using the representation of a penalized spline model as a mixed model. Figure 3 shows the smooth fitted and predicted trend for diameter and height, the $95 %$ confidence interval (grey line) and the $95 %$ prediction interval (dashed lines). Note that the prediction is done backward and forward; the proposed methodology allows us to obtain the prediction for any range of the independent variable. This could not be done by using methods developed in the series temporal framework.

Figure 3:

Fit, forecast, $95 %$ confidence interval (grey lines) and $95 %$ prediction interval (dashed lines) of the additive smooth term for diameter (left panel) and for height (right panel), result of applying the methodology of a dataset on the stem biomass

The estimation of the covariance parameters was carried out by the REML through the SOP algorithm, proposed in Rodríguez-Álvarez et al. (2018).

4.2 Forecasting

{SO}_{2}

concentration levels

Now, we analyse data where out-of-sample prediction is needed and missing observations are present in the data. In this situation, the one-stage approach offers an elegant solution to solve simultaneously both problems. We consider measurements on sulphur dioxide ( ${SO}_{2}$ ) concentration levels (in $μ g / m^{3}$ ) over station AT02 from January 1990 to December 2001. Figure 4 shows the dataset in which we can see that there are some missing observations between October 1995 and March 1999.

Figure 4:

Time series plot of $\log ({SO}_{2})$ data for station AT02

The data were collected through the ‘European Monitoring and Evaluation Programme’ (EMEP) under the Co-operative Programme for Monitoring and Evaluation of the Long-range Transmission of Air Pollutants in Europe (see further information available at htpp:\\www.emep.int). Note that there is a clear evidence of temporal trends and seasonal effects and that there are some missing observations (mostly due to equipment failure, replacement or calibration); there is a big gap between October 1995 and March 1999.

First of all, let us introduce the smooth modulation model based on P-splines suggested by Eilers et al. (2008):

[y_{i}| x_{i}] = f (x_{i}) + \sum_{j = 1}^{J} \{g_{j} (x_{i}) \cos (j ω x_{i}) + h_{j} (x_{i}) \sin (j ω x_{i})\},

(4.2)

where $f (\cdot)$ accounts for the smooth trend, $g (\cdot)$ and $h (\cdot)$ are smooth functions that describe the local amplitudes of cosine and sine waves, and $ω = 2 π / p$ , with $p$ being the period (e.g., $p = 12$ for monthly data). The number of harmonics functions, $J$ , required for the seasonal component is usually taken as $1$ or $2$ to reduce the number of parameters to be estimated. If $J = 1$ , the model (4.2) can be written in matrix form as:

y = \overset{}{B} \overset{}{θ} + ε, ε \sim N (0, R),

(4.3)

where $R$ is the covariance matrix of the error—we work with uncorrelated i.i.d. errors—that is, $R = σ_{ε}^{2} I$ ; and $\overset{}{B}$ is the regression matrix, $\overset{}{B} = [B | C B | S B],$ where $B$ is a B-spline basis, $C = diag \{\cos (ω x)\}$ and $S = diag \{\sin (ω x)\}$ .

If we want to obtain the coefficients $\overset{}{θ} = (θ, θ_{c}, θ_{s})$ in (4.3) with the P-splines method, we have to minimize the following function of $\overset{}{θ}$ :

S = (y - \overset{}{B} \overset{}{θ})^{'} (y - \overset{}{B} \overset{}{θ}) + {\overset{}{θ}}^{'} \overset{}{P} \overset{}{θ},

(4.4)

where $\overset{}{P} = \overset{}{λ} \otimes {\overset{}{D}}^{'} \overset{}{D}$ , with $\overset{}{λ} = diag (λ, λ_{c}, λ_{s})$ and $\overset{}{D} = blockdiag (D_{q}^{'} D_{q}^{'}, I_{2} \otimes D_{q_{cs}}^{'} D_{q_{cs}})$ , $q$ and $q_{cs}$ are usually $2$ and $1$ , respectively.

To obtain the fit and the forecast simultaneously, we just have to extend the B-spline basis for the trend and the modulation components ${\overset{}{B}}_{+} = [B_{+} | C_{+} B_{+} | S_{+} B_{+}],$ where $C_{+} = diag (\cos (ω x_{+}))$ and $S_{+} = diag (\sin (ω x_{+}))$ , for the additive modulation blocks $C_{+} B_{+}$ and $S_{+} B_{+}$ , and consider the following penalty matrix:

{\overset{}{P}}_{+} = blockdiag (λ {D_{+}}_{q}^{'} {D_{+}}_{q}, λ_{c} {D_{+}}_{q_{c}}^{'} {D_{+}}_{q_{c}}, λ_{s} {D_{+}}_{q_{s}}^{'} {D_{+}}_{q_{s}}) .

Once ${\overset{}{B}}_{+}$ and ${\overset{}{P}}_{+}$ are computed, ${\overset{}{θ}}_{+}$ can be easily computed through the formula in (2.7) by using ${\overset{}{B}}_{+}$ instead of $B_{+}$ and ${\overset{}{P}}_{+}$ instead of $λ P_{+}$ .

Keeping in mind that in a penalized spline modulation model, we have the trend and the seasonality components, it is straightforward to represent it as a mixed model.

As we have seen the order of the penalty, $q$ denotes the $q - 1$ grade polynomial of the fixed part when the smoothing parameters are very large (the null model), and hence, the random part can be considered as smooth deviates from the null model. In model (4.3), it is considered a first-order penalty for the modulation terms, then the null terms for the modulation are $\cos (ω x)$ and $\sin (ω x)$ . For instance, for $J = 1$ in (4.2), the fixed effect matrix for the smooth modulation model (4.3) is a design matrix of a harmonic regression model: $\overset{}{X} = [1_{n} | x | x^{2} | . . . | x^{q - 1} | \cos (ω x) | \sin (ω x)] .$

The matrix of the random component, $\overset{}{Z}$ , is a block matrix: $\overset{}{Z} = [Z | C Z_{cs} | S Z_{cs}],$ where $Z = B Ω_{r}$ with $Ω_{r} = U_{r} {\tilde{Σ}}^{- 1 / 2}$ ( $U_{r}$ and $\tilde{Σ}$ are obtained from the eigenvalue decomposition of $D_{q}^{'} D_{q}$ , with $D_{q}$ a penalty matrix of order usually $q = 2$ ) and $Z_{cs} = B \overset{}{Ω}$ , with $\overset{}{Ω} = {\overset{}{U}}_{r} {\overset{}{\tilde{Σ}}}^{- 1 / 2}$ ( ${\overset{}{U}}_{r}$ and $\overset{}{\tilde{Σ}}$ are obtained from the eigenvalue decomposition of $D_{q_{cs}}^{'} D_{q_{cs}}$ , with $D_{q_{cs}}$ a penalty matrix of order usually $q_{cs} = 1$ ). The smoothing parameters $λ$ , $λ_{c}$ , $λ_{s}$ in (4.4) become the ratio between the error variance term and the random effect variances $σ^{2}$ , $σ_{c}^{2}$ and $σ_{s}^{2}$ for the trend and modulation terms respectively. The covariance matrix is a block diagonal matrix, $\overset{}{G} = blockdiag (G, G_{c}, G_{s})$ , where $G = σ^{2} I_{c - q}$ , $G_{c} = σ_{c}^{2} I_{c - 1}$ and $G_{s} = σ_{s}^{2} I_{c - 1}$ .

The one-stage approach can deal with missing observations within and out of sample just by setting to zero the corresponding diagonal entries of $R_{+}^{- 1}$ . Figure 5 shows the forecasted trend and final predictions (including the seasonal projections) for AT02 station with second and third penalty orders for the trend component. Note that, as it is stated in Corollary 1, for second-order penalty ( $q = 2$ ) the trend forecast is linear, and for third-order penalty ( $q = 3$ ) the trend forecast is quadratic.

Figure 5:

Forecast for AT02 station. Top and bottom left figures show the data (points), the fitted and forecasted trend (black line), 95% confidence interval (grey lines) and 95% prediction interval (dashed lines) for second and third order penalties, respectively. Top and bottom right figures show the data (points), the fit and the forecast in the modulation model (black line), 95% confidence interval (grey lines) and 95% prediction interval (dashed lines) for second and third order penalties, respectively

5 Conclusions

Smoothing techniques have become a very popular tool for describing patterns in noisy data. However, prediction is still an open area of research. In this article we have proposed a general framework for prediction of new observations in penalized regression. The methodology proposed can be accommodated to the different frameworks in which smoothing is carried out, regardless of the basis and quadratic penalty used:

1. Extend the basis used for regression and the penalty to control the smoothness in the framework of penalized regression based on quadratic penalties. 2. Extend the fixed and random components in the context of mixed models and carry out the prediction as a two-stage procedure or using a new approach that allows us to fit and predict simultaneously.

The approach given in Section 3.1.2 is a particular case of prediction in the context of Gaussian process regression, under the assumption that the regression function is a realization of a stochastic process. Then, prediction of new values are based on Gaussian process prior and a P-spline covariance matrix. Hence, we can use the flexibility of the Gaussian process and the choice of suitable covariance matrix to model any non-linear model nonparametrically. In this context, prediction is straightforward due to the properties of Gaussian processes.

In the context of penalties based on differences between adjacent coefficients, we have proved that all the approaches compared in this work are equivalent, but the one-stage approach can be used regardless of the mixed model reparameterization used. Furthermore, when predicting out-of-sample observations, the horizon of prediction has no effect on the fit or on the selection of the smoothing parameter.

These results are based on the basis constructed from equally spaced knot. However, the results apply also to the non-equal spaced knots case, defining the appropriate scaled penalty matrices.

To illustrate the methodology and the proved results, we have analysed two real datasets and showed the performance of the extended mixed model approach with B-spline basis and penalties based on differences. These examples illustrate how to obtain out-of-sample prediction in the case of additive models. In the first example, we show how the method can do backward and forward predictions, and in the second example, we show how the one-stage approach can deal with out-of-sample and missing value predictions simultaneously.

The results obtained here have opened new areas of research; in particular, we have already obtained similar results in the case of multidimensional smoothing and preliminary work is being done to extend this work to the case of non-Gaussian data.

Supplementary material

Supplementary material for this article are available from http://www.statmod.org/smij/archive.html.

Footnotes

Declaration of conflicting interests

Funding

The first and second authors acknowledge financial support from the Spanish Ministry of Economy and Competitiveness MTM2014-52184-P. The fourth author acknowledges the financial support by the Basque Government through the BERC 2018-2021 programme and by Spanish Ministry of Economy and Competitiveness MINECO through BCAM Severo Ochoa excellence accreditation SEV-2013-0323 and through project MTM2017-82379-R funded by (AEI/FEDER, UE) and acronym AFTERAM.

References

Sinn

Goude

Pompey

(2012) Adaptive learning of smoothing functions: Application to electricity load forecasting. In Advances in Neural Information Process- ing Systems 25, edited by Pereira

Burges

Bottou

Weinberger

2519–27.

Besag

Green

Higdon

Mengersen

(1995) Bayesian computation and stochastic system. Statistical Science , 10, 3–66.

Cressie NAC (1993) Statistics for Spatial Data . New York, NY: Wiley.

Currie

Durban

(2002) Flexible smoo- thing with P-splines: A unified approach. Statistical Modelling , 2, 333–49.

Currie

Durban

Eilers

(2004). Smoothing and forecasting mortality rates. Statistical Modelling . doi: 10.1191/1471082X04st080oa

De Boor

(1972) On calculating with b-splines. Journal of Approximation Theory , 6, 50–62.

Durban

Currie

(2003) A note on p-spline additive models with correlated errors. Computational Statistics , 18, 251–62.

Eilers

Gampe

Marx

Rau

(2008) Modulation models for seasonal time series and incidence tables. Statistics in Medicine , 27, 3430–41.

Eilers

PHC

Marx

(1996) Flexible smoothing with b-splines and penalties. Statistical Science , 11, 89–121.

10.

Eilers

PHC

Marx

(2010) Splines, knots, and penalties. Computational Statistics , 2, 637–53.

11.

Gilmour

Cullis

Welham

Gogel

Thompson

(2004) An efficient computing strategy for prediction in mixed linear models. Computational Statistics and Data Analysis , 44, 571–86.

12.

Green

(1987) Penalized likelihod for general semi-parametric regression models. Inter- national Statistical Review , 55, 245–59.

13.

Harville

(2000) Matrix Algebra from a Statistician's Perspective . New York: Springer.

14.

Henderson

(1975) Best linear unbiased estimation and prediction under a selection model. Biometrics , 31, 423–47.

15.

Hyndman

Koehler

Ord

Snyder

(2008) Forecasting with Exponential Smoothing . Berlin Heidelberg: Springer.

16.

Jones

Schonlau

William

(1998) Efficient global optimization of expensive black-box functions. Journal of Global Optimization , 13, 455–92.

17.

Kauermann

(2005) A note on smoothing parameter selection for penalised spline smoothing. Journal of Statistical Planing and Inference , 127, 53–69.

18.

Krivobokoa

Kauermann

(2007) A note on penalized spline smoothing with correlated errors. Journal of the American Statistical Association , 102, 1328–37.

19.

Lang

Brezger

(2004) Bayesian p-splines. Journal of Computational and Graphical Statistics , 13, 183–212.

20.

Patterson

Thompson

(1971) Recovery of inter-block information when block sizes are unequal. Biometrika , 58, 545–54.

21.

Rivas-Martínez

Díaz

Fernández-González

Izco

Loidi

Lousã

Penas

(2002). Vascular plant communities of Spain and Portugal. Addenda to the Syntaxonomical checklist of 2001. Itinera Geobotanica , 15, 5–922.

22.

Rodrguez-Alvarez

Lee

D-J

Kneib

Durban

Eilers

(2018) On the estimation of variance parameters in non-standard generalised linear mixed models: Application to penalised smoothing. Statistics and Computing , 29, 1–18.

23.

Rue

Held

(2005) Gaussian Markov Random Fields . Boca Raton: Chapman & Hall.

24.

Ruppert

Wand

Carroll

(2003) Semiparametric Regression (Cambridge Series in Statistical and Probabilistic Mathematics) . New York: Cambridge University Press.

25.

Sacks

Welch

Mitchell

Wynn

(1989) Design and analysis of computer experiments. Statistical Science , 4, 409–35.

26.

Sánchez-González

Durbán

Lee

Cañellas

Sixto

(2016) Smooth additive mixed models for predicting aboveground biomass. Journal of Agricultural, Biological and Environmental Statistics , 22, 23–41.

27.

Wand

(2003) Smoothing and mixed models. Computational Statistics , 18, 223–49.

A general framework for prediction in penalized regression

Abstract

Keywords

1 Introduction

2 Prediction with smooth models and quadratic penalties

Example of an extended basis to the right of the data (forward)

3.1.1 Standard methodology for prediction

Equivalent results would be obtained using ML.

4.1 Prediction of aboveground biomass

Figure 2:

Plot of weight versus diameter (left panel) and plot of weight versus height (right panel)

Fit, forecast, 95 % confidence interval (grey lines) and 95 % prediction interval (dashed lines) of the additive smooth term for diameter (left panel) and for height (right panel), result of applying the methodology of a dataset on the stem biomass

Figure 4:

Time series plot of log ( SO 2 ) data for station AT02

Supplementary material

Footnotes

Declaration of conflicting interests

Funding

References

Fit, forecast, $95 %$ confidence interval (grey lines) and $95 %$ prediction interval (dashed lines) of the additive smooth term for diameter (left panel) and for height (right panel), result of applying the methodology of a dataset on the stem biomass

Time series plot of $\log ({SO}_{2})$ data for station AT02