Penalty parameter selection and asymmetry corrections to Laplace approximations in Bayesian P-splines models

Abstract

Laplace P-splines (LPS) combine the P-splines smoother and the Laplace approximation in a unifying framework for fast and flexible inference under the Bayesian paradigm. The Gaussian Markov random field prior assumed for penalized parameters and the Bernstein-von Mises theorem typically ensure a razor-sharp accuracy of the Laplace approximation to the posterior distribution of these quantities. This accuracy can be seriously compromised for some unpenalized parameters, especially when the information synthesized by the prior and the likelihood is sparse. Therefore, we propose a refined version of the LPS methodology by splitting the parameter space in two subsets. The first set involves parameters for which the joint posterior distribution is approached from a non-Gaussian perspective with an approximation scheme tailored to capture asymmetric patterns, while the posterior distribution for the penalized parameters in the complementary set undergoes the LPS treatment with Laplace approximations. As such, the dichotomization of the parameter space provides the necessary structure for a separate treatment of model parameters, yielding improved estimation accuracy as compared to a setting where posterior quantities are uniformly handled with Laplace. In addition, the proposed enriched version of LPS remains entirely sampling-free, so that it operates at a computing speed that is far from reach to any existing Markov chain Monte Carlo approach. The methodology is illustrated on the additive proportional odds model with an application on ordinal survey data.

Keywords

Additive model Bayesian P-splines Laplace approximation Skewness

1 Introduction

By publishing his Mémoire sur la probabilité des causes par les événements (Laplace, 1774), the young French polymath Pierre-Simon de Laplace (1749–1827) seeded an idea today known as the Laplace approximation. At that time, Laplace probably could not have imagined that almost two centuries later, his approximation technique would be resurrected (see e.g., Leonard, 1982; Tierney and Kadane, 1986; Rue et al., 2009) to play a pivotal role in the modern Bayesian literature. Essentially, the Laplace approximation is a Gaussian distribution centered at the maximum a posteriori (MAP) of the target distribution with a variance–covariance matrix that coincides with the inverse of the negative Hessian of the log-posterior target evaluated at the MAP estimate. Recently, the Laplace approximation has crossed the path of P-splines, the brainchild of Paul Eilers and Brian Marx (Eilers and Marx, 1996), to inaugurate a new approximate Bayesian methodology labelled as ‘Laplace P-splines’ (LPS) with promising applications in survival analysis (Gressani and Lambert, 2018, Gressani et al., 2022b; Lambert and Kreyenfeld, 2023), generalized additive models (Gressani and Lambert, 2021), nonparametric double additive location-scale models for censored data (Lambert, 2021) and infectious disease epidemiology (Gressani et al., 2022a,c). The sampling-free inference scheme delivered by Laplace approximations combined with the possibility of smoothing different model components with P-splines in a flexible fashion paves the way for a robust and much faster alternative to existing simulation-based methods.

Although LPS shares some methodological aspects with the popular integrated nested Laplace approximations (INLA) approach (Rue et al., 2009), there are fundamental points of divergence. First, the tools in INLA and its associated R-INLA software are originally built to compute approximate posteriors of univariate latent variables, contrary to LPS that natively delivers approximations to the (multivariate) joint posterior distribution of the latent vector. The key benefit of working with an approximate version of the joint posterior is that pointwise estimators and credible intervals for subsets of the latent vector (and functions thereof) can be straightforwardly constructed. Second, by working with closed-form expressions for the gradients and Hessians involved in the model, LPS is computationally more efficient than the numerical differentiation proposed in INLA. Third, while INLA can be combined with various techniques for smoothing nonlinear model components, LPS is entirely devoted to P-splines smoothers with the key advantage of having full control over the penalization scheme (as the approximate posterior distribution of the penalty parameter(s) is analytically available). In this regard, LPS has more in common with the work of Wood and Fasiolo (2017) than with INLA, especially in the class of (generalized) additive models (Wood, 2017).

The success of Laplace approximations in Bayesian statistics owes much to a central limit type argument. Under certain regularity conditions, the Bernstein-von Mises theorem (see e.g., Van der Vaart, 1998) ensures that posterior distributions in differentiable models converge to a Gaussian distribution under large samples. In situations involving small to medium sample sizes, the suitability of the Laplace approximation can be questioned as it does not take into account the potential skewness or kurtosis of posterior distributions (Ruli et al., 2016). Even under relatively large samples, the Laplace approximation might fail in scenarios involving binary data as the latter are poorly informative for the model parameters and can result in a flat log-likelihood function, thus complicating inference (Ferkingstad and Rue, 2015; Gressani and Lambert, 2021).

Laplace P-splines belong to the class of latent Gaussian models, where model parameters are dichotomized between a vector of latent variables $ξ$ (including penalized B-spline coefficients, regression coefficients and other parameters of interest) that are assigned a Gaussian prior and another vector of hyperparameters $η$ that involves nuisance parameters, such as the smoothing parameter inherent to P-splines, and for which prior assumptions need not be Gaussian. Combining Bayes’ rule and a simplified Laplace approximation, the conditional posterior distribution of $ξ$ under the LPS framework is approximated by a Gaussian distribution denoted by ${\tilde{p}}_{G} (ξ | \hat{η}, D)$ , where $\overset{⌢}{η}$ is a summary statistic of the posterior hyperparameter vector (e.g., the MAP estimate, the posterior mean or median) and $D$ denotes the observed data. Although the latter approximation is typically accurate for penalized B-spline coefficients, it might be less appropriate for other candidates in $ξ$ with large prior variance. In that case, the misfit between the Laplace approximation and a potentially asymmetric (or heavy-tailed) target posterior distribution for a parameter can have a detrimental effect on posterior summary statistics and on any results relying on the generated approximation for the posterior distribution of the model parameters. This motivates us to develop an approach that corrects for potential posterior misfits provided by the Laplace approximation.

