Penalized complexity priors for degrees of freedom in Bayesian P-splines

1 Introduction

Penalized spline (P-spline) regression is a well-established and numerically stable approach for smoothing (Eilers and Marx, 1996; Ruppert et al., 2003). Typically, P-spline components are implemented in Bayesian additive regression models (Fahrmeir et al., 2013) to fit non-linear covariate effects or higher dimensional effects such as spatial and spatio-temporal smooth trends. The P-spline approach proposed by Eilers and Marx (1996) uses equally spaced B-splines and constructs a smooth function as the sum of these B-splines scaled by spline coefficients. A regularizing penalty is assumed on these coefficients to control the degree of smoothness of the fitted function. A common approach is to penalize the sum of second-order-squared differences between adjacent spline coefficients, but specific penalties can be designed to drive the fit towards desired features (Eilers and Marx, 2010). This is a very useful strategy in presence of prior information about the shape, or degree of smoothness, of the function to be estimated.

The Bayesian approach to P-spline by Lang and Brezger (2004) assumes an intrinsic Gaussian Markov random field (IGMRF) prior on the spline coefficients. An IGMRF is a multivariate normal distribution with rank deficient precision matrix Q ( τ ) , depending on a precision hyper-parameter τ . Similarly to a regularizing penalty, the IGMRF forces the spline coefficients to be shrunk towards an infinite smooth model, which we will denote as the ‘base model’. The degree of smoothness of the base model depends on the order of the IGMRF; for instance, an IGMRF of order 2 forces shrinkage towards a linear trend, that is, a polynomial of degree one (Rue and Held, 2005).

The amount of shrinkage towards the corresponding base model depends on the IGMRF precision τ . The prior π ( τ ) can have a substantial impact on the posterior distribution of the spline coefficients and hence, to some extent, on the shape of the fit. A common strategy in Bayesian P-splines is to adopt the conjugate gamma family, that is, gamma ( a , b ) , with shape a and rate b (Fahrmeir and Kneib, 2009; Lang and Brezger, 2004). Lang and Brezger (2004) suggest to choose a = 1 and small b , for example, b = 5 · 10 − 4 , leading to a diffuse prior for τ − 1 . Jullion and Lambert (2007) note that the choice of b clearly affects the smoothness of the fitted curve, when sample size or signal-to-noise ratio is small, and propose a mixture of gamma distributions with different b values. Another popular choice is the gamma ( ε , ε ) , with small ε (e.g., ε = 0.001 , which is the default option in the software BayesX (Belitz et al., 2000)) as an attempt of vagueness on τ − 1 . The suitability of the gamma family as a non-informative prior for the scale parameters in hierarchical models has been debated in the literature (Gelman, 2006); overfitting due to gamma priors has been demonstrated in Frühwirth-Schnatter and Wagner (2010, 2011) and Simpson et al. (2014). In particular, in Bayesian P-splines, the main difficulty with using a gamma prior on τ is that τ scales differently according to the amount of noise present in the data and the number (and location) of knots selected by the user.

The present work proposes a new prior for τ which is informative about model complexity and implicitly accounts for different choices about number (and location) of knots. A suitable measure of complexity of the P-spline model is the number of ‘effective degrees of freedom’ denoted as d in the following text, calculated as the trace of the hat matrix (Hastie and Tibshirani, 1990)). The value d relates to the degree of a polynomial equivalent to the smooth function to be estimated. An expert user who has a prior guess about the shape of this function may find easy to elicit d . As an example, for a monotonic cubic trend, one may elicit d in a range between 3 and 5 and assign very low prior probability to d > 5 . The key point is that, in presence of this prior information, elicitation of a range for d is intuitive and immediate, whereas elicitation of a distribution for τ , directly, is very difficult.

The challenge is to design a prior distribution on a model property (i.e., d ) rather than on a parameter of the model (i.e., τ ). To achieve this, we follow the penalized complexity (PC) prior approach proposed by Simpson et al. (2014). Within this framework, we derive the PC prior for d and calibrate it by two intuitive parameters: U , an upper bound for d , and α , the prior probability assigned to d > U . In the aforementioned example, the user would only need to set U = 5 and α = 0.01 , or some other small value. As a further challenging point, note that d depends on the noise variance characterizing the dataset. Thus, implementing the proposed PC prior for degrees of freedom in real datasets, where the noise variance is typically unknown, implies defining a joint prior on two quantities, the IGMRF precision and the noise precision.

