Statsmodels: how to run and interpret a Gamma regression?

Question

I have an endogenous variable that is continuous and non-negative. From what I can gather, a Poisson regression is not appropriate because the values of the response variable are not natural number, therefore I should use a Gamma regression.

This is the code that I am using:

gamma_model = sm.GLM(cpl.total, sm.add_constant(cpl.old), family=sm.families.Gamma())
gamma_results = gamma_model.fit()
gamma_results.summary()

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                  total   No. Observations:                  752
Model:                            GLM   Df Residuals:                      750
Model Family:                   Gamma   Df Model:                            1
Link Function:          inverse_power   Scale:                         0.80069
Method:                          IRLS   Log-Likelihood:                -5349.3
Date:                Fri, 09 Oct 2020   Deviance:                       554.79
Time:                        15:39:03   Pearson chi2:                     601.
No. Iterations:                     8                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.0061      0.000     19.275      0.000       0.005       0.007
old        -1.506e-07      1e-08    -15.028      0.000    -1.7e-07   -1.31e-07
==============================================================================

There are a couple of things that tell me this isn't a reliable result:

• I get an error about the link function: DomainWarning: The inverse_power link function does not respect the domain of the Gamma family.
• The coefficient of interest (old) is negative. I have reason to believe it should be positive.
  • A plot seems to indicate that as old increases, so does total.
    • 
  • A standard linear regression gives a positive coefficient

    • linear_model = ols('total ~ old', data=cpl)
      linear_results = linear_model.fit()
      linear_results.summary() 
      
                                  OLS Regression Results                            
      ==============================================================================
      Dep. Variable:                  total   R-squared:                       0.201
      Model:                            OLS   Adj. R-squared:                  0.200
      Method:                 Least Squares   F-statistic:                     188.4
      Date:                Fri, 09 Oct 2020   Prob (F-statistic):           2.09e-38
      Time:                        16:02:55   Log-Likelihood:                -5780.3
      No. Observations:                 752   AIC:                         1.156e+04
      Df Residuals:                     750   BIC:                         1.157e+04
      Df Model:                           1                                         
      Covariance Type:            nonrobust                                         
      ==============================================================================
                       coef    std err          t      P>|t|      [0.025      0.975]
      ------------------------------------------------------------------------------
      Intercept   -257.6356     61.257     -4.206      0.000    -377.890    -137.381
      old            0.0330      0.002     13.724      0.000       0.028       0.038
      ==============================================================================
      Omnibus:                      305.139   Durbin-Watson:                   1.940
      Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1367.182
      Skew:                           1.846   Prob(JB):                    1.32e-297
      Kurtosis:                       8.478   Cond. No.                     8.10e+04
      ==============================================================================
      
      Warnings:
      [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
      [2] The condition number is large, 8.1e+04. This might indicate that there are
      strong multicollinearity or other numerical problems.
      

I'd like to know whether I am using the correct specification. Perhaps I'm using the wrong link function. Ultimately, if the model is correct, I'd also like to know how to interpret the coefficient of a gamma regression.

Answer 1

Whether to use Poisson or Gamma regression shouldn't depend on whether the data are integer-valued, that is a minor consideration. In the quasi-GLM framework you can use Poisson regression with non-integer data. The key difference between Gamma and Poisson regression is how the mean/variance relationship is encoded in the model. The Poisson approach models the variance as being proportional to the mean, the Gamma approach models the standard deviation as being proportional to the mean. This is a major difference. Also note that the Tweedie GLM interpolates between Poisson and Gamma, in that the variance is proportional to the mean raised to power p; so p=1 correponds to Poisson and p=2 corresponds to Gamma but any p in [1, 2] is possible (p outside this range is also possible to but harder to fit).

The code below shows how to fit Poisson and Gamma GLMs to simulated data in statsmodels. Note that the mean structure parameters are estimated well even if the family is not correct. This is theoretically expected due to robustness of quasi-GLMs to variance misspecification (although p-values and CI's are wrong unless using robust inference). Also, I am using Gamma with the non-canonical log link here, which makes it directly comparable to Poisson (and is probably the best way to go in practice too).

import numpy as np
import statsmodels.api as sm
import pandas as pd                                                                                                                    

x1 = np.random.normal(size=5000)
x2 = np.random.normal(size=5000)
lp = x1/2 - x2
ev = np.exp(lp)
fml = "y ~ x1 + x2"

for k in range(2):
                                                                                                                                         
    if k == 0:        
        y = np.random.poisson(ev)                                                                                                             
        df = pd.DataFrame({"y": y, "x1": x1, "x2": x2})                                                                                  
    else:                                                                                             
        # Coefficient of variation for Gamma, scale is cv^2                                                                               
        cv = 3                                                                                                                          
        y = np.random.gamma(1/cv**2, ev*cv**2)                                                                                        
    df = pd.DataFrame({"y": y, "x1": x1, "x2": x2})                                                                                                                                                                                 
    mp = sm.GLM.from_formula(fml, family=sm.families.Poisson(), data=df)
    rp = mp.fit(scale="X2")                                                                                                               
    print(rp.summary())                                                                                                                   
    print(rp.scale)                                                                                                                  
    mg = sm.GLM.from_formula(fml,  
             family=sm.families.Gamma(link=sm.families.links.log()), data=df)
    rg = mg.fit(scale="X2")
    print(rg.summary())
    print(rg.scale)