A recent technique proposed by Chiuchiolo et al. (2022) in the INLA framework consists in using a skew Gaussian copula to correct for skewness when posterior latent variables have a non-negligible deviation from Gaussianity. Our proposal in models involving P-splines consists in splitting the latent vector $ξ$ into a set of parameters $γ$ for which the posterior distribution (conditional on the hyperparameters) is approximated in a non-Gaussian fashion with an emphasis on capturing asymmetries, and a set of parameters $θ$ (that typically involves penalized B-spline coefficients) for which the conditional posterior is approached with Laplace approximations. Our refined LPS approach thus allows to obtain an approximation to the joint posterior distribution of $ξ$ (given $η$ ) together with an approximation to the posterior of the hyperparameters $η$ without relying on an MCMC sampling scheme.

A simple motivating example inspired by the infectious disease model of Gressani et al. (2022c) helps framing the problem. Let $D = {(x_{i}, y_{i}) : i = 1, \dots, n}$ be a sample of $n = 120$ independent pairs where $x_{i} = i$ and $y_{i}$ has a negative binomial distribution $NB (μ (x_{i}), ϕ)$ following the parametrization of Piegorsch (1990) with mean $E (y | x) = μ (x)$ , variance $V (y | x) = μ (x) + ϕ μ {(x)}^{2}$ and overdispersion parameter $ϕ > 0$ . We model the conditional mean with P-splines $\log μ (x) = θ^{⊤} b (x)$ , where $b (\cdot)$ is a cubic B-spline basis on the interval $[1,120]$ and $θ$ is a vector of B-spline coefficients. Following Lang and Brezger (2004), a global smoothness prior is assumed for the B-spline parameters, $p (θ | λ) \propto \exp (- \frac{λ}{2} θ^{⊤} P θ)$ , where $P$ is a penalty matrix and $λ$ is the penalty parameter to which we assign a weakly informative Gamma prior (with mean $a / b$ and variance $a / b^{2}$ ), denoted by $λ \sim G (a, b)$ . To complete the model specification, a Gamma prior with large variance is assumed for the overdispersion parameter $ϕ^{- 1}$ . The left panel of Figure 1 shows a data set of size $n = 120$ simulated from the above negative binomial model with a nonlinear function for $μ (x)$ and $ϕ^{- 1} = 6$ . The histogram for $\log (ϕ)$ on the right panel of Figure 1 is obtained from a long MCMC chain with a Metropolis-within-Gibbs algorithm. There is an important misfit between the Laplace approximation (solid curve) and the MCMC output, so that quantities like the posterior standard deviation or selected posterior quantiles for $\log ϕ$ will be poorly estimated using a straightforward Laplace approximation to $p (ϕ | D)$ . The dashed curve represents the skewed distribution that we propose as an alternative candidate to the Laplace approximation and that will be thoroughly discussed in the next section within a Bayesian P-splines context. The fitted distribution is able to capture the asymmetry that is apparent in the MCMC sample, improving the precision of posterior estimates for $\log ϕ$ as compared to Laplace and at a much lower computational cost than MCMC. The article is organized as follows. Section 2 presents the Bayesian Laplace P-spline model and gives a detailed description of the proposed asymmetric posterior approximation methodology for non-penalized parameters. In Section 3, we illustrate the method in an additive proportional odds model for ordinal data. Finally, Section 4 concludes with a discussion.

Figure 1

Left panel: Count data ( $n = 120$ ) generated using a negative binomial $NB (μ (x), ϕ)$ with $ϕ^{- 1} = 6$ . Right panel: Histogram of a MCMC sample for $\log ϕ$ compared to Laplace (solid) and skew- $t$ (dashed) approximations.

2 Laplace approximation and Bayesian P-splines

2.1 Model specification

Consider a regression model describing the conditional distribution of a response $y$ for given covariates $x$ . Denote by $ξ$ the model parameters: it includes the regression and spline parameters, plus possibly the (log of the) scale and (unconstrained transformed) shape parameters. Denote by $p (ξ | η)$ the joint prior density of $ξ$ conditionally on a vector of hyperparameters $η$ . In the context of a P-spline model, the latter parameters can come from $J$ unknown smooth functions specified as $f_{j} (\cdot) = \sum_{l = 1}^{L} θ_{l j} b_{j l} (\cdot)$ ( $j = 1, \dots, J$ ) where $B_{j} = {b_{j l} (\cdot) : l = 1, \dots, L}$ denotes a B-spline basis with equidistant knots spanning the argument range. Vector $η$ would typically include positive roughness penalty parameters $λ_{j}$ with prior density $p (λ_{j})$ for the $j$ th function. The frequentist penalty on changes in differences of neighbour spline parameters (Eilers and Marx, 1996) can be translated in a Bayesian context using a conditional prior on $θ_{j} = (θ_{j 1}, \dots, θ_{j L})^{⊤}$ (Lang and Brezger, 2004), $p (θ_{j} | λ_{j}) \propto \exp (- 12 θ_{j}^{⊤} (λ_{j} P) θ_{j})$ , with $P = D_{r}^{⊤} D_{r}$ denoting a penalty matrix corresponding to a finite difference penalty matrix $D_{r}$ of order $r$ . For example, when $r = 2$ , one has $θ_{j}^{⊤} P θ_{j} = | | D_{r} θ_{j} {| |}_{2}^{2} = \sum_{l = 1}^{L - 2} {(θ_{l + 2, j} - 2 θ_{l + 1, j} + θ_{l, j})}^{2}$ . The penalty parameter $λ_{j} > 0$ is used to tune the smoothness of the associated additive term with, at the limit when $λ_{j} \to + \infty$ , a polynomial of order $r - 1$ for $f_{j} (\cdot)$ . Different prior distributions could be chosen for $λ_{j}$ with Brezger and Lang (2006) suggesting to take Gamma priors $λ_{j} \sim G (a_{j}, b_{j})$ (with mean $a_{j} / b_{j}$ and variance $a_{j} / b_{j}^{2}$ ). A small value for $b_{j}$ ( ${= 10}^{- 4}$ , say) combined with $b_{j} = a_{j}$ or $a_{j} = 1$ ensures a large prior variance with some more weight set on small or large values of $λ_{j}$ , respectively. Mixtures of Gamma densities were also investigated in Jullion and Lambert (2007) with $(λ | δ) \sim G (ν / 2, ν δ / 2)$ and $δ \sim G (a_{δ}, b_{δ})$ yielding, in the special case $a_{δ} = b_{δ} = .5$ and $ν = 1$ , a Beta prime distribution $B' (.5,.5)$ for $λ$ or equivalently a half-Cauchy prior for $\sqrt{λ}$ with $p (λ) \propto λ^{- .5} {(1 + λ)}^{- 1}$ (Lambert and Bremhorst, 2019).

