22
$\begingroup$

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.

Improve this question
$\endgroup$
0

2 Answers 2

Reset to default
15
$\begingroup$

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.

Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ Thanks. In R, is there a way to specify what $ \phi = K $ is? From this link, 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)? $\endgroup$
    – Jon Claus
    Commented May 9, 2013 at 15:58
  • $\begingroup$ 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? $\endgroup$
    – Glen_b
    Commented May 9, 2013 at 17:48
  • $\begingroup$ 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. $\endgroup$
    – Jon Claus
    Commented May 9, 2013 at 18:10
  • 2
    $\begingroup$ 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? $\endgroup$
    – Glen_b
    Commented May 9, 2013 at 18:31
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jon Claus
    Commented May 9, 2013 at 19:11
18
$\begingroup$

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)
Improve this answer
$\endgroup$
1
  • 3
    $\begingroup$ 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" $\endgroup$
    – Maverick Meerkat
    Commented Sep 6, 2020 at 18:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.