How do I obtain the parameters for rgamma() from a specified gamma model

[JD Young]

I would like to obtains input parameters for the function gamma(a,b), where a is the gamma shape parameter and b is the scale parameter. I would like to obtain these parameters from a glm model, e.g.,

glm y1 x1 x2 x3, family(gamma) link(log)

It seems the scale parameter is reported in the standard output if I am understanding correctly. But I am not sure how to go about calculated the shape parameter. Any help is greatly appreciated.

[Jared Greathouse]

Is it in the help file?

[JD Young]

Lord all mighty...

[Jared Greathouse]

Yeah the help usually lists what's kept and what isn't! Super helpful when you need it!

[JD Young]

Yes, I am familiar. It contains the scale parameter. However, I am not convinced it is that straight forward. I was able to find these equation from another source that may be of some use. shape = (m^2)/v, scale = v/m, where m is the median and v is the variance. It results in the attached graph that compares the raw data to the modeled (in my case). Any help on verifying these equations would be great. Seems to make sense visually... but inquiring to make sure

[Leonardo Guizzetti]

Just some playing around here. You'll need to be careful as well with how the gamma distribution is parameterized. The two main forms are described on the Wiki article. Stata's -rgamma()- uses the first notation from the Wiki article, where Stata's a and b correspond to k and theta in the Wiki notation. However, the -glm- model reparameterizes the model in terms of the mean, and so it's yet different.

Below I show two ways to estimate the gamma distribution paramters. One based on the glm, the other based on Nick Cox's -gammafit- (SSC).

. set seed 17
. set obs 100000
. gen y = rgamma(3, .5) // rgamma(a, b)

. summ y, detail // show mean and variance of y to check parameterization of rgamma()

                              y
-------------------------------------------------------------
      Percentiles      Smallest
 1%     .2151081       .0064835
 5%     .4091736       .0164218
10%     .5508196       .0374375       Obs             100,000
25%     .8646331       .0379513       Sum of wgt.     100,000

50%     1.340371                      Mean           1.500639
                        Largest       Std. dev.      .8641069
75%     1.962899        7.33879
90%     2.659449       7.341051       Variance       .7466807
95%        3.145       7.628283       Skewness       1.132031
99%      4.18282       9.040347       Kurtosis       4.890583

. * estimate based on regression
. glm y, fam(gamma) link(id) nolog

Generalized linear models                         Number of obs   =    100,000
Optimization     : ML                             Residual df     =     99,999
                                                  Scale parameter =   .3315757
Deviance         =  35162.47525                   (1/df) Deviance =   .3516283
Pearson          =  33157.24187                   (1/df) Pearson  =   .3315757

Variance function: V(u) = u^2                     [Gamma]
Link function    : g(u) = u                       [Identity]

                                                  AIC             =   2.811801
Log likelihood   = -140589.0698                   BIC             =   -1116119

------------------------------------------------------------------------------
             |                 OIM
           y | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   1.500639   .0027326   549.17   0.000     1.495283    1.505994
------------------------------------------------------------------------------

. di "mean = " _b[_cons]
mean = 1.5006385

. di "scale, a = " 1/e(phi)
scale, a = 3.0159022

. di "shape, b = " _b[_cons] * e(phi)
shape, b = .49757532

. * ML fit of gamma distribution
. gammafit y

< output omittted >

ML fit of two-parameter gamma distribution        Number of obs   =     100000
                                                  Wald chi2(0)    =          .
Log likelihood = -115482.58                       Prob > chi2     =          .

------------------------------------------------------------------------------
           y | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
alpha        |
       _cons |   3.000254   .0127422   235.46   0.000      2.97528    3.025228
-------------+----------------------------------------------------------------
beta         |
       _cons |   .5001704   .0023122   216.32   0.000     .4956386    .5047023
------------------------------------------------------------------------------

[JD Young]

Thanks Leonardo! Makes perfect sense. I see where I was going wrong before.