Penalties can also be combined and extended in multiple ways, see, for example, the book by Eilers and Marx (2021) for inspiring examples. More generally, we assume that the joint conditional prior for the vector $θ$ , stacking all the vectors of penalized B-spline coefficients in the model, can be written as

p (θ | λ) \propto \exp (- \frac{1}{2} θ^{⊤} 𝒫_{λ} θ),

(2.1)

where $𝒫_{λ}$ is a positive semi-definite matrix. The vector of model parameters can be reorganized as follows, $ξ = (γ^{⊤}, θ^{⊤})^{⊤} \in I R^{k_{1} + k_{2}}$ , where $γ \in I R^{k_{1}}$ denotes the vector of non-penalized parameters. If $D$ generically denotes the available data and if $λ$ stands for the vector of hyperparameters $η$ in the specific context of P-spline models, then the joint posterior for $ξ$ directly follows from Bayes’ theorem,

p (ξ, λ | D) \propto L (ξ | D) p (γ) p (θ | λ) p (λ),

where $L (ξ | D)$ denotes the likelihood. It is typically explored using Markov chain Monte Carlo methods (MCMC). In this article, we build up on the methodology described in Gressani and Lambert (2018, 2021) and in Lambert (2021), where Laplace approximations to the conditional posterior of $(ξ | λ, D)$ and an additional approximation to the marginal posterior of $(λ | D)$ enable to bypass sampling algorithms, see Section 2.2.

2.2 Laplace approximation and penalty parameter selection

Assume that closed form expressions can be derived for the gradient and Hessian of $\log p (ξ | λ, D)$ ,

U_{λ} = U_{λ} (ξ) = \partial \log p (ξ | λ, D) / \partial ξ; H_{λ} = H_{λ} (ξ) = \partial^{2} \log p (ξ | λ, D) / \partial ξ \partial ξ^{⊤} .

The conditional posterior mode ${\hat{ξ}}_{λ}$ of $ξ$ can be quickly obtained using the Newton-Raphson (NR) algorithm with the substitution, $ξ \leftarrow ξ - H_{λ}^{- 1} U_{λ}$ , repeated until convergence. The Levenberg-Marquardt algorithm (Marquardt, 1963) could be preferred if good initial conditions are not easily found to ensure convergence. A Laplace approximation to the conditional posterior distribution of $ξ$ directly follows: $(ξ | λ, D) \overset{\cdot}{\sim} N_{k_{1} + k_{2}} ({\hat{ξ}}_{λ}, Σ_{λ})$ where $Σ_{λ} = - H_{λ}^{- 1}$ . Thanks to the Gaussian Markov random field (GMRF) prior (Rue and Held, 2005), $p (θ | λ)$ , assumed in (2.1) for the penalized parameters, the Gaussian approximation to the conditional posterior of $(θ | λ, D)$ is usually excellent, see Rue et al. (2009) for the same argument in latent Gaussian models. However this might not be true for some non-penalized parameters in $ξ$ , especially when the combined information coming from their prior and the likelihood is sparse, see Section 2.3 for a specific handling.

The preceding Laplace approximation can be used to approximate the marginal posterior distribution of the penalty parameters $λ$ with the Normal approximation substituted in the denominator of the following identity, $p_{λ} (λ | D) = p (ξ, λ | D) / p (ξ | λ, D)$ , yielding

\begin{matrix} {\tilde{p}}_{λ} (λ | D) = \frac{p (ξ, λ | D)}{{\tilde{p}}_{G} (ξ | λ, D)} \propto p ({\hat{ξ}}_{λ}, λ | D) | {\hat{Σ}}_{λ} |^{\frac{1}{2}} \\ \propto \underset{Marginal likelihood}{\underset{ุ}{L ({\hat{ξ}}_{λ} D) p ({\hat{ξ}}_{λ} | λ, D) | {\hat{Σ}}_{λ} |^{\frac{1}{2}}}} \times p (λ), \end{matrix}

(2.2)

see Tierney and Kadane (1986) for the same strategy in the approximation of a marginal distribution. One might prefer to work with $υ = \log λ$ and its approximate marginal posterior,

{\tilde{p}}_{υ} (υ | D) = {\tilde{p}}_{λ} (e^{υ} | D) \prod_{j} e^{υ_{j}} .

(2.3)

The maximization of (2.2) or of the marginal likelihood (as with ‘empirical Bayes’ methods) can be used to select a specific value for $λ$ . Alternatively, it could be derived from the log-penalty using (2.3), yielding larger penalty values when the selection is made from the marginal posterior instead of the (parametrization invariant) marginal likelihood.

2.3 Asymmetric posterior for non-penalized parameters

Assume that $ξ = (γ^{⊤}, θ^{⊤})^{⊤} \in I R^{k_{1} + k_{2}}$ with $γ \in I R^{k_{1}}$ is suspected to have a non-symmetric marginal posterior distribution. Let ${\hat{ξ}}_{λ} = ({\hat{γ}}_{λ}^{⊤}, {\hat{θ}}_{λ}^{⊤})^{⊤}$ denote the posterior mode of $p (ξ | λ, D)$ and ${\hat{Σ}}_{λ}^{- 1}$ the observed information matrix structured in blocks as follows,

{\hat{Σ}}_{λ} = [\begin{matrix} {\hat{Σ}}_{λ}^{γ γ} & {\hat{Σ}}_{λ}^{γ θ} \\ {\hat{Σ}}_{λ}^{θ γ} & {\hat{Σ}}_{λ}^{θ θ} \end{matrix}] .

The conditional posterior of $(θ | γ, λ, D)$ has an approximate Normal distribution resulting from the GMRF prior in (2.1) for $(θ | λ)$ . One has

