# Comparing log likelihood of "nested" GAM models (mgcv)

I am using the mgcv package to fit a GAM model to my data set to determine whether the trajectories along two (or more) different processes are different. More specifically I fit the following two models:

null = gam(y ~ s(time, bs='cr', k=6), family = "nb", knots = knotList, weights = weights, method="ML")
alt = gam(y ~ process + s(time, by=process, bs='cr', k=6), family= "nb", knots=knotList, weights = weights, method="ML")


where process is a factor that defines which process each observation belongs to (edit of this post contains additional details about the set-up, but my question is really just related to the above).

I now want to compare how well each of these models fit my data and have examined the values of logLik(null) and logLik(alt). My expectation was that the null model was nested inside of the alternative model and hence logLik(null) $$\leq$$ logLik(alt). However, I find that in some cases logLik(null) $$>$$ logLik(alt). Is this just likely the result of some small numerical/optimization inaccuracy and I should just replace cases where logLik(alt) $$-$$ logLik(null) $$<0$$ with zero? Or are these models not nested as I thought? Does the penalization applied to the likelihood during fitting complicate this comparison?

> summary(null)

Family: Negative Binomial(7.948)

Formula:
y ~ s(time, bs = "cr", k = 6)

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.39277    0.03935   60.81   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(time) 3.907  4.422  90.24  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.429   Deviance explained = 37.3%
-ML = 488.53  Scale est. = 1         n = 240

> summary(alt)

Family: Negative Binomial(7.662)

Formula:
y ~ process + s(time, by = process, bs = "cr",
k = 6)

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)     2.37912    0.05695  41.775   <2e-16 ***
process2  0.02401    0.07957   0.302    0.763
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(time):process1 3.322  3.868  44.20  <2e-16 ***
s(time):process2 3.512  4.049  42.16  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.421   Deviance explained = 36.8%
-ML = 492.37  Scale est. = 1         n = 240
$$$$


You should probably be comparing

gam(y ~ process + s(time, bs = 'cr', k = 6), family = "nb", knots = knotList)


with

gam(y ~ process + s(time, by = process, bs = 'cr', k = 6), family = "nb",
knots = knotList)


If the aim is to determine if you need process-specific smooths rather than a single common smooth.

It's difficult to give concrete advice as to what is happening — because you haven't provided any information on the fits, like summary output etc — but I'd be looking at how wiggly the smooths in the two models are; in the first model the smooth could be modelling differences in group means and fit the data better by using many more EDF than are used by the model with two smooths. Just because your null model has fewer model terms doesn't mean that it is less complex than a model with more terms.

Some notes:

1. why are you passing formulas as character vectors?
2. you almost surely want method = "REML" or method = "ML"` when fitting models, and when comparing models with different fixed effects (as here) where those terms aren't fully penalized, then you shouldn't use REML so ML is the way to go (at least for the comparisons).
• Thanks @Gavin! First, my apologies about the typo in specifying the model formulas, I’ve edited it. I’ve also edited my question to include an example of the model summaries where logLik(null) > logLik(alt). I’ve set the method to “ML” for both models; could you elaborate on why this method is more appropriate? It does seem that the edf of the smooth term in the null model is greater than the edf of the individual smooth terms in the alt model. Could I perhaps specify that the edf of the smooth term of the null model should not exceed the max edf of the smooth terms of the alt model? Jun 8, 2022 at 22:14
• I also did not include “process” as a term in my null model as I wanted the null model to consider no difference between the two processes, so no difference in mean along with no difference in smooths. Including “process” as a term in the null model also generates cases where logLik(null) > logLik(alt). Jun 8, 2022 at 22:15