The plan of the article is as follows. In Section 2.1, the Bayesian P-spline approach is revised with focus on the challenges to be addressed in defining a prior for τ . The principles behind the construction of a PC prior for τ are revised in Section 3. In Section 4, the PC prior for d is derived and its properties are described in the case of known noise variance. In Section 5, we show how to implement the PC prior for d when the noise variance is unknown, focusing on an additive P-splines model framework. Results from a simulation study assessing the impact of the proposed prior compared to standard gamma priors and other PC priors proposed in the recent literature are illustrated in Section 6. An application of these new priors in a P-spline model for nitrate concentrations observed in river Oglio, Lombardy, Italy, is illustrated in Section 7. The article closes with a discussion in Section 8.

2 Background on P-splines

Let y = ( y 1 , . . . , y n ) T be observations of a response variable, x be a continuous covariate, f be a smooth function describing the effect of the covariate on the response and ε be independent errors with zero mean and variance τ ε − 1 . The P-spline model (Eilers and Marx, 1996) is y = f ( x ) + ε , f ( x ) = B β , where B is a n × K basis matrix containing K B-spline functions built on a set of equally spaced (for simplicity) knots within the covariate domain, while β is a K × 1 vector of unknown spline coefficients. The method requires to select a generous number of knots to over fit the data and then add a penalty on β which smoothes adjacent spline coefficients. In the frequentist approach, β is estimated via penalized maximum likelihood, conditional on a tuning parameter regulating the degree of smoothness of f . The optimal tuning can be found via cross-validation (Wood, 2006) or estimated via restricted maximum likelihood in a mixed model representation (Ruppert et al., 2003 chapter 4). P-splines are widely used in generalized additive models (Hastie and Tibshirani, 1990; Wood, 2006) or structured additive regression models (Fahrmeir et al., 2004; 2013). Higher dimensional smooth functions can also be represented as P-splines, using the tensor product of marginal B-spline bases (Currie et al., 2006; Eilers et al., 2006). For a systematic presentation of the different approaches to p-spline regression, see Ruppert et al. (2003); for an excellent review of spline methods and their applications in statistical modelling, see Hastie et al. (2009), Wakefield (2013) and Wood (2006).

2.1 Bayesian P-splines

The Bayesian approach to P-splines (Lang and Brezger, 2004) assumes an IGMRF prior on the spline coefficients,

π ( β | τ β ) = ( 2 π ) − rank ( R ) / 2 ( | τ β R | * ) 1 / 2 exp − τ β 2 β T R β

(2.1)

where the precision Q ( τ β ) is given by τ β R . Matrix R is denoted as the structure of the IGMRF, that is, a K × K sparse matrix with non-zero entries indicating conditional dependencies among the spline coefficients, τ β is a scalar precision hyper-parameter and | τ β R | * is the generalized determinant. Throughout the article, we will assume R = D r T D r , where D r is a ( K − r ) × K matrix such that D r β = Δ r β (Eilers et al., 2006), with Δ r the r th -order difference operator. In this form, R is the structure of an r th -order random walk on β (Rue and Held, 2005 Chapter 3), with rank ( R ) = K − r , where r indicates the order of the IGMRF (2.1).

The IGMRF (2.1) describes deviation from a base model, which is a polynomial of degree r − 1 . The amount of deviation depends on τ β . A fully Bayesian specification requires priors on τ β and τ ε . Since we usually have enough information in the data to estimate τ ε , the prior π ( τ ε ) has less impact on the fit. The hyper-parameter τ β enters at a lower level in the hierarchy, the data bring little information on it and the prior π ( τ β ) can have a substantial impact on the posterior distribution of β and, as a consequence, on the smoothness of f . Therefore, specification of π ( τ β ) should be as consistent as possible with the prior information actually available about the smoothness of the function to be estimated.

The marginal variance of the IGMRF (2.1), given by the diagonal elements of the generalized inverse matrix Σ * = τ β − 1 R − 1 , depends on K . We denote this as the ‘scaling issue’ (Sørbye and Rue, 2014), meaning that the amount of deviation from the base model depends on the number of knots. This is illustrated in Figure 1, where the two panels report the marginal standard deviation of the smooth f ( x ) = B β , for K = { 50 , 100 } . On the other hand, results (not shown here) show that the degree of the B-splines has little or no impact on the marginal variance of B β , especially when K is large enough, say K > 50 .

