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.

Histogram of the response (log transformed)

The other three plots are the fitted vs response plots I am referring to:

  1. the first one (from the left) is the plot from the model (without transforming the variable).
  2. 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.

  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?

fitted vs response plots

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).

enter image description here

  • 6
    $\begingroup$ You need to give more details! What does your response variable represent in real life? Since you use negbin, it is probably a count variable, what does it count? Sample size? Maybe show us some plots? What are your predictor variables? $\endgroup$ Jun 24, 2023 at 1:29
  • $\begingroup$ It is a count variable of number of messages for different events (sum of messages of each event) which is a highly right skewed variable, sample size about 800. I use predictors of number of messages, which I had transformed (square root) to make output more understandable because they are highly skewed themselves. They are strictly positive continuous variables and specifically the IVs are averaged features of these messages (e.g., the average sentiment of all messages for an event) or associated features like the area of the region where events happened $\endgroup$
    – kk68
    Jun 24, 2023 at 8:36
  • $\begingroup$ This is messier than your post admits as your graphs imply using log(count + 1). That's one step even more dubious than the text shows. It's going to be hard to explain and defend what you did, or as you state to relate this to existing literature. Backing up, there is much to be said for trying Poisson regression here: the propaganda about using negative binomial regression instead is exaggerated. $\endgroup$
    – Nick Cox
    Jun 24, 2023 at 9:39
  • 1
    $\begingroup$ Thanks for clarification. Labeling in terms of the original scale … 0.1 1 10 … would avoid the need for such roundabout complications. $\endgroup$
    – Nick Cox
    Jun 24, 2023 at 11:30
  • 3
    $\begingroup$ If the response variable excludes zero values it might be worth trying a truncated negative binomial distribution as the response; you can fit this in glmmTMB using family = truncnbinom[12] (unfortunately, at present you can't combine this with mgcv-style smooths, although you can use something like splines::ns(x, 8) to denote a spline with a pre-specified number of knots/degrees of freedom ... $\endgroup$
    – Ben Bolker
    Jun 24, 2023 at 18:22

1 Answer 1


Transforming the response and then fitting the NB model is not a sensible approach. For one, I assume that you have a count response if you're considering an NB model. But the log of a count variable isn't a count any more. So the NB model would not be appropriate.

What do you mean by a "fitted vs. predicted plot"? Usually, fitted and predicted values are considered synonymous in a regression setting.

EDIT: (Added to address OP's edit)

If your response has a minimum value of 1, then it can't be NB (or Poisson) distributed. Those distributions place a positive mass at 0. Perhaps your response minus 1 could be reasonably described by the NB distribution, though.

The heteroskedasticity that you noticed is normal for NB (or Poisson) distributed responses. Their variance increases with their mean (see, e.g., this reference). You can't interpret Figures 1 and 2 the way you would if you were assessing a linear regression model.

  • $\begingroup$ Some more detail as asked by @kjetil b halvorsen and clarification. You’re right about the fact I have count data (count of messages, highly right skewed variable). You’re right that the log transformation makes them unsuitable to negative binomial. I have overlooked this simple basic fact and the package let me fit the model without giving me an error. About plot, it is indeed the fitted vs response plot (edited). Still, results are empirically good when applying the transformation. I can get a similar plot if I log-transform the fitted vs response values after fitting a standard negbin model $\endgroup$
    – kk68
    Jun 24, 2023 at 8:23
  • $\begingroup$ I also added a few plots, as asked by @kjetil b halvorsen, which might make the situation and the question clearer $\endgroup$
    – kk68
    Jun 24, 2023 at 9:19
  • $\begingroup$ I'll guess this answer would be even more sceptical if you noted that the transformation is log(count + 1), but you must speak for yourself. $\endgroup$
    – Nick Cox
    Jun 24, 2023 at 9:41
  • $\begingroup$ Please check my answer and edited question (the response minimum is 1 and I didn't add 1 for fitting the negative binomial). Sorry for the confusion. $\endgroup$
    – kk68
    Jun 24, 2023 at 9:51
  • 1
    $\begingroup$ Poisson regression doesn't require a marginal Poisson regression. At root the functional form is $Y = \exp(Xb)$ and it could work well across a variety of set-ups. The main concern beyond overall goodness of fit is getting plausible standard errors. See e.g. blog.stata.com/2011/08/22/… $\endgroup$
    – Nick Cox
    Jun 25, 2023 at 5:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.