Likelihood ratio statistic differs when computed from logLik, deviance, and anova functions in R

Question

I'm confused about getting different results when trying to perform the likelihood ratio test (LRT) with R in different ways that are supposed to be equivalent. Below is the R code with some simulated data. I use GAM model with one predictor, whose significance I want to evaluate by comparing the models with and without that predictor.

As far as I know, the theory says the LRT statistic can be computed either from the log likelihoods of the original and the reduced models, or from their residual deviances, as both ways lead to identical calculations:

$$LRstatistic = Deviance_{reduced} - Deviance_{orig} \\ = -2×(logLik_{reduced}-logLik_{saturated}) + 2x(logLik_{orig}-logLik_{saturated}) \\ = -2×(logLik_{reduced} - logLik_{orig})$$

But in the example below, the LRT statistic computed from the log likelihoods (as provided by the logLik function) is 17.7, while the same statistic computed from the residual deviances (as provided by the deviance function) is 38555.7. It is way far from identical.

Moreover, when I perform the LRT test using the anova function, it reports the same deviances and residual DF, but the final p value differs from either of the two manual calculations.

Could anyone give me some insight why I see those differences and which of the three ways is actually correct? Thanks in advance!

# data
set.seed(100)
pred <- runif(100)
resp <- 0.5*pred - 5*pred^2 - 50*pred^3 + rnorm(100, sd=50)
plot(pred, resp)

# model fit
m <- gam(resp~s(pred))         # original model
m.red <- update(m, ~.-s(pred)) # reduced model

# LRT using logLik
ll <- logLik(m)                 # log likelihood of the orig. model
npar <- attr(ll, "df")          # number of estimated parameters in the orig. model
ll.red <- logLik(m.red)         # log likelihood of the reduced model
npar.red <- attr(ll.red, "df")  # number of estimated parameters in the reduced model
lrt <- -2*(ll.red[1] - ll[1])   # likelihood ratio statistic (17.7)
df <- npar - npar.red           # degrees of freedom (3.2)
1-pchisq(lrt, df)               # p value (0.0007)

# LRT using deviance
dev <- deviance(m)            # residual scaled deviance of the orig. model
dfr <- df.residual(m)         # residual degrees of freedom of the orig. model
dev.red <- deviance(m.red)    # residual scaled deviance of the reduced model
dfr.red <- df.residual(m.red) # residual degrees of freedom of the reduced model
lrt2 <- dev.red - dev         # likelihood ratio statistics (38555.7)
df2 <- dfr.red - dfr          # degrees of freedom (3.2)
1-pchisq(lrt2, df2)           # p value (0)

# LRT using anova
anova(m, m.red, test="LRT") # p value (0.001)
```

Answer 1

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:

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