top of page

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida

you are kindly asked to include the complete citation if you used this material in a publication

â€‹

â€‹

Code 5.2 GLM logistic regression in Python

==================================================

import numpy as np
import statsmodels.api as sm

â€‹

# Data
x = np.array([13, 10, 15, 9, 18, 22, 29, 13, 17, 11, 27, 21, 16, 14, 18, 8])
y = np.array([1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0])

â€‹

X = np.transpose(x)

# Fit
mu = (y + 0.5) / 2                                                        # initialize mu
eta = np.log(mu/(1-mu))                                             # initialize eta with the Bernoulli link

for i in range(8):
w = mu * (1 - mu);                                                 # variance function
z = eta + (y - mu)/(mu * (1 - mu))                         # working response
mod = sm.WLS(z, X, weights=w).fit()                 # weigthed regression
eta = mod.predict()                                                # linear predictor
mu = 1/(1 + np.exp(-eta))                                      # fitted value
print(mod.params)                                                # print iteration log

â€‹

# Output
print(mod.summary())

â€‹

# Write data as dictionary
mydata = {}
mydata['x'] = x
mydata['y'] = y

â€‹

# fit using glm package

import statsmodels.formula.api as smf

mylogit = smf.glm(formula='y ~ x', data=mydata, family=sm.families.Binomial())
res = mylogit.fit()
print(res.summary())

==================================================

Output on screen:

â€‹

WLS Regression Results
==============================================================================
Dep. Variable:                                                            y   R-squared:                                                            0.045
Method:                                                 Least Squares   F-statistic:                                                           0.6658
Date:                                                Sun, 18 Dec 2016   Prob (F-statistic):                                                 0.428
Time:                                                              20:51:45   Log-Likelihood:                                               -34.192
No. Observations:                                                     16   AIC:                                                                    72.38
Df Residuals:                                                            14   BIC:                                                                    73.93
Df Model:                                                                   1
Covariance Type:                                         nonrobust
==============================================================================
coef                    std err                         t                         P>|t|                      [95.0% Conf. Int.]
------------------------------------------------------------------------------------------------------------------------------------
const                 1.2861                    1.657                  0.776                        0.451                      -2.269         4.841
x1                    -0.0792                    0.097                 -0.816                        0.428                      -0.287         0.129
==============================================================================
Omnibus:                                                          15.189   Durbin-Watson:                                                      1.774
Prob(Omnibus):                                                  0.001   Jarque-Bera (JB):                                                   2.191
Skew:                                                                -0.089   Prob(JB):                                                                0.334
Kurtosis:                                                             1.196   Cond. No.                                                                51.9
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

â€‹

â€‹

â€‹

Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:                                                         y   No. Observations:                                                         16
Model:                                                              GLM   Df Residuals:                                                                14
Model Family:                                            Binomial   Df Model:                                                                       1
Method:                                                             IRLS   Log-Likelihood:                                                   -10.684
Date:                                              Sun, 18 Dec 2016   Deviance:                                                              21.369
Time:                                                            20:51:45   Pearson chi2:                                                            15.9
No. Iterations:                                                          6
==============================================================================
coef                        std err                        z                        P>|z|                    [95.0% Conf. Int.]
------------------------------------------------------------------------------------------------------------------------------------
Intercept          1.2861                         1.553                 0.828                       0.408                        -1.758     4.330
x                     -0.0792                         0.091                -0.871                      0.384                         -0.257     0.099
==============================================================================

bottom of page