dgamma and dgamma2 producing different results

[jebyrnes]

So, I've been working with some Gamma GLMs lately and, while I really appreciate the dgamma2 density for making life easy, I find that it produces different results than dgamma (when dgamma is used properly) and am trying to struggle with why and/or how to create a model that will produce equivalent results. Let's start with an example data set from Witman et al. 2016 Ecological Monographs (happy to provide ref) that is looking at the duration of ecological experiments over the recent history of the discipline.

library(rethinking)
#load gamma_sample
source("https://gist.githubusercontent.com/jebyrnes/27fdece3e28dbcec7684546caddf68e8/raw/c83cc4d6dfcffd74a533a788356a70eaa3023b5e/gamma_sample.R")

We can model this with a Gamma distribution (durations!) and a log link.

g_glm <- glm(Duration ~ Year.c, data = gamma_sample, family = Gamma(link = log))
coef(g_glm)

which leads to

 (Intercept)       Year.c 
 5.798676770 -0.002591553 

OK - now, with dgamma2

gamma2_mod <- alist(
  #likelihood
  Duration ~ dgamma2(mu, scale),
  
  log(mu) <- a + b*Year.c,
  
  #priors
  a ~ dnorm(0,100),
  b ~ dnorm(0,1),
  scale ~ dexp(2)
  
)

gamma2_fit <- map(gamma2_mod, gamma_sample)

This leads to somewhat different coefficient estimates (and is quite sensitive to the scale prior) (which is where I suspect the trouble might be coming in)

> coef(gamma2_fit)
            a             b         scale 
  5.186479426  -0.001695557 153.327060740 

Now, with the dgamma formulation, you have to do a little but more work in model construction, making rate = log_mu/shape

gamma_mod <- alist(
  #likelihood
  #using dgamma(shape, rate) formulation
  Duration ~ dgamma(shape, rate),
  
  rate <- shape/exp(log_mu),
  log_mu <- a + b*Year.c,
  
  #priors
  a ~ dnorm(0,100),
  b ~ dnorm(0,1),
  shape ~ dexp(2)
  
)

gamma_fit <- map(gamma_mod, gamma_sample)

This works pretty well, though, as compared to the original.

> coef(g_glm)
 (Intercept)       Year.c 
 5.798676770 -0.002591553 
> coef(gamma_fit)
           a            b        shape 
 5.798670099 -0.002591604  0.564214438 

So... what's going on here?

Note, I've found the same problem with the dzagamma2. Using some simulated data

source("https://gist.githubusercontent.com/jebyrnes/195b32f7d5d9c44f0cd3bb5ec17ac2f1/raw/318a60383e9ad6cf28b39931e768d9437e976950/sim_catches.R")

I built the corresponding dzaggama2 model and then wrote

dzagamma <- function (x, prob, shape, rate, log = FALSE) 
{
  K <- as.data.frame(cbind(x = x, prob = prob, shape = shape, rate = rate))
  llg <- dgamma(x, shape = shape, rate = rate, log = TRUE)
  ll <- ifelse(K$x == 0, log(K$prob), log(1 - K$prob) + llg)
  if (log == FALSE) 
    ll <- exp(ll)
  ll
}

dzagamma2 was close but no cigar to recovering the coefficients that I used to generate the data (although it got the probability right). dzagamma was right on the money, as it were.

Advice or thoughts?

[jswesner]

+1. I had a similar issue with dgamma2 and dzagamma2.

[rmcelreath]

Different parameterizations impose different priors. There doesn't appear to be any bug here.

[jebyrnes]

OK, I get this - but, so, for your parameterization of a gamma, what is the correct way to create a flat prior, if not an exponential? Would be helpful to show students this.