(θ | γ, λ, D) \overset{\cdot}{\sim} N_{k_{2}} (E (θ | γ, λ, D), {\hat{Σ}}_{λ}^{θ | γ}),

where

\begin{matrix} E (θ | γ, λ, D) = {\hat{θ}}_{λ} + {\hat{Σ}}_{λ}^{θ γ} {({\hat{Σ}}_{λ}^{γ γ})}^{- 1} (γ - {\hat{γ}}_{λ}); \\ {\hat{Σ}}_{λ}^{θ | γ} = {\hat{Σ}}_{λ}^{θ θ} - {\hat{Σ}}_{λ}^{θ γ} {({\hat{Σ}}_{λ}^{γ γ})}^{- 1} {\hat{Σ}}_{λ}^{γ θ} . \end{matrix}

(2.4)

Hence, starting from the following identity, $p_{γ} (γ | λ, D) = p (γ, θ | λ, D) / p (θ | γ, λ, D)$ , one gets the approximation

p_{γ} (γ | λ, D) \approx \frac{p (γ, θ | λ, D)}{{\tilde{p}}_{G} (θ | γ, λ, D)} \propto p (γ, E (θ | γ, λ, D) | λ, D) | {\hat{Σ}}_{λ}^{θ | γ} |^{\frac{1}{2}},

(2.5)

see Eq. (2) in Tierney et al. (1989) for a similar expression. We propose to reparametrize $γ$ by projecting it on the eigenvectors of the singular value decomposition (SVD) of ${\hat{Σ}}_{λ}^{γ γ} = V Z V^{⊤}$ where $V = [v_{1} \dots v_{k_{1}}]$ denotes the matrix of orthonormal eigenvectors, $ζ$ the eigenvalues and $Z = diag (ζ)$ . It yields $\tilde{γ} = Z^{- \frac{1}{2}} V^{⊤} (γ - {\hat{γ}}_{λ})$ and $γ = {\hat{γ}}_{λ} + V Z^{\frac{1}{2}} \tilde{γ}$ , with

p_{\tilde{γ}} (\tilde{γ} | λ, D) \propto p_{γ} ({\hat{γ}}_{λ} + V Z^{\frac{1}{2}} \tilde{γ} | λ, D).

(2.6)

The posterior dependence between the components of $\tilde{γ}$ is expected to be milder than under the original $γ$ parametrization. Therefore, conditionally on $λ$ , we propose to approximate the joint posterior density of $(\tilde{γ} | λ, D)$ by the product of the marginal densities of its components:

p_{\tilde{γ}} (\tilde{γ} | λ, D) \approx \prod_{s = 1}^{k_{1}} p_{{\tilde{γ}}_{s}} ({\tilde{γ}}_{s} | λ, D).

(2.7)

Under that working independence hypothesis, each univariate marginal in the product in (2.7) is equal to its conditional with the other components set equal to an arbitrary value. Combined with (2.6), it implies that

\begin{matrix} p_{{\tilde{γ}}_{s}} ({\tilde{γ}}_{s} | λ, D) = p_{{\tilde{γ}}_{s} | {\tilde{γ}}_{- s}} ({\tilde{γ}}_{s} | {\tilde{γ}}_{- s} = 0, λ, D) \\ \propto p_{\tilde{γ}} ({\tilde{γ}}_{s} e_{s} | λ, D) \\ \propto p_{γ} ({\hat{γ}}_{λ} + {\tilde{γ}}_{s} \sqrt{ζ_{s}} v_{s} | λ, D), \end{matrix}

(2.8)

where $e_{s}$ denotes the $s$ th unit vector in $I R^{k_{1}}$ such that ${[e_{s}]}_{k} = δ_{k s}$ . The univariate marginal posterior density in (2.8) can be evaluated for any value of ${\tilde{γ}}_{s}$ using (2.5). We suggest to approximate it using a skew-normal (SN) or a skew- $t$ (ST) distribution to also handle kurtosis. Here, we provide details for the ST distribution. By definition, $X \sim ST (ψ, ω^{2}, α, υ)$ if $X \in IR$ and has density

t (x | ψ, ω^{2}, α, υ) = \frac{2}{ω} t_{υ} (\frac{x - ψ}{ω}) T_{υ} (α \frac{x - ψ}{ω}),

with location parameter $ψ \in IR$ , scale parameter $ω > 0$ , slant (or skewness) parameter $α \in IR$ and $υ > 0$ degrees of freedom (d.f.), where $t_{υ}$ , $T_{υ}$ respectively denote the density and c.d.f. of a standard Student distribution with $υ$ d.f., see, for example, Azzalini and Capitanio (2014) for more details on the ST definition and its properties. An approximation to the target distribution can for example be derived by minimizing its Jensen-Shannon divergence from a skew- $t$ distribution, giving $({\tilde{γ}}_{s} | λ, D) \overset{\cdot}{\sim} ST ({\tilde{ψ}}_{s}, {\tilde{ω}}_{s}^{2}, {\tilde{α}}_{s}, {\tilde{υ}}_{s})$ . Substituting these ST densities in (2.7) provides an approximation to the joint posterior of $(\tilde{γ} | λ, D)$ and an efficient method to sample from it using the independence of its components. An analytic form for the joint posterior density of $(γ | λ, D)$ follows from the combination of (2.7) with $\tilde{γ} = Z^{- \frac{1}{2}} V^{⊤} (γ - {\hat{γ}}_{λ})$ , giving

p_{γ} (γ | λ, D) = \prod_{s = 1}^{k_{1}} \frac{1}{\sqrt{ζ_{s}}} t ({\tilde{γ}}_{s} | {\tilde{ψ}}_{s}, {\tilde{ω}}_{s}^{2}, {\tilde{α}}_{s}, {\tilde{υ}}_{s}).