2.2 Degrees of freedom

The scaling issue can be avoided if we consider building priors on the number of effective degrees of freedom (Hastie and Tibshirani, 1990), d = tr B T B + τ β τ ε R − 1 B T B . If we think of the smooth f ( x ) = B β as a polynomial, then d can be thought of as the degree of this polynomial. In presence of prior information on the degree of an equivalent polynomial, it seems to be a sensible approach to design a prior for d , π ( d ) , instead of τ β .

A fundamental issue when building π ( d ) is that d depends on both precisions τ β and τ ε . The former regulates the number of effective degrees of freedom, conditionally on the latter. When τ ε is known, the construction of π ( d ) can be based on the prior π ( τ β ) (see Section 4). When τ ε is unknown, the prior π ( d ) will be specified in terms of the joint π ( τ β | τ ε ) π ( τ ε ) , following a fully Bayesian approach (see Section 5).

The degrees of freedom can be reduced to

d = tr I + τ β τ ε R ( B T B ) − 1 − 1 = ∑ k = 1 K 1 1 + τ β τ ε v k ,

(2.2)

where v 1 , . . . , v K are the eigenvalues of R ( B T B ) − 1 , whose null space has dimension r (the rank deficiency of R ). When the factor τ β / τ ε goes to infinity, we obtain the minimum number of degrees of freedom, d = r . When τ β / τ ε goes to zero, we obtain the maximum number of degrees of freedom, K , corresponding to the most flexible model under the assumed IGMRF.

The prior π ( d ) depends on the eigenvalues of R ( B T B ) − 1 , hence on the choice of B . Hereafter, B will be referred to simply as ‘design’, because it is determined by both the assumptions made by the user (location and number of knots, order of the B-splines) and the assumptions purely made by design (location and number of observations along the covariate domain). Since the degrees of freedom depend on B , it follows that π ( d ) automatically accounts for the design. This will be discussed in detail in Section 4.2.

OPEN IN VIEWER Figure 1:

The scaling issue. The two panels show the marginal standard deviation of f ( x ) = B β for varying dimension of the basis B , K = { 50 , 100 } , where B is a matrix of cubic B-splines (coloured lines) defined over the interval ( 0 , 1 ) and β is an IGMRF of order 2. The standard deviation (black line) is calculated as the squared root diagonal entries of B Σ * B T , with τ β = 1 .

The PC prior framework by Simpson et al. (2014) introduces a new concept for building priors in hierarchical additive models, where the latent structure is given by the sum of a number of model components described by a small number of flexibility parameters. Each model component is seen as a flexible extension of a base model. For instance, τ β is a flexibility parameter for the IGMRF component π ( β | τ β ) and a natural base model corresponds to π ( β | τ β − 1 = 0 ) . Below, the four principles underpinning the construction of a PC prior for τ β are reviewed, following Simpson et al. (2014).

Parsimony: A simple model should be preferred unless there is enough evidence for a more flexible one. Under this principle, the prior probability mass assigned to models of increasing complexity should decay as their distance from the base model (measured in terms of model complexity) increases. In Bayesian P-splines, the IGMRF prior operates on β (the object of inference), but we can extend the notion of base and flexible model to the ‘spline-modelled’ function; we denote with f 0 = B β 0 the base model, which is a polynomial of degree r − 1 , and with f = B β the flexible model, which reflects any deviation from such polynomial.

Constant rate penalization: Flexible models are penalized by a constant decay rate for increasing δ . Following this principle, the PC prior is defined as an exponential distribution on the distance scale, π PC ( δ ) = λ exp ( − λ δ ) , with constant rate λ . It follows that the mode of a PC prior is always at the base model. By a change of variable and setting the rate λ = θ / K τ β 0 , Simpson et al. (2014) obtain the PC prior for τ β as,

π PC ( τ β ) = λ exp − λ τ β 0 K / τ β ∂ τ β 0 K / τ β ∂ τ β = θ 2 τ β − 3 / 2 exp − θ / τ β ,

(3.1)

which is a Gumbel ( 1 / 2 , θ ) type 2 distribution, θ > 0 .

