Using R for GLM with Gamma distribution

Question

I currently have a problem understanding the syntax for R for fitting a GLM using the Gamma distribution.

I have a set of data, where each row contains 3 co-variates ($X_1, X_2, X_3$), a response variable ($Y$), and a shape parameter ($K$). I want to model the scale of the Gamma distribution as a linear function of the 3 covariates, but I don't understand how to set the shape of the distribution to $K$ for each row of data.

A situation that I think is analogous is that for a binomial distribution, the GLM requires that the number of trials ($N$) is known for each data entry.

score 18 · Answer 1 · edited Jun 09 '21 at 03:17

18

I used the gamma.shape function of MASS package as described by Balajari (2013) in order to estimate the shape parameter afterwards and then adjust coefficients estimations and predictions in the GLM. I advised you to read the lecture as it is, in my opinion, very clear and interesting concerning the use of gamma distribution in GLMs.

    glmGamma <- glm(response ~ x1, family = Gamma(link = 
                    "identity")
    library(MASS)
    myshape <- gamma.shape(glmGamma)
    gampred <- predict(glmGamma , type = "response", se = TRUE, 
                      dispersion = 1/myshape$alpha) 
    summary(glmGamma, dispersion = 1/myshape$alpha)

edited Jun 09 '21 at 03:17

kjetil b halvorsen

63,378
26
142
467

answered Jul 24 '14 at 11:27

Xochitl C.

381
3
9

3

This is preferred because "the dispersion offered by the default GLM summary command does not take into account the special information about the dispersion that can be calculated by using the Gamma distribution" – Maverick Meerkat Sep 06 '20 at 18:39

Glen_b · Accepted Answer · 2013-05-09T22:42:54.013

14

The usual gamma GLM contains the assumption that the shape parameter is constant, in the same way that the normal linear model assumes constant variance.

In GLM parlance the dispersion parameter, $\phi$ in $\text{Var}(Y_i)=\phi\text{V}(\mu_i)$ is normally constant.

More generally, you have $a(\phi)$, but that doesn't help.

It might perhaps be possible to use a weighted Gamma GLM to incorporate this effect of a specified shape parameter, but I haven't investigated this possibility yet (if it works it is probably the easiest way to do it, but I am not at all sure that it will).

If you had a double GLM you could estimate that parameter as a function of covariates... and if the double glm software let you specify an offset in the variance term you could do this. It looks like the function dglm in the package dglm let you specify an offset. I don't know if it will let you specify a variance model like (say) ~ offset(<something>) + 0 though.

Another alternative would be to maximize the likelihood directly.

> y <- rgamma(100,10,.1)

> summary(glm(y~1,family=Gamma))

Call:
glm(formula = y ~ 1, family = Gamma)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.93768  -0.25371  -0.05188   0.16078   0.81347  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.0103660  0.0003486   29.74   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.1130783) 

    Null deviance: 11.223  on 99  degrees of freedom
Residual deviance: 11.223  on 99  degrees of freedom
AIC: 973.56

Number of Fisher Scoring iterations: 5

The line where it says:

   (Dispersion parameter for Gamma family taken to be 0.1130783)

is the one you want.

That $\hat\phi$ is related to the shape parameter of the Gamma.

edited May 09 '13 at 22:42

answered May 09 '13 at 07:45

Glen_b

257,508
32
553
939

1

Thanks. In R, is there a way to specify what $ \phi = K $ is? From this [link](https://stat.ethz.ch/pipermail/r-help/2001-February/011366.html), it seems that I don't have to decide on a given $ K $ until I print the results. Am I correct in saying that if there is a fixed $ K $, then it does not affect the result for $ \beta $, the coefficient vector? If so, how do I decide on the best $ K $ to fit the data manually (without using R)? – Jon Claus May 09 '13 at 15:58
If there is a fixed shape parameter for the Gamma, it does not affect the estimate of $\mu$, and hence not the coefficient vector either. You can compute *an* estimate from the GLM output, but it's not maximum likelihood. If I wanted to identify the shape parameter, I'd use the relevant functions in the package `MASS`. Why is it important to avoid using R, and why would you try to do it manually rather than use a computer? – Glen_b May 09 '13 at 17:48
I misspoke. By manually, I meant I wanted a decently simple algorithm I could implement myself outside of R. Also, when I try testing `glm(V4 ~ V3 + V2 + V1, family=Gamma)`, where $ V_1, V_2, V_3 $ are predictor co-variates and $ V_4 $ is the response, it improperly determines $ \boldsymbol \beta $, the coefficient vector. I know it improperly determines it because I generated sample data with known coefficients to determine the scale and a constant shape of 5. – Jon Claus May 09 '13 at 18:10
2

Well you can implement anything outside R that could be implemented within it; you could maximize the likelihood, for example, or you could use the estimate based off $\hat\phi$. Can you explain in more detail what you mean by "improper" here? – Glen_b May 09 '13 at 18:31
1

For the purpose of testing my own code, I generated a set of data with 10,000 tuples. To generate it, I fixed $ \boldsymbol \beta $, generated sample $ {\bf V} $, calculated $ \theta = (\boldsymbol \beta^T {\bf V})^{-1} $ (the scale parameter with the inverse link function), and generated a random variable from the distribution $ Y \sim \text{Gamma}(5, \theta) $. When I run R on the data set, its predicted $ \hat{\boldsymbol \beta} $ is nowhere near $ \boldsymbol \beta $. When I have done this for other distributions, R's prediction has been nearly exactly correctly. – Jon Claus May 09 '13 at 19:11
Looks to me like you're confused about what the glm is modelling there ... because that's NOT how you do it. Note that a GLM is a model for the mean, not the $\theta$ parameter of a Gamma. If you create a model for $\theta$ in a Gamma$(5,\theta)$... you don't have the mean of $Y$. Your estimates of the parameters in $\beta$ will be out by a factor of that same size as the factor by which your mean is out. – Glen_b May 09 '13 at 19:43
My mistake. I actually did have that $ \mu = (\boldsymbol \beta^T {\bf V})^{-1} $ and $ \theta = \frac{\mu}{K} $. As stated above, $ K = 5 $ for all tuples and hence $ Y \sim \text{Gamma}(5 , \frac{\mu}{5}) $. However, I had originally done the error you pointed out and failed to regenerate the data after fixing it. It works correctly now. Also, you mentioned I could estimate $\phi $ off of $ \hat{\phi} $. How do I compute the latter? Basically, I am just looking for a reference that elaborates on how to predict $ \phi $ more easily than directly solving the likelihood function. – Jon Claus May 09 '13 at 20:55
$\hat\phi$ is produced in the glm output. Let me add an example to my answer. – Glen_b May 09 '13 at 22:40
Alright, what I had meant was that I wanted a very basic algorithm for how R arrives at its estimate for the dispersion parameter. I decided to go with $ \phi = \frac{1}{n - p} \sum_{i = 1}^n \frac{\left(y_i - \mu_i\right)^2}{V\left(\mu_i\right)} $, where $ p = \log_\mu (V(\mu)) $. – Jon Claus May 10 '13 at 18:11
Aside from the value of $p$ (unless I missed something), that looks like the usual estimate. – Glen_b May 15 '13 at 04:53

Using R for GLM with Gamma distribution

2 Answers2

Linked