An approximation to the marginal distribution of its $s$ th component $(γ_{s} | λ, D)$ can be obtained using fast Monte Carlo techniques. Indeed, an approximate large random sample from $(\tilde{γ} | λ, D)$ can first be generated by sampling its independent skew- $t$ components $ST ({\tilde{ψ}}_{s}, {\tilde{ω}}_{s}^{2}, {\tilde{α}}_{s}, {\tilde{υ}}_{s})$ , yielding ${{\tilde{γ}}^{(m)} : m = 1, \dots, M}$ . The associated random sample for $(γ | λ, D)$ is given by ${γ^{(m)} = {\hat{γ}}_{λ} + V Z^{\frac{1}{2}} {\tilde{γ}}^{(m)} : m = 1, \dots, M}$ . Then, a skew- $t$ approximation to the marginal posterior of $γ_{s}$ can be fitted to ${γ_{s}^{(m)} : m = 1, \dots, M}$ , yielding $(γ_{s} | λ, D) \overset{\cdot}{\sim} ST (ψ_{s}, ω_{s}^{2}, α_{s}, υ_{s})$ . Point estimates or credible regions for $γ_{s}$ can be computed from it. These different steps provide a convenient approximation to the joint posterior of the model parameters. Indeed, based on the factorization of the joint posterior density, $p (ξ, λ | D) = p (θ | γ, λ, D) p (γ | λ, D) p (λ | D)$ , one has the following stochastic representation for $(ξ, λ | D)$ ,

(ξ, λ | D) \overset{\cdot}{\sim} N_{k_{2}} (E (θ | γ, λ, D), {\hat{Σ}}_{λ}^{θ | γ}) \times \prod_{s = 1}^{k_{1}} ST ({\tilde{γ}}_{s} | {\tilde{ψ}}_{s}, {\tilde{ω}}_{s}^{2}, {\tilde{α}}_{s}, {\tilde{υ}}_{s}) \times (λ | D),

(2.9)

with mean and variance–covariance matrix in the first factor given in (2.4). It can be used to generate an arbitrarily large number of independent copies from the joint posterior much faster than with MCMC. This is for example particularly useful to make inference on complicated functions of the model parameters or for predictive purposes.

3 Illustration

3.1 The additive proportional odds model for ordinal data

Denote by $ε_{m} = {1, \dots, m}$ the set containing the first $m$ positive integers. Assume that $n$ independent units are observed with data $D = {(x_{i}, y_{i}) : i = 1, \dots, n}$ where $y$ is an ordinal random variable taking values in $ε_{R}$ and $x$ a vector of covariates. The proportional odds (PO) model is a popular choice when the response is ordinal (Agresti, 2010). It assumes that

logit [P (Y \leq r | x)] = η_{r} = γ_{r} + x^{⊤} θ (r \in ε_{R - 1}),

with a specific intercept $γ_{r}$ for each cumulative logit, but a shared vector of regression parameters $θ$ . Consequently,

\log \frac{\Pr (Y \leq r | x_{1}) / \Pr (Y > r | x_{1})}{\Pr (Y \leq r | x_{2}) / \Pr (Y > r | x_{2})} = θ (x_{1} - x_{2})

is independent of $r$ and provides a clear interpretation to the regression parameters $θ$ with a change $θ_{k}$ in the log-odds of $Y$ taking values in the lower end of the ordinal scale for every unit increase in the $k$ th component of $x$ . Let $F_{i r} = P (Y_{i} \leq r | x_{i}) = e^{η_{i r}} / (1 + e^{η_{i r}})$ for $r$ in $ε_{R - 1}$ and let $F_{i 0} = 0$ , $F_{i R} = 1$ . Then, $π_{i r} = P (Y_{i} = r | x_{i}) = F_{i r} - F_{i, r - 1}$ for $r$ in $ε_{R}$ . Let $ξ = (γ^{⊤}, θ^{⊤})^{⊤}$ and assume a prior of the following form,

p (ξ | λ) \propto \exp (-_{2}^{1} (ξ - e)^{⊤} K_{λ} (ξ - e)),

(3.1)

conditionally on a vector of parameters $λ$ and for a positive semi-definite matrix $K_{λ}$ . The log-likelihood takes a simple form, $l (ξ | D) = \sum_{i = 1}^{n} \log π_{i y_{i}} (ξ)$ , with the resulting conditional posterior for $ξ$ ,

\log p (ξ | λ, D) = l (ξ | D) -_{2}^{1} (ξ - e)^{⊤} K_{λ} (ξ - e).

(3.2)

Explicit analytical forms can be derived for the associated gradient and Hessian matrix, see Appendix 1. Assume now for simplicity a model with $J$ continuous covariates $x_{1}, \dots, x_{J}$ with smooth additive terms $f_{j} (x_{j})$ ( $j = 1, \dots, J$ ) describing their effects on the conditional log-odds,

logit [P (Y \leq r | x)] = η_{r} = γ_{r} + f_{1} (x_{1}) + \dots + f_{J} (x_{J}) .

Following Eilers and Marx (1996), consider now a basis of $(L + 1)$ cubic B-splines ${s_{j l}^{*} (\cdot)}_{l = 1}^{L + 1}$ associated to a generous number of equally spaced knots on the range $(x_{j}^{\min}, x_{j}^{\max})$ of values for $x_{j}$ (Marx and Eilers, 1998). They are recentered for identification purposes in the additive model using $s_{j l} (\cdot) = s_{j l}^{*} (\cdot) - 1 x_{j}^{\max} - x_{j}^{\min} \int_{x_{j}^{\min}}^{x_{j}^{\max}} s_{j l}^{*} (u) d u (l = 1, \dots, L)$ . Then, the additive terms in the conditional model can be approximated using linear combinations of these recentered B-splines, $f_{j} (x_{j}) = \sum_{l = 1}^{L} s_{j l} (x_{i j}) θ_{l j}$ . In a Bayesian framework, as reminded in Section 2.1, smoothness can be forced on these additive terms by taking GMRF priors for the vectors of spline coefficients

p (θ_{j} | λ_{j}) \propto \exp (-_{2}^{1} θ_{j}^{⊤} (λ_{j} P) θ_{j}),

