What I have tried and done so far:

I am running GLMs and GAMs on my positive, continuous response variable. I have determined that the Gamma distribution would fit my data best by plotting various QQ plots with EnvStats::qqPlot() function against possible known distributions (Gaussian, Lognorm, Weibull, Gamma).
Then, I started using the Gamma family with it's canonical link function (which is "inverse" in R) but models were bad. Following several statistics books (e.g. Dunn&Smyth 2018) I tried other links like log link, but that was even worse. Only when I tried the "identity" link within Gamma family argument I finally started to get nicer models where diagnostic plots from gratia::appraise() are acceptable and in general the model results make sense. Also, plotting the models with gratia::draw() finally yielded partial effect plots where the residuals actually fall along the smoother lines and mostly inside the confidence bands.
I have learned here, here, here and there that the combination is not ideal as the identity link could potentially yield negative predictions. Yet, it seems to be used in some cases and be regularly tested as one potential option (e.g. Dunn&Smyth 2018). Similar to the cited CV posts I receive these warnings when I run the models:

In log(ifelse(y == 0, 1, y/mu)) : NaNs produced

Unlike the other CV posts, I never get the other part, i.e. the error message:

Error: no valid set of coefficients has been found: please supply starting values

My questions now are as follows:

  1. Does the missing error message mean I am still kind of okay?
  2. Considering that all alternatives with other links always perform worse in all regards (explained Deviance, diagnostic plots, plotted models), no matter which other predictors I use and how I try changing the model formula: Can I still keep the identity link, if I am cautious?
  3. If yes, are there indicators to look out for that would tell me the prediction of negative values turned from a theoretical to an actual problem?
    • I have read that negative intercepts are a bad sign: can that be an indicator?
    • Would I actually see negative fitted values in the fitted vs residual plot of gratia::appraise() output if the negative predictions became a real issue?
  4. If all the above is not an option: What else can I do to get similarly nice results but avoiding the identity link problem?

Please excuse me for not providing a reprex: I can't reproduce my problem with dummy data. But I know that several of my observations have zeros for some of the predictor values. I also know that my predictors are often not normally distributed, therefore have often used transformed predictor variables, but that does not seem to improve the models with log link or inverse link in any way; actually, the models with untransformed predictor vars tend to have higher explained Deviance and adj. R².

Update: Details

For illustrating what I mean when I say that in comparison to log link and inverse link models the identity models look better I have created two graphs derived from 3 versions of one and the same gamma model, fitted with this formula mgcv::gam(y ~ s(var1) + s(var2) + s(var3) + s(var4) +grouping_factor, family=Gamma(link="xyz"),method="REML",data = data) (whereby they only differ wrt the link specified instead of "xyz")
The first graph has diagnostic plots of all 3 models generated by gratia::appraise(). I have marked in orange, what looked suspicious to me and in blue what I am starting to suspect looks wrong for you: enter image description here The second graph shows the gratia::draw() outputs for the same 3 models. enter image description here

Still grateful for any help and suggestions, Thanks in advance!

  • $\begingroup$ What do you think is gamma-distributed from QQ plots? Do you mean the entirety of your $y$ variable? $\endgroup$
    – Dave
    Commented Aug 15, 2023 at 18:52
  • $\begingroup$ @ Dave: Yes, I guess that is what the QQ plot tells me. I also did a lot of histograms of the response variable and even tried what happens if I subject the response variable to various possible transformations including BoxCox, but have read in stats books that manually transforming the response variable is kind "of a last resort" if nothing else helps. Also I'd rather take the option of using an existing link function to ease back-transformation of predictions etc. PS: Code looked like this: EnvStats::qqPlot(data$RespVar, distribution = 'gamma', estimate.params = T, add.line = T) $\endgroup$ Commented Aug 15, 2023 at 19:00
  • $\begingroup$ Have you looked at the residuals and other model output? Especially wrt your out-of-sample data? $\endgroup$
    – user78229
    Commented Aug 15, 2023 at 19:05
  • $\begingroup$ @Mike Hunter: Yes, I have looked at the residuals a lot as they are being part of the model diagnostic plots, e.g. produced by gratia::appraise() function: it plots QQ, residuals vs linear predictor, observed vs fitted and a histogram of residuals. Regarding out-of-sample-data: I only have a small dataset with <100 observations, hence I was under the impression that it's too small to do out-of-sample. (sorry, had mixed up appraise() & draw() in the previous post, corrected it here and clarified above) $\endgroup$ Commented Aug 16, 2023 at 9:47
  • $\begingroup$ It's useful to note that, without oos data, these are calibration plots. 100 data points? That's not too few to do a bootstrap or random forest-type iterative, approximating approach to oos. Some of the circled residuals indicate a lack of fit while others show pretty good fit. Is time a feature? If so, plots of residuals by time would be helpful. Also, predicted vs actual data. $\endgroup$
    – user78229
    Commented Aug 16, 2023 at 10:51

1 Answer 1

  1. You can't just look at a histogram of the untransformed response and decide what the distribution is. This is a regression problem, so you need to see the residuals from a correctly specified model to decide what distribution best fits the conditional response, making the whole thing a complete paradox.

  2. The two errors couldn't be more different. If you look at the code for Gamma you'll find that the step dev.resids <- function(y, mu, wt) -2 * wt * (log(ifelse(y == 0, 1, y/mu)) - (y - mu)/mu) contains the issue. The mean vector mu and response y have both evaluated to 0, since in R 0/0=NaN. Debug by setting trace=T in the glm so you can inspect the iterations and you can recreate the problem. The most obvious solution is to drop the identity link because of the challenges fitting this type of model.

  3. The lack of error message does not mean that the fit is well behaved, as you see from the calculation of the deviance residuals, the identity link will not place much emphasis on observations very close to 0, but the deviance residuals based on the ratio can be quite large.

  • $\begingroup$ @ AdamO: Thank you for your answers! I am not sure I totally understand what you are saying, though: 1. "complete paradox" means I should do what exactly? Did I not specify correct models when using the log or inverse link and looking at the residuals was not happy, therefore turned somewhere else? & 3. Are yousaying a large explained Deviance is no indication of a useful model? What if all the residual and diagnostic plots and plotted model curves point at the same direction and imply the high-Dev-explained-models are also better in all other regards? $\endgroup$ Commented Aug 16, 2023 at 8:27
  • $\begingroup$ @user_20201213 there are other posts on the network about how to pick a GLM link. When the regression parameters are any non-trivial value (i.e. not 0), then the actual response is a convolution of whatever distribution you are modeling - in Poisson regression the Y is a convolution of Poissons. You can't just look at a convolved distribution and say, "oh yes these are a bunch of Poissons". You can certainly use a model selection approach, although the validity (as with all multiple testing scenarios) is compromised. $\endgroup$
    – AdamO
    Commented Aug 18, 2023 at 14:56
  • $\begingroup$ @user_20201213 my preference is to robustify GLMs by using robust standard errors when $n>40$ so that the choice of link is less important. I previously did a plethora of simulations to convince myself the inference is usually valid even when the link is wrong (but not very wrong). $\endgroup$
    – AdamO
    Commented Aug 18, 2023 at 14:57

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.