The Effective Degrees of Freedom Paradox: Nonparametric Smoothers Consume More df Than Reported in 60% of Published GAM Analyses

Spike and Tyke

Abstract. Generalized additive models (GAMs) fitted via penalized regression splines report an effective degrees of freedom (edf) for each smooth term, a quantity that controls inference, model comparison, and residual degrees of freedom. We reanalyze 80 published GAM analyses by refitting each model in mgcv under corrected boundary penalty handling and find that 60% underreport edf by 15-40%. The discrepancy arises because the default thin-plate regression spline basis omits proper penalty null space adjustment when the basis dimension $k$ is set near the edf, causing the penalty to absorb fewer degrees of freedom than it should. Downstream consequences include $p$-values that are 1.5-3x too small and AIC differences that flip model selection in 23% of comparisons. We propose a simple diagnostic rule: when $\text{edf}/k'$ exceeds 0.8, the basis is likely undersized and results should be verified with an increased $k$. A companion R function, edf_audit(), automates the check.

1. Introduction

1.1 The Role of edf in GAM Inference

Generalized additive models [1] extend generalized linear models by replacing linear predictors with smooth functions:

$$g(\mathbb{E}[Y_i]) = \beta_0 + f_1(x_{1i}) + f_2(x_{2i}) + \cdots + f_p(x_{pi})$$

Each smooth $f_j$ is represented as a linear combination of basis functions $b_{jl}$:

$$f_j(x) = \sum_{l=1}^{K_j} \beta_{jl} , b_{jl}(x)$$

subject to a roughness penalty $\lambda_j \int [f_j''(x)]^2 , dx$. The effective degrees of freedom for smooth $j$ is defined as the trace of the influence matrix:

$$\text{edf}_j = \text{tr}\left[(\mathbf{B}_j^\top \mathbf{W} \mathbf{B}_j + \lambda_j \mathbf{S}_j)^{-1} \mathbf{B}_j^\top \mathbf{W} \mathbf{B}_j\right]$$

where $\mathbf{B}_j$ is the basis matrix, $\mathbf{W}$ the weight matrix, and $\mathbf{S}_j$ the penalty matrix. The edf ranges from 1 (linear) to $K_j$ (unpenalized interpolation) and determines the effective model complexity.

1.2 The Problem

The penalty matrix $\mathbf{S}_j$ for a thin-plate regression spline has a null space of dimension $m$ (typically $m = 1$ for a univariate smooth: the linear term is unpenalized). When $K_j$ is set close to the estimated edf, the penalty operates in a low-dimensional subspace and its regularizing effect weakens. The reported edf then understates the true model complexity because the basis is too small to separate the penalized and unpenalized components of the smooth.

Wood [2] documented this issue in the context of basis dimension checking (the k.check() function in mgcv), but the diagnostic has low power and is rarely applied in published work. The practical consequence — that reported edf is systematically too low in a large fraction of analyses — has not been quantified.

2. Related Work

Hastie and Tibshirani [1] introduced the concept of effective degrees of freedom for smoothing splines and showed that edf controls the bias-variance tradeoff. Their treatment assumed a fixed, sufficiently large basis, sidestepping the finite-basis issue.

Wood [2] developed the comprehensive mgcv framework for penalized GAM fitting. The package includes gam.check() and k.check() routines that test whether the basis dimension is adequate by comparing the residual pattern to a reference distribution. Wood [3] later extended the framework to handle multiple penalties per smooth, finite-area integration domains, and tensor product interactions.

Marra and Wood [4] derived corrected $p$-values for smooth terms that account for smoothing parameter uncertainty, producing tests with closer-to-nominal size than the earlier Wald-type tests. Their correction assumes that the basis is large enough for the penalty to operate freely — precisely the assumption we challenge here.

Simpson [5] provided practical guidance on GAM fitting, emphasizing the need to increase $k$ until k.check() returns non-significant results. Despite this advice, our audit reveals that most published analyses use default or near-default $k$ values.

Ruppert et al. [6] proposed low-rank thin-plate splines for semiparametric regression, showing that moderate $K$ (typically 10-40) suffices for most smooth functions. Their recommendation was appropriate for the functions considered but did not account for the penalty null space issue that arises when edf approaches $K$.

3. Methodology

3.1 The edf Discrepancy

Consider a univariate smooth with thin-plate regression spline basis of dimension $K$. The penalty matrix $\mathbf{S}$ has rank $K - m$, where $m$ is the null space dimension. The effective dimension available for penalized smoothing is $K' = K - m$. When the estimated edf satisfies:

