Model misfit with DHARMa - What needs/can be done?

Question

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:

1. How bad is the fit? Are the discrepancies too severe or am I still able to report on the fitted models results?
2. 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.

Answer 1

Interesting problem! In addition to what Florian has suggested, here are my thoughts:

1. The mixed effects models you fitted may not be the best for teasing out the effect of person-level predictors (just recently I came upon a reference discussing the challenge with interpreting such effects - I'll see if I can find it and add it to this post);

2. If you stick with mixed effects models (as opposed to, say, GEE style models), it may be worth trying to fit them using the bam() function from the mgcv package, which is designed to accommodate large datasets. This may cut back on computational time. See https://stat.ethz.ch/pipermail/r-help/2016-April/438227.html for how to specify random effects in a bam() call.

3. When you are trying various refinments of your model, can you fit them on a smaller random sample of subjects just to get a sense of how reasonable they would be? Then you can fit them for all available subjects.

4. Within the context of using mixed effects models, I think you may need to be more careful with how you model the effect of time in your model. It seems that currently you model the effect of time as being linear (?). If you plot the residuals for each of your model against month (1 through 12) and year (e.g., 2018, 2019, etc.), you will have an opportunity to see if there is any systematic temporal structure in those residuals that was left unaccounted for by the model. If necessary, you could try replacing the time predictor in your model with a month predictor and a year predictor. Then perhaps include something like month + year or even month + year + month:year with month coded as a factor and year as a numeric; or go via the GAM route and include smooth additive effects of month (coded as numeric in R) and year (coded as numeric) or a smooth interaction between month and year.

Answer 2

I think what's pretty clear is that the chosen distributions don't fit the data very well. I don't find this particular surprising. For example active days gambling per month is not really a count variable, as the month has a strict max (30/31) days, so if you count how many of those days someone gambles, this is more like a k/n binomial. The other two cases are less obviously not count data, but if you look at the data-generating process, there is nothing that really resembles a classical poisson respectively gamma process. In particular

• For the second model, I would also consider fitting this as a k/n where you fit the probability that someone hits the loss limit, given that they play (k = losses, n = plays)

• For the gamma, I would simply try other models or transformations of the response to fit the data better.

Regarding your questions of "how important" these misfits are: as a rule of thumb, the most important thing for p-values is to get the dispersion roughly correct. As you fit variable dispersion models, this is out of the way. Beyond this, the impact of residual problems is typically smaller, but what I see here seems large enough for me to have a meaningful impact on p-values or estimates, so I would try to change the models to get a better fit.