5
$\begingroup$

I understand that Gamma distribution generates only positive values. And this is reflected in R gamma family glm function which does not run when the dependent variable contains zeros or negative values. However, in Stata, the glm procedure for Gamma family runs and gives output with estimated coefficents even when there are zeros in the dependent variable. 

glm cost salary [pw=wg], family(gamma) link(log) exposure(exposure_time) vce(cluster unique_key)

Iteration 0:   log pseudolikelihood = -1602.3771  
Iteration 1:   log pseudolikelihood = -1579.4773  
Iteration 2:   log pseudolikelihood = -1579.3417  
Iteration 3:   log pseudolikelihood = -1579.3417  

Generalized linear models                         No. of obs      =     22,624
Optimization     : ML                             Residual df     =     22,622
                                                  Scale parameter =   .0233436
Deviance         =   195.440412                   (1/df) Deviance =    .009039
Pearson          =  504.7344288                   (1/df) Pearson  =   .0233436

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

                                                  AIC             =    .146258
Log pseudolikelihood = -1579.341674               BIC             =  -215625.8

                        (Std. Err. adjusted for 16,223 clusters in unique_key)
------------------------------------------------------------------------------
             |               Robust
        cost |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      salary |   .0000121   2.87e-06     4.21   0.000     6.45e-06    .0000177
       _cons |   4.112714   .1696957    24.24   0.000     3.780116    4.445311
ln(exposu~e) |          1  (exposure)
------------------------------------------------------------------------------
count if cost==0
  9,265

It does not look like Stata drops the zero values by default.

glm cost salary [pw=wg] if cost!=0, family(gamma) link(log) exposure(exposure_time) vce(cluster unique_key)

Iteration 0:   log pseudolikelihood = -1299.6486  
Iteration 1:   log pseudolikelihood = -1286.4665  
Iteration 2:   log pseudolikelihood = -1286.3631  
Iteration 3:   log pseudolikelihood = -1286.3631  

Generalized linear models                         No. of obs      =     13,359
Optimization     : ML                             Residual df     =     13,357
                                                  Scale parameter =   .0232079
Deviance         =  188.1527815                   (1/df) Deviance =   .0152264
Pearson          =   286.780419                   (1/df) Pearson  =   .0232079

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

                                                  AIC             =   .2084899
Log pseudolikelihood = -1286.363084               BIC             =  -116241.2

                         (Std. Err. adjusted for 8,402 clusters in unique_key)
------------------------------------------------------------------------------
             |               Robust
        cost |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      salary |     .00001   2.48e-06     4.03   0.000     5.14e-06    .0000149
       _cons |    4.44726   .1527983    29.11   0.000     4.147781    4.746739
ln(exposu~e) |          1  (exposure)
------------------------------------------------------------------------------

Anyone has any idea why this is so? How exactly was the model generated in former case?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

The Stata output tells you that it is fitting a pseudo-likelihood rather than a true likelihood. Pseudo-likelihood is not well documented in Stata, but I think it is equivalent here to what is usually called quasi-likelihood in glm theory. That would mean that it is assuming a mean-variance relationship but not assuming any specific response distribution. For quasi-Gamma glms, the assumed variance function is $V(\mu)=\mu^2$.

You can fit the same sort of quasi-Gamma glm in R using:

fit <- glm(y ~ x, family = quasi(link="log", variance="mu^2"))

This code will not complain if y contains exact zeros, but will otherwise give identical results to a gamma glm in R.

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ This works with zeros in y, but produces errors for values very close to zero. I think the problem is how glm chooses starting values, so it might be safest to manually specify start. $\endgroup$
    – Paul Harrison
    Commented Mar 28, 2020 at 1:25

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.