multinomial logit python

1In [1]: import statsmodels.api as st
2 
3In [2]: iris = st.datasets.get_rdataset('iris', 'datasets')
4 
5In [3]: ### get the y 
6 
7In [4]: y = iris.data.Species
8 
9In [5]: print y.head(3)
100    setosa
111    setosa
122    setosa
13Name: Species, dtype: object
14 
15In [6]: ### get the x 
16 
17In [7]: x = iris.data.ix[:, 0]
18 
19In [8]: x = st.add_constant(x, prepend = False)
20 
21In [9]: print x.head(3)
22   Sepal.Length  const
230           5.1      1
241           4.9      1
252           4.7      1
26 
27In [10]: ### specify the model
28 
29In [11]: mdl = st.MNLogit(y, x)
30 
31In [12]: mdl_fit = mdl.fit()
32Optimization terminated successfully.
33         Current function value: 0.606893
34         Iterations 8
35 
36In [13]: ### print model summary ###
37 
38In [14]: print mdl_fit.summary()
39                          MNLogit Regression Results                          
40==============================================================================
41Dep. Variable:                Species   No. Observations:                  150
42Model:                        MNLogit   Df Residuals:                      146
43Method:                           MLE   Df Model:                            2
44Date:                Fri, 23 Aug 2013   Pseudo R-squ.:                  0.4476
45Time:                        22:22:58   Log-Likelihood:                -91.034
46converged:                       True   LL-Null:                       -164.79
47                                        LLR p-value:                 9.276e-33
48=====================================================================================
49Species=versicolor       coef    std err          z      P>|z|      [95.0% Conf. Int.]
50--------------------------------------------------------------------------------------
51Sepal.Length           4.8157      0.907      5.310      0.000         3.038     6.593
52const                -26.0819      4.889     -5.335      0.000       -35.665   -16.499
53--------------------------------------------------------------------------------------
54Species=virginica       coef    std err          z      P>|z|      [95.0% Conf. Int.]
55-------------------------------------------------------------------------------------
56Sepal.Length          6.8464      1.022      6.698      0.000         4.843     8.850
57const               -38.7590      5.691     -6.811      0.000       -49.913   -27.605
58=====================================================================================
59 
60In [15]: ### marginal effects ###
61 
62In [16]: mdl_margeff = mdl_fit.get_margeff()
63 
64In [17]: print mdl_margeff.summary()
65       MNLogit Marginal Effects      
66=====================================
67Dep. Variable:                Species
68Method:                          dydx
69At:                           overall
70=====================================================================================
71    Species=setosa      dy/dx    std err          z      P>|z|      [95.0% Conf. Int.]
72--------------------------------------------------------------------------------------
73Sepal.Length          -0.3785      0.003   -116.793      0.000        -0.385    -0.372
74--------------------------------------------------------------------------------------
75Species=versicolor      dy/dx    std err          z      P>|z|      [95.0% Conf. Int.]
76--------------------------------------------------------------------------------------
77Sepal.Length           0.0611      0.022      2.778      0.005         0.018     0.104
78--------------------------------------------------------------------------------------
79Species=virginica      dy/dx    std err          z      P>|z|      [95.0% Conf. Int.]
80-------------------------------------------------------------------------------------
81Sepal.Length          0.3173      0.022     14.444      0.000         0.274     0.360
82=====================================================================================
83 
84In [18]: ### aic and bic ###
85 
86In [19]: print mdl_fit.aic
87190.06793279
88 
89In [20]: print mdl_fit.bic
90202.110473966
91<div class="open_grepper_editor" title="Edit & Save To Grepper"></div>