# Why do I get such a large divergence between Gaussian and Gamma family GLMs?

From what I understand the Gamma distribution is a good choice for positive right skew data. My data is a magnitude quantity with heavy outliers as you can see here from a histogram of the pooled set:

The data I am analyzing comes from a 2X2 factorial design (Region and Cond) with metric response variable (called here y). In R, the default (Gaussian) GLM fit for the multiplicative model comes in about as I'd expect given that the pooled mean of the data is around 990.

Call: glm(formula = y ~ Region * Cond, data = myData)

Coefficients:
(Intercept)         RegionPOST             CondON  RegionPOST:CondON
920.70             127.14              23.42             -19.90


And I get similar results using a Bayesian ANOVA with homogeneous variance. However I am uncomfortable with using these tests, because my data is so far from normal. I thought that I could change my models to handle Gamma distributions instead by calling the same command with the "family" argument set to Gamma but I get:

Call:  glm(formula = y ~ Region * Cond, family = Gamma, data = myData)

Coefficients:
(Intercept)         RegionPOST             CondON  RegionPOST:CondON
1.076e-03         -1.445e-04         -1.251e-05          4.766e-05


These values bear no relation to my data that I can see, neither to means nor to the shape and rate of any likely Gamma distributions.

Why does changing the family of the GLM to Gamma change the results so dramatically, and how can I interpret these results in the context of the original data?

You need to take into account the link function that is being employed. Checking the documentation, we see that the default link (your code did not change this) is the inverse. If you take the inverse of your intercept, you will see that it is approximately equal to the estimate of the intercept in the normal model (it shouldn't necessarily be expected that they be identical):

1/1.076e-03
# [1] 929.368


R uses the inverse by default because it is the canonical link (cf., Difference between 'link function' and 'canonical link function' for GLM). However, there are various ways in which people don't like the inverse for Gamma models (e.g., they are less interpretable); many people prefer to use a log link instead (that's what I would use). If you chose to do that, you would use:

glm(y~Region*Cond, family=Gamma(link="log"), data=myData)


The predicted values on the linear scale would then be exponentiated to get the values on the original scale. For more on this, see: How to interpret parameters in GLM with family=Gamma.

• Thanks and one follow-up question: You can see that in the Gaussian GLM my contrasts are much smaller than my intercept. However in the gamma result, they are also smaller for example the inverse of CondON is around 83,000 and that is far outside the range of my data. So is your answer only applicable to the intercept? Jun 30, 2017 at 16:40
• @barnhillec, w/ an inverse link, the interpretation of the betas will be complicated. (Again, I think the log is more intuitive, which is why I would use that.) The most foolproof way to work w/ what you have is to solve for the predicted value for whatever covariates you are interested in, then take the inverse to get the number on the original scale. For the intercept, that's just taking the inverse of the estimate, for the rest, you need to do some multiplication & addition, that's why I used the intercept here. Jun 30, 2017 at 16:44
• @barnhillec 1. Inverse links are sometimes readily interpretable (e.g. if your data are times the inverse will be in terms of rates/speeds); 2, on top of the effect of the link function the relative sizes of coefficients depends on the scales of the variables; a variable with large values will tend to have smaller coefficients. 3. You don't invert the coefficients when trying to consider the original scale effects of a variable in an inverse link model. That doesn't work. Jul 2, 2017 at 3:35