About Multivariate Generalized Regression Mixed Models with negative binomial distribution

Question

I have a research question about the intra-individual variations of contact patterns before and during the lockdown. I will explain the definition of contact patterns and the model I built with test data. My question is I am not sure whether the model I construct can answer the research question.

Contact patterns refer to the age-specific, location-stratified contact rates. Here is a plot that may give you a hint.

. I take the number of contact and the location as the outcome, since this is count data and its distribution is right skewed. So my options are three: 1) Multivariate Generalized Linear Mixed Models with Poisson distribution, 2) Multivariate Generalized Regression Mixed Models with negative binomial distribution, and 3) Multivariate Generalized Regression Mixed Models with negative binomial distribution. After comparing the AIC, the model I choose is the second one.

Here is the test data and the model I built using R.

data <- data.frame(
  contacts_education = rpois(100, lambda = 5),
  contacts_work = rpois(100, lambda = 10),
  contacts_household = rpois(100, lambda = 8),
  contacts_others = rpois(100, lambda = 3),
  token = rep(as.character(c("A", "B", "C", "D", "E")), each = 20), 
  age = rep(c(56, 77, 15, 32, 27), each = 20),
  sex = factor(sample(c("Male", "Female"), 100, replace = TRUE)),
  hh_size = rnorm(100, mean = 3, sd = 1),
  city = factor(sample(c("CityA", "CityB", "CityC"), 100, 
                       replace = TRUE)),
  period = factor(sample(c("1", "2", "3"), 100, replace = TRUE))
)

data_long <- data %>%
  pivot_longer(cols = starts_with("contacts"), 
               names_to = "location", 
               values_to = "num_contacts")

nb_model <- glmmTMB(num_contacts ~ period * location * (I(age/100) + 
    sex + I(hh_size) + city) + (1 | location) + (1|token ) + (1|city), 
                    family = nbinom2, 
                    data = data_long)

Here, token is the participant id, age is the participant age, sex is the participant gender, hh_size is the household size. Because the research question is about the intraindividual variation, so I only include participants who reports at least contact information at two time points, even they are at the same period. period refers to before lockdown, during lockdown with seasonal effect, and during lockdown.