where $P$ stands for the penalty matrix. A multivariate Normal prior could be taken for $γ$ to complete the model specification, $γ \sim N (\tilde{e}, (Q)^{- 1})$ . Under the general formulation in (3.1), one has $ξ = (γ^{⊤}, θ^{⊤})^{⊤}$ , $θ = (θ_{1}^{⊤}, \dots, θ_{J}^{⊤})^{⊤}$ , $e = ({\tilde{e}}^{⊤}, 0_{J L}^{⊤})^{⊤}$ , and a block-diagonal penalty matrix $K_{λ} = diag (Q, P_{λ})$ where $P_{λ} = Λ \otimes P$ with ${[Λ]}_{j j^{'}} = δ_{j j^{'}} λ_{j}$ . The conditional posterior for $ξ$ is given by (3.2). With a Gamma prior for the penalty parameters, $λ_{j} \sim G (a, b)$ , one has $(λ_{j} | ξ, D) \sim G (a + .5 ρ (P), b + .5 θ_{j}^{⊤} P θ_{j})$ . Starting from these conditional posterior distributions, a Metropolis-within-Gibbs algorithm can be set up to generate a random sample from the joint posterior for $(ξ, λ)$ , with Gibbs steps for the penalty parameters $λ$ and Metropolis steps for the regression and splines parameters $ξ$ . Alternatively, proposals for $ξ$ could be made using the modified Langevin (Roberts and Tweedie, 1996; Lambert and Eilers, 2009) or the Metropolis-Hastings algorithm with proposals based on the local topological information provided by the explicit analytic forms for the gradient and Hessian matrix (Gamerman, 1997). Such a sampling approach based on MCMC will be compared to the strategy proposed in Section 2.

3.2 Application on survey data

Consider now an illustration of the proposed methodology on data coming from the European Social Survey (ESS Round 9, 2018) with a specific focus on the French speaking respondents from Wallonia, one of the three regions in Belgium. Each of the participants (aged at least 15) was asked to react to the following statement, Gay men and lesbians should be free to live their own life as they wish, with a positioning on a Likert scale going from 1 (= Agree strongly) to 5 (= Disagree strongly), with 3 labelled as Neither agree nor disagree (with relative frequencies 1: 54.9% ; 2: 30.4% ; 3: 8.2% ; 4: 5.4% ; 5: 1.1%). That ordinal response effectively recorded on $n = 552$ respondents was analysed using the proportional odds model described above with the number of completed years of education ( $14.1 \pm 4.4$ years) and age ( $47.3 \pm 18.5$ years) entering as additive terms with $L = 10$ recentered B-splines spanning each covariate range. Starting from Gamma priors, $λ_{j} \sim G ({1,10}^{- 4})$ ( $j = 1,2$ ), the penalty parameters $λ_{1}$ and $λ_{2}$ associated to $f_{1} (eduyrs)$ and $f_{2} (age)$ , respectively, were selected by maximizing $p (λ | D)$ in (2.2) using the Levenberg-Marquardt algorithm, yielding ${\hat{λ}}_{1} = 191.8$ (e.d.f.=1.24), ${\hat{λ}}_{2} = 18.4$ (e.d.f.=2.55), the value in brackets standing for the effective degrees of freedom. Alternatively the maximization of the marginal likelihood in (2.2) would yield a very large value for ${\hat{λ}}_{1}$ (e.d.f.=1.00), suggesting linearity for $f_{1} (eduyrs)$ , and ${\hat{λ}}_{2} = 18.5$ (e.d.f.=2.55) practically unchanged. The fitted additive terms are visible on Figure 2 with their pointwise 95% credible intervals, suggesting a statistically non-significant effect of eduyrs. but a tolerant perception of homosexuality tending to decrease with age, with a marked change in attitude revealed beyond age 60.

Figure 2

ESS dataset: fitted additive terms for eduyrs and age with pointwise 95% credible intervals: this suggests a growing hostility to homosexuality beyond the age of 60, while it appears that the number of completed years of education does not play a statistically significant role.

To compare the merits of our proposal, a MCMC algorithm was run to explore $p (ξ | \hat{λ}, D)$ , the generated samples and their properties being compared to the analytical approximations suggested in Section 2.3. The estimated additive terms and their 95% credible intervals based on MCMC are practically identical to our estimates in Figure 2, confirming the excellent quality of the Laplace approximation to the conditional posterior distribution of $(θ | λ, D)$ underlying our calculations. Let us now focus on the non-penalized regression parameter in $γ$ standing for the four intercepts in the proportional odds model. The scatterplot of the MCMC sample ${γ^{(m)} : m = 1, \dots, M}$ generated from $p (γ | \hat{λ}, D)$ using the modified Langevin algorithm can be found in the left panel of Figure 3 where the posterior dependence between the vector components clearly stands out. The reparametrization suggested in Section 2.3 along the principal axes corresponding to the eigenvectors of the SVD decomposition of ${\hat{Σ}}_{λ}^{γ γ}$ yields $\tilde{γ}$ , with the scatterplot of the associated MCMC sample ${{\tilde{γ}}^{(m)} = Z^{- \frac{1}{2}} V^{⊤} (γ^{(m)} - {\hat{γ}}_{λ}) : m = 1, \dots, M}$ , visible in the right panel of Figure 3 confirming that the posterior dependence between the vector components of $\tilde{γ}$ is very mild and probably negligible for most practical purposes.

Figure 3

ESS dataset: scatterplots of the MCMC samples for $(γ | λ, D)$ and $(\tilde{γ} | λ, D)$ when $λ = \hat{λ}$ .

The suggested analytical approximation to $p_{{\tilde{γ}}_{s}} ({\tilde{γ}}_{s} | λ, D)$ in (2.8) was evaluated and added to the MCMC sample taken as a trustful proxy of the true marginal posterior distribution of $({\tilde{γ}}_{s} | λ, D)$ , see Figure 4. The quality of the analytical approximations is excellent with a noticeable left asymmetry for $({\tilde{γ}}_{1} | λ, D)$ in particular. When transformed back to the $γ$ -parametrization, the approximating skew-Normal densities shown in Figure 5 closely match the distribution of the MCMC samples for $(γ_{s} | λ, D)$ . The positive skewness is non-negligible for the marginal posterior distribution of $γ_{4}$ : that asymmetry would not be captured by a simple Laplace approximation. It is caused by the small proportion of respondents in the survey expressing a strong disagreement with the submitted statement on the freedom of gays and lesbians.

