# Using R for GLM with Gamma distribution

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.

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.

• 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)? – Jon Claus May 9 '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 -Reinstate Monica May 9 '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 9 '13 at 18:10
• 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 -Reinstate Monica May 9 '13 at 18:31
• 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 9 '13 at 19:11

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 = T, dispersion = 1/myshape$alpha) summary(glmGamma, dispersion = 1/myshape$alpha)