I fit a negative binomial GAM model with the R mgcv package. I noticed some heteroskedasticity in the fitted vs response plots. Then I noticed that if I log transform the dependent variable before fitting the model, it becomes far better in terms of the fitted vs response values: it does not show any heteroskedasticity and also the simulateResiduals test from the DHARMa package looks good (differently from the model with the non-log-transformed variable).
However, to my knowledge, log transforming the dependent variable before fitting a negative binomial is uncommon, as the negative binomial already uses a log link. Hence this would be like performing a consecutive log transformation. Is that acceptable? Or have you any idea of an alternative and more straightforward strategy to apply in such a situation?
Edit: I add a few plots as I was asked to. The first is the histogram of the response.
The other three plots are the fitted vs response plots I am referring to:
- the first one (from the left) is the plot from the model (without transforming the variable).
- The second one is the same plot by applying a log transformation to both the fitted and the response variables: I applied this transformation simply because the first plot was unreadable, and I then realized that there was this good but heteroskedastic relation that became evident when taking the log.
EDIT: I only added 1 to avoid unimportant negative values in the plot, as the minimum value of the fitted is 0.08 (equal to -2.525729 in log scale), but the response variable minimum is 1.
- The third plot is the fitted vs response plot of the negative binomial model fitted to the log-transformed variable (I did not add anything to its values, as the minimum is 1). As correctly said by @Rachel Altman, this is not a proper way to proceed. Still, the output is interesting, and it strikes me that it looks so good. Why does this happen? What am I missing? Is there anything I can do to improve the model by including the information derived from these exploratory analyses?
EDIT: I also tried to fit a Poisson model. The following are the fitted vs response plots. In this case, the DHARMa simulateResiduals plots are totally out. However, if I run a correlation between fitted vs response, the correlation of the Poisson model is very high (Pearson 0.9880688, Spearman 0.5639613) and even higher than the one of the (wrongly) log-transformed negative binomial (Pearson 0.8398934, Spearman 0.771207) – whose plot looks, however, better to me (and the DHARMa simulateResiduals plots do not identify relevant issues).