Generalized additive model

Data

Response values: \(y_1, \ldots, y_N\) with \(y_i \in \mathcal{Y} \subset \Re\)

Covariate values: \(x_1, \ldots, x_N\) with \(x_i \in \Re^p\)

Model¹

\(Y_1, \ldots, Y_N\) are independent conditionally on the covariate values

\(Y_i \mid x_i\) has an exponential family distribution with mean \(\mu_i = E(Y_i \mid x_i)\) and dispersion parameter \(\phi\), with

\[ g(\mu_i) = \alpha^\top w_i + \sum_{j = 1}^q f_j(z_{ij}) \]

• \(w_i\) and \(z_i\) are sub-vectors of \(x_i\)
• \(\alpha\) is a vector of unknown scalar parameters, \(f_1, \ldots, f_q\) are unknown smooth functions

Spline basis expansions

Typically, \(f_1, \ldots f_q\) are defined via spline basis expansions, resulting in the linear predictor structure \[ g(\mu_i) = \alpha^\top w_i+ \sum_{j = 1}^q \beta_j^\top b_{ij} \] where \(b_{ij}\) is determined solely by the basis expansion and the covariates \(z_{ij}\).

Maximum penalized likelihood estimation

Estimation typically proceeds by maximizing the penalized log-likelihood \[ l_P = l - \frac{1}{2\phi} \sum_{j=1}^q \lambda_j \beta_j^\top S_j \beta_j \]

• \(\beta_j^\top S_j \beta_j\) is a quadratic roughness penalty

• \(\lambda_1, \ldots, \lambda_q\) are positive smoothing parameters

• \(l\) denotes the log-likelihood as a function of \(\theta = (\alpha^\top, \beta^\top)^\top\)

• \(\phi\) is treated as fixed

The penalty prevents overfitting by balancing function fit to the data and smoothness.

Penalized (left) and unpenalized (right) smooth fits.

Maximum penalized likelihood estimation

Regression parameters

For fixed \(\lambda\), a P-IRLS iteration can be used to estimate \(\alpha\) and \(\beta\), under identifiability constraints for the functions¹.

Smoothing parameters

\(\lambda_j\) can be estimated using REML, GCV, UBRE via an outer iteration around P-IRLS.

See, for example, ?mgcv::gam for details, and other options.

Dispersion parameter

\(\phi\) can be estimated using the adjusted Pearson statistic estimator of Fletcher (2012).

Bias

As is the case for maximum likelihood (ML) estimation for GLMs, maximum penalized likelihood (MPL) estimators for GAMs can also suffer from finite-sample bias.

Beyond standard sampling bias, bias may also be amplified by the smoothing penalty¹.

Bias reduction for ML estimation

Explicit methods

\(\tilde\theta = \hat\theta - \hat{B}\)

e.g. Asymptotic bias correction¹, Jackknife², Bootstrap³

Implicit methods

e.g. Bias-reducing adjusted score equations⁴

The solution \(\tilde\theta\) of \(s + A = 0\), with \(s = \nabla l\), has \(E(\tilde\theta - \theta) = o(n^{-1})\) if \[ A_t = \frac{1}{2} {\rm trace}\left\{ i^{-1}\left(P_t + Q_t\right) \right\} \]

\(i = E(s s{}^\top)\), \(P_t = E( s_t s s{}^\top)\), \(Q_t = E(s_t \nabla \nabla^\top l)\)

Bias-reducing adjusted score equations

Bias reduction using adjusted score equations is becoming a standard approach, because it does not require repeated fitting and has useful side-effects.

e.g. in logistic regression¹

• If the model matrix is of full rank, \(\tilde\theta\) has all of its components finite.

• \(\tilde\theta\) shrinks the model towards equi-probability across observations relative to the ML estimator, which typically leads to smaller estimated variances.

• \(\tilde\theta\) has the same asymptotic distribution expected by standard ML estimators.

Adjusted score equations in GAM estimation

The adjusted score equations approach cannot be directly applied to GAMs, because:

• Even when using spline basis expansions, a GAM may have more parameters than observations, and, hence, the information matrix is rank deficient.

• The existence of the smoothing penalty results in bias in the estimating equations.

Bias-reducing adjustments for GAMs

Adjusted penalized score equations

Let \(\tilde\theta\) be the solution of \(\nabla l - \frac{1}{\phi} S_\lambda \theta + A_P = 0\), where \(A_P = O(1)\).

Bias expansion

If \(\lambda_j = O(1)\) and \(\phi\) is treated as fixed, then the bias of \(\tilde\theta\) is

\[ E(\tilde\theta - \theta) = -\underbrace{i_P^{-1} A}_{\text{sampling}} - \overbrace{\frac{1}{\phi} i_P^{-1} S_\lambda \theta}^{\text{smoothing}} + \underbrace{i_P^{-1} A_P}_{\text{adjustment}} + O(n^{-3/2}) \]

• \(i_P = i + S_\lambda / \phi\), assumed invertible for \(\lambda_j > 0\) under identifiability constraints

