Using GAM to investigate the within-subject variation between several periods

Question

my research question is to analyse intra-individual (within-subject) heterogeneity in contact patterns and identify determinants of stable contact behaviour before and during the lockdown. Since the contact patterns are associated with the age of the participants and the location of the contact, the contact patterns of the study population are visualized in the plot below.

As shown in the plot, since the association between contact number and age stratified by location seems to follow different patterns, the contact data is stratified by location. And models are built according to different locations.

First I build a preliminary model to estimate the theta

prelim_model <- gam(n ~ s(age, by = period, k = 5) +  
                      s(token, bs = "re") + s(Bundesland, bs = "re"), 
                    data = data_education, 
                    family = nb())

theta_est <- prelim_model$family$getTheta(TRUE) 
theta_est 

Because I want to examine the effect of variables on the number of contacts according to period and within-subject, after comparing the AIC and BIC, the final model is like this:

gam(n ~ s(age, by = period, k = 3) + period * sex + s(hh_size, by = period, k = 3)  + 
               period * employment_2cat + period * education_2cat + period * marital_2cat + 
               period * weekend +  period * region + period * transportation_2cat + period + 
               s(token, bs = "re") + s(Bundesland, bs = "re"), 
             data = data_education, 
             family = nb(theta = theta_est)) # using the theta I got from the preliminary model

From the results:

• marital status is marginal statistically siginificant
• neither of the interaction terms is statistically siginificant
• age has significant non-linear effect on the number of contacts
• household size has a significant non-lineat effect on the number of contacts only for period 3
• Personal behaviour greatly influences contact numbers.

I use DHARMa residuals and gam.check to check the model

My questions are as follows:

1. Is this method reasonable that I divided the dataset according to the location and build model seperately?
2. I have already tried other methods, including GLMM with poisson distribution, GLMM with ng distribution, Zero-inflated model with ng distribution. When I check the model with DHARMa residuals, the results are much worse.
3. Because the high Std.Errors for the Intercept, period2, period3 from the analysis result, and the siginificant divviation from the QQ plot residuals, and as suggested by the result of gam.check, how should I improve the model?