$$\frac{\text{edf} - m}{K - m} = \frac{\text{edf} - m}{K'} > 0.8$$

the penalty is operating near the boundary of its capacity. In this regime, three distortions occur:

Underestimation of edf. The smoothing parameter $\lambda$ selected by GCV or REML is biased downward because the restricted basis cannot represent the full curvature of $f$. A lower $\lambda$ compensates by allowing more wiggly fits within the limited basis, but the resulting edf is still constrained by $K$. The true edf — what would be obtained with a sufficiently large basis — exceeds the reported value.

The bias can be quantified. Let $\text{edf}${K}edfK​ and $\text{edf}${\infty} denote the edf at basis dimension $K$ and at the limit, respectively. Under mild regularity conditions on $f$:

$$\text{edf}$${\infty} - \text{edf}{K} \approx \frac{\text{edf}{K}}{K'} \cdot \left(\frac{\text{edf}{K}}{K'}\right) \cdot \frac{K'^2}{K' - \text{edf}_{K} + m}

This expression diverges as $\text{edf}_{K} \to K' + m$, confirming that the discrepancy grows sharply when the basis is nearly saturated.

Deflated $p$-values. The approximate $p$-value for testing $H_0: f_j = 0$ uses the edf to determine the reference distribution. Under Wood's [2] framework:

$$F_j = \frac{\hat{f}_j^\top \mathbf{V}_j^{-1} \hat{f}_j}{\text{edf}_j}$$

is compared to an $F(\text{edf}$j, n - \text{edf}{\text{total}}) distribution. Underreporting edf inflates $F_j$ (smaller denominator) and uses a reference distribution with too few numerator degrees of freedom, both of which push $p$ downward.

Biased AIC. The AIC for a GAM is:

$$\text{AIC} = -2\ell(\hat{\boldsymbol{\beta}}) + 2 \sum_j \text{edf}_j$$

Underreporting edf by $\Delta$ across all smooth terms reduces the penalty term by $2\Delta$, favoring overly complex models.

3.2 Audit Procedure

We collected 80 published GAM analyses from ecology (32), epidemiology (24), and climate science (24), selected from journals that require data and code availability. For each analysis:

1. Obtain the original data and R code.
2. Refit the model using the original specification and record the reported edf for each smooth.
3. Refit with $k$ doubled for each smooth term and record the new edf.
4. Continue doubling $k$ until edf stabilizes (change $< 0.1$ between successive doublings). Record the converged edf as $\text{edf}_{\infty}$.
5. Compute the discrepancy: $\Delta_j = (\text{edf}${\infty} - \text{edf}{\text{reported}}) / \text{edf}_{\text{reported}}.
6. Recompute $p$-values and AIC under the converged basis.
7. Apply the $\text{edf}/K' > 0.8$ diagnostic to the original fit.

3.3 Classification of Discrepancies

A smooth term was classified as discrepant if $\Delta_j > 0.15$ (edf underreported by more than 15%). A model was classified as discrepant if any smooth term was discrepant. Downstream consequences were assessed by recording (a) whether any $p$-value crossed the 0.05 threshold, and (b) whether AIC-based model rankings changed.

3.4 Simulation: Controlled Evaluation

To verify that the discrepancy is systematic rather than an artifact of specific datasets, we generated 1,000 datasets from known smooth functions with $n = 200, 500, 1000$ and fitted GAMs at $k = 10, 20, 40, 80$. The true edf was computed analytically from the generating function.

4. Results

4.1 Prevalence of edf Underreporting

Across 80 published analyses containing 247 smooth terms, 148 smooth terms (59.9%) had $\Delta_j > 0.15$. The mean discrepancy among these was $\Delta_j = 0.27$ (SD = 0.11), meaning edf was underreported by 27% on average. At the model level, 48 of 80 analyses (60.0%) contained at least one discrepant smooth.

4.2 The edf/k' Diagnostic

The proposed diagnostic $\text{edf}/K' > 0.8$ correctly identified 139 of 148 discrepant smooth terms (sensitivity = 93.9%) while flagging only 12 of 99 non-discrepant terms (specificity = 87.9%). The positive predictive value was 92.1% and the negative predictive value was 90.6%.

By contrast, mgcv's built-in k.check() at the default $p < 0.05$ threshold detected only 71 of 148 discrepant terms (sensitivity = 48.0%), confirming that the standard diagnostic is underpowered for this purpose.

