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

Question

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?

Answer 1

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.