PC priors can be helpful as ‘default’ priors in complex hierarchical models where, typically, ‘it is difficult to elicit information about structural parameters that are further down the model hierarchy’(Simpson et al., 2014). In addition, the user-defined scaling approach enables to build ‘informative’ priors for the original parameter τ β or for a property of the associated model component, by tuning two intuitive parameters U and α . In the next section, we introduce a new scaling approach to derive the PC prior for the degrees of freedom of a P-spline model component B β . Other approaches might be possible though. For instance, recently Klein and Kneib (2015) proposed PC priors for the scale (or range of variation) of B β and showed via simulation that these outperformed the gamma family in cases where the data are weakly informative and/or the size of the effects is close to the base model.

So far we have worked under the assumption of known noise precision τ ε . The number of effective degrees of freedom is a function of τ β , which scales differently according to the level of noise present in the data (see the right panel in Figure 2). Knowing τ ε is then crucial to scale the PC prior correctly, in order to guarantee the upper bound for degrees of freedom specified by the user.

The noise precision is typically unknown in applications. One could estimate τ ε from the data and then specify the PC prior for d conditional on this estimate; this strategy has been proposed by Fong et al. (2010) to define gamma-induced priors for degrees of freedom. We, instead, adopt a fully Bayesian model,

y | β , τ ε ∼ N ( B β , τ ε − 1 I ) ,

(5.1)

β | τ β ∼ N ( 0 , τ β − 1 R − 1 ) ,

(5.2)

τ β | τ ε ∼ Gumbel ( 1 / 2 , θ ( τ ε ) ) ,

(5.3)

τ ε ∼ π ( τ ε ) ∝ 1 / τ ε ,

(5.4)

We developed a Markov chain Monte Carlo (MCMC) algorithm to fit the model in equations (5.1)–(5.4); see pseudo-code reported in Appendix A.2 algorithm 1. The algorithm includes a Metropolis–Hastings step to jointly update ( τ ε , τ β , β ) . In our experience, block updating ensures good mixing and fast convergence of the proposed MCMC algorithms; see Rue and Held (2005) for details on block updating in hierarchical models with GMRF components. A brief description on how algorithm 1 works follows. At iteration j , both precision parameters (here denoted simply as τ ) are sampled from the proposal distribution adopted in Knorr-Held and Rue (2002): q ( τ * | τ ( j − 1 ) ) = t τ ( j − 1 ) , where τ * and τ ( j − 1 ) are, respectively, the proposed and current values at iteration j , t is random with density π ( t ) ∝ 1 + 1 / t , for t ∈ 1 / T , T and T > 1 is a tuning parameter; in our experience, setting T approximately equal to 1.5 works well in most applications. The proposal q ( · ) has the advantage that the ratio q ( τ * | τ ( j − 1 ) ) / q ( τ ( j − 1 ) | τ * ) equals one; for more details on this, see Rue and Held (2005), Chapter 4.2. Given τ ε * and τ β * , we draw β * from the full conditional π ( β | τ ε * , τ β * , y ) . Within this scheme, the acceptance probability for ( τ ε * , τ β * , β * ) simplifies to (dropping the superscript ( j − 1 ) to ease the notation)

a = min 1 , π ( τ ε * , τ β * , β * | y ) π ( τ ε , τ β , β | y ) π ( β | τ ε , τ β , y ) π ( β * | τ ε * , τ β * , y ) q ( τ ε | τ ε * ) q ( τ ε * | τ ε ) q ( τ β | τ β * ) q ( τ β * | τ β ) = min 1 , π ( τ ε * , τ β * | y ) π ( τ ε , τ β | y ) ,

(5.5)

5.1 Additive P-splines

We now focus on an additive P-spline modelling framework, where the linear predictor is the sum of a number of smooth functions. Let y be a Gaussian response and x j , j = 1 , . . . , J , be a set of J continuous covariates; the model is

y = ∑ j = 1 J f j ( x j ) + ε ; ε ∼ N ( 0 , τ ε − 1 ) f j ( x j ) = B j β j , j = 1 , . . . , J β j | τ β j ∼ N ( 0 , τ β j − 1 R − 1 ) ,

where B j is the n × K j B-spline basis matrix and β j the vector of spline coefficients associated to the smooth function f j ; with no loss of generality, we consider the same number of knots ∀ j , yielding K = K j , j = 1 , . . . , J . We assume the joint prior ∏ j = 1 J π PC ( τ β j | τ ε ) π ( τ ε ) , where π ( τ ε ) ∝ 1 / τ ε and π PC ( τ β j | τ ε ) = Gumbel ( 1 / 2 , θ j ( τ ε ) ) , j = 1 , . . . , J . The scaling parameter θ j ( τ ε ) = − log ( α j ) d j − 1 ( U j | τ ε ) , where ( U j , α j ) are the parameters calibrating the induced PC prior for the degrees of freedom of f j . The mapping for the degrees of freedom of f j is given by d j ( τ β j | τ ε ) = tr B j T B j + τ β j τ ε R − 1 B j T B j .