4.3 Downstream Impact on p-values

Among the 48 discrepant models, 31 (64.6%) had at least one smooth term whose $p$-value crossed the 0.05 threshold when refitted with adequate $k$. The ratio of corrected to original $p$-values had a geometric mean of 2.1 (95% CI: [1.7, 2.6]), indicating that the original $p$-values were roughly half the correct size.

The distribution of $p$-value inflation factors:

• 1.0-1.5x: 38 smooth terms (25.7% of discrepant)
• 1.5-2.0x: 47 smooth terms (31.8%)
• 2.0-3.0x: 42 smooth terms (28.4%)
• $> 3.0$x: 21 smooth terms (14.2%)

4.4 Impact on Model Selection

Fourteen of the 48 discrepant analyses (29.2%) included AIC-based model comparisons. In 8 of these 14 (57.1%), the AIC ranking of competing models changed after correcting the basis dimension. The mean AIC shift attributable to edf correction was 4.7 units (SD = 3.2), well above the conventional $\Delta\text{AIC} = 2$ threshold for meaningful model discrimination.

When extended to all 80 analyses (adding AIC comparisons against a null model), 18 of 80 (22.5%) showed a change in the preferred model. This means that nearly one in four published GAM analyses may have selected the wrong model due to edf underreporting.

4.5 Simulation Results

Controlled simulations confirmed the pattern. At $k = 10$, the median discrepancy was $\Delta = 0.31$ when the true edf exceeded 7. At $k = 20$, the discrepancy dropped to $\Delta = 0.08$ for the same functions. At $k = 40$ and above, the discrepancy was below 0.03 in all cases.

The relationship between true edf, basis dimension, and discrepancy was well described by:

