GAM summary F values

Does anyone know how the F values are computed in summary(gam.object). I've looked through the help page, and the gam package, and the book "Statistical models in S", but I've yet to come across where this is explicitly laid out. For example, I would really like to know how to understand the tables in

library(gam)
data(gam.data)

fit=gam(y~z+s(x), data=gam.data)
summary(fit)


I would really to know how the F values are computed.

Here is the paper on them:

http://biomet.oxfordjournals.org/content/early/2012/10/18/biomet.ass048.short

Basically, the F-tests on smooth terms (rather than for the full model) are joint tests for equality to zero for all of the coefficients making up a single spline term. So, a model that you specify as

gam(y~s(x))

is actually
$$y = \alpha + \displaystyle\sum_k\gamma_kb(x) + \epsilon$$ (in the gaussian case). $b(x)$ is your basis function and the $\gamma$'s, along with the intercept, are estimated. Given that the $\gamma$'s jointly represent a single term, a t-test of equality to zero for any one of them doesn't help you much if you want to know if the whole term is significant. So you need to test them jointly.

This is straightforward when there is no penalization, the typical Wald statistic is $\hat\beta^T(\hat{V}_{\hat\beta})^{-1}\hat\beta$. You could do this to test any subset of the parameters of an unpenalized regression model, if you just subset your coefficients and your parameter covariance matrix appropriately for the test you want. (In case it isn't obvious to you, a T-test is what results when $\beta$ and $V$ are scalars).

The problem when coefficients are penalized is that the typical F-test is overpowered. You will get values that are too-high (over-rejecting the null). This is because the penalized coefficients are being scaled by a variance matrix that doesn't take their penalization into account. Hence Simon Wood uses $\hat\beta^T(\hat{V}_{\hat\beta})^{-r}\hat\beta$, which is a "rank $r$ pseudoinverse". Sort of like taking a scalar $s^{-r}$ rather than $s^{-1}$ -- if $r<1$ (and greater than zero), you get a larger number back.

Here is a simple example. First I fit a gam to a straight-line DGP and everything is penalized to OLS. Both F-stats are the same:

> set.seed(1)
> N <- 500
> x <- rnorm(N)
> y <- x + rnorm(N, sd = .1)
> m <- gam(y~s(x)-1)
> summary(m)

Family: gaussian

Formula:
y ~ s(x) - 1

Approximate significance of smooth terms:
edf Ref.df     F p-value
s(x)   1      1 44103  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.989   Deviance explained = 98.9%
GCV = 0.011509  Scale est. = 0.011486  n = 500
> Fnaive <- coef(m) %*% solve(m$Vp) %*% coef(m) > Fnaive [,1] [1,] 44102.9  Now I make the DGP sinusoidal, and the F-stats start to diverge: > y <- sin(x) + rnorm(N, sd = .1) > m <- gam(y~s(x)-1) > summary(m) Family: gaussian Link function: identity Formula: y ~ s(x) - 1 Approximate significance of smooth terms: edf Ref.df F p-value s(x) 7.943 8.721 2304 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.976 Deviance explained = 97.6% GCV = 0.010555 Scale est. = 0.010387 n = 500 > Fnaive <- coef(m) %*% solve(m$Vp) %*% coef(m)
> Fnaive
[,1]
[1,] 20089.11

• Can we compare the F statistics between GAM terms and say they represent the strength of evidence for including each term in the model? Sep 3, 2021 at 2:56

If I'm understanding your question correctly, it has more to do with Anova than GAMs. For the sake of illustration, here's the output you're referencing:

# > summary(fit)
#
# Call: gam(formula = y ~ z + s(x), data = gam.data)
# Deviance Residuals:
#   Min       1Q   Median       3Q      Max
# -0.68109 -0.21431  0.02907  0.24551  0.53192
#
# (Dispersion Parameter for gaussian family taken to be 0.0841)
#
# Null Deviance: 57.7496 on 99 degrees of freedom
# Residual Deviance: 7.9077 on 93.9999 degrees of freedom
# AIC: 44.0543
#
# Number of Local Scoring Iterations: 2
#
# Anova for Parametric Effects
# Df Sum Sq Mean Sq  F value Pr(>F)
# z          1  0.003   0.003   0.0343 0.8535
# s(x)       1 38.313  38.313 455.4321 <2e-16 ***
#   Residuals 94  7.908   0.084
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Anova for Nonparametric Effects
# Npar Df Npar F     Pr(F)
# (Intercept)
# z
# s(x)              3 45.484 < 2.2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


There are two sets of "F-values" that you might be referring to, but I'm going to assume you mean the ones in the Anova table for parametric effects (for example, an F-value of 455.4321 for s(x)). Correct me if I'm wrong here.

The F-test is basically the ratio of the explained and unexplained variance. Using s(x) from your model, it's calculated as

summary(fit)$parametric.anova[2,2] / summary(fit)$parametric.anova[3,3]
# [1] 455.4321


which is the sum of squares divided by the residuals, which you can also calculate directly from the R output

$$\frac{38.313}{0.084} = 455.4321$$

• Yes, thank you. I was really curious about both. But, from what I understand the "F-value" is computed by For each nonparametric term in the model the nonlinear component is set to zero and the parametric part of the model is refit by weighted least squares, holding the other nonlinear components fixed." which seems like it's not so easy to just calculate like you did above. Jul 3, 2015 at 18:36