Full article: Smoothing Parameter and Model Selection for General Smooth Models

ABSTRACT
Formulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

ABSTRACT

This article discusses a general framework for smoothing parameter estimation for models with regular likelihoods constructed in terms of unknown smooth functions of covariates. Gaussian random effects and parametric terms may also be present. By construction the method is numerically stable and convergent, and enables smoothing parameter uncertainty to be quantified. The latter enables us to fix a well known problem with AIC for such models, thereby improving the range of model selection tools available. The smooth functions are represented by reduced rank spline like smoothers, with associated quadratic penalties measuring function smoothness. Model estimation is by penalized likelihood maximization, where the smoothing parameters controlling the extent of penalization are estimated by Laplace approximate marginal likelihood. The methods cover, for example, generalized additive models for nonexponential family responses (e.g., beta, ordered categorical, scaled t distribution, negative binomial and Tweedie distributions), generalized additive models for location scale and shape (e.g., two stage zero inflation models, and Gaussian location-scale models), Cox proportional hazards models and multivariate additive models. The framework reduces the implementation of new model classes to the coding of some standard derivatives of the log-likelihood. Supplementary materials for this article are available online.

KEYWORDS: Additive model, AIC, Distributional regression, GAM, Location scale and shape model, Ordered categorical regression, Penalized regression spline, REML, Smooth Cox model, Smoothing parameter uncertainty, Statistical algorithm, Tweedie distribution.

1. Introduction

This article is about smoothing parameter estimation and model selection in statistical models with a smooth regular likelihood, where the likelihood depends on smooth functions of covariates and these smooth functions are the targets of inference. Simple Gaussian random effects and parametric dependencies may also be present. When the likelihood (or a quasi-likelihood) decomposes into a sum of independent terms each contributed by a response variable from a single parameter exponential family distribution, then such a model is a generalized additive model (GAM, Hastie and Tibshirani 1986, 1990). GAMs are widely used in practice (see, e.g., Ruppert, Wand, and Carroll 2003; Fahrmeir et al. 2013) with their popularity resting in part on the availability of statistically well founded smoothing parameter estimation methods that are numerically efficient and robust (Wood 2000, 2011) and perform the important task of estimating how smooth the component functions of a model should be.

The purpose of this article is to provide a general method for smoothing parameter estimation when the model likelihood does not have the convenient exponential family (or quasi-likelihood) form. For the most part we have in mind regression models of some sort, but the proposed methods are not limited to this setting. The simplest examples of the extension are generalized additive models where the response distribution is not in the single parameter exponential family. For example, when the response has a Tweedie, negative binomial, beta, scaled t, or some sort of ordered categorical or zero inflated distribution. Examples of models with a less GAM like likelihood structure are Cox proportional hazard and Cox process models, scale-location models, such as the GAMLSS class of Rigby and Stasinopoulos (2005), and multivariate additive models (e.g., Yee and Wild 1996). Smooth function estimation for such models is not new: what is new here is the general approach to smoothing parameter estimation, and the wide variety of smooth model components that it admits.

The proposed method broadly follows the strategy of Wood (2011) that has proved successful for the GAM class. The smooth functions will be represented using reduced rank spline bases with associated smoothing penalties that are quadratic in the spline coefficients. There is now a substantial literature showing that the reduced rank approach is well-founded, and the basic issues are covered in an online Supplementary Appendix A (henceforth online “SA A”). More importantly, from an applied perspective, a wide range of spline and Gaussian process terms can be included as model components by adopting this approach (Figure 1). We propose to estimate smoothing parameters by Newton optimization of a Laplace approximate marginal likelihood criterion, with each Newton step requiring an inner Newton iteration to find maximum penalized likelihood estimates of the model coefficients. Implicit differentiation is used to obtain derivatives of the coefficients with respect to the smoothing parameters. This basic strategy works well in the GAM setting, but is substantially more complex when the simplifications of a GLM type likelihood no longer apply.

Figure 1. Examples of the rich variety of smooth model components that can be represented as reduced rank basis smoothers, with quadratic penalties and therefore can routinely be incorporated as components of a GAM. This article develops methods to allow their routine use in a much wider class of models. (a) One dimensional smooths such as cubic, P- and adaptive splines. (b) isotropic smooths of several variables, such as thin plate splines and Duchon splines. (c) Nonisotropic tensor product splines used to model smooth interactions. (d) Gaussian Markov random fields for data on discrete geographies. (e) Finite area smoothers, such as soap film smoothers. (f) Splines on the sphere. Another important class are simple Gaussian random effects.

Display full size

Our aim is to provide a general method that is as numerically efficient and robust as the GAM methods, such that (i) implementation of a model class requires only the coding of some standard derivatives of the log-likelihood for that class and (ii) much of the inferential machinery for working with such models can reuse GAM methods (e.g., interval estimation or p-value computations). An important consequence of our approach is that we are able to compute a simple correction to the conditional AIC for the models considered, which corrects for smoothing parameter estimation uncertainty and the consequent deficiencies in a conventionally computed conditional AIC (see Greven and Kneib 2010). This facilitates the part of model selection distinct from smoothing parameter estimation.

The article is structured as follows. Section 2 introduces the general modeling framework. Section 3 then covers smoothness selection methods for this framework, with Section 3.1 developing a general method, Section 3.2 illustrating its use for the special case of distributional regression, and Section 3.3 covering the simplified methods that can be used in the even more restricted case of models with a similar structure to generalized additive models. Section 4 then develops approximate distributional results accounting for smoothing parameter uncertainty which are applied in Section 5 to propose a corrected AIC suitable for the general model class. The remaining sections present simulation results and examples, while extensive further background, and details for particular models, are given in the supplementary appendices (referred to as online “SA A,” “SA B,” etc., below).

2. The General Framework

Consider a model for an n-vector of data, $y$ , constructed in terms of unknown parameters, $θ$ , and some unknown functions, g_j, of covariates, x_j. Suppose that the log-likelihood for this model satisfies the Fisher regularity conditions, has four continuous derivatives, and can be written $l (θ, g_{1}, g_{2}, ..., g_{M}) = log f (y | θ, g_{1}, g_{2}, ..., g_{M})$ . In contrast to the usual GAM case, the likelihood need not be based on a single parameter exponential family distribution, and we do not assume that the log-likelihood can be written in terms of a single additive linear predictor. Now let the g_j(x_j) be represented via basis expansions of modest rank (k_j), $g_{j} (x) = \sum_{i = 1}^{k_{j}} β_{j i} b_{j i} (x),$ where the β_ji are unknown coefficients and the b_ji(x) are known basis functions such as splines, usually chosen to have good approximation theoretical properties. With each g_j is associated a smoothing penalty, which is quadratic in the basis coefficients and measures the complexity of g_j. Writing all the basis coefficients and $θ$ in one p-vector $β$ , then the $j th$ smoothing penalty can be written as $β^{T} S^{j} β$ , where S^j is a matrix of known coefficients, but generally has only a small nonzero block. The estimated model coefficients are then (1) $\hat{β} = \underset{β}{argmax} \{l (β) - \frac{1}{2} \sum_{j}^{M} λ_{j} β^{T} S^{j} β\}$ (1) given M smoothing parameters, λ_j, controlling the extent of penalization. A slight extension is that the smoothing penalties may be such that several $λ_{i} β^{T} S^{i} β$ are associated with one g_j, for example when g_j is a nonisotropic function of several variables. Note also that the framework can incorporate Gaussian random effects, provided the corresponding precision matrices can be written as $\sum λ_{i} β^{T} S^{i} β$ (where the Sⁱ are known).

From a Bayesian viewpoint $\hat{β}$ is a posterior mode for $β$ . The Bayesian approach views the smooth functions as intrinsic Gaussian random fields with prior f_λ given by N(0, S^{λ −}) where S^{λ −} is a Moore–Penrose (or other suitable) pseudoinverse of ∑_jλ_jS^j. Then the posterior modes are $\hat{β}$ from (1), and in the large sample limit, assuming fixed smoothing parameter vector, $λ$ , we have $β | y \sim N (\hat{β}, {(I + S^{λ})}^{- 1})$ , where $I$ is the expected negative Hessian of the log-likelihood (or its observed version) at $\hat{β}$ . An empirical Bayesian approach is appealing here as it gives well calibrated inference for the g_j (Wahba 1983; Silverman 1985; Nychka 1988; Marra and Wood 2012) in a GAM context. Appropriate summations of the elements of $diag {{(I + S^{λ})}^{- 1} I}$ provide estimates of the “effective degrees of freedom” of the whole model, or of individual smooths.

Under this Bayesian view, smoothing parameters can be estimated to maximize the log marginal likelihood (2) $V_{r} (λ) = log \int f (y | β) f_{λ} (β) d β,$ (2) or a Laplace approximate version of this (e.g., Wood 2011). In practice optimization is with respect to $ρ$ where ρ_i = log λ_i. Marginal likelihood estimation of smoothing parameters in a Gaussian context goes back to Anderssen and Bloomfield (1974) and Wahba (1985), while Shun and McCullagh (1995) showed that Laplace approximation of more general likelihoods is theoretically well founded. That marginal likelihood is equivalent to REML (in the sense of Laird and Ware 1982) supports its use when the model contains Gaussian random effects. Theoretical work by Reiss and Ogden (2009) also suggests practical advantages at finite sample sizes, in that marginal likelihood is less prone to multiple local minima than GCV (or AIC). Supplementary Appendix B (SA B) also demonstrates how Laplace approximate marginal likelihood (LAML) estimation of smoothing parameters maintains statistical consistency of reduced rank spline estimates. The use of Laplace approximation and demonstration of statistical consistency requires the assumption that $dim (β) = O (n^{α})$ where α < 1/3.

3. Smoothness Selection Methods

This section describes the general smoothness selection method, and a simplified method for the special case in which the likelihood is a simple sum of terms for each observation of a univariate response, and there is a single GAM like linear predictor.

The nonlinear dependencies implied by employing a general smooth likelihood result in unwieldy expressions unless some care is taken to establish a compact notation. In the rest of this article, Greek subscripts denote partial differentiation with respect to the given variable, while Roman superscripts are indices associated with the derivatives. Hence, D^ij_βθ = ∂²D/∂β_i∂θ_j. Similarly $D_{β θ}^{i j} = \partial^{2} D / \partial β_{i} \partial θ_{j} |_{\hat{β}}$ . Roman subscripts denote vector or array element indices. For matrices the first Roman sub- or superscript denotes rows, the second columns. Roman superscripts without a corresponding Greek subscript are labels, for example $β^{1}$ and $β^{2}$ denote two separate vectors $β$ . For Hessian matrices only, D^β_i_j^θ is element i, j of the inverse of the matrix with elements Dⁱ_β^j_θ. If any Roman index appears in two or more multiplied terms, but the index is absent on the other side of the equation, then a summation over the product of the corresponding terms is indicated (the usual Einstein summation convention being somewhat unwieldy in this context). To aid readability, in this article summation indices will be highlighted in bold. For example, the equation $a_{i j} b_{i k} c^{i l} + d_{j k l} = 0$ is equivalent to ∑_ia_ijb_ikc^il + d_jkl = 0. An indexed expression not in an equation is treated like an equation with no indices on the other side (so $a_{i j} b_{j}$ is interpreted as ∑_ja_ijb_j).

3.1. General Model Estimation