Figure 4

ESS dataset: approximated marginal posterior density for $(\tilde{γ} | λ, D)$ compared to MCMC samples when $λ = \hat{λ}$ .

Figure 5

ESS dataset: approximated marginal posterior density for $(γ | λ, D)$ compared to MCMC samples when $λ = \hat{λ}$ .

4 Discussion

In this article, the Laplace P-spline (LPS) approach has been extended to improve the accuracy of inference in a Bayesian framework. Indeed, when information is sparse, the posterior distribution of non-penalized parameters may exhibit a non-negligible skewness that can have adverse effects on inference or predictions when ignored. The proposed approximation to the joint posterior density in (2.9) takes a simple form that can be used in a much faster way than MCMC to make predictions or inference on functions of the model parameters.

An approximation to the marginal posterior distribution of the penalty parameters $λ$ keeps playing an important role in the procedure. Point estimates for $λ$ can be derived from it with a subsequent empirical Bayes approach (Carlin and Louis, 2000) to handle these hyperparameters, see Section 3 for an illustration. Alternatively, the uncertainty in the selection of $λ$ could be accounted for by marginalizing over it with a Monte Carlo or a grid-based integration in $p (ξ | D) = \int_{λ} p (ξ | λ, D) p (λ | D) d λ .$ However, in the context of generalized additive models (Gressani and Lambert, 2021) and nonparametric double additive location-scale models (Lambert, 2021), simulation studies suggest that coverages of credible intervals resulting from an empirical Bayes approach for model parameter estimation are already close to their nominal values, even with moderate sample sizes.

The proposed methodology diverges from the proposal made by Rue et al. (2009) and underlying INLA where the size of the latent vector $ξ$ increases with sample size and where the marginal distribution of the scalar components of $ξ$ are the research focus. The joint distribution of $(ξ | λ, D)$ with asymmetric forms for the non-penalized components is here available for the whole vector $ξ$ and completed by an approximation to the marginal posterior distribution of the hyperparameters $λ$ .

The code necessary to reproduce the results in the article can be downloaded from https://github.com/plambertULiege/ordgam.

This article and, more broadly, our research on smoothing methods, owe much to Brian Marx, who left us too soon. His joint work with Paul Eilers will continue to endure and shape the field for years to come.

A Gradient and Hessian in the PO model

Consider the proportional odds model defined in Section 3 and the notations therein. Let $v_{i r} = F_{i r} (1 - F_{i r})$ , $w_{i r} = (1 + π_{i r} - 2 F_{i r})$ , $z_{i r} = (1 - 2 F_{i r}) v_{i r}$ for $r$ in $ε_{R}$ and take $s, t \in ε_{R - 1}$ . One has:

\frac{\partial l}{\partial γ_{s}} = \frac{\sum_{i} \partial \log π_{i y_{i}}}{\partial γ_{s}};_{\partial γ_{s}}^{\partial \log π_{i r}} =_{π_{i r}}^{1} (δ_{r s} v_{i r} - δ_{r - 1, s} v_{i, r - 1})

\frac{\partial l}{\partial θ_{k}} = \sum_{i} \frac{\partial \log π_{i y_{i}}}{\partial θ_{k}}; \frac{\partial \log π_{i r}}{\partial θ_{k}} = x_{i k} w_{i r}

and

\frac{\partial^{2} l}{\partial γ_{s} \partial γ_{t}} = \sum_{i} \frac{1}{π_{i y_{i}}} \{δ_{y_{i}, s, t} z_{i, y_{i}} - δ_{y_{i} - 1, s, t} z_{i, y_{i} - 1}\} - \sum_{i} \frac{\partial \log π_{i y_{i}}}{\partial γ_{s}} \frac{\partial \log π_{i y_{i}}}{\partial γ_{t}};

\frac{\partial^{2} l}{\partial^{2} θ_{k} θ_{l}} = \sum_{i} x_{i k} x_{i l} (π_{i y_{i}} w_{i y_{i}} - 2 \sum_{j = 1}^{y_{i}} π_{i j} w_{i j});

\frac{\partial^{2} l}{\partial θ_{k} \partial γ_{s}} = - \sum_{i} x_{i k} (δ_{y_{i} s} v_{i y_{i}} + δ_{y_{i} - 1, s} v_{i, y_{i} - 1}) .

Therefore, given $λ$ , one has

U_{λ} = \frac{\partial \log p (ξ | λ, D)}{\partial ξ} = \frac{\partial l}{\partial ξ} - K_{λ} (ξ - e); H_{λ} = \frac{\partial^{2} \log p (ξ | λ, D)}{\partial ξ \partial ξ^{⊤}} = \frac{\partial^{2} l}{\partial ξ \partial ξ^{⊤}} - K_{λ} .

Footnotes

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

Philippe Lambert acknowledges the support of the ARC project IMAL (grant 20/25-107) financed by the Wallonia-Brussels Federation and granted by the Académie Universitaire Louvain.

References

Agresti

(2010) Analysis of ordinal categorical data . Hoboken, New Jersey:

John Wiley & Sons, Inc.

Azzalini

and Capitanio

(2014) The Skew-Normal and Related Families . New-York:

Cambridge University Press.

Brezger

and Lang

(2006) Generalized structured additive regression based on Bayesian P-splines. Computational Statistics and Data Analysis , 50, 967–991.

Carlin

and Louis

(2000) Bayes and Empirical Bayes Methods for Data Analysis . New-York: Chapman and Hall / CRC Press.

Chiuchiolo

, van Niekerk

and Rue

(2022) Joint posterior inference for latent Gaussian models with R-INLA. Journal of Statistical Computation and Simulation . doi: 10.1080/00949655.2022.2117813.

Eilers

PHC

and Marx

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

Eilers

PHC

and Marx

(2021) Practical Smoothing: The Joys of P-splines . Cambridge University Press.

ESS Round 9 (2018). European Social Survey Round 9. Data file edition 3.1. Sikt - Norwegian Agency for shared services in education and research, Norway - Data Archive and distributor of ESS data for ESS ERIC. doi: 10.21338/NSD-ESS9-2018.

Ferkingstad

and Rue

