Modelos lineales generalizados con Python y statsmodels

Avatar Tutor | octubre 28, 2018

El modelo lineal generalizado amplía los modelos lineales, de manera que las variables dependientes están relacionadas linealmente con factores y las covariables mediante una determinada función.

Una de las ventajas es que este modelo permite que la variable dependiente tenga una distribución no normal.

Como ejemplo, una empresa de transporte puede utilizar modelos lineales generalizados para ajustar una regresión de Poisson a las frecuencias de rechazo de paso de sus barcos en el canal de Panamá en varios períodos de tiempo. El modelo resultante puede ayudar a determinar cuales son las fechas de mayor rechazo de paso de sus barcos (por atochamiento, reparaciones del canal etc.).

Los modelos lineales generalizados actualmente admiten la estimación utilizando las familias exponenciales de un solo parámetro.

Usando un modelo Gamma

import statsmodels.api as sm
data = sm.datasets.scotland.load()
data.exog = sm.add_constant(data.exog)
gamma_model = sm.GLM(data.endog, data.exog, family=sm.families.Gamma())
gamma_results = gamma_model.fit()
print(gamma_results.summary())

Resultados

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                   32
Model:                            GLM   Df Residuals:                       24
Model Family:                   Gamma   Df Model:                            7
Link Function:          inverse_power   Scale:                       0.0035843
Method:                          IRLS   Log-Likelihood:                -83.017
Date:                Mon, 14 May 2018   Deviance:                     0.087389
Time:                        21:48:07   Pearson chi2:                   0.0860
No. Iterations:                     6   Covariance Type:             nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0178      0.011     -1.548      0.122      -0.040       0.005
x1          4.962e-05   1.62e-05      3.060      0.002    1.78e-05    8.14e-05
x2             0.0020      0.001      3.824      0.000       0.001       0.003
x3         -7.181e-05   2.71e-05     -2.648      0.008      -0.000   -1.87e-05
x4             0.0001   4.06e-05      2.757      0.006    3.23e-05       0.000
x5         -1.468e-07   1.24e-07     -1.187      0.235   -3.89e-07    9.56e-08
x6            -0.0005      0.000     -2.159      0.031      -0.001   -4.78e-05
x7         -2.427e-06   7.46e-07     -3.253      0.001   -3.89e-06   -9.65e-07
==============================================================================

Datos binomiales de respuesta

%matplotlib inline

from __future__ import print_function
import numpy as np
import statsmodels.api as sm
from scipy import stats
from matplotlib import pyplot as plt

data = sm.datasets.star98.load()
data.exog = sm.add_constant(data.exog, prepend=False)

# La variable dependiente es N por 2 (Exito: NABOVE, Fallo: NBELOW):

print(data.endog[:5,:])

Resultado

[[452. 355.]
 [144.  40.]
 [337. 234.]
 [395. 178.]
 [  8.  57.]]

Ajuste y Resumen

glm_binom = sm.GLM(data.endog, data.exog, family=sm.families.Binomial())
res = glm_binom.fit()
print(res.summary())

Resultado

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:           ['y1', 'y2']   No. Observations:                  303
Model:                            GLM   Df Residuals:                      282
Model Family:                Binomial   Df Model:                           20
Link Function:                  logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -2998.6
Date:                Mon, 14 May 2018   Deviance:                       4078.8
Time:                        21:45:26   Pearson chi2:                 4.05e+03
No. Iterations:                     5   Covariance Type:             nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
x1            -0.0168      0.000    -38.749      0.000      -0.018      -0.016
x2             0.0099      0.001     16.505      0.000       0.009       0.011
x3            -0.0187      0.001    -25.182      0.000      -0.020      -0.017
x4            -0.0142      0.000    -32.818      0.000      -0.015      -0.013
x5             0.2545      0.030      8.498      0.000       0.196       0.313
x6             0.2407      0.057      4.212      0.000       0.129       0.353
x7             0.0804      0.014      5.775      0.000       0.053       0.108
x8            -1.9522      0.317     -6.162      0.000      -2.573      -1.331
x9            -0.3341      0.061     -5.453      0.000      -0.454      -0.214
x10           -0.1690      0.033     -5.169      0.000      -0.233      -0.105
x11            0.0049      0.001      3.921      0.000       0.002       0.007
x12           -0.0036      0.000    -15.878      0.000      -0.004      -0.003
x13           -0.0141      0.002     -7.391      0.000      -0.018      -0.010
x14           -0.0040      0.000     -8.450      0.000      -0.005      -0.003
x15           -0.0039      0.001     -4.059      0.000      -0.006      -0.002
x16            0.0917      0.015      6.321      0.000       0.063       0.120
x17            0.0490      0.007      6.574      0.000       0.034       0.064
x18            0.0080      0.001      5.362      0.000       0.005       0.011
x19            0.0002   2.99e-05      7.428      0.000       0.000       0.000
x20           -0.0022      0.000     -6.445      0.000      -0.003      -0.002
const          2.9589      1.547      1.913      0.056      -0.073       5.990
==============================================================================

Valores de Interés

print('Numero de pruebas:',  data.endog[0].sum())
print('Parametros: ', res.params)
print('valores-T: ', res.tvalues)

Resultado

numero de pruebas: 807.0
Parametros:  [-1.68150366e-02  9.92547661e-03 -1.87242148e-02 -1.42385609e-02
  2.54487173e-01  2.40693664e-01  8.04086739e-02 -1.95216050e+00
 -3.34086475e-01 -1.69022168e-01  4.91670212e-03 -3.57996435e-03
 -1.40765648e-02 -4.00499176e-03 -3.90639579e-03  9.17143006e-02
  4.89898381e-02  8.04073890e-03  2.22009503e-04 -2.24924861e-03
  2.95887793e+00]
valores-T:  [-38.74908321  16.50473627 -25.1821894  -32.81791308   8.49827113
   4.21247925   5.7749976   -6.16191078  -5.45321673  -5.16865445
   3.92119964 -15.87825999  -7.39093058  -8.44963886  -4.05916246
   6.3210987    6.57434662   5.36229044   7.42806363  -6.44513698
   1.91301155]

Primeras diferencias: mantenemos todas las variables explicativas constantes en sus medios y manipulamos el porcentaje de hogares de bajos ingresos para evaluar su impacto en las variables de respuesta:

means = data.exog.mean(axis=0)
means25 = means.copy()
means25[0] = stats.scoreatpercentile(data.exog[:,0], 25)
means75 = means.copy()
means75[0] = lowinc_75per = stats.scoreatpercentile(data.exog[:,0], 75)
resp_25 = res.predict(means25)
resp_75 = res.predict(means75)
diff = resp_75 - resp_25

La primera diferencia intercuartil para el porcentaje de hogares de bajos ingresos en un distrito escolar es:

print("%2.4f%%" % (diff*100))

Resultado

-11.8753%

Written by Tutor

02:18

Related View More