Consider the general case in which the log-likelihood depends on several smooth functions of predictor variables, each represented via a basis expansion and each with one or more associated penalties. The likelihood may also depend on some strictly parametric model components. The log-likelihood is assumed to satisfy the Fisher regularity conditions and in addition we usually assume that it has 4 bounded continuous derivatives with respect to the parameters (with respect to g_j(x) for any relevant fixed x in the case of a smooth, g_j). Let the model coefficients be $β$ (recalling that this includes the vector $θ$ of parametric coefficients and nuisance parameters). The penalized log-likelihood is then $L (β) = l (β) - \frac{1}{2} λ_{j} β^{T} S^{j} β,$ and we assume that the model is well enough posed that this has a positive definite maximum (at least after dealing with any parameter redundancy issues that can be addressed by linear constraint). Let $\hat{β}$ be the maximizer of $L$ and let $H$ be the negative Hessian, with elements $- L_{\hat{β}}^{i}_{\hat{β}}^{j}$ . The log LAML (see online SA C) is $V (λ) = L (\hat{β}) + \frac{1}{2} log | S^{λ} |_{+} - \frac{1}{2} log | H | + \frac{M_{p}}{2} log (2 π),$ where $S^{λ} = λ_{j} S^{j}$ and |S^λ|₊ is the product of the positive eigenvalues of S^λ. M_p is the number of zero eigenvalues of S^λ, when all λ_j are strictly positive. The basic strategy is to optimize $V$ with respect to $ρ = log (λ)$ via Newton’s method. This requires $\hat{β}$ to be obtained for each trial $ρ$ via an inner Newton iteration, and derivatives of $\hat{β}$ must be obtained by implicit differentiation. The log determinant computations have the potential to be computationally unstable, and reparameterization is needed to deal with this. The full Newton method based on computationally exact derivatives has the substantial practical advantage that it can readily be detected when $V$ is indefinite with respect to a particular ρ_i, since then $\partial V / \partial ρ_{i} = \partial^{2} V / \partial ρ_{i}^{2} ≃ 0$ . Such indefiniteness occurs when a smoothing parameter, λ_i, →∞ or a variance component tends to zero, both of which are perfectly legitimate. Dropping a ρ_i from Newton update when such indefiniteness is detected ensures that it takes a value which can be treated as “working infinity” without overflowing. Methods which use an approximate Hessian, or none, do not have this advantage.

The proposed general method consists of outer and inner iterations, as follows.

Outer algorithm for $ρ$

Obtain initial values for $ρ = log (λ)$ , to ensure that the effective degrees of freedom of each smooth lies away from its maximum or minimum possible values.

Find initial $\hat{β}$ guesstimates (model specific).

Perform the initial reparameterizations required in Section 3.1.1 to facilitate stable computation of log |S^λ|₊.

Repeat the following standard Newton iteration until convergence is detected at Step (c).

	(a)	Find $\hat{β}$ , $V_{ρ}^{i}$ and $V_{ρ ρ}^{i j}$ by the inner algorithm.
	(b)	Drop any $V_{ρ}^{i}$ , $V_{ρ ρ}^{i j}$ and $V_{ρ ρ}^{j i}$ for which $V_{ρ}^{i} ≃ V_{ρ ρ}^{i i} ≃ 0$ . Let $I$ denote the indices of the retained terms.
	(c)	Test for convergence, that is, all $V_{ρ}^{i} ≃ 0$ and the Hessian (elements $- V_{ρ ρ}^{j i}$ ) is positive semidefinite.
	(d)	If necessary perturb the Hessian (elements $- V_{ρ ρ}^{j i}$ ) to make it positive definite (guaranteeing that the Newton step will be a descent direction).
	(e)	Define $Δ_{I}$ as the subvector of $Δ$ indexed by $I$ , with elements $- V_{i j}^{ρ ρ} V_{ρ}^{j}$ , and set Δ_j = 0 ∀ $j \notin I$ .
	(f)	While $V (ρ + Δ) < V (ρ)$ set $Δ \leftarrow Δ / 2$ .
	(g)	Set $ρ \leftarrow ρ + Δ$ .

Reverse the Step 3 reparameterization.

The method for evaluating $V$ and its gradient and Hessian with respect to $ρ$ is as follows, where $L_{k}^{\hat{β}}_{j}^{\hat{β}}$ denotes the inverse of $L_{\hat{β}}^{k}_{\hat{β}}^{j}$ .

Inner algorithm for $β$

1.	Reparameterize to deal with any “type 3” penalty blocks as described in Section 3.1.1 so that computation of log \|S^λ\|₊ is stable, and evaluate the derivatives of log \|S^λ\|₊.
2.	Use Newton’s method to find $\hat{β}$ , regularizing the Hessian, and applying step length control, to ensure convergence even when the Hessian is indefinite and/or $\hat{β}$ is not identifiable, as described in Section 3.1.2.
3.	Test for identifiability of $\hat{β}$ at convergence by examining the rank of the $H$ as described in Section 3.1.2. Drop unidentifiable coefficients.
4.	If coefficients were dropped, find the reduced $\hat{β}$ by further steps of Newton’s method (Section 3.1.2).
5.	Compute $d {\hat{β}}_{i} / d ρ_{k} = L_{i}^{\hat{β}}_{j}^{\hat{β}} λ_{k} S_{j l}^{k} {\hat{β}}_{l}$ and hence $l_{\hat{β}}^{i}_{\hat{β}}^{j}_{ρ}^{l} = l_{\hat{β}}^{i}_{\hat{β}}^{j}_{\hat{β}}^{k} d {\hat{β}}_{k} / d ρ_{l}$ (Section 3.1.3).
6.	Compute $d^{2} {\hat{β}}_{i} / d ρ_{k} d ρ_{l} = L_{i}^{\hat{β}}_{j}^{\hat{β}} {(- l * *_{\hat{β}}^{j} {\hat{β}}_{ρ}^{p} l + λ_{l} S_{j p}^{l}) d {\hat{β}}_{p} / d ρ_{k} + λ_{k} S_{j p}^{k} d {\hat{β}}_{p} / d ρ_{l}} + δ_{k}^{l} d {\hat{β}}_{i} / d ρ_{k},$ (Section 3.1.3).
7.	Compute $L * _{k}^{\hat{β}} {\hat{β}}_{j} l *_{\hat{β}}^{j} k_{\hat{β}} ρ_{ρ}^{p} v$ (model specific). (3.1.3)
8.	The derivatives of $V$ can now be computed according to Section 3.1.4.
9.	For each parameter dropped from $\hat{β}$ during fitting, zeroes must be inserted in $\hat{β}$ , $\partial \hat{β} / \partial ρ_{j}$ and the corresponding rows and columns of $L_{k}^{\hat{β}}_{j}^{\hat{β}}$ . The Step 1 reparameterization is then reversed.

The following subsections fill in the method details, but note that to implement a particular model in this class it is necessary to be able to compute, l, lⁱ_β and lⁱ_β^j_β, given $β$ , along with $l * *_{\hat{β}}^{i} j_{\hat{β}}^{k} ρ$ given $d \hat{β} / d ρ_{k}$ , and $L * *_{k}^{\hat{β}} {\hat{β}}_{j} l * *_{\hat{β}}^{j} k_{\hat{β}} ρ_{ρ}^{p} v$ given $d^{2} \hat{β} / d ρ_{k} d ρ_{l}$ . The last of these is usually computable much more efficiently than if $l * *_{\hat{β}}^{j} k_{\hat{β}} ρ_{ρ}^{p} v$ was computed explicitly.

3.1.1. Derivatives and Stable Evaluation of log |S^λ|₊

This section covers the details for outer Step 3 and inner Step 1. Stable evaluation of the log determinant terms is the key to stable computation with the LAML. The online SA C explains the issue. Wood (2011) proposed a solution which involves orthogonal transformation of the whole parameter vector $β$ , but in the general case the likelihood may depend on each smooth function separately and such a transformation is therefore untenable. It is necessary to develop a reparameterization strategy which does not combine coefficients from different smooths. This is possible if we recognize that S^λ is block diagonal, with different blocks relating to different smooths. For example, if $S^{j}$ denotes the nonzero sub-block of S^j, $S^{λ} = (\begin{matrix} λ_{1} S^{1} & . & . & . & . \\ . & λ_{2} S^{2} & . & . & . \\ . & . & λ_{j} S^{j} & . & . \\ . & . & . & . & . \\ . & . & . & . & . \end{matrix}) .$ That is, there are some blocks with single smoothing parameters, and others with a more complicated additive structure. There are usually also some zero blocks on the diagonal. The block structure means that the generalized determinant, its derivatives with respect to ρ_k = log λ_k and the matrix square root of S^λ can all be computed blockwise. So for the above example, $\begin{matrix} log | S^{λ} |_{+} & = & rank (S^{1}) log (λ_{1}) + log {| S^{1} |}_{+} + rank (S^{2}) log (λ_{2}) \\ + log | S^{2} |_{+} + log {| λ_{j} S^{j} |}_{+} + \dots \end{matrix}$ For any ρ_k relating to a single parameter block we have $\frac{\partial log | S^{λ} |_{+}}{\partial ρ_{k}} = rank (S^{k})$ and zero second derivatives. For multi-λ blocks there will generally be first and second derivatives to compute. There are no second derivatives “between-blocks.”

In general, there are three block types, each requiring different preprocessing.

1.	Single parameter diagonal blocks. A reparameterization can be used so that all nonzero elements are one, and the rank precomputed.
2.	Single parameter dense blocks. An orthogonal reparameterization, based on the eigenvectors of the symmetric eigen-decomposition of the block, can be used to make these blocks look like the previous type (by similarity transform). Again the rank is computed.
3.	Multi-λ blocks will require the reparameterization method of Wood (2011) appendix B to be applied for each new $ρ$ proposal, since the numerical problem that the reparameterization avoids is $ρ$ dependent (see online SA C). Initially, before the smoothing parameter selection iteration, it is necessary to reparameterize to separate the parameters corresponding to the block into penalized and unpenalized subvectors. This initial reparameterization can be based on the eigenvectors of the symmetric eigen decomposition of the “balanced” version of the block penalty matrix, $\sum_{j} S^{j} / {∥ S^{j} ∥}_{F}$ , where ‖ · ‖_F is the Frobenious norm. The balanced penalty is used for maximal numerical stability, and is usable because formally the spaces for the penalized and unpenalized components do not change with the smoothing parameters.

The reparameterizations from each block type are applied to the model, usually to the model matrices X^j of the individual smooth terms. The reparameterization information must be stored so that we can return to the original parameterization at the end.

After the one off initial reparameterization just described, then step one of the inner algorithm requires only that the reparameterization method of Wood (2011) Appendix B be applied to the parameters corresponding to type 3 blocks, for each new set of smoothing parameters.

3.1.2. Newton Iteration for $\hat{β}$

This section provides details for inner Steps 2–4. Newton iteration for $\hat{β}$ requires the gradient vector, $G$ , with elements $L_{β}^{i} = l_{β}^{i} - λ_{k} S_{i j}^{k} β_{j}$ and negative Hessian matrix $H$ with elements $- L_{β}^{i}_{β}^{j} = - l_{β}^{i}_{β}^{j} + λ_{k} S_{i j}^{k}$ (we will also use $H$ to denote the Hessian of the negative unpenalized log-likelihood with elements − lⁱ_β^j_β). In principle Newton iteration proceeds by repeatedly setting $β$ to $β + Δ$ , where $Δ = H^{- 1} G$ . In practice, Newton’s method is only guaranteed to converge to a maximum of $L$ , provided (i) that the Hessian is perturbed to be positive definite if it is not, guaranteeing that the Newton direction is an ascent direction, (ii) that step reduction is used to ensure that the step taken actually increases $L$ and (iii) that the computation of the step is numerically stable (see Nocedal and Wright 2006).

