There are quite a few things wrong in your derivation or code. I can't check your derivation of the LRT from the likelihoods just now, but somewhere you have gone wrong as your code yields
> lrt
[1] 13.98288
> lrt2
[1] 35780.14
and the only latter is the correct value (up to the sign).
Another problem is exactly how the EDF of the models is computed and thence their difference. Both your methods yield the same value:
> df
[1] 3.620239
> df2
[1] 3.620239
This turns out to be the difference of the standard EDF calculation for the smooths in the model:
> sum(m$edf) - sum(m.red$edf)
[1] 3.620239
That's how logLik.gam() computes the EDF (technically, you want the order reversed
> sum(m.red$edf) - sum(m$edf)
[1] -3.620239
but that doesn't really matter for the purposes of computing the test statistic and the p values.)
In anova.gam(), the EDFs of the models are computed using somewhat different definitions of the EDF. If it's available we'd prefer the values in $edf2, which include a correction to the EDF of each basis function for the fact that we selected the smoothing parameters. The EDF you computed and what is in $edf assumes the value(s) of the smoothing parameter(s) were known before you fitted the model. However, the smoothing corrected EDFs $edf2 are only available when you do smoothness selection using method = "REML" or method = "ML". If $edf2 is not available, $edf1 is used instead of $edf. ?gamObject somewhat unhelpfully has this to say about $edf1
edf: estimated degrees of freedom for each model parameter.
Penalization means that many of these are less than 1.
edf1: similar, but using alternative estimate of EDF. Useful for
testing.
edf2: if estimation is by ML or REML then an edf that accounts for
smoothing parameter uncertainty can be computed, this is it.
‘edf1’ is a heuristic upper bound for ‘edf2’.
If you want the gory details they are in sections 6.1.2 and 6.12.1 of Simon's GAM book (Wood, 2017). The main difference between $edf and $edf1 is that the latter accounts for smoothing bias, and simulation studied by Simon and people in his lab have shown that this smoothing bias-corrected version of the EDF is preferred over the normal EDF. If we have $edf2 we'd ideally use that, because it is corrected for smoothing bias like edf1 but also includes a correction due to selecting the smoothing parameters.
So, the correct EDF for the test in this case of models fitted with GCV is:
> sum(m.red$edf1) - sum(m$edf1)
[1] -4.50227
You wouldn't have gone too far wrong not knowing about the EDF thing, but the final problem is that your manual calculation of the test statistic for the LRT is incorrect as it doesn't scale the difference in deviance of the two models (the LR) by the dispersion parameter of the full model. anova.gam() ultimately ends up calling stat.anova(..., scale = m$sig2) and the test statistic that it computes is:
# EDF diff
edf_diff <- sum(m.red$edf1) - sum(m$edf1)
edf_diff
#: [1] -4.50227
# LR
d_dif <- deviance(m) - deviance(m.red)
d_dif
#: [1] -35780.14
# test statistic
scale <- m$sig2 # dispersion parameter of the reference model
test_stat <- d_dif / scale * sign(edf_diff)
test_stat
#: [1] 14.3143
and then the p value of the test is computed as
pchisq(test_stat, abs(edf_diff), lower.tail = FALSE)
#: [1] 0.009519332
which matches the output from anova.gam():
> anova(m, m.red, test = "LRT")
Analysis of Deviance Table
Model 1: resp ~ s(pred)
Model 2: resp ~ 1
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 94.498 238412
2 99.000 274192 -4.5023 -35780 0.009519 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
References
Wood, S.N., 2017. Generalized Additive Models: An Introduction with R, Second Edition. CRC Press.