Identifiability constraints are important in additive P-splines. The IGMRF prior on β j controls deviations of the smooth term B j β j from a polynomial base model; in the following, for simplicity, we consider an IGMRF of order 2 (which forces shrinkage towards a linear base model). All smooths B j β j include the linear base model, thus they all compete to capture the mean of the data. To ensure identifiability, we adopt the following reparametrization,

y = μ + ∑ j = 1 J x j γ j + ∑ j = 1 J B j β j ULC + ε ,

(5.6)

where μ is the intercept, γ j is the slope coefficient for covariate x j and β j ULC are the spline coefficients β j under the two following linear constraints: c T B j β j = 0 and l T B j β j = 0 , with ‘constant’ vector c = 1 n and ‘line’ vector l = 1 , 2 , . . . , n T . In this way, the smooth term B j β j ULC captures residual variations from the linear base model μ + γ j x j . In other words, the constrained model (5.6) allows each smooth component to be identified, by separating the linear and flexible terms which coexist in each B j β j . Model (5.6) can be expressed in compact form as

y = X γ + B β ULC + ε ,

(5.7)

where B = B 1 : . . . : B j is the n × ( KJ ) joint basis matrix, β ULC is the joint vector of spline coefficients subject to the linear constraints, X = 1 n : x 1 : . . . : x J is the n × ( J + 1 ) matrix of covariates with an additional column of ones for the intercept term and γ is the vector of fixed effects. We assume γ ∼ N ( 0 , τ γ I J + 1 ) with a small precision, for example, τ γ = 10 − 4 , as a prior for the fixed effects. Other covariates can be added to X in model (5.7), if we assume them to have a simple linear effect.

We wrote a block updating MCMC algorithm implementing the joint prior in the model described above; pseudo-code is given in Appendix A.2 algorithm 2. This includes Metropolis–Hastings steps to jointly update the blocks ( τ ε , γ ) and ( τ β j , β j ULC ) , for each j = 1 , . . . , J ; in our experience, this scheme gives good mixing (and convergence) properties. Algorithm 2 presents two main changes with respect to algorithm 1. First, rescaling the conditional PC prior π PC ( τ β j | τ ε ) is no longer necessary at each iteration; θ ( τ ε * ) must be recomputed only when ( τ ε * , γ * ) is accepted. Second, the spline coefficients are sampled under linear constraints. To do this, we use the algorithm proposed in Rue and Held (2005), Chapter 2, which samples first the unconstrained coefficients and then ‘corrects’ them for the constraints. To compute the acceptance probability for the candidate ( τ β j * , β j * ULC ) in step (14) of algorithm 2, Appendix 10, the full conditional density π ( β * | τ β * , y ) is evaluated at the constrained β j * ULC ; for computational details, see Rue and Held (2005), formula 2.31.

The scaling parameter θ of the PC prior equation (3.1) can be tuned in several ways. Our proposal is to select θ through assumptions on the number of effective degrees of freedom, d , of the function f ( x ) = B β . This approach seems intuitive in the Gaussian case where quantity d relates immediately to the degree of an equivalent polynomial (which an expert user might have prior information about). Moreover, literature on smoothing often refers to degrees of freedom as a way to summarize model complexity. In a recent paper, Klein and Kneib (2015) specify θ through assumptions on the ‘scale’, or range of variation, of f ( x ) = B β and denote this PC prior as ‘scale-dependent prior’ (SD prior). Precisely, θ is derived numerically by requiring that Pr ( | f ( x ) | ≤ c ) ≥ 1 − α for each x in the covariate domain, where α ∈ ( 0 , 1 ) and c indicates an upper bound for the scale of f . Both scaling approaches lead to a PC prior which is invariant under design, in the sense that the computation of θ accounts for the adopted B-spline design B . However, the two PC priors differ regarding conditioning on the noise variance. Our PC prior, π PC ( d ) , is a joint prior π ( τ β | τ ε ) = π PC ( τ β | τ ε ) π ( τ ε ) , while the SD prior is defined unconditionally on τ ε as the scale of f does not depend on τ ε .