$L$ may be indefinite away from a maximum, but even near the maximum there are two basic impediments to stability and positive definiteness. First, some elements of $β$ may be unidentifiable. This issue will be dealt with by dropping parameters at convergence, as described shortly. The second issue is that some smoothing parameters may legitimately become very large during fitting, resulting in very large λ_jS^j components, poor scaling, poor conditioning and, hence, computational singularity. However, given the initial and Step 1 reparameterizations, such large elements can be dealt with by diagonal preconditioning of $H$ . That is define diagonal matrix $D$ such that $D_{i i} = {| H_{i i} |}^{- 1 / 2}$ , and preconditioned Hessian $H^{'} = D H D$ . Then $H^{- 1} = D H^{' - 1} D$ , with the right-hand side resulting in much better scaled computation. In the work reported here the pivoted Cholesky decomposition of the perturbed Hessian $R^{T} R = H^{'} + ϵ I$ is repeated with increasing ε, starting from zero, until positive definiteness is obtained. The Newton step is then computed as $Δ = D R^{- 1} R^{- T} D G .$ If the step to $β + Δ$ fails to increase the likelihood then $Δ$ is repeatedly halved until it does. Note that the perturbation of the Hessian does not change the converged state of a Newton algorithm (although varying the perturbation strength can change the algorithm convergence rate).

At convergence $H$ can at worst be positive semi- definite, but it is necessary to test for the possibility that some parameters are unidentifiable. The test should not depend on the particular values of the smoothing parameters. This can be achieved by constructing the balanced penalty S = ∑_jS^j/‖S^j‖_F (‖ · ‖_F is the Frobenius norm, but another norm could equally well be used), and then forming the pivoted Cholesky decomposition $P^{T} P = H / {∥ H ∥}_{F} + S / {∥ S ∥}_{F}$ . The rank of $P$ can then be estimated by making use of Cline et al. (1979). If this reveals rank deficiency of order q then the coefficients corresponding to the matrix rows and columns pivoted to the last q positions should be dropped from the analysis. The balanced penalty is used to avoid dropping parameters simply because some smoothing parameters are very large. Given the nonlinear setting it is necessary to repeat the Newton iteration to convergence with the reduced parameter set, in order that the remaining parameters adjust to the omission of those dropped.

3.1.3. Implicit Differentiation

This section provides the details for inner Steps 5–7. We obtain the derivatives of the identifiable elements of $\hat{β}$ with respect to $ρ$ . All computations here are in the reduced parameter space, if parameters were dropped. At the maximum penalized likelihood estimate we have $L_{\hat{β}}^{i} = l_{\hat{β}}^{i} - λ_{k} S_{i j}^{k} {\hat{β}}_{j} = 0$ and differentiating with respect to ρ_k = log λ_k yields $\begin{matrix} L_{\hat{β}}^{i}_{ρ}^{k} & = & l_{\hat{β}}^{i}_{\hat{β}}^{j} \frac{d {\hat{β}}_{j}}{d ρ_{k}} - λ_{k} S_{i j}^{k} {\hat{β}}_{j} - λ_{l} S_{i j}^{l} \frac{d {\hat{β}}_{j}}{d ρ_{k}} = 0 and rearranging, \\ \frac{d {\hat{β}}_{i}}{d ρ_{k}} & = & L_{i}^{\hat{β}}_{j}^{\hat{β}} λ_{k} S_{j l}^{k} {\hat{β}}_{l}, \end{matrix}$ given which we can compute $l * *_{\hat{β}}^{i} {\hat{β}}_{ρ}^{j} l = l * *_{\hat{β}}^{i} {\hat{β}}_{\hat{β}}^{j} k d {\hat{β}}_{k} / d ρ_{l}$ from the model specification. $- l * *_{\hat{β}}^{i} {\hat{β}}_{ρ}^{j} l + δ_{k}^{l} λ_{k} S_{i j}^{k}$ are the elements of $\partial H / \partial ρ_{l}$ , required in the next section (δ^l_k is 1 for l = k and 0 otherwise). Then $\frac{d^{2} {\hat{β}}_{i}}{d ρ_{k} d ρ_{l}} = L_{i}^{\hat{β}}_{j}^{\hat{β}} \{(- l_{\hat{β}}^{j}_{\hat{β}}^{p}_{ρ}^{l} + λ_{l} S_{j p}^{l}) \frac{d {\hat{β}}_{p}}{d ρ_{k}} + λ_{k} S_{j p}^{k} \frac{d {\hat{β}}_{p}}{d ρ_{l}}\} + δ_{k}^{l} \frac{d {\hat{β}}_{i}}{d ρ_{k}},$ which enables computations involving $\partial^{2} H / \partial ρ_{k} \partial ρ_{l}$ , with elements $- l * *_{\hat{β}}^{i} j_{\hat{β}} ρ_{ρ}^{k} l + δ_{k}^{l} λ_{k} S_{i j}^{k}$ , and $l * *_{\hat{β}}^{i} j_{\hat{β}} ρ_{ρ}^{k} l = l_{\hat{β} \hat{β} \hat{β} \hat{β}}^{i j r t} \frac{d {\hat{β}}_{r}}{d ρ_{k}} \frac{d {\hat{β}}_{t}}{d ρ_{l}} + L_{\hat{β} \hat{β} \hat{β}}^{i j r} \frac{d^{2} {\hat{β}}_{r}}{d ρ_{k} d ρ_{l}} .$ As mentioned in Section 3.1, it will generally be inefficient to form this last quantity explicitly, as it occurs only in the summations involved in computing the final trace in (3).

3.1.4. The Remaining Derivatives

Recalling that $H$ is the matrix with elements $- L_{β}^{i}_{β}^{j} = - l_{β}^{i}_{β}^{j} + λ_{k} S_{i j}^{k}$ , we require (inner Step 8) $\frac{\partial V}{\partial ρ_{k}} = - \frac{λ_{k}}{2} {\hat{β}}^{T} S^{k} \hat{β} + \frac{1}{2} \frac{\partial log | S^{λ} |_{+}}{\partial ρ_{k}} - \frac{1}{2} \frac{\partial log | H |}{\partial ρ_{k}}$ and $\begin{matrix} \frac{\partial^{2} V}{\partial ρ_{k} \partial ρ_{l}} & = & - δ_{k}^{l} \frac{λ_{k}}{2} {\hat{β}}^{T} S^{k} \hat{β} - \frac{d {\hat{β}}^{T}}{d ρ_{l}} H \frac{d \hat{β}}{d ρ_{k}} + \frac{1}{2} \frac{\partial^{2} log {| S^{λ} |}_{+}}{\partial ρ_{k} \partial ρ_{l}} \\ - \frac{1}{2} \frac{\partial^{2} log | H |}{\partial ρ_{k} \partial ρ_{l}}, \end{matrix}$ where components involving $L_{\hat{β}}^{j}$ are zero by definition of $\hat{β}$ . The components not covered so far are (3) $\begin{matrix} \frac{\partial log | H |}{\partial ρ_{k}} & = & tr (H^{- 1} \frac{\partial H}{\partial ρ_{k}}) and \\ \frac{\partial^{2} log | H |}{\partial ρ_{k} \partial ρ_{l}} & = & - tr (H^{- 1} \frac{\partial H}{\partial ρ_{k}} H^{- 1} \frac{\partial H}{\partial ρ_{l}}) + tr (H^{- 1} \frac{\partial^{2} H}{\partial ρ_{k} \partial ρ_{l}}) . \end{matrix}$ (3) The final term above is expensive if computed naively by explicitly computing each term $\partial^{2} H / \partial ρ_{k} \partial ρ_{l}$ , but this is unnecessary and the computation of $tr (H^{- 1} \partial^{2} H / \partial ρ_{k} \partial ρ_{l})$ can usually be performed efficiently as the final part of the model specification, keeping the total cost to O(Mnp²): see online SA G and Section 3.2 for illustrative examples.

The Cox (1972) proportional hazards model provides a straightforward application of the general method, and the requisite computations are set out in online SA G in a manner that maintains O(Mnp²) computational cost. Another example is the multivariate additive model, in which the means of a multivariate Gaussian response are given by separate linear predictors, which may optionally share terms. This model is covered in the online SA H and Section 8. Section 3.2 considers how another class of models falls into the general framework.

3.2. A Special Case: GAMLSS Models

The GAMLSS (or “distributional regression”) models discussed by Rigby and Stasinopoulos (2005) (and also Yee and Wild 1996; Klein et al. 2014, 2015) fall within the scope of the general method. The idea is that we have independent univariate response observations, y_i, whose distributions depend on several unknown parameters, each of which is determined by its own linear predictor. The log-likelihood is a straightforward sum of contributions from each y_i (unlike the Cox models, e.g.), and the special structure can be exploited so that implementation of new models in this class requires only the supply of some derivatives of the log-likelihood terms with respect to the distribution parameters. Given the notational conventions established previously, the expressions facilitating this are rather compact (without such a notation they can easily become intractably complex).

Let the log-likelihood for the ith observation be l(y_i, η¹_i, η²_i, …) where the $η^{k} = X^{k} β^{k}$ are K linear predictors. The Newton iteration for estimating $β = {(β^{1 T}, β^{2 T}, ...)}^{T}$ requires $l_{β^{l}}^{j} = l_{η^{l}}^{i} X_{i j}^{l}$ and $l_{β^{l}}^{j}_{β^{m}}^{k} = l_{η^{l}}^{i}_{η^{m}}^{i} X_{i j}^{l} X_{i k}^{m},$ which are also sufficient for first-order implicit differentiation.

LAML optimization also requires $\begin{matrix} l_{{\hat{β}}^{l}}^{j}_{{\hat{β}}^{m}}^{k}_{ρ}^{p} & = & l_{{\hat{β}}^{l}}^{j}_{{\hat{β}}^{m}}^{k}_{{\hat{β}}^{q}}^{r} \frac{d {\hat{β}}_{r}^{q}}{d ρ_{p}} = l_{{\hat{η}}^{l}}^{i}_{{\hat{η}}^{m}}^{i}_{{\hat{η}}^{q}}^{i} X_{i j}^{l} X_{i k}^{m} X_{i r}^{q} \frac{d {\hat{β}}_{r}^{q}}{d ρ_{p}} \\ = & l_{{\hat{η}}^{l}}^{i}_{{\hat{η}}^{m}}^{i}_{{\hat{η}}^{q}}^{i} X_{i j}^{l} X_{i k}^{m} \frac{d {\hat{η}}_{i}^{q}}{d ρ_{p}} . \end{matrix}$ Notice how this is just an inner product $X^{T} V X$ , where the diagonal matrix V is the sum over q of some diagonal matrices. At this stage the second derivatives of $\hat{β}$ with respect to $ρ$ can be computed, after which we require only $\begin{matrix} l_{{\hat{β}}^{l}}^{j}_{{\hat{β}}^{m}}^{k}_{ρ}^{p}_{ρ}^{v} & = & l_{{\hat{β}}^{l}}^{j}_{{\hat{β}}^{m}}^{k}_{{\hat{β}}^{q}}^{r}_{{\hat{β}}^{s}}^{t} \frac{d {\hat{β}}_{r}^{q}}{d ρ_{p}} \frac{d {\hat{β}}_{t}^{s}}{d ρ_{v}} + l_{{\hat{β}}^{l}}^{j}_{{\hat{β}}^{m}}^{k}_{{\hat{β}}^{q}}^{r} \frac{d^{2} {\hat{β}}_{r}^{q}}{d ρ_{p} d ρ_{v}} \\ = & l_{{\hat{η}}^{l}}^{i}_{{\hat{η}}^{m}}^{i}_{{\hat{η}}^{q}}^{i}_{{\hat{η}}^{s}}^{i} X_{i j}^{l} X_{i k}^{m} \frac{d {\hat{η}}_{i}^{q}}{d ρ_{p}} \frac{d {\hat{η}}_{i}^{s}}{d ρ_{v}} + l_{{\hat{η}}^{l}}^{i}_{{\hat{η}}^{m}}^{i}_{{\hat{η}}^{q}}^{i} X_{i j}^{l} X_{i k}^{m} \frac{d^{2} {\hat{η}}_{i}^{q}}{d ρ_{p} d ρ_{v}} . \end{matrix}$ So to implement a new family for GAMLSS estimation requires mixed derivatives up to fourth order with respect to the parameters of the likelihood. In most cases what would be conveniently available is, for example, $l_{{\hat{μ}}^{l}}^{i}_{{\hat{μ}}^{m}}^{i}_{{\hat{μ}}^{q}}^{i}_{{\hat{μ}}^{s}}^{i}$ rather than $l_{{\hat{η}}^{l}}^{i}_{{\hat{η}}^{m}}^{i}_{{\hat{η}}^{q}}^{i}_{{\hat{η}}^{s}}^{i}$ , where μ^k is the $k th$ parameter of the likelihood and is given by h^k(μ^k) = η^k, h^k being a link function.

To get from the μ derivatives to the η derivatives, the rules (A.1)--(A.4) from Appendix A are used. This is straightforward for any derivative that is not mixed. For mixed derivatives containing at least one first-order derivative the transformation rule applying to the highest order derivative is applied first, followed by the transformations for the first-order derivatives. This leaves only the transformation of $l_{{\hat{μ}}^{j}}^{i}_{{\hat{μ}}^{j}}^{i}_{{\hat{μ}}^{k}}^{i}_{{\hat{μ}}^{k}}^{i}$ as at all awkward, but we have $\begin{matrix} l * *_{{\hat{η}}^{j}}^{i} i_{{\hat{η}}^{j}}^{i} {\hat{η}}^{k}_{{\hat{η}}^{k}}^{i} & = & (l * *_{{\hat{μ}}^{j}}^{i} i_{{\hat{μ}}^{j}}^{i} {\hat{μ}}^{k}_{{\hat{μ}}^{k}}^{i} / h_{i}^{j' 2} - l * *_{{\hat{μ}}^{j}}^{i} i_{{\hat{μ}}^{k}}^{i} {\hat{μ}}^{k} h_{i}^{j''} / h_{i}^{j' 3}) / \\ h_{i}^{k' 2} - (l * *_{{\hat{μ}}^{j}}^{i} i_{{\hat{μ}}^{j}}^{i} {\hat{μ}}^{k} / h_{i}^{j' 2} - l * *_{{\hat{μ}}^{j}}^{i} i_{{\hat{μ}}^{k}} h_{i}^{j''} / h_{i}^{j' 3}) h_{i}^{k''} / h_{i}^{k' 3} . \end{matrix}$ ^{Figure 2}

The general method requires $L_{k j}^{\hat{β} \hat{β}} l_{\hat{β} \hat{β} ρ ρ}^{j k p v}$ to be computed, which would have $O {M (M + 1) n P^{2} / 2}$ cost if the terms $l_{\hat{β} \hat{β} ρ ρ}^{j k p v}$ were computed explicitly for this purpose (where $P$ is the dimension of combined $β$ ). However, this can be reduced to $O (n P^{2})$ using a trick most easily explained by switching to a matrix representation. For simplicity of presentation assume K = 2, and define matrix B to be the inverse of the penalized Hessian, so that $B_{i j} = L_{i}^{\hat{β}}_{j}^{\hat{β}}$ . Defining $v_{i}^{l m} = l_{{\hat{η}}^{l}}^{i}_{{\hat{η}}^{m}}^{i}_{{\hat{η}}^{q}}^{i}_{{\hat{η}}^{s}}^{i} \frac{d {\hat{η}}_{i}^{q}}{d ρ_{p}} \frac{d {\hat{η}}_{i}^{s}}{d ρ_{v}} + l_{{\hat{η}}^{l}}^{i}_{{\hat{η}}^{m}}^{i}_{{\hat{η}}^{q}}^{i} \frac{d^{2} {\hat{η}}_{i}^{q}}{d ρ_{p} d ρ_{v}} and [4 p t] V^{l m} = diag (v_{i}^{l m}) we have [4 p t]$ (4) $L_{k j}^{\hat{β} \hat{β}} l_{\hat{β} \hat{β} ρ ρ}^{j k p v} = tr \{B (\begin{matrix} X^{1 T} V^{11} X^{1} X^{1 T} V^{12} X^{2} \\ X^{2 T} V^{12} X^{1} X^{2 T} V^{22} X^{2} \end{matrix})\} [4 p t] = tr \{B {(\begin{matrix} X^{1} 0 \\ 0 X^{2} \end{matrix})}^{T} (\begin{matrix} V^{11} X^{1} V^{12} X^{2} \\ V^{12} X^{1} V^{22} X^{2} \end{matrix})\} .$ (4) Hence, following the one off formation of $B {(\begin{matrix} X^{1} & 0 \\ 0 & X^{2} \end{matrix})}^{T}$ (which need only have $O (n P^{2})$ cost), each trace computation has $O (M n P)$ cost (since $tr (C^{T} D) = D_{i j} C_{i j}$ ).

See online SA I where a zero inflated Poisson model provides an example of the details. Figure 2 shows estimates for the model ${accel}_{i} \sim N (f_{1} (t_{i}), σ_{i}^{2})$ where log σ_i = f₂(t_i), f₁ is an adaptive P-spline and f₂ a cubic regression spline, while SA F.2 provides another application. Package mgcv also includes multinomial logistic regression implemented this way and further examples are under development. An interesting possibility with any model which has multiple linear predictors is that one or more of those predictors should depend on some of the same terms, and online SA H shows how this can be handled.

Figure 2. A smooth Gaussian location scale model fit to the motorcycle data from Silverman (1985), using the methods developed in Section 3.2. The left plot shows the raw data as open circles and an adaptive p-spline smoother for the mean overlaid. The right plot shows the simultaneous estimate of the standard deviation in the acceleration measurements, with the absolute values of the residuals as circles. Dotted curves are approximate 95% confidence intervals. The effective degrees of freedom of the smooths are 12.5 and 7.3 respectively.

Display full size

3.3. A More Special Case: Extended Generalized Additive Models

For models with a single linear predictor, in which the log-likelihood is a sum of contributions per y_i, it is possible to perform fitting by iterative weighted least squares, enabling profitable reuse of some components of standard GAM fitting methods, including the exploitation of very stable orthogonal methods for solving least squares problems. Specifically, consider observations y_i, and let the corresponding log-likelihood be of the form $l = \sum_{i} l_{i} (y_{i}, μ_{i}, θ, ϕ),$ where the terms in the summation may also be written as l_i for short, and μ_i is often $E (y_{i})$ , but may also be a latent variable (as in the ordered categorical model of SA K). Given h, a known link function, h(μ_i) = η_i where $η = X β + o$ , $X$ is a model matrix, $β$ is a parameter vector and $o$ is an offset (often simply 0). $θ$ is a parameter vector, containing the extra parameters of the likelihood, such as the p parameter of a Tweedie density (see online SA J), or the cut points of an ordered categorical model (see online SA K). Notice that in this case $θ$ is not treated as part of $β$ , since $θ$ can not always be estimated by straightforward iterative regression. Instead $θ$ will be estimated alongside the smoothing parameters. φ is a scale parameter, often fixed at one. Let ${\tilde{l}}_{i} = {max}_{μ_{i}} l_{i} (y_{i}, μ_{i}, θ, ϕ)$ denote the saturated log-likelihood. Define the deviance corresponding to y_i as $D_{i} = 2 ({\tilde{l}}_{i} - l_{i}) ϕ$ , where φ is the scale parameter on which D_i does not depend. Working in terms of the deviance is convenient in a regression setting, where deviance residuals are a preferred method for model checking and the proportion deviance explained is a natural substitute for the r² statistic as a measure of goodness of fit (but see the final comment in online SA I).

In general the estimates of $β$ will depend on some log smoothing parameter ρ_j = log λ_j, and it is notationally expedient to consider these to be part of the vector $θ$ , although it is to be understood that l does not actually depend on these elements of $θ$ . Given $θ$ , estimation of $β$ is by minimization of the penalized deviance $D (β, θ) = \sum_{i} D_{i} (β, θ) + \sum_{j} λ_{j} β^{T} S^{j} β,$ with respect to $β$ . This can be achieved by penalized iteratively reweighted least squares (PIRLS), which consists of iterative minimization of $\sum_{i} w_{i} {(z_{i} - X_{i} β)}^{2} + \sum_{j} λ_{j} β^{T} S^{j} β,$ where the pseudodata and weights are given by $z_{i} = η_{i} - o_{i} - \frac{1}{2 w_{i}} \frac{\partial D_{i}}{\partial η_{i}}, w_{i} = \frac{1}{2} \frac{d^{2} D_{i}}{d η_{i}^{2}} .$ Note that if w_i = 0 (or w_i is too close to 0), the penalized least squares estimate can be computed using only w_iz_i, which is then well defined and finite when z_i is not.

Estimation of $θ$ , and possibly φ, is by LAML. Writing $W$ as the diagonal matrix of w_i values, the log LAML is given by $\begin{matrix} V (θ, ϕ) & = & - \frac{D (\hat{β}, θ)}{2 ϕ} + \tilde{l} (θ, ϕ) - \frac{log | X^{T} WX + S^{λ} | - log | S^{λ} |_{+}}{2} \\ + \frac{M_{p}}{2} log (2 π ϕ), \end{matrix}$ where $W$ is evaluated at the $\hat{β}$ implied by $θ$ . To compute the derivatives of $V$ with respect to $θ$ the derivatives of $\hat{β}$ with respect to $θ$ are required. Note that $V$ is a full Laplace approximation, rather than the “approximate” Laplace approximation used to justify PQL (Breslow and Clayton 1993), so that PQL’s well known problems with binary and low count data are much reduced. In particular: (i) most PQL implementations estimate φ when fitting the working linear mixed model, even in the binomial and Poisson cases, where it is fixed at 1. For binary and low count data this can give very poor results. (ii) PQL uses the expected Hessian rather than the Hessian, and these only coincide for the canonical link case. (iii) PQL is justified by an assumption that the iterative fitting weights only vary slowly with the smoothing parameters, an assumption that is not needed here.

The parameters $θ$ and φ can be estimated by maximizing $V$ using Newton’s method, or a quasi-Newton method. Notice that $V$ depends directly on the elements of $θ$ via $D$ , $\tilde{l}$ and S^λ, but also indirectly via the dependence of $\hat{μ}$ and $W$ on $\hat{β}$ and hence on $θ$ . Hence, each trial $θ, ϕ$ requires a PIRLS iteration to find the corresponding $\hat{β}$ , followed by implicit differentiation to find the derivatives of $\hat{β}$ with respect to $θ$ . Once these are obtained, the chain rule can be applied to find the derivatives of $V$ with respect to $θ$ and φ.

As illustrated in SA C, there is scope for serious numerical instability in the evaluation of the determinant terms in $V$ , but for this case we can reuse the stabilization strategy from Wood (2011), namely for each trial $θ$ and φ:

1.	Use the orthogonal reparameterization from Appendix B of Wood (2011) to ensure that log \|S^λ\|₊ can be computed in a stable manner.
2.	Estimate $\hat{β}$ by PIRLS using the stable least squares method for negatively weighted problems from Section 3.3 of Wood (2011), setting structurally unidentifiable coefficients to zero.
3.	Using implicit differentiation, obtain the derivatives of $V$ required for a Newton update.

Step 3 is substantially more complicated than in Wood (2011), and is covered in Appendix A.

3.3.1. Extended GAM New Model Implementation

The general formulation above assumes that various standard information is available for each distribution and link. What is needed depends on whether quasi-Newton or full Newton is used to find $\hat{θ}$ . Here is a summary of what is needed for each distribution

1.	For finding $\hat{β}$ . Dⁱ_μ, Dⁱ_μⁱ_μ, h′, and h′′.
2.	For $\hat{ρ}$ via quasi-Newton. h′^′′, Dⁱ_μ^j_θ, Dⁱ_θ, Dⁱ_μⁱ_μⁱ_μ, and Dⁱ_μⁱ_μ^j_θ.
3.	For $\hat{ρ}$ via full Newton. h′^′′′, Dⁱ_θ^j_θ, Dⁱ_μ^j_θ^k_θ, Dⁱ_μⁱ_μⁱ_μⁱ_μ, Dⁱ_μⁱ_μⁱ_μ^j_θ, and Dⁱ_μⁱ_μ^j_θ^k_θ.

In addition, first and second derivatives of $\tilde{l}$ with respect to its arguments are needed. All of these quantities can be obtained automatically using a computer algebra package. $E D_{μ}^{i}_{μ}^{i}$ is also useful for further inference. If it is not readily computed then we can substitute Dⁱ_μⁱ_μ, but a complication of penalized modeling is that Dⁱ_μⁱ_μ can fail to be positive definite at $\hat{β}$ . When this happens $E D_{μ}^{i}_{μ}^{i}$ can be estimated as the nearest positive definite matrix to Dⁱ_μⁱ_μ.

We have implemented beta, negative binomial, scaled t models for heavy tailed data, simple zero inflated Poisson, ordered categorical and Tweedie additive models in this way. The first three were essentially automatic: the derivatives were computed by a symbolic algebra package and coded from the results. Some care is required in doing this, to avoid excessive cancellation error, underflow or overflow in the computations. Overly naive coding of derivatives can often lead to numerical problems: The online SA I on the zero inflated Poisson provides an example of the sort of issues that can be encountered. The ordered categorical and Tweedie models are slightly more complicated and details are therefore provided in the online SA J and K (including further examples of the need to avoid cancellation error).

4. Smoothing Parameter Uncertainty

Conventionally in a GAM context smoothing parameters have been treated as fixed when computing interval estimates for functions, or for other inferential tasks. In reality smoothing parameters must be estimated, and the uncertainty associated with this has generally been ignored except in fully Bayesian simulation approaches. Kass and Steffey (1989) proposed a simple first-order correction for this sort of uncertainty in the context of iid Gaussian random effects in a one way ANOVA type design. Some extra work is required to understand how their method works when applied to smooths. It turns out that the estimation methods described above provide the quantities required to correct for smoothing parameter uncertainty.

Assume we have several smooth model components, let ρ_i = log λ_i and S^λ = ∑_jλ_jS^j. Writing ${\hat{β}}_{ρ}$ for $\hat{β}$ , to emphasize the dependence of $\hat{β}$ on the smoothing parameters, we use the Bayesian large sample approximation (see SB.4) (5) $β | y, ρ \sim N ({\hat{β}}_{ρ}, V_{β}) where V_{β} = {(\hat{I} + S^{λ})}^{- 1}$ (5) which is exact in the Gaussian case, along with the large sample approximation (6) $ρ | y \sim N (\hat{ρ}, V_{ρ}),$ (6) where V_ρ is the inverse of the Hessian of the negative log marginal likelihood with respect to $ρ$ . Since the approximation (6) applies in the interior of the parameter space, it is necessary to substitute a Moore-Penrose pseudoinverse of the Hessian if a smoothing parameter is effectively infinite, or otherwise to regularize the inversion (which is equivalent to placing a Gaussian prior on $ρ$ ). Conventionally (5) is used with $\hat{ρ}$ plugged in and the uncertainty in $ρ$ neglected. To improve on this note that if (5) and (6) are correct, while $z \sim N (0, I)$ and independently $ρ^{*} \sim N (\hat{ρ}, V_{ρ})$ , then $β | y \overset{d}{=} {\hat{β}}_{ρ^{*}} + R_{ρ^{*}}^{T} z$ where $R_{ρ^{*}}^{T} R_{ρ^{*}} = V_{β}$ (and V_β depends on $ρ^{*}$ ). This provides a way of simulating from $β | y$ , but it is computationally expensive as ${\hat{β}}_{ρ^{*}}$ and $R_{ρ^{*}}$ must be computed afresh for each sample. (The conventional approximation would simply set $ρ^{*} = \hat{ρ}$ .) Alternatively consider a first-order Taylor expansion $β | y \overset{d}{=} {\hat{β}}_{\hat{ρ}} + J (ρ - \hat{ρ}) + R_{\hat{ρ}}^{T} z + \sum_{k} {\frac{\partial R_{ρ}^{T} z}{\partial ρ_{k}}|}_{\hat{ρ}} (ρ_{k} - {\hat{ρ}}_{k}) + r,$ where r is a lower order remainder term and $J = d \hat{β} {/ d ρ |}_{\hat{ρ}}$ . Dropping r, the expectation of the right-hand side is ${\hat{β}}_{\hat{ρ}}$ . Denoting the elements of R_ρ by R_ij, tedious but routine calculation shows that the three remaining random terms are uncorrelated with covariance matrix (7) $\begin{matrix} V_{β}^{'} & = & V_{β} + V^{'} + V^{''}, where V^{'} = {JV}_{ρ} J^{T} and \\ V_{j m}^{''} & = & \sum_{i}^{p} \sum_{l}^{M} \sum_{k}^{M} \frac{\partial R_{i j}}{\partial ρ_{k}} V_{ρ, k l} \frac{\partial R_{i m}}{\partial ρ_{l}}, \end{matrix}$ (7) which is computable at O(Mp³) cost (see online SA D). Dropping V′′ we have the Kass and Steffey (1989) approximation $β | y \sim N ({\hat{β}}_{\hat{ρ}}, V_{β}^{*}) where V_{β}^{*} = V_{β} + J V_{ρ} J^{T} .$ (A first-order Taylor expansion of $\hat{β}$ about $ρ$ yields a similar correction for the frequentist covariance matrix of $\hat{β}$ : $V_{\hat{β}}^{*} = {(\hat{I} + S^{λ})}^{- 1} \hat{I} {(\hat{I} + S^{λ})}^{- 1} + J V_{ρ} J^{T}$ , where $\hat{I}$ is the negative Hessian of the log-likelihood).

The online SA D shows that in a Demmler-Reinsch like parameterization, for any penalized parameter β_i with posterior standard deviation $σ_{β_{i}}$ , $\frac{d \hat{β_{i}} / d ρ_{j}}{d {(R^{T} z)}_{i} / d ρ_{j}} ≃ \frac{{\hat{β}}_{i}}{z_{i} σ_{β_{i}}} .$ So the $J (ρ - \hat{ρ})$ correction is dominant for components that are strongly nonzero. This offers some justification for using the Kass and Steffey (1989) approximation, but not in a model selection context, where near zero model components are those of most interest: hence, in what follows we will use (7) without dropping V′′.

5. An Information Criterion for Smooth Model Selection

When viewing smoothing from a Bayesian perspective, the smooths have improper priors (or alternatively vague priors of convenience) corresponding to the null space of the smoothing penalties. This invalidates model selection via marginal likelihood comparison. An alternative is a frequentist AIC (Akaike 1973), based on the conditional likelihood of the model coefficients, rather than the marginal likelihood. In the exponential family GAM context, Hastie and Tibshirani (1990, §6.8.3) proposed a widely used version of this conditional AIC in which the effective degrees of freedom of the model, τ₀, is used in place of the number of model parameters (in the general setting $τ_{0} = tr {V_{β} \hat{I}}$ is equivalent to the Hastie and Tibshirani (1990) proposal). But Greven and Kneib (2010) showed that this is overly likely to select complex models, especially when the model contains random effects: the difficulty arises because τ₀ neglects the fact that the smoothing parameters have been estimated and are, therefore, uncertain (a marginal AIC based on the frequentist marginal likelihood, in which unpenalized effects are not integrated out, is equally problematic, partly because of underestimation of variance components and consequent bias toward simple models). A heuristic alternative is to use $τ_{1} = tr (2 \hat{I} V_{β} - \hat{I} V_{β} \hat{I} V_{β})$ as the effective degrees of freedom, motivated by considering the number of unpenalized parameters required to optimally approximate a bias corrected version of the model, but the resulting AIC is too conservative (see, Section 6, e.g.). Greven and Kneib (2010) show how to exactly compute an effective modified AIC for the Gaussian additive model case based on defining the effective degrees of freedom as $\sum_{i} \partial {\hat{y}}_{i} / \partial y_{i}$ (as proposed by Liang et al. 2008). Yu and Yau (2012) and Säfken et al. (2014) considered extensions to generalized linear mixed models. The novel contribution of this section is to use the results of the previous section to avoid the problematic neglect of smoothing parameter uncertainty in the conditional AIC computation in a manner that is easily computed and applicable to the general model class considered in this article.

The derivation of AIC (see, e.g., Davison 2003, sec. 4.7) with the MLE replaced by the penalized MLE is identical up to the point at which the AIC score is represented as (8) $AIC = - 2 l (\hat{β}) + 2 E \{{(\hat{β} - β_{d})}^{T} I_{d} (\hat{β} - β_{d})\} [4 p t]$ (8) (9) $= - 2 l (\hat{β}) + 2 tr [E {(\hat{β} - β_{d}) {(\hat{β} - β_{d})}^{T}} I_{d}],$ (9) where $β_{d}$ is the coefficient vector minimizing the KL divergence and $I_{d}$ is the corresponding expected negative Hessian of the log-likelihood. In an unpenalized setting $E {(\hat{β} - β_{d}) {(\hat{β} - β_{d})}^{T}}$ is estimated as the observed inverse information matrix ${\hat{I}}^{- 1}$ and $τ^{'} = tr {E (\hat{β} - β_{d}) {(\hat{β} - β_{d})}^{T} I_{d}}$ is estimated as $tr ({\hat{I}}^{- 1} \hat{I}) = k$ . Penalization means that the expected inverse covariance matrix of $\hat{β}$ is no longer well approximated by $\hat{I}$ , and there are then two ways of proceeding.

The first is to view $β$ as a frequentist random effect, with predicted values $\hat{β}$ . In that case the covariance matrix for the predictions, $\hat{β}$ , corresponds to the posterior covariance matrix obtained when taking the Bayesian view of the smoothing process, so we have the conventional estimate $τ = tr {V_{β} \hat{I}}$ if we neglect smoothing parameter uncertainty, or $τ = tr (V_{β}^{'} \hat{I})$ accounting for it using (7).^{Figure 3}

The frequentist random effects formulation is not a completely natural way to view smooths, since we do not usually expect the smooth components of a model to be resampled from the prior with each replication of the data. However in the smoothing context V_β has the interpretation of being the frequentist covariance matrix for $\hat{β}$ plus an estimate of the prior expectation of the squared smoothing bias (matrix), which offers some justification for using the same τ estimate as in the strict random effects case. To see this consider the decomposition $E {(\hat{β} - β_{d}) {(\hat{β} - β_{d})}^{T}} = E {(\hat{β} - E \hat{β}) {(\hat{β} - E \hat{β})}^{T}} + Δ_{β} Δ_{β}^{T},$ where $Δ_{β}$ is the smoothing bias in $\hat{β}$ . The first term on the right-hand side, above, can be replaced by the standard frequentist estimate $V_{\hat{β}} = {(\hat{I} + S^{λ})}^{- 1} \hat{I} {(\hat{I} + S^{λ})}^{- 1}$ . Now expand the penalized log-likelihood around $β_{d}$ : $\begin{matrix} l_{p} (β^{'}) & ≃ & l (β_{d}) + \frac{\partial l}{\partial β^{T}} (β^{'} - β_{d}) - \frac{1}{2} {(β^{'} - β_{d})}^{T} I_{d} (β^{'} - β_{d}) \\ - \frac{1}{2} β^{' T} S^{λ} β^{'} . \end{matrix}$ Differentiating with respect to $β^{'}$ and equating to zero we obtain the approximation $\hat{β} ≃ {(I_{d} + S^{λ})}^{- 1} (I_{d} β_{d} + {\frac{\partial l}{\partial β}|}_{β_{d}}) .$ ${E d l / d β |}_{β_{d}} = 0$ by definition of $β_{d}$ , so taking expectations of both sides we have $E (\hat{β}) ≃ {(I_{d} + S^{λ})}^{- 1} I_{d} β_{d}$ . Hence estimating $I_{d}$ by $\hat{I}$ we have ${\tilde{Δ}}_{β} ≃ {{(\hat{I} + S^{λ})}^{- 1} \hat{I} - I} β_{d}$ . Considering the expected value of ${\tilde{Δ}}_{β} {\tilde{Δ}}_{β}^{T}$ according to the prior mean and variance assumptions of the model, we have the following.

Lemma 1.

Let the setup be as above and let $E_{π}$ denote expectation assuming the prior mean and covariance for $β$ . Treating $\hat{I}$ as fixed, then $V_{\hat{β}} + E_{π} ({\tilde{Δ}}_{β} {\tilde{Δ}}_{β}^{T}) = V_{β} .$

For proof see online SA D. This offers some justification for again using $τ = tr {V_{β} \hat{I}}$ , or $τ = tr (V_{β}^{'} \hat{I})$ accounting for $ρ$ uncertainty. So both the frequentist random effects perspective and the prior expected smoothing bias approach result in (10) $AIC = - 2 l (\hat{β}) + 2 tr (\hat{I} V_{β}^{'}) .$ (10) This is the conventional Hastie and Tibshirani (1990) conditional AIC with an additive correction $2 tr {\hat{I} (V^{'} + V^{''})}$ , accounting for smoothing parameter uncertainty. The correction is readily computed for any model considered here, provided only that the derivatives of $\hat{β}$ and V_β can be computed: the methods of Section 3 provide these. Section 6 provides an illustration of the efficacy of (10).

6. Simulation Results

The improvement resulting from using the corrected AIC of Section 5 can be illustrated by simulation. Simulations were conducted for additive models with true expected values given by η = f₀(x₀) + f₁(x₁) + f₂(x₂) + f₃(x₃), where the f_j are shown in the online SA E, and the x covariates are all independent U(0, 1) deviates. Two model comparisons were considered. In the first a 40 level Gaussian random effect was added to η, with the random effect standard deviation being varied from 0 (no effect) to 1. AIC was then used to select between models with or without the random effect included, but where all smooth terms were modeled using penalized regression splines. In the second case models with and without f₀ were compared, with the true model being based on cf₀ in place of f₀, where the effect strength c was varied from 0 (no effect) to 1. Model selection was based on (i) conventional conditional generalized AIC using τ₀ from Section 5, (ii) the corrected AIC of Section 5, (iii) a version of AIC in which the degrees of freedom penalty is based on τ₁ from Section 5, (iv) AIC based on the marginal likelihood with the number of parameters given by the number of smoothing parameters and variance components plus the number of unpenalized coefficients in the model, and (v) The Greven and Kneib (2010) corrected AIC for the Gaussian response case. The marginal likelihood in case (iv) is a version in which unpenalized coefficients are not integrated out (to avoid the usual problems with fixed effect differences and REML, or improper priors and marginal likelihood).

Results are shown in the top row of Figure 3 for a sample size of 500 with Gaussian sampling error and standard deviation of 2. For the random effect comparison, conventional conditional AIC is heavily biased toward the more complex model, selecting it on over 70% of occasions. The ML based AIC is too conservative for an AIC criterion with 3.5% selection of the larger model when it is not correct, as against the roughly 16% one might expect from AIC comparison of models differing in 1 parameter. The known underestimation of variance components estimated by this sort of marginal likelihood is partly to blame. The AIC based on τ₁ from Section 5 also lacks power, performing even less well than the ML based version. By contrast, the new corrected AIC performs well, and in this example is a slight improvement on Greven and Kneib (2010).^{Figure 4}For the smooth comparison the different calculations differ much less, although the alternatives are slightly less biased toward the more complex model than the conventional conditional generalized AIC, with the corrected Section 5 version showing the smallest bias. The lower row of Figure 3 shows equivalent power plots for the same Gaussian random effect and linear predictor η, but with Bernoulli, beta and Cox proportional hazard (partial) likelihoods (the first two using logit links).

Figure 3. Simulation based illustration of the problems with previous AIC type model selection criteria and the relatively good performance of the Section 5 version. In all panels: (i) the solid curves are for conventional conditional AIC, (ii) the dotted curves are for the Section 5 version, (iii) the middle length dashed curves are for AIC based on the heuristic upper bound degrees of freedom, (iv) the dashed dot curves are for the marginal likelihood based AIC and (v) the long dashed curves are for the Greven and Kneib (2010) corrected AIC (top row only). (a) Observed probability of selecting the larger model as the effect strength of the differing term is increased from zero, for a 40 level random effect and Gaussian likelihood. (b) whole model effective degrees of freedom used in the alternative conditional AIC scores for the left hand panel as effect size increases. (c) Same as (a), but where the term differing between the two models was a smooth curve. (d) As (a) but for a Bernoulli likelihood. (e) As (a) for a beta likelihood. (f) As (a) for a Cox proportional hazards partial likelihood.

Display full size

The purpose of this article is to develop methods to allow the rich variety of smoothers illustrated in Figure 1 to be used in models beyond the exponential family, a task for which general methods were not previously available. However, for the special case of univariate P-splines (Eilers and Marx 1996; Marx and Eilers 1998) some comparison with existing methods is possible, in particular using R package gamlss (Rigby and Stasinopoulos 2005, 2014) and the BayesX package (Fahrmeir and Lang 2001; Fahrmeir, Kneib, and Lang 2004; Brezger and Lang 2006; Umlauf et al. 2015; Belitz et al. 2015, www.bayesx.org). For this special case both packages implement models using essentially the same penalized likelihoods used by the new method, but they optimize localized marginal likelihood scores within the penalized likelihood optimization algorithm to estimate the smoothing parameters.

The comparison was performed using data simulated from models with the linear predictor given above (but without any random effect terms). Comparison of the new method with GAMLSS was only possible for negative binomial, beta, scaled t and simple zero inflated Poisson families, and with BayesX was only possible for the ordered categorical model (BayesX has a negative binomial family, but it is currently insufficiently stable for a sensible comparison to be made). Simulations with both uncorrelated and correlated covariates were considered. Three hundred replicates of the sample size 400 were produced for each considered family at three levels of noise (see SA E for further details). Models were estimated using the correct link and additive structure, and using P-splines with basis dimensions of 10, 10, 15, and 8 which were chosen to avoid any possibility of forced oversmoothing, while keeping down computational time.

Model performance for the negative binomial (nb), beta, scaled t (scat), and zero inflated Poisson (zip) families was compared via MSE, $n^{- 1} \sum_{i = 1}^{n} {\{\hat{η} (x_{i}) - η_{t} (x_{i})\}}^{2},$ on the additive predictor scale. The Brier score, $\frac{1}{n} \sum_{i = 1}^{n} \sum_{j = 1}^{R} {(p_{i j} - {\hat{p}}_{i j})}^{2},$ was used to measure the performance for the ordered categorical (ocat) family, where R is a number of categories, p_ij are true category probabilities and ${\hat{p}}_{i j}$ their estimated values. In addition, the computational performance (CPU time) of the alternative methods was recorded. Figure 4 summarizes the results. In general, the new method provides a small improvement in statistical performance, which is slightly larger when covariates are correlated. The correlated covariate setting is the one in which local approximate smoothness selection methods would be expected to perform less well, relative to “whole model” criteria. In terms of speed and reliability the new method is an improvement, especially for correlated covariates, which tend to lead to reduced numerical stability, leading the alternative methods to fail in up to 4% of cases.

Figure 4. Results of simulation comparison with gamlss (beta, nb, scat, zip) and BayesX (ocat) packages for one dimensional P-spline models. The two plots at lower right show comparisons of log ₁₀ computing times for the case with the smallest time advantage for the new method — Beta regression. The remaining panels show boxplots of replicate by replicate difference in MSE/Brier’s score each standardized by the average MSE or Brier’s score for the particular simulation comparison. Each panel shows three box plots, one for each noise to signal level. Positive values indicate that the new method is doing better than the alternative. Boxplots are shaded grey when the difference is significant at the 5% level (all three for nb correlated should be gray). In all cases where the difference is significant at 5% the new method is better than the alternative, except for the zero inflated Poisson with uncorrelated data, where the alternative method is better at all noise levels.

Display full size

7. Example: Predicting Prostate Cancer

This section and^{Figure 5, Figure 6}the next provide example applications of the new methods, while the online SA F provides further examples in survival analysis and animal distribution modeling. Figure 5 shows representative protein mass spectra from serum taken from patients with a healthy prostate, relatively benign prostate enlargement and prostate cancer (see Adam et al. 2002). To avoid the need for intrusive biopsy there is substantial interest in developing noninvasive screening tests to distinguish cancer, healthy and more benign conditions. One possible model is an ordered categorical signal regression in which the mean of a logistically distributed latent variable z is given by $μ_{i} = α + \int f (D) ν_{i} (D) d D,$ where f(D) is an unknown smooth function of mass D (in Daltons) and ν_i(D) is the $i th$ spectrum. The probability of the patient lying in category 1, 2, or 3 corresponding to “healthy,” “benign enlargement” and “cancer” is then given by the probability of z_i lying in the range ( − ∞, −1], ( − 1, θ] or (θ, ∞), respectively (see online SA K).

Figure 5. Three representative protein mass spectra (centered and normalized) from serum taken from patients with apparently healthy prostate, enlarged prostate, and prostate cancer. It would be useful to be able to predict disease status from the spectra. The red and blue spectra have been shifted upward by 5 and 10 units, respectively.

Display full size

Given the methods developed in this article, estimation of this model is routine, as is the exploration of whether an adaptive smooth should be used for f, given the irregularity of the spectra. Figure 6 shows some results of model fitting. The estimated f(D) is based on a rank 100 thin plate regression spline. Its effective degrees of freedom is 29. An adaptive smooth gives almost identical results. The right panel shows a QQ-plot of ordered deviance residuals against simulated theoretical quantiles (Augustin, Sauleau, and Wood 2012). There is modest deviation in the lower tail. The middle panel shows boxplots of the probability of cancer according to the model for the three observed categories. Cancer and healthy are quite well separated, but cancer and benign enlargement less so. For cases with cancer, the model gave cancer a higher probability than normal prostate in 92% of cases, and a higher probability that either other category in 83% of cases. For healthy patients the model gave the normal category higher probability than cancer in 85% of cases and the highest probability in 77% of cases. These results are somewhat worse than those reported by Adam et al. (2002) for a relatively complex machine learning method which involved first preprocessing the spectra to identify peaks believed to be discriminating. On the other hand the signal regression model here would allow the straightforward inclusion of further covariates, and does automatically supply uncertainty estimates.

Figure 6. Results from the ordered categorical prostate model fit. (a) The estimated coefficient function f(D) with 95% confidence interval. (b) Boxplots of the model probability of cancer, for the 3 observed states (1, healthy, 2, enlarged and 3, cancer). (c) QQ-plot of ordered deviance residuals against simulated theoretical quantiles, indicating some mismatch in the lower tail.

Display full size

8. Multivariate Additive Modeling of Fuel Efficiency

Figure 7 shows part of a dataset on the fuel efficiency of 207 U.S. car models, along with their characteristics (Bache and Lichman 2013). Two efficiency measures were taken: miles per gallon (MPG) in city driving, and the same for highway driving. One possible model might be a bivariate additive model, as detailed in the online SA H, where the two mpg measurements are modeled as bivariate Gaussian, with means given by separate linear predictors for the two components. A priori, it might be expected that city efficiency would be highly influenced by weight and highway efficiency by air resistance and, hence, by frontal area or some other combination of height and width of the car.

Figure 7. Part of a dataset from the USA on fuel efficiency of cars.

Display full size

The linear predictors for the two components were based on the additive fixed effects of factors “fuel type” (petrol or diesel), “style” of car (hatchback, sedan, etc.) and “drive” (all-, front- or rear-wheel). In addition i.i.d. Gaussian random effects of the 22 car manufacturers were included, as well as smooth additive effects of car weight and horsepower. Additive and tensor product smooths of height and width were tried as well as a smooth of the product of height and width, but there was no evidence to justify their inclusion-term selection penalties (Marra and Wood 2011) remove them, p-values indicate they are not significant and AIC suggests that they are better dropped.

The possibility of smooth interactions between weight and horsepower were also considered, using smooth main effects plus smooth interaction formulations of the form f₁(h) + f₂(w) + f₃(h, w). The smooth interaction term f₃ can readily be constructed in a way that excludes the main effects of w and h, by constructing its basis using the usual tensor product construction (e.g., Wood 2006), but based on marginal bases into which the constraints ∑_if₁(h_i) = 0 and ∑_if₂(w_i) = 0 have already been absorbed by linear reparameterization. The marginal smoothing penalties and, hence, the induced tensor product smoothing penalties are unaffected by the marginal constraint absorption. This construction is the obvious generalization of the construction of parametric interactions in linear models, and is simpler than the various schemes proposed in the literature.

The interactions again appear to add nothing useful to the model fit, and we end up with a model in which the important smooth effects are horse power (hp) and weight, while the important fixed effects are fuel type and drive, with diesel giving lower fuel consumption than petrol and all wheel drive giving higher consumption than the two-wheel drives. These effects were important for both city and highway, whereas the random effect of manufacturer was only important for the city. Figure 8 shows the smooth and random effects for the city and highway linear predictors. Notice the surprising similarity between the effects although the city smooth effects are generally slightly less pronounced than those for the highway. The overall r² for the model is 85% but with the city and highway error MPG standard deviation estimated as 1.9 and 2.3 MPG respectively. The estimated correlation coefficient is 0.88.

Figure 8. Fitted smooth and random effects for final car fuel efficiency model. Panels (a)–(c) relate to the city fuel consumption, while (d)–(f) are for the highway. (c) and (f) are normal QQ-plots of the predicted random effects for manufacturer, which in the case of highway MPG are effectively zero.

Display full size

9. Discussion

This article has outlined a practical framework for smooth regression modeling with reduced rank smoothers, for likelihoods beyond the exponential family. The methods build seamlessly on the existing framework for generalized additive modeling, so that practical application of any of the models implemented as part of this work is immediately accessible to anyone familiar with GAMs via penalized regression splines. The key novel components contributed here are (i) general, reliable and efficient smoothing parameter estimation methods based on maximized Laplace approximate marginal likelihood, (ii) a corrected AIC and distributional results incorporating smoothing parameter uncertainty to aid model selection and further inference, and (iii) demonstration of the framework’s practical utility by provision of the details for some practically important models. The proposed methods should be widely applicable in situations in which effects are really smooth, and the methods scale well with the number of smooth model terms. In situations in which some component functions are high rank random fields, then the INLA approach of Rue, Martino, and Chopin (2009) will be much more efficient; however, there are trade-offs between efficiency and stability in this case, since pivoting, used by our method to preserve stability, has instead to be employed to preserve sparsity in the INLA method (see online SA K).

The methods are implemented in R package mgcv from version 1.8 (see online SA M).

Supplemental material

Supplementary Materials

Download PDF (795 KB)

Related Research Data

Smoothing Parameter and Model Selection for General Smooth Models

Source: Figshare

Smoothing Parameter and Model Selection for General Smooth Models

Source: Figshare

Smoothing parameter and model selection for general smooth models

Source: Figshare

Smoothing parameter and model selection for general smooth models

Source: Figshare

Smoothing parameter and model selection for general smooth models

Source: Figshare

Linking provided by

Acknowledgment

We thank the anonymous referees for a large number of very helpful comments that substantially improved the paper and Phil Reiss for spotting an embarrassing error in Supplementary Appendix A.

Supplementary Materials

The online supplementary materials contain additional appendices for the article.

Funding

SNW and NP were funded by EPSRC grant EP/K005251/1 and NP was also funded by EPSRC grant EP/I000917/1. BS was funded by the German Research Association (DFG) Research Training Group “Scaling Problems in Statistics” (RTG 1644). SNW is grateful to Carsten Dorman and his research group at the University of Freiburg, where the extended GAM part of this work was carried out.

References

Adam, B.-L., Qu, Y., Davis, J. W., Ward, M. D., Clements, M. A., Cazares, L. H., Semmes, O. J., Schellhammer, P. F., Yasui, Y., Feng, Z., et al. (2002), “Serum Protein Fingerprinting Coupled With a Pattern-Matching Algorithm Distinguishes Prostate Cancer From Benign Prostate Hyperplasia and Healthy Men,” Cancer Research, 62, 3609–3614. [PubMed], [Web of Science ®], [Google Scholar]
Akaike, H. (1973), “Information Theory and an Extension of the Maximum Likelihood Principle,” in International Symposium on Information Theory, eds. B. Petran, and F. Csaaki, Budapest: Akadeemiai Kiadi, pp. 267–281. [Google Scholar]
Anderssen, R., and Bloomfield, P. (1974), “A Time Series Approach to Numerical Differentiation,” Technometrics, 16, 69–75. [Taylor & Francis Online], [Web of Science ®], [Google Scholar]
Augustin, N. H., Sauleau, E.-A., and Wood, S. N. (2012), “On Quantile Quantile Plots for Generalized Linear Models,” Computational Statistics & Data Analysis, 56, 2404–2409. [Crossref], [Web of Science ®], [Google Scholar]
Bache, K., and Lichman, M. (2013), “UCI Machine Learning Repository,” available at http://archive.ics.uci.edu/ml/. [Google Scholar]
Belitz, C., Brezger, A., Kneib, T., Lang, S., and Umlauf, N. (2015), “Bayesx: Software for Bayesian Inference in Structured Additive Regression Models,” available at http://www.statistik.lmu.de/~bayesx/bayesx.html. [Google Scholar]
Breslow, N. E., and Clayton, D. G. (1993), “Approximate Inference in Generalized Linear Mixed Models,” Journal of the American Statistical Association, 88, 9–25. [Taylor & Francis Online], [Web of Science ®], [Google Scholar]
Brezger, A., and Lang, S. (2006), “Generalized Structured Additive Regression Based on Bayesian p-Splines,” Computational Statistics & Data Analysis, 50, 967–991. [Crossref], [Web of Science ®], [Google Scholar]
Cline, A. K., Moler, C. B., Stewart, G. W., and Wilkinson, J. H. (1979), “An Estimate for the Condition Number of a Matrix,” SIAM Journal on Numerical Analysis, 16, 368–375. [Crossref], [Web of Science ®], [Google Scholar]
Cox, D. (1972), “Regression Models and Life Tables,” Journal of the Royal Statistical Society, Series B, 34, 187–220. [Google Scholar]
Davison, A. C. (2003), Statistical Models, Cambridge, UK: Cambridge University Press. [Crossref], [Google Scholar]
Eilers, P. H. C., and Marx, B. D. (1996), “Flexible Smoothing With B-Splines and Penalties,” Statistical Science 11, 89–121. [Crossref], [Web of Science ®], [Google Scholar]
Fahrmeir, L., Kneib, T., and Lang, S. (2004), “Penalized Structured Additive Regression for Space-Time Data: A Bayesian Perspective,” Statistica Sinica, 14, 731–761. [Web of Science ®], [Google Scholar]
Fahrmeir, L., Kneib, T., Lang, S., and Marx, B. (2013), “Regression Models, New York: Springer. [Crossref], [Google Scholar]
Fahrmeir, L., and Lang, S. (2001), “Bayesian Inference for Generalized Additive Mixed Models Based on Markov Random Field Priors,” Applied Statistics, 50, 201–220. [Crossref], [Web of Science ®], [Google Scholar]
Greven, S., and Kneib, T. (2010), “On the Behaviour of Marginal and Conditional AIC in Linear Mixed Models,” Biometrika, 97, 773–789. [Crossref], [Web of Science ®], [Google Scholar]
Hastie, T., and Tibshirani, R. (1986), “Generalized Additive Models” (with discussion), Statistical Science, 1, 297–318. [Crossref], [Google Scholar]
Hastie, T., and Tibshirani, R. (1990), Generalized Additive Models, London: Chapman & Hall. [Google Scholar]
Kass, R. E., and Steffey, D. (1989), “Approximate Bayesian Inference in Conditionally Independent Hierarchical Models (Parametric Empirical Bayes Models),” Journal of the American Statistical Association, 84, 717–726. [Taylor & Francis Online], [Web of Science ®], [Google Scholar]
Klein, N., Kneib, T., Klasen, S., and Lang, S. (2014), “Bayesian Structured Additive Distributional Regression for Multivariate Responses,” Journal of the Royal Statistical Society, Series C, 64, 569–591. [Google Scholar]
Klein, N., Kneib, T., Lang, S., and Sohn, A. (2015), “Bayesian Structured Additive Distributional Regression With an Application to Regional Income Inequality in Germany,” Annals of Applied Statistics, 9, 1024–1052. [Crossref], [Web of Science ®], [Google Scholar]
Laird, N. M., and Ware, J. H. (1982), “Random-Effects Models for Longitudinal Data,” Biometrics, 38, 963–974. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
Liang, H., Wu, H., and Zou, G. (2008), “A Note on Conditional AIC for Linear Mixed-Effects Models,” Biometrika, 95, 773–778. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
Marra, G., and Wood, S. N. (2011), “Practical Variable Selection for Generalized Additive Models,” Computational Statistics & Data Analysis, 55, 2372–2387. [Crossref], [Web of Science ®], [Google Scholar]
Marra, G., and Wood, S. N. (2012), “Coverage Properties of Confidence Intervals for Generalized Additive Model Components,” Scandinavian Journal of Statistics, 39, 53–74. [Crossref], [Web of Science ®], [Google Scholar]
Marx, B. D., and Eilers, P. H. (1998), “Direct Generalized Additive Modeling With Penalized Likelihood,” Computational Statistics and Data Analysis, 28, 193–209. [Crossref], [Web of Science ®], [Google Scholar]
Nocedal, J., and Wright, S. (2006), “Numerical Optimization (2nd ed.), New York: Springer Verlag. [Google Scholar]
Nychka, D. (1988), “Bayesian Confidence Intervals for Smoothing Splines,” Journal of the American Statistical Association, 83, 1134–1143. [Taylor & Francis Online], [Web of Science ®], [Google Scholar]
Reiss, P. T., and Ogden, T. R. (2009), “Smoothing Parameter Selection for a Class of Semiparametric Linear Models,” Journal of the Royal Statistical Society, Series B, 71, 505–523. [Crossref], [Google Scholar]
Rigby, R., and Stasinopoulos, D. M. (2005), “Generalized Additive Models for Location, Scale and Shape,” Journal of the Royal Statistical Society, Series C, 54, 507–554. [Crossref], [Google Scholar]
Rigby, R. A., and Stasinopoulos, D. M. (2014), “Automatic Smoothing Parameter Selection in GAMLSS With an Application to Centile Estimation,” Statistical Methods in Medical Research, 23, 318–332. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
Rue, H., Martino, S., and Chopin, N. (2009), “Approximate Bayesian Inference for Latent Gaussian Models by Using Integrated Nested Laplace Approximations,” Journal of the Royal Statistical Society, Series B, 71, 319–392. [Crossref], [Google Scholar]
Ruppert, D., Wand, M. P., and Carroll, R. J. (2003), Semiparametric Regression, New York: Cambridge University Press. [Crossref], [Google Scholar]
Säfken, B., Kneib, T., van Waveren, C.-S., and Greven, S. (2014), “A Unifying Approach to the Estimation of the Conditional Akaike Information in Generalized Linear Mixed Models,” Electronic Journal of Statistics, 8, 201–225. [Crossref], [Web of Science ®], [Google Scholar]
Shun, Z., and McCullagh, P. (1995), “Laplace Approximation of High Dimensional Integrals,” Journal of the Royal Statistical Society, Series B, 57, 749–760. [Google Scholar]
Silverman, B. W. (1985), “Some Aspects of the Spline Smoothing Approach to Non-Parametric Regression Curve Fitting,” Journal of the Royal Statistical Society, Series B, 47, 1–53. [Google Scholar]
Umlauf, N., Adler, D., Kneib, T., Lang, S., and Zeileis, A. (2015), “Structured Additive Regression Models: An r Interface to Bayesx,” Journal of Statistical Software, 63, 1–46. [Crossref], [Web of Science ®], [Google Scholar]
Wahba, G. (1983), “Bayesian Confidence Intervals for the Cross Validated Smoothing Spline,” Journal of the Royal Statistical Society, Series B, 45, 133–150. [Google Scholar]
Wahba, G. (1985), “A Comparison of GCV and GML for Choosing the Smoothing Parameter in the Generalized Spline Smoothing Problem,” The Annals of Statistics, pp. 1378–1402. [Crossref], [Web of Science ®], [Google Scholar]
Wood, S. N. (2000), “Modelling and Smoothing Parameter Estimation With Multiple Quadratic Penalties,” Journal of the Royal Statistical Society, Series B, 62, 413–428. [Crossref], [Google Scholar]
Wood, S. N. (2006), “Low-Rank Scale-Invariant Tensor Product Smooths for Generalized Additive Mixed Models,” Biometrics, 62, 1025–1036. [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
Wood, S. N. (2011), “Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation of Semiparametric Generalized Linear Models,” Journal of the Royal Statistical Society, Series B, 73, 3–36. [Crossref], [Google Scholar]
Yee, T. W., and Wild, C. (1996), “Vector Generalized Additive Models,” Journal of the Royal Statistical Society, Series B, 481–493. [Google Scholar]
Yu, D., and Yau, K. K. (2012), “Conditional Akaike Information Criterion for Generalized Linear Mixed Models,” Computational Statistics & Data Analysis, 56, 629–644. [Crossref], [Web of Science ®], [Google Scholar]

Appendix A: Implicit Differentiation in the Extended Gam Case

Let $D_{i}^{\hat{β}}_{j}^{\hat{β}}$ denote elements of the inverse of the Hessian matrix ( $X^{T} WX + S^{λ}$ ) with elements $D_{\hat{β}}^{i}_{\hat{β}}^{j}$ , and note that $\hat{β}$ is the solution of $D_{\hat{β}}^{i} = 0 .$ Finding the total derivative with respect to $θ$ of both sides of this we have $D_{\hat{β}}^{i}_{\hat{β}}^{k} \frac{d {\hat{β}}_{k}}{d θ_{j}} + D_{\hat{β} θ}^{i j} = 0, implying that \frac{d {\hat{β}}_{k}}{d θ_{j}} = - D_{k}^{\hat{β}}_{i}^{\hat{β}} D_{\hat{β} θ}^{i j}$ Differentiating once more yields $\begin{matrix} \frac{d^{2} {\hat{β}}_{i}}{d θ_{j} d θ_{k}} & = & - D_{i}^{\hat{β}}_{l}^{\hat{β}} (D_{\hat{β}}^{l}_{\hat{β}}^{p}_{\hat{β}}^{q} \frac{d {\hat{β}}_{q}}{d θ_{j}} \frac{d {\hat{β}}_{p}}{d θ_{k}} + D_{\hat{β}}^{l}_{\hat{β} θ}^{p}^{j} \frac{d {\hat{β}}_{p}}{d θ_{k}} \\ + D_{\hat{β}}^{l}_{\hat{β} θ}^{p k} \frac{d {\hat{β}}_{p}}{d θ_{j}} + D_{\hat{β} θ θ}^{l j k}) . \end{matrix}$ The required partials are obtained from those generically available for the distribution and link used and by differentiation of the penalty. Generically we can obtain derivatives of D_i w.r.t μ_i and $θ$ .

The preceding expressions hold whether θ_j is a parameter of the likelihood or a log smoothing parameter. Suppose Λ denotes the set of log smoothing parameters, then $D_{β}^{i}_{θ}^{j} = \{\begin{matrix} 2 exp (θ_{j}) S_{i k}^{j} β_{k} & θ_{j} \in Λ \\ D_{β}^{i}_{θ}^{j} & otherwise, \end{matrix}$ where S^j here denotes the penalty matrix associated with θ_j. Similarly $\begin{matrix} D_{β}^{l}_{β}^{p}_{θ}^{j} & = & \{\begin{matrix} 2 exp (θ_{j}) S_{l p}^{j} & θ_{j} \in Λ \\ D_{β}^{l}_{β}^{p}_{θ}^{j} & otherwise \end{matrix} while \\ D_{β}^{l}_{θ}^{j}_{θ}^{k} & = & \{\begin{matrix} 2 exp (θ_{j}) S_{l m}^{j} β_{m} & j = k; θ_{j}, θ_{k} \in Λ \\ D_{β}^{l}_{θ}^{j}_{θ}^{k} & θ_{j}, θ_{k} \notin Λ \\ 0 & otherwise. \end{matrix} \end{matrix}$ Derivatives with respect to $η$ are obtained by standard transformations (A.1) $D_{η}^{i} = D_{μ}^{i} / h_{i}^{'},$ (A.1) where h′_i = h′(μ_i) and more primes indicate higher derivatives. Furthermore, (A.2) $D_{η}^{i}_{η}^{i} = D_{μ}^{i}_{μ}^{i} / h_{i}^{' 2} - D_{μ}^{i} h_{i}^{''} / h_{i}^{' 3},$ (A.2) where the expectation of the second term on the right-hand side is zero at the true parameter values. (A.3) $Also D_{η}^{i}_{η}^{i}_{η}^{i} = D_{μ}^{i}_{μ}^{i}_{μ}^{i} / h_{i}^{' 3} - 3 D_{μ}^{i}_{μ}^{i} h_{i}^{''} / h_{i}^{' 4} [4 p t] + D_{μ}^{i} (3 h_{i}^{'' 2} / h_{i}^{' 5} - h_{i}^{'''} / h_{i}^{' 4}), and [4 p t]$ (A.3) (A.4) $D_{η}^{i}_{η}^{i}_{η}^{i}_{η}^{i} = D_{μ}^{i}_{μ}^{i}_{μ}^{i}_{μ}^{i} / h_{i}^{' 4} - 6 D_{μ}^{i}_{μ}^{i}_{μ}^{i} h_{i}^{''} / h_{i}^{' 5} + D_{μ}^{i}_{μ}^{i} (15 h_{i}^{'' 2} / h_{i}^{' 6} [4 p t] - 4 h^{'''} / h_{i}^{' 5}) - D_{μ}^{i} (15 h_{i}^{'' 3} / h_{i}^{' 7} - 10 h_{i}^{''} h_{i}^{'''} / h_{i}^{' 6} + h_{i}^{''''} / h_{i}^{' 5}) .$ (A.4) Mixed partial derivatives with respect to $η$ / $μ$ and $θ$ transform in the same way, the formula to use depending on the number of η subscripts. The rules relating the derivatives w.r.t $η$ to those with respect to $β$ are much easier: $D_{β}^{i} = D_{η}^{k} X_{k i}, D_{β}^{i}_{β}^{j} = D_{η η}^{k k} X_{k i} X_{k j}, D_{β}^{i}_{β}^{j}_{β}^{k} = D_{η}^{l}_{η}^{l}_{η}^{l} X_{l i} X_{l j} X_{l k} .$ Again mixed partials follow the rule appropriate for the number of β subscripts present. It is usually more efficient to compute using the definitions, rather than forming the arrays explicitly.

The ingredients so far are sufficient to compute $\hat{β}$ and its derivatives with respect to $θ$ . We now need to consider the derivatives of $V$ with respect to $θ$ . Considering $D$ first, the components relating to the penalties are straightforward. The deviance components are then $\begin{matrix} \frac{d D}{d θ_{i}} & = & D_{\hat{η}}^{j} \frac{d {\hat{η}}_{j}}{d θ_{i}} + D_{\hat{θ}}^{i} and \frac{d^{2} D}{d θ_{i} d θ_{j}} = D_{\hat{η} \hat{η}}^{k k} \frac{d {\hat{η}}_{k}}{d θ_{i}} \frac{d {\hat{η}}_{k}}{d θ_{j}} + D_{\hat{η}}^{k} \frac{d^{2} {\hat{η}}_{k}}{d θ_{i} d θ_{j}} \\ + D_{\hat{η} θ}^{k j} \frac{d {\hat{η}}_{k}}{d θ_{i}} + D_{\hat{η} θ}^{k i} \frac{d {\hat{η}}_{k}}{d θ_{j}} + D_{θ θ}^{i j}, \end{matrix}$ where the derivatives of $\hat{η}$ are simply $X$ multiplied by the derivatives of $\hat{β}$ . The partials of $\tilde{l}$ are distribution specific. The derivatives of the determinant terms are obtainable using Wood (2011) once derivatives of w_i with respect to $θ$ have been obtained. These are $\frac{d w_{i}}{d θ_{j}} = \frac{1}{2} D_{\hat{η}}^{i}_{\hat{η}}^{i}_{\hat{η}}^{i} \frac{d {\hat{η}}_{i}}{d θ_{j}} + \frac{1}{2} D_{η \hat{η} θ}^{i i j}, [4 p t]$ $\frac{d^{2} w_{i}}{d θ_{j} d θ_{k}} = \frac{1}{2} D_{\hat{η}}^{i}_{\hat{η}}^{i}_{\hat{η}}^{i}_{\hat{η}}^{i} \frac{d {\hat{η}}_{i}}{d θ_{j}} \frac{d {\hat{η}}_{i}}{d θ_{k}} + \frac{1}{2} D_{\hat{η}}^{i}_{\hat{η}}^{i}_{\hat{η}}^{i} \frac{d^{2} {\hat{η}}_{i}}{d θ_{j} d θ_{k}} + \frac{1}{2} D_{\hat{η}}^{i}_{\hat{η}}^{i}_{\hat{η}}^{i}_{θ}^{k} \frac{d {\hat{η}}_{i}}{d θ_{j}} [4 p t] + \frac{1}{2} D_{\hat{η} \hat{η} \hat{η} θ}^{i i i j} \frac{d {\hat{η}}_{i}}{d θ_{k}} + \frac{1}{2} D_{\hat{η} \hat{η} θ θ}^{i i j k} .$

Theory and Methods

Smoothing Parameter and Model Selection for General Smooth Models

Theory and Methods

Smoothing Parameter and Model Selection for General Smooth Models

ABSTRACTFormulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

1. Introduction

Published online:

2. The General Framework

3. Smoothness Selection Methods

3.1. General Model Estimation

3.1.1. Derivatives and Stable Evaluation of log |Sλ|+

3.1.2. Newton Iteration for <img src="/na101/home/literatum/publisher/tandf/journals/content/uasa20/2016/uasa20.v111.i516/01621459.2016.1180986/20170309/images/uasa_a_1180986_ilm0092.gif" alt="" />β^

3.1.3. Implicit Differentiation

3.1.4. The Remaining Derivatives

3.2. A Special Case: GAMLSS Models

Published online:

3.3. A More Special Case: Extended Generalized Additive Models

3.3.1. Extended GAM New Model Implementation

4. Smoothing Parameter Uncertainty

5. An Information Criterion for Smooth Model Selection

6. Simulation Results

Published online:

Published online:

7. Example: Predicting Prostate Cancer

Published online:

Published online:

8. Multivariate Additive Modeling of Fuel Efficiency

Published online:

Published online:

9. Discussion

Supplementary Materials

Related Research Data

Acknowledgment

Supplementary Materials

Funding

References

Appendix A: Implicit Differentiation in the Extended Gam Case

Article Metrics

Browse journals by subject

Information for

Open access

Help and info

Keep up to date

ABSTRACT
Formulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

3.1.1. Derivatives and Stable Evaluation of log |S^λ|₊

3.1.2. Newton Iteration for $\hat{β}$