4
$\begingroup$

What are the assumptions when doing hypothesis testing using a Gamma GLM or GLMM? Are the residuals suppose to be normally distributed and is heteroscedasticity a concern like the Gaussian (normal) distribution? Do you test the assumptions in the same fashion Levene Tests of individual fixed effects, residual and qqplots, and Shapiro-Wilks test?

Discussion of when to use Gamma Glm and Glmm is found here:

When to use gamma GLMs?

However, I found little to no discussion of the assumptions and how to test them?

Edit 1: Requested Example residual plots- are any of these a concern. Levene's test, comes back significant.

Residuals are plotted using this code. Time (Obs) is the ordinal factor levels from the data

modelresiduals = residuals(Gama_model)
plot(modelresiduals)
plot(y = modelresiduals, x = Model_data2$Obs)

1. All resids-

enter image description here

2. By time

enter image description here

Edit 2 Based on comment suggestions- Patterns look good

ypred = predict(gama_model)
res = residuals(gama_model, type = 'pearson')
plot(ypred,res)
hist(res)
plot(res)
plot(y = res, x = Model_data$Obs)

3. Pearson histograms- based on answer code enter image description here

4. Pearson vs pred (link scale)- based on answer code enter image description here

5. All Pearson residual- based on answer code enter image description here

6. Pearson residual by Time- based on answer code enter image description here

$\endgroup$

3 Answers 3

3
$\begingroup$

Are the residuals suppose to be normally distributed...

No. The distribution of $y\vert x$ is assumed to be gamma, and I don't think that $X - \mathbb{E}(X)$ is gamma for $X\sim \text{Gamma}(k,\theta)$. The deviance residuals, on the other hand, should be normal.

is heteroscedasticity a concern like the Gaussian (normal) distribution?

Heteroscedasticity should be observed since the variance changes with the mean.

Here is an example of Gamma regression in R.

library(tidyverse)

#Simulate
set.seed(0)
N <- 250
x <- runif(N, -1, 1)
a <- 0.5
b <- 1.2
y_true <- exp(a + b * x)
shape <- 2
y <- rgamma(N, scale = y_true/shape, shape = shape)

#Model
plot(x,y)
model <- glm(y ~ x, family = Gamma(link = "log"))
summary(model)

#Deviance residuals
ypred = predict(model)
res = residuals(model, type = 'deviance')
plot(ypred,res)
hist(res)

#Deviance GOF

deviance = model$deviance
p.value = pchisq(deviance, df = model$df.residual, lower.tail = F)
# Fail to reject; looks like a good fit.

For more on GLMs, I would recommend Extending the Linear Model by Faraway.

$\endgroup$
14
  • $\begingroup$ so am I correct in saying that if my data is hetersocedastic in time (ie. Cone shaped), if I want to do hypothesis tests with pairwise comparisons and I'm doing a gamma glmm, it's not a concern like it would be in a linear mixed model(Gaussian)? $\endgroup$ Jan 31, 2019 at 1:23
  • $\begingroup$ That should be correct. The variance should change with the mean for a Gamma Glm $\endgroup$ Jan 31, 2019 at 1:37
  • 1
    $\begingroup$ "Heteroscedasticity should be observed since the variance changes with the mean." This is true for the modeled response variable but not for the residuals vs. the predicted values as part of the model validation process? Your plot(ypred,res) looks good and there is no residual pattern present. $\endgroup$
    – Stefan
    Jan 31, 2019 at 2:44
  • 1
    $\begingroup$ @Stefan, yes, thank you for clarifying. The response should be heteroskedastic, but the deviance residuals should not. $\endgroup$ Jan 31, 2019 at 2:45
  • 1
    $\begingroup$ I agree with Demetri. You shouldn't use those tests to assess whether assumptions are met. I mostly use Pearson residuals and not Deviance residuals for model validation - although both (Pearson and Deviance residuals) are defined in Faraway's book to do diagnostic plots. This paper might be helpful too (see Step 7): besjournals.onlinelibrary.wiley.com/doi/full/10.1111/… $\endgroup$
    – Stefan
    Jan 31, 2019 at 3:10
4
$\begingroup$

A nice approach for checking the fit of your assumed model to the data, accounting for features, such as, over-dispersion, non-normality, zero-inflation is the simulated scaled residuals provided by the DHARMa package. If your assumed model is correct, these residuals should have a uniform distribution. You can find more details on the procedure they are defined and used in the vignette of the package.

As an example, using the simulated example above, I compare below the fit of the Gamma model to the fit of the wrong normal model:

set.seed(0)
N <- 250
x <- runif(N, -1, 1)
a <- 0.5
b <- 1.2
y_true <- exp(a + b * x)
shape <- 2
y <- rgamma(N, scale = y_true / shape, shape = shape)

gamma_model <- glm(y ~ x, family = Gamma(link = "log"))
normal_model <- glm(y ~ x, family = gaussian())

library("DHARMa")
check_gamma_model <- simulateResiduals(fittedModel = gamma_model, n = 500)
plot(check_gamma_model)

check_normal_model <- simulateResiduals(fittedModel = normal_model, n = 500)
plot(check_normal_model)

$\endgroup$
2
  • 1
    $\begingroup$ note if using glmmtmb for mixed models...random effects can make residual vs predicted line plots go haywire....despite good fit and lack of overdispersion..thanks for addition $\endgroup$ Jan 31, 2019 at 19:59
  • 2
    $\begingroup$ In my GLMMadaptive package (drizopoulos.github.io/GLMMadaptive) that fits similar models as the glmmTMB but using the adaptive Gaussian quadrature instead of the Laplace approximation it works. For an example you may check: drizopoulos.github.io/GLMMadaptive/articles/… $\endgroup$ Jan 31, 2019 at 20:01
1
$\begingroup$

For checking the deviance goodness-of-fit statistic for a Generalized Linear Model like a Gamma Regression Model I recommend to follow the instructions on slide 18/19 of Technical University Munich respectively on slide 104 of Technical University Graz.

That is, if dispersion parameter disp can not be assumed to be '1', the residual deviance test, i.e. calculating the p-value, should be applied on the scaled deviance, since the scaled statistic follows the chi square distribution (as outlined in the cited slides), e.g.:

p.value = pchisq(myModel$deviance/summary(myModel$disp), df = myModel$df.residual, lower.tail = F)

$\endgroup$
1
  • $\begingroup$ Actually you should use summary(myModel)$disp $\endgroup$
    – igorkf
    Mar 21, 2023 at 23:01

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.