In this section, we present a simulation study which investigates further the relevance of degrees of freedom for designing priors for Gaussian P-splines. The objective of our study is to assess the behaviour of our joint prior in scenarios with different noise levels, and to compare this with two alternative priors which are defined unconditionally on τ ε , namely the conjugate gamma prior and the SD prior. Therefore, our simulation study does not aim to generically assess the behaviour of PC priors compared to standard priors. For an extended simulation study evaluating the performance of PC priors (particularly SD priors) compared to several alternative hyper-priors for variance parameters, in both Gaussian and non-gaussian contexts, the reader is referred to Klein and Kneib (2015). Furthermore, our simulation is restricted to the Gaussian case.

We consider two different models regarding the shape of the true f ( x ) :

f 1 ( x ) = sin ( x ) ; x ∈ ( − 1 , 1 )

Different scenarios are generated by setting: n = { 20 , 50 } (small and moderate sample sizes), K = { 20 , 30 } , τ ε = { 0.25 , 1 , 5 } (high, moderate and low noise). We aim to assess the model fit obtained by the following priors:

conjugate gamma prior on τ β , with two specifications widely used in applications: gamma ( 10 − 3 , 10 − 3 ) and gamma ( 1 , 5 · 10 − 4 ) .

6.1 Results

In each scenario, goodness of fit was assessed for each of 1 000 simulated datasets by the mean squared error (MSE) of f ̂ ( x ) = B β ̂ , where β ̂ is the posterior mean, as MSE = n − 1 ∑ i = 1 n ( f ̂ ( x i ) − f ( x i ) ) 2 . The posterior for β was computed using R-INLA (Rue et al., 2009, www.r-inla.org) for the gamma and SD priors, and using MCMC algorithm 1 for the joint prior. For the sake of comparison between the three classes of priors, we assume π ( τ ε ) = gamma ( 1 , 5 ⋅ 10 − 4 ) throughout the simulation study (hence, our joint prior is π PC ( τ β | τ ε ) gamma ( 1 , 5 · 10 − 4 ) ).

Figure 5 reports log ( MSE ) for small sample size n = 20 (which is when the hyper-prior is expected to be most influential on the posterior) and K = 20 . We do not see much change for increasing n and K , thus results for other scenarios are only reported in the supplemental material. Our main findings are as follows:

The gamma priors are generally outperformed by the two PC priors (both for degrees of freedom and scale), unless the data are very informative about the true model (e.g., in model f 2 with low noise, see bottom-right panel in Figure 5). As expected, the gamma ( 10 − 3 , 10 − 3 ) overfits when model f 1 is the true one (see left panels in Figure 5), while the gamma ( 1 , 5 · 10 − 4 ) performs poorly in scenario f 2 (especially with high noise, see top-right panel in Figure 5).

For sensible choices of the upper bound U , the joint prior performs better than, or at least as good as, the SD prior. The main difference is noticed in scenario f 2 with high noise (top-right panel in Figure 5) and scenario f 1 with low noise (bottom-left panel in Figure 5). In the former case, setting the joint prior with U = { 5 , 7 , 10 } outperforms most SD prior specifications; in particular, the SD prior with c = { 1.5 , 2 } achieves poor performance because the implied upper bound for d is too small (i.e., below 3 , see Table 1) for this case, where the true effect is highly non-linear and the data are very noisy. The SD prior with c = 3 gives similar performance because it implies a larger upper bound for degrees of freedom. Instead, in scenario f 1 with low noise, SD priors are outperformed by the joint prior with U = { 2 , 3 } , because any choice c = { 1.5 , 2 , 3 } implies an upper bound for d clearly larger than needed (i.e., above 4 , see Table 1) in this case, where data are very informative and the true effect is close to the base model.

In cases where the choice of the upper bound U is inappropriate to describe the complexity of the true function f , the joint prior performs poorly. For instance, when the true curve is close to the base model (i.e., scenario f 1 ), the joint prior with U ≥ 5 is outperformed by both the SD priors and the gamma ( 1 , 5 · 10 − 4 ) (see left panels in Figure 5); in other words, when data are generated under a linear model, a joint prior with large upper bound U will result in a bad choice, as it allows ‘far more’ degrees of freedom than needed. For the same reasons, when the true curve is highly non-linear (i.e., scenario f 2 ), the joint prior with U = { 2 , 3 } generally achieves poor performance (see the right panels in Figure 5), because it assigns almost all weight to the base model or close to it.

