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) 
X = sm.add_constant(X)                                             # add intercept

 

# 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
Model:                                                                 WLS   Adj. R-squared:                                                   -0.023
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
Link Function:                                                    logit   Scale:                                                                           1.0
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
==============================================================================