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 only the 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:

```r
> sum(m$edf) - sum(m.red$edf)                                          
[1] 3.620239
```

That's how `logLik.gam()` computes the EDF (technically, you want the order reversed

```r
> 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:

```r
> 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:

```r
# 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

```r
pchisq(test_stat, abs(edf_diff), lower.tail = FALSE)
#: [1] 0.009519332
```

which matches the output from `anova.gam()`:

```r
> 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.