In summary, the simulation study showed that both PC priors (i.e., the SD prior and our joint prior) provide potentially better performance than the standard gamma priors for sensible choices of the scaling parameter θ . None of these PC priors shows uniformly better performance than the other one. However, we think that the particular scaling approach used to tune θ plays an important role. Scaling the PC prior in terms of degrees of freedom automatically accounts for the level of noise present in the data. We were able to show via simulation that our joint prior has the potential to perform well in general, that is, both in high and low noise cases, provided that the elicited U is an appropriate upper bound for the degrees of freedom of the function f underlying the data. Therefore, we conclude that the number of effective degrees of freedom is a relevant quantity for building PC priors for smoothing Gaussian data, more suitable than other transformations of τ β that do not depend on τ ε .

OPEN IN VIEWER

Figure 5:

Simulation results: log ( MSE ) for f 1 (left panels) and f 2 (right panels), in presence of high noise ( τ ε = 0.25 , top panels) and low noise ( τ ε = 5 , bottom panels), sample size n = 20 , K = 20 . In the legend on the right, label ‘G’ indicates the Gamma prior; ‘PC_d’ indicates our PC prior for degrees of freedom (joint prior), with α = 0.01 and U = { 2 , 3 , 5 , 7 , 10 } ; ‘SD’ denotes scale dependent prior with α = 0.01 and c = { 1.5 , 2 , 3 } .

We demonstrate the use of the joint prior within an additive P-spline framework for modelling nitrate concentration in river Oglio, Lombardy, Italy. A total of n = 576 observations of NO 3 − concentration were collected during 2010–2012 by taking one sample at each season (spring, summer, autumn and winter) in 48 gauging stations located along the river catchment. The response variable is log ( NO 3 − ij ) , measured at station i = 1 , . . . , 48 and season j = 1 , . . . , 4 . Covariate stream_i is the distance from each station i to the Iseo lake (i.e., the river source) measured in km along the stream network; the first station in proximity to the lake is at stream ≈ 0 , while the last station downstream the river is at stream ≈ 150 . The goal is to understand ‘river enrichment’ in terms of nitrates, by studying the behaviour of log ( NO 3 − ) as the stream distance increases. A substantial amount of information comes from previous studies (see Bartoli et al., 2012; Delconte et al., 2014 and references therein), suggesting that river enrichment in terms of nitrates may vary non-linearly as the stream distance increases. Due to different characteristics in terms of groundwater interactions and irrigation practices, the river catchment can be divided into upstream, middle and downstream reach. Different processes are expected within the three reaches and between seasons; hence, the enrichment curve may show different shapes in the three river segments and seasons.

In order to investigate possible seasonal effects on river enrichment, we adopt the model: log ( NO 3 − ij ) = μ + γ j + f j ( stream i ) + ε ij , i = 1 , . . . , 48 , j = 1 , . . . , 4 , where μ is the overall intercept, γ j is the season-specific intercept and f j ( stream i ) is the season-specific smooth function of the stream distance, modelled with a P-spline with joint prior on ( τ β j , τ ε ) as described in Section 5,

f j ( stream ) = B j β j j = 1 , . . . , 4 β j | τ β j ∼ N ( 0 , τ β j − 1 R − 1 ) τ β j | τ ε ∼ Gumbel ( 1 / 2 , θ j ( τ ε ) ) τ ε ∼ π ( τ ε ) ∝ 1 / τ ε .

(7.1)

To construct the matrices X and B of model (5.5), we need to create suitable dummy vectors of length n : dummy j , j = 2 , . . . , 4 , taking value 1 when the observation is from season j and 0 elsewhere (these are associated to the season-specific intercepts); stream × dummy j , j = 1 , . . . , 4 , taking the actual stream distance when the observation is from season j and 0 elsewhere (these are associated to the season-specific slopes). The n × 8 fixed effect design matrix is X = [ 1 , D , S ] , with D = [ dummy 2 , . . . , dummy 4 ] and S = [ stream × dummy 1 , . . . , stream × dummy 4 ] . Each basis B j contains K = 30 cubic B-splines, evaluated on the season-specific stream distances stream × dummy j . The n × ( 4 K ) B-spline design matrix is B = [ B 1 , B 2 , B 3 , B 4 ] . To separate the season-specific slopes and intercepts from the season-specific smooth variation captured by the B-splines, suitable linear constraints must be applied to each f j , as discussed in Section 5.1. We fit this model via MCMC based on the pseudo-code reported in Appendix A.2 algorithm 2; in particular, we use a two blocks-updating scheme, which ensured good mixing: first block is ( τ ε , γ ) and the other one contains all the four sets of spline coefficients and associated precisions, ( τ β 1 , . . . , τ β 4 , β 1 ULC , . . . , β 4 ULC ) . During MCMC iterations, the PC prior for each precision τ β j needs to be rescaled according to the current τ ε , but this can be done at a negligible computational cost. The R code is provided as supplemental material.

