Different Results of the Same GAM model depends on "discrete = TRUE"

Question

I am using bam function fitting GAMs.

The code for the model is:

    mod_1 <- bam(n ~ s(age, by = period, k = 15) + 
               s(hh_size, by = period,  k = 9) +
               period +
               s(token, bs = "re") + s(Bundesland, bs = "re") + s(period, bs = "re"),
             data = halle_data_household,
             method = "fREML", discrete = TRUE,
             family = nb(),
             nthreads = c(4,1))
    mod_2 <- bam(n ~ s(age, by = period, k = 15) + 
               s(hh_size, by = period,  k = 9) +
               period +
               s(token, bs = "re") + s(Bundesland, bs = "re") + s(period, bs = "re"),
             data = halle_data_household,
             method = "fREML", 
             family = nb(),
             nthreads = c(4,1))

The difference between the two models is "discrete = TRUE" to faster the computational time. But I got two different results.

The result of the mod_1

    Family: Negative Binomial(13.366) 
    Link function: log 
    
    Formula:
    n ~ s(age, by = period, k = 15) + s(hh_size, by = period, k = 9) + 
        period + s(token, bs = "re") + s(Bundesland, bs = "re") + 
        s(period, bs = "re")
    
    Parametric coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept)   0.5409     0.1354   3.996 6.74e-05 ***
    period2      -0.1253     0.1674  -0.749    0.454    
    period3      -0.1319     0.1600  -0.825    0.410    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Approximate significance of smooth terms:
                             edf  Ref.df     F  p-value    
    s(age):period1     1.000e+00   1.000 0.641 0.423337    
    s(age):period2     4.683e+00   5.656 3.799 0.001034 ** 
    s(age):period3     4.569e+00   5.444 4.928 0.000116 ***
    s(hh_size):period1 1.960e+00   2.412 1.186 0.303077    
    s(hh_size):period2 1.000e+00   1.000 1.798 0.180227    
    s(hh_size):period3 1.000e+00   1.000 5.930 0.014997 *  
    s(token)           4.406e+02 859.000 1.242  < 2e-16 ***
    s(Bundesland)      5.429e-05  15.000 0.000 0.801883    
    s(period)          1.082e-14   3.000 0.000 0.999933    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    R-sq.(adj) =  0.323   Deviance explained = 42.4%
    fREML = 3657.9  Scale est. = 1         n = 2115

The result of mod_2

    Family: Negative Binomial(13.32) 
    Link function: log 
    
    Formula:
    n ~ s(age, by = period, k = 15) + s(hh_size, by = period, k = 9) + 
        period + s(token, bs = "re") + s(Bundesland, bs = "re") + 
        s(period, bs = "re")
    
    Parametric coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept)  0.32541    0.05293   6.148 9.83e-10 ***
    period2     -0.12288    0.07204  -1.706  0.08824 .  
    period3     -0.20909    0.06826  -3.063  0.00222 ** 
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Approximate significance of smooth terms:
                              edf  Ref.df     F  p-value    
    s(age):period1      1.000e+00   1.000 0.641 0.423576    
    s(age):period2      4.685e+00   5.661 3.801 0.001061 ** 
    s(age):period3      4.568e+00   5.442 4.963 0.000117 ***
    s(hh_size):period1  1.961e+00   2.417 1.036 0.311454    
    s(hh_size):period2  1.000e+00   1.000 1.811 0.178523    
    s(hh_size):period3  1.000e+00   1.000 5.930 0.014986 *  
    s(token)            4.405e+02 857.000 1.226  < 2e-16 ***
    s(Bundesland)       6.066e-05  15.000 0.000 0.802106    
    s(period)          -9.299e-16   3.000 0.000 0.031014 *  
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    R-sq.(adj) =  0.322   Deviance explained = 42.4%
    fREML = 3406.3  Scale est. = 1         n = 2115

Answer 1

I see where your comment about negative EDF came from! That's just floating point precision for you. You can treat that s(period) as being irrelevant given the other terms in the model and ignore the p value. You might also consider sending this as a bug report to Simon Wood, mgcv's maintainer as a negative EDF doesn't make sense (and here it isn't really negative it's just so small that floating point issues have cropped up and made the number that is numerically 0 [or very close to it] less than 0.)

Anyway, the bigger problem and the cause of your issues is that you are including the period means twice:

1. through the period parametric term, and also
2. through the s(period, bs = "re") random effect term.

Those two terms are not uniquely identifiable from the data. What you are seeing a numerical differences putting more (or less) weight on one of those two things that likely depend on the data used in the model; with discrete = TRUE the predictor data is slightly different - it has been discretized slightly and that could easily see the weight shift from the parametric version to the random effect version of the group mean effects you included in the model.

Basically, use one or the other, not both. You don't need to include the group means via a parametric effect if you fit a factor by; you just need to account for the group means in the model and you can do that either by a parametric term, a random effect term, or some other term that also includes the group means (and "fs" or a "sz" smooth will include the group means for example)