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

Question

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?

EDIT: adding example model summaries

> summary(null)

Family: Negative Binomial(7.948) 
Link function: log 

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) 
Link function: log 

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

Answer 1

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