• \(S_\lambda\) is such that \(\theta^\top S_\lambda \theta = \sum_{j=1}^q \lambda_j \beta_j^\top S_j \beta_j\)

• \(A_t = {\rm trace}\{ i_P^{-1} (P_t + Q_t)\} / 2\)

Bias-reducing adjustments for GAMs

Bias-reducing adjustments to the penalized score equations

If \(\lambda_j = O(1)\) and \(A_P = A + S_\lambda \theta / \phi\), then \(E(\tilde\theta - \theta) = o(n^{-1})\).

No smoothing

For \(\lambda = 0\), \(i_P = i\) and \(S_\lambda = 0\). Hence, we recover the bias-reducing adjusted score equations for standard GLMs.

Implementation

For GAMs, implementations of \(i\), \(P_t\) and \(Q_t\) as functions of \(\theta\) can be computed by the enrichwith R package (using the corresponding glm object).

P-IRLS can be extended to solve the adjusted score equations. The caveat is that as for GLMs¹, \(\alpha\), \(\beta\) and \(\phi\) need to be updated at each iteration.

The smoothing parameters can be estimated in the same way as for maximum penalized likelihood.

Simulation study

\(Y_1, \ldots, Y_N\) are independent conditionally on covariates \(w_1, \ldots w_N, z_1, \ldots, z_N\).

\(Y_i \mid w_i, z_i \sim \operatorname{Binomial}(10, \pi_i)\) with

\[ \log\frac{\pi_i}{1 - \pi_i} = \alpha_0 + \sum_{k=1}^{10} \alpha_k w_{ik} + f_1(z_{i1}) + f_2(z_{i2}) \]

• \(f_1(x) = \sin(17.5 x^2)\)

• \(f_2(x) = \begin{cases} 2\sin(\pi x/1.5), & 0\le x<0.75,\\ 2, & 0.75\le x\le1. \end{cases}\)

\(z_{i1}\) and \(z_{i2}\) are realisations of independent \(\operatorname{Unif}(0,1)\) random variables.

Scalar covariates and effects

\(g_{i1}, \ldots, g_{i10}\) are realisations of independent standard normal variables, \(I(\cdot)\) is the indicator function, and \([a]\) is the floor of \(a\).¹

\(\alpha_0\) is set so that \(E(Y_i | w_i, z_i) \approx 4\).

Setup

\(N = 100\)

\(30\) knots per function with a second-order P-spline¹

\(1000\) samples

\(\lambda_1 = \lambda_2 = \lambda \in \{0.1, 0.2, 0.4, 0.8, 1.6, 3.2\}\) along with the data-driven mgcv smoothing parameter estimates

Estimation of scalar effects

Bias for the ten scalar effects

Estimation of scalar effects

MSE for the ten scalar effects

Estimation of \(f_1(\cdot)\)

Estimation of \(f_2(\cdot)\)

Smoothing parameter estimators

The bias expansion is derived under the assumption that the smoothing parameters are \(O(1)\), but usual smoothing parameter estimators (e.g. using the methods implemented in mgcv) need not have \(O(1)\) expectation.

\(\log\hat\lambda\) against \(\log N\); blue line has slope \(1/2\)

Nevertheless, experiments suggest the method still delivers improved estimation with data-driven smoothing parameter selection.

Abalone data

Data

Attributes of \(N=4177\) abalones, a type of shellfish¹:

• Response: an age proxy, with value \(1\) if the abalone has more than \(10\) rings, and \(0\) otherwise

• Covariates: sex (Male, Female, Infant), and seven weight and length measurements

Model

Independent Bernoulli responses conditional on the covariates with probability \(\pi_i\), such that

\[ \log\frac{\pi_i}{1 - \pi_i} = \alpha_0+\alpha_1 w_{i1}+\alpha_2 w_{i2}+ \sum_{j=1}^7f_j(z_{ij}) \]

\(w_{i1}, w_{i2}\): two dummy variables encoding Sex with infant as baseline

\(z_{i1}, \ldots, z_{i7}\): seven weight and length measurements

Abalone-based simulation design

1. Represent each functional component using a P-spline basis constructed from \(10\) knots

2. Take covariate subsamples for \(N \in \{50, 100, 200, 400\}\)

3. Compute the estimates of the scalar and functional effects using mgcv on the full data set

4. Generate \(500\) response vectors for each covariate sub-sample using the estimates from step 3 as truth

Estimation of scalar effects

Bias of the estimators of the scalar effects

Estimation of scalar effects

Variance of the estimators of the scalar effects; A value of \(3\) indicates variance \(>3\).

Estimation of functional effects

Average estimates of the seven functional effects for \(n = 400\)

Remarks

Adjusted penalized score equations provide a practical route to improved GAM estimation. They

• Target both sampling and smoothing-induced bias

• Reduce bias in scalar effects

• Improve bias and boundary behaviour for smooth functions

• Induce extra regularization that appears to stabilize fitting with small samples