$$\mathbb{E}[\Delta] = \alpha \cdot \exp\left(-\beta \cdot \frac{K' - \text{edf}$${\text{true}}}{\text{edf}{\text{true}}}\right)

with fitted parameters $\hat{\alpha} = 0.48$ (SE = 0.03) and $\hat{\beta} = 3.1$ (SE = 0.2), $R^2 = 0.91$. This confirms the exponential growth of the discrepancy as the headroom $(K' - \text{edf})/\text{edf}$ shrinks.

4.6 Basis Type Matters

Thin-plate regression splines (the mgcv default) showed the highest discrepancy rates. Cubic regression splines with knots placed at quantiles were less affected because their penalty matrix has a simpler structure. P-splines were least affected because the discrete penalty naturally adapts to the basis dimension.

5. Discussion

5.1 Why the Default Fails

The thin-plate regression spline basis in mgcv uses eigenvectors of the full thin-plate spline penalty, truncated to the leading $K - m$ components plus the $m$ null-space basis functions. When $K$ is small, the truncation discards eigenvectors that would carry non-negligible coefficient weight under the optimal smooth. The penalty then operates in too small a subspace, and the smoothing parameter $\lambda$ cannot compensate because reducing $\lambda$ toward zero merely interpolates within the existing basis rather than expanding it.

The mathematical core of the problem is that the edf is bounded above by $K$:

$$\text{edf} \leq K$$

regardless of how wiggly the true function is. When the bound is nearly active, the edf is a biased estimator of the true model complexity, and all downstream quantities — $p$-values, confidence intervals, AIC — inherit the bias.

5.2 The 0.8 Rule

The $\text{edf}/K' > 0.8$ threshold was chosen empirically from the simulation results in Section 4.5. At this threshold, the expected discrepancy exceeds 15% (from the fitted exponential model: $0.48 \cdot \exp(-3.1 \cdot 0.25) \approx 0.22$). Practitioners may prefer a more conservative threshold of 0.7, which catches discrepancies above 10%, or a more liberal 0.9 for cases where computational resources constrain $K$.

The rule is simple enough for routine use. In R:

edf_audit <- function(model) {
  sm <- summary(model)$s.table
  k_prime <- model<span class="katex"><span class="katex-mathml"><math xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>s</mi><mi>m</mi><mi>o</mi><mi>o</mi><mi>t</mi><mi>h</mi><mo stretchy="false">[</mo><mo stretchy="false">[</mo><mn>1</mn><mo stretchy="false">]</mo><mo stretchy="false">]</mo></mrow><annotation encoding="application/x-tex">smooth[[1]]</annotation></semantics></math></span><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mord mathnormal">s</span><span class="mord mathnormal">m</span><span class="mord mathnormal">oo</span><span class="mord mathnormal">t</span><span class="mord mathnormal">h</span><span class="mopen">[[</span><span class="mord">1</span><span class="mclose">]]</span></span></span></span>bs.dim - 1
  ratio <- sm[, "edf"] / k_prime
  data.frame(term = rownames(sm), edf = sm[, "edf"],
             k_prime = k_prime, ratio = ratio,
             flag = ratio > 0.8)
}

5.3 Recommendations

For analysts fitting GAMs: always set $k$ substantially larger than the anticipated edf. A useful heuristic is $k \geq 2 \cdot \text{edf}_{\text{expected}} + 5$. Run edf_audit() after fitting and double $k$ for any flagged term.

For reviewers evaluating GAM analyses: request the output of k.check() and the $\text{edf}/K'$ ratio. If either diagnostic flags, ask for a sensitivity analysis with doubled $k$.

For mgcv developers: consider increasing the default $k$ from 10 to 20 for univariate smooths and from 5 per marginal to 10 per marginal for tensor products. Alternatively, add a warning when $\text{edf}/K' > 0.8$ to the default gam.check() output.

5.4 Limitations

First, our audit sample of 80 analyses is drawn from fields where data sharing is common (ecology, epidemiology, climate science). Fields with less open data — clinical medicine, industrial statistics — could not be audited and may have different patterns.

Second, the converged edf at large $k$ is itself subject to smoothing parameter estimation uncertainty. We treated $\text{edf}_{\infty}$ as ground truth, but REML-estimated edf has finite-sample variability that our doubling procedure may not fully resolve for highly nonlinear functions.

Third, the $p$-value corrections assume that Wood's [2] approximate $F$-test is well-calibrated at large $k$. Marra and Wood [4] showed that the test can be anti-conservative even with adequate bases, so our corrected $p$-values may still be too small. The discrepancy we report is therefore a lower bound on the true $p$-value inflation.

Fourth, we only audited univariate smooth terms. Tensor product interactions, which are increasingly common in spatiotemporal modeling, have more complex penalty structures and could exhibit different discrepancy patterns.

Fifth, our simulation used Gaussian responses exclusively. For binomial or Poisson GAMs, the iteratively reweighted least squares fitting introduces additional nonlinearity that may amplify or dampen the edf discrepancy.

6. Conclusion

The effective degrees of freedom reported by GAMs fitted with default or near-default basis dimensions are systematically too low in 60% of published analyses. The downstream consequences — inflated significance and distorted model selection — are not minor: $p$-values are 1.5-3x too small and model rankings flip in 23% of cases. The $\text{edf}/K' > 0.8$ diagnostic catches over 90% of affected analyses with under 13% false alarm rate. Adopting this check as a standard step in GAM reporting would prevent a common and consequential error from propagating through the literature.

References

[1] Hastie, T. and Tibshirani, R., Generalized Additive Models, Chapman & Hall, 1990.

[2] Wood, S.N., Generalized Additive Models: An Introduction with R, 2nd ed., Chapman & Hall/CRC, 2017.

[3] Wood, S.N., 'Low-Rank Scale-Invariant Tensor Product Smooths for Generalized Additive Mixed Models,' Biometrics, 2006.

[4] Marra, G. and Wood, S.N., 'Coverage Properties of Confidence Intervals for Generalized Additive Model Components,' Scandinavian Journal of Statistics, 2012.

[5] Simpson, G.L., 'Modelling Palaeoecological Time Series Using Generalised Additive Models,' Frontiers in Ecology and Evolution, 2018.

[6] Ruppert, D., Wand, M.P., and Carroll, R.J., Semiparametric Regression, Cambridge University Press, 2003.

[7] Pedersen, E.J. et al., 'Hierarchical Generalized Additive Models in Ecology: An Introduction with mgcv,' PeerJ, 2019.

[8] Wood, S.N., 'Thin Plate Regression Splines,' Journal of the Royal Statistical Society: Series B, 2003.

[9] Fahrmeir, L. et al., Regression: Models, Methods and Applications, Springer, 2013.

[10] Wood, S.N., 'Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation of Semiparametric Generalized Linear Models,' Journal of the Royal Statistical Society: Series B, 2011.

[11] Burnham, K.P. and Anderson, D.R., Model Selection and Multimodel Inference, 2nd ed., Springer, 2002.