The results displayed in Figure 6 reveal different river enrichment curves ( f ̂ j ( stream ) ) in different seasons: A distinctive pattern is observed in summer, where the fitted curve shows a fast increase upstream ( stream < 50 km) and tendency to decrease downstream the river ( stream > 80 km). This pattern supports the argument given in Delconte et al. (2014) about an upstream reach (from 0 to 25 km) where nitrates are stable, reflecting the chemistry of the lake; a middle reach (from 25 to 80 km) showing an increase of NO 3 − concentration, probably due to groundwater inputs as replacement of river water abstracted for irrigation; a downstream reach (from 80 to 150 km) where nitrates should remain constant or even decrease, mainly due to the dilution of river water with NO 3 − -deprived inputs.

OPEN IN VIEWER Figure 6

Estimated river enrichment curves ( f ̂ j ( stream ) , black line), for each season, and 95 % credible bands (grey). The curve for summer shows clearly a distinctive pattern with respect to the other seasons. Partial residuals (Wood, 2006) are plotted in each panel (dots), indicating larger variability, than what assumed by the model, for the log NO 3 − concentrations observed in spring and autumn (with possible outliers in the upstream river segment in autumn).

PC priors are defined to penalize complexity with respect to a given base model, the magnitude of the penalty being elicited by the user using an intuitive scaling approach. The scaling tool allows the user to derive the PC prior on an interpretable scale, different from the scale of the original parameter, provided that the link between the two is known. We took advantage of this nice feature and derived PC priors for the number of effective degrees of freedom d of a P-spline model for Gaussian data.

For non-Gaussian responses, the idea presented in this article follows straightforwardly by assuming the definition of the degrees of freedom of a generalized P-spline model (Hastie and Tibshirani, 1990), d = tr B T W B + τ β τ ε R − 1 B T W B , where W is a diagonal matrix, with entries depending on the linear predictor of the model (i.e., on B β ) and the adopted link function. MCMC methods similar to those proposed in this article can then be developed being careful about implementing the mapping d ( τ β | τ ε , W ) , which is conditional on W in the generalized case (note that for most distributions in the exponential family, τ ε is known, e.g., for Poisson we have τ ε = 1 ). For Poisson and binomial responses, a convenient approach is to use auxiliary variable methods (Frühwirth-Schnatter and Frühwirth, 2007; Frühwirth-Schnatter et al., 2008) and work with an equivalent ‘augmented’ P-spline model for Gaussian (pseudo) data. The PC prior for d can then be implemented in the same way as described in this article, with algorithms 1 and 2 needing only the inclusion of a Gibbs step to update the augmented parameters.

The potential advantages of using PC priors for degrees of freedom are twofold. First, they are ‘easy-to-elicit’ by the user, who has to define two intuitive scaling parameters: U , an upper bound for d , and α , the prior mass assigned to d > U . This scaling tool can be handled flexibly. For instance, elicitation of the median M for the degrees of freedom results from fixing α = 0.5 . In this case, the PC prior density is bimodal: one mode is set at the base model (by definition) and another mode is set around M degrees of freedom. This bimodal behaviour is due to the attraction to the base model implicit in PC priors.

As a second advantage, these PC priors avoid overfitting and are invariant under design, which means that the parameters U and α do not need to be rescaled if the design changes. In other words, the PC prior is able to code into the model the prior knowledge on the complexity of the curve, or its degrees of freedom, in a design-adaptive way. The ability to adapt to design and avoid overfitting by construction makes π PC ( d ) an appealing default choice in additive models where the latent structure includes several smooth functions (built on a basis of B-splines, e.g., P-splines) and other types of structures, such as individual, spatial and spatio-temporal random effects.