(2015) Improving the INLA approach for approximate Bayesian inference for latent Gaussian models. Electronic Journal of Statistics , 9, 2706–2731.

10.

Gamerman

(1997) Sampling from the posterior distribution in generalized linear mixed models. Statistics and Computing , 7, 57–68.

11.

Gressani

and Lambert

(2018) Fast Bayesian inference using Laplace approximations in a flexible promotion time cure model based on P-splines. Computational Statistics & Data Analysis , 124, 151–167.

12.

Gressani

and Lambert

(2021) Laplace approximations for fast Bayesian inference in generalized additive models based on P-splines. Computational Statistics & Data Analysis , 154, 107088.

13.

Gressani

, Faes

and Hens

(2022a) An approximate Bayesian approach for estimation of the instantaneous reproduction number under misreported epidemic data. Biometrical Journal . doi: 10.1002/bimj.202200024.

14.

Gressani

, Faes

and Hens

(2022b) Laplacian-P-splines for Bayesian inference in the mixture cure model. Statistics in Medicine , 41, 2602–2626.

15.

Gressani

, Wallinga

, Althaus

, Hens

and Faes

(2022c) EpiLPS: A fast and flexible Bayesian tool for estimation of the time-varying reproduction number. PLoS Computational Biology , 18, 1010618.

16.

Jullion

and Lambert

(2007) Robust specification of the roughness penalty prior distribution in spatially adaptive Bayesian P-splines models. Computational Statistics and Data Analysis , 51, 2542–2558.

17.

Lambert

(2021) Fast Bayesian inference using Laplace approximations in nonparametric double additive location-scale models with right- and interval-censored data. Computational Statistics & Data Analysis , 161, 107250.

18.

Lambert

and Bremhorst

(2019) Estimation and identification issues in the promotion time cure model when the same covariates influence long- and short-term survival. Biometrical Journal , 61, 275–289.

19.

Lambert

and Eilers

PHC

(2009) Bayesian density estimation from grouped continuous data. Computational Statistics and Data Analysis , 53, 1388–1399.

20.

Lambert

and Kreyenfeld

(2023) Exogenous time-varying covariates in double additive cure survival model with application to fertility. arXiv:2302.00331.

21.

Lang

and Brezger

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

22.

Laplace

(1774) Mémoire sur la probabilité des causes par les événements. Mémoires de l’Académie Royale des Sciences de Paris (Savants étrangers) , 6, 621–656.

23.

Leonard

(1982) A simple predictive density function: Comment. Journal of the American Statistical Association , 77, 657–658.

24.

Marquardt

(1963) An algorithm for least-squares estimation of nonlinear parameters. Applied Mathematics , 11, 431–441.

25.

Marx

and Eilers

PHC

(1998) Direct generalized additive modeling with penalized likelihood. Computational Statistics and Data Analysis , 28, 193–209.

26.

Piegorsch

(1990) Maximum likelihood estimation for the negative binomial dispersion parameter. Biometrics , 46, xfp –867.

27.

Roberts

and Tweedie

(1996) Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli , 2, 341–363.

28.

Rue

and Held

(2005) Gaussian Markov Random Fields: Theory and Applications . Boca Raton: CRC Press.

29.

Rue

, Martino

and Chopin

(2009) Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology) , 71, 319–392.

30.

Ruli

, Sartori

and Ventura

(2016) Improved Laplace approximation for marginal likelihoods. Electronic Journal of Statistics , 10, 3986–4009.

31.

Tierney

and Kadane

(1986) Accurate approximations for posterior moments and marginal densities. Journal of the American Statistical Association , 81, 82–86.

32.

Tierney

, Kass

and Kadane

(1989) Fully exponential Laplace approximations to expectations and variances of nonpositive functions. Journal of the American Statistical Association , 84, 710–716.

33.

Van der Vaart

(1998) Asymptotic statistics . Cambridge University Press.

34.

Wood

(2017) Generalized Additive Models: An Introduction with R (2nd ed.). Boca Raton: CRC Press.

35.

Wood

and Fasiolo

(2017) A generalized Fellner-Schall method for smoothing parameter optimization with application to Tweedie location, scale and shape models. Biometrics , 73, 1071–1081.

Penalty parameter selection and asymmetry corrections to Laplace approximations in Bayesian P-splines models

Abstract

Keywords

1 Introduction

Figure 1

Left panel: Count data ( n = 120 ) generated using a negative binomial NB ( μ ( x ), ϕ ) with ϕ − 1 = 6 . Right panel: Histogram of a MCMC sample for log ϕ compared to Laplace (solid) and skew- t (dashed) approximations.

2.1 Model specification

3.1 The additive proportional odds model for ordinal data

Figure 2

ESS dataset: fitted additive terms for eduyrs and age with pointwise 95% credible intervals: this suggests a growing hostility to homosexuality beyond the age of 60, while it appears that the number of completed years of education does not play a statistically significant role.

ESS dataset: scatterplots of the MCMC samples for ( γ | λ , D ) and ( γ ˜ | λ , D ) when λ = λ ^ .

ESS dataset: approximated marginal posterior density for ( γ ˜ | λ , D ) compared to MCMC samples when λ = λ ^ .

ESS dataset: approximated marginal posterior density for ( γ | λ , D ) compared to MCMC samples when λ = λ ^ .

A Gradient and Hessian in the PO model

Footnotes

Declaration of Conflicting Interests

Funding

References

Left panel: Count data ( $n = 120$ ) generated using a negative binomial $NB (μ (x), ϕ)$ with $ϕ^{- 1} = 6$ . Right panel: Histogram of a MCMC sample for $\log ϕ$ compared to Laplace (solid) and skew- $t$ (dashed) approximations.

ESS dataset: scatterplots of the MCMC samples for $(γ | λ, D)$ and $(\tilde{γ} | λ, D)$ when $λ = \hat{λ}$ .

ESS dataset: approximated marginal posterior density for $(\tilde{γ} | λ, D)$ compared to MCMC samples when $λ = \hat{λ}$ .

ESS dataset: approximated marginal posterior density for $(γ | λ, D)$ compared to MCMC samples when $λ = \hat{λ}$ .