Background
I am working with a large dataset that contains longitudinal data on gambling behavior of 184,113 participants. The data is based on complete tracking of electronic gambling behavior within a gambling operator. Gambling behavior data is aggregated on a monthly level, a total of 70 months. I have an ID variable separating participants, a time variable (months, 0-69), as well as numerous gambling behavior variables. I am especially interested in investigating the role of age and gender in predicting 3 gambling behavior outcomes. Age is handled categorically with 6 categories, 18-29, 30-39, 40-49, 50-59, 60-69, 70+ years. Gender is dichotomous male/female. Gambling behavior outcomes are active days gambling per month (count variable), number of times a person meets a set loss limit for each month (count variable), and total bet size divided on number of active days in a month (continuous variable, termed “gambling intensity).
Question
I have run several models for each outcome variable and selected the model with best fit according to AIC scores. However, when I use the DHARMa package in R to do further model diagnostics on chosen models it reveals that the simulated residuals from the fitted models do not fit well with the observed residuals. See below for plots and glmmTMB code for models. I have some questions regarding the DHARMa results and how to proceed:
- How bad is the fit? Are the discrepancies too severe or am I still able to report on the fitted models results?
- Are there any simple modifications I can do to improve the models?
Some comments regarding the second question: I was thinking about attempting some different distributions/likelihood functions. So, for the negative binomial models I have originally used the traditional negative binomial 2 model, but I could try NB1 as well. For the gamma model I could try to specify the same model with GLMMadaptive which I understand uses a different approach than Laplace approximation (which glmmTMB uses). I could also try another link function, but all the other models have log links which makes the overall results interpretations harmonize better. Regarding adding new predictors to improve fit: Unfortunately, there are no more non-gambling behavior predictors available to me in the data set. Age and gender are the only person-related variables.
The reason that I have not already tried all of these possible modifications before posting is that the computation time for each new model is >4 days so I want new analyses choices to be well-founded. Computation time is due to the big data size, there are 3,231,544 observations.
I have already used several months on data analysis, so a practical concern is to avoid too much more time being used on tweaking the models. Naturally, I want to reconcile this with reporting models that are good enough.
DHARMa plots for each model
Active gambling days model
Loss limit model
Gambling intensity model
More information on models
Active gambling days in a month can take values between 1-31 depending on specific month. Zero active days are treated as missing, thus a zero-truncated negative binomial model was used. glmmTMB code for model:
DaysPlayedConditionalAgeGenderTruncated <- glmmTMB(daysPlayed ~ 1 + time + ageCategory * gender + (time | id), dfLong, family = truncated_nbinom2)
Loss limits met per month can take values between 0-31 depending on specific months, meaning loss limits are specific to each day played. Most people won’t reach the loss limits, thus a zero-inflated negative binomial model was used. glmmTMB code for model:
LossLimitConditionalAgeGenderTruncated <- glmmTMB(lossLimit ~ 1 + time + ageCategory * gender + (time | id), ziformula =~1, dfLong, family = nbinom2)
Gambling intensity, total bet size (in currency) divided on number of active days in a month, the variable can take on continuous positive values above 0. It is also skewed. Thus, a gamma distribution model was used. glmmTMB code for model:
gamblingIntensityConditionalAgeGender <- glmmTMB(gamblingIntensity ~ 1 + time + ageCategory * gender + (time | id), dfLong, family = Gamma(link = “log”))
PS: For those interested in more information on the model comparisons that were done before the DHARMa simulations of the selected models, please see this other post (the same predictors were compared for the gamma model): Do conventional thresholds for global fit indices (e.g. AIC) hold for models based on very large data sets?
UPDATE New attempts and some follow-up questions
Thank you Florian and Isabella for your suggestions! I am a PhD student and our research group lacks experience with these type of models so I appreciate the assistance. I have done some more analyses on a random sub-sample of 50,000 participants. Here is an update on my attempts so far and some follow-up questions. New attempts:
- New time predictors: I tried the same models with month as factor and year as numeric and with a random intercept. DHARMa simulations for model 1 (days gambling) and model 2 (loss limits met) perform almost identical as before. Model 3 did not converge (gambling intensity, still a gamma for this analysis).
- Variations of model 3 (gambling intensity): I tried to model this as a linear mixed model with log(gambling intensity). DHARMa simulations appear almost identical to the generalized linear mixed gamma model. I also tried to use the GLMMadaptive package which I understand uses a different estimation technique. The model did not converge.
Does the fact that the DHARMa simulations vary very little when using different time predictors, going from random slope to random intercept, and (in the case of model 3) going from gamma distribution to normal distribution tell us anything? Florian also suggested to use a k/n binomial for models 1 and 2. I am not familiar with this. Are there any accessible sources to recommend on the topic? I am not too interested in the models in themselves but rather the estimates being adequate. Florian says the model fit likely has a meaningful impact on p-values and estimates. Does this mean that there is a higher probability of type 1 errors? Could the big N offset some of this weakness? Predicted effects from the ggeffects package are all highly significant.
I tried to look into General Additive Models, but I am unsure if I can still account for the zero-inflation (in loss limit data) and zero-truncation (in days gambled data). We briefly considered General Estimating Equations early in the project but did not go with these because we could not find a solution/r package to deal with the zero issue within them. In terms of research aims however, GEE would fit nicely as we are primarily interested in population effects. That is, showing differences in high-risk gambling between different groups of age (hence the 6-category approach to age) and gender.