Gamma shape and scale seemingly swapped - where am I going wrong?

[JockLawrie]

Hi there,

I am using a Gamma GLM with LogLink. I am interested in estimating the coefficients as well as the scale parameter so that, for a given input vector x, I can reconstruct the whole Gamma distribution and not just the mean.

The example below illustrates my current attempt. I get reasonable estimates, but the shape and scale parameters seem to be swapped (see the last line). I've been banging my head on this for a while - any ideas where I'm going wrong?

Cheers

using GLM

# True parameter values: shape=5, scale=2. Aim is to estimate these from data
d = Gamma(5.0, 2.0)

# Data
y = rand(d, 1000)
X = fill(1.0, 1000, 1)

# Model
dist  = Gamma()
lnk   = LogLink()
model = glm(y, X, dist, lnk)

# Estimate mean and variance at x=1
b  = coef(model)
lp = b[1]  # linear_predictor = intercept
mu = GLM.linkinv(lnk, lp)
v  = GLM.dispersion(model) * GLM.glmvar(dist, mu)

# From the mean and variance, extract shape and scale for the Gamma distribution at x=1
dist_params(d::Gamma, mu, v) = (shape = mu*mu/v, scale = v/mu)

dist_params(dist, mu, v)  # Approx (shape=2, scale=5) instead of (shape=5, scale=2)

[andreasnoack]

If v the same as McCullagh and Nelder's ν then it's the scale parameter of the Gamma and it's reciprocal to the dispersion parameter, i.e.

julia> v = inv(GLM.dispersion(model))
2.349964824015671

and therefore

julia> (shape = mu/v, scale = v)
(shape = 4.331704011939433, scale = 2.349964824015671)

[JockLawrie]

Thanks Andreas, most helpful.