In [1]:

    

import pandas as pd
import glob
from os import path
import numpy as np
from __future__ import division 
import sqlalchemy as sqla
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels
from sklearn import linear_model
import matplotlib
import rpy2.robjects as robjects
#import rpy2.robjects.packages as rpackages
from rpy2.robjects import pandas2ri


    

In [2]:

    

#%matplotlib qt
%matplotlib inline


    

In [3]:

    

engine = sqla.create_engine('postgresql://postgres:postgres@localhost:5432/TaxiData',echo=False)

sqlmycols={'MOGE001':'total_commuters',\
           'MOGE011':'car_commuters',\
          'MRUE001':'pc_income',\
          'totalpopulation':'population',\
          'MOGM001':'me_total_commuters',\
          'MOGM011':'me_car_commuters',\
          'MRUM001':'me_pc_income',\
          'time_dif_derived_approxcount':'approxcount',\
          'time_dif_derived_approxcount_percerror':'perc_error_approxcount',\
           'time_dif_derived_approxcount_error':'error_approxcount',\
           'twentythirteen_full_count':'count',\
           'fipscodes':'fipscodes',\
           'perc_white':'white_ratio'
          }
columnstring=', '.join(['"%s" AS %s' % (i,j) for i,j in sqlmycols.iteritems()])


    

In [4]:

    

full=pd.read_sql_query('SELECT %s FROM lotsofdata WHERE totalpopulation>=1000' % (columnstring) ,engine).set_index('fipscodes')


    

In [5]:

    

full['driver_ratio']=full['car_commuters']/full['total_commuters']
full['approxcount_pc']=full['approxcount']/full['population']
full['count_pc']=full['count']/full['population']


    

In [6]:

    

borocods=pd.read_sql_table('testborocodedict',engine)
borocods.drop('county_name_test',inplace=True,axis=1)
borocods.set_index('code',inplace=True)
borocods


    

    
Out[6]:






  

    

      

      
int
      
name
      
county_name
    
    

      
code
      

      

      

    
  
  

    

      
005
      
2
      
Bronx
      
Bronx
    
    

      
047
      
3
      
Brooklyn
      
Kings
    
    

      
061
      
1
      
Manhattan
      
New York
    
    

      
081
      
4
      
Queens
      
Queens
    
    

      
085
      
5
      
Staten Island
      
Richmond
    
  

In [7]:

    

full['borocode']=borocods.loc[full.index.to_series().str[2:5].values]['int'].values


    

In [8]:

    

psub=full[(full['borocode']!=5)]
#psub.drop(outliers,axis=0,inplace=True)
model=statsmodels.discrete.discrete_model.Poisson((psub['count']),\
                                                  sm.add_constant(np.log(psub['pc_income']),prepend=False),\
                                                 exposure=psub['population'].values)
fitp=model.fit(start_params=[2.2948,-21.5754])
fitp.summary()


    

    




Optimization terminated successfully.
         Current function value: 50432.844615
         Iterations 4






    
Out[8]:





Poisson Regression Results

  
Dep. Variable:
       
count
      
  No. Observations:  
   
  1979
   


  
Model:
              
Poisson
     
  Df Residuals:      
   
  1977
   


  
Method:
               
MLE
       
  Df Model:          
   
     1
   


  
Date:
          
Fri, 04 Mar 2016
 
  Pseudo R-squ.:     
   
0.6460
   


  
Time:
              
17:21:35
     
  Log-Likelihood:    
 
-9.9807e+07


  
converged:
           
True
       
  LL-Null:           
 
-2.8196e+08


  
 
                      
 
        
  LLR p-value:       
   
 0.000
   




      
         
coef
     
std err
      
z
      
P>|z|
 
[95.0% Conf. Int.]
 


  
pc_income
 
    2.1338
 
    0.000
 
 1.75e+04
 
 0.000
 
    2.134     2.134


  
const
     
  -19.7852
 
    0.001
 
-1.45e+04
 
 0.000
 
  -19.788   -19.783

In [9]:

    

R=robjects.r
pandas2ri.activate()


    

In [10]:

    

r_model=R.glm('count~1+log(pc_income)', family=R.quasipoisson(),data=psub,offset=np.log(psub.population))


    

In [11]:

    

deviance=R.summary(r_model).rx2('deviance')[0]
'%.4g' % deviance


    

    
Out[11]:





'1.996e+08'

In [12]:

    

dispersion=R.summary(r_model).rx2('dispersion')[0]


    

In [13]:

    

'%.4g' % dispersion


    

    
Out[13]:





'5.176e+05'

In [14]:

    

deviance/dispersion


    

    
Out[14]:





385.57723365839183

In [15]:

    

SIsub=full[(full['borocode']==5)]


    

In [16]:

    

# Plotting the data:
#plt.rc('text', **{'usetex':True,'latex.preamble':[
#       r'\usepackage{siunitx}',   
#       r'\sisetup{detect-all}',   
#       r'\usepackage{helvet}',    
#       r'\usepackage{sansmath}',  
#       r'\sansmath'               
#]  })
comd=0
buff=100

matplotlib.rc('legend', fontsize=10)

border=0.2
boro=1

#asub=full[(full.drivercommuterratio<=border)]
#bsub=full[(full.drivercommuterratio>border)]
XX=np.linspace(psub.pc_income.min(),psub.pc_income.max())

scatter_kwargs = {"zorder":100}
error_kwargs = {"lw":.5, "zorder":0}

a=6
gr=(1+np.sqrt(5))/2
plt.figure(figsize=[a*gr,a])
plt.autoscale(tight=True)
plt.xlabel('Per-capita income')
plt.ylabel('2013 yellow cab dropoffs per-capita')
plt.title("Yellow cab dropoffs vs income of census tract")
colorlist=plt.cm.gist_rainbow(np.linspace(0,2.8/4,4)).tolist()
colorlist.reverse()
for i in borocods[borocods['int']<5]['int'].values:
    #print(i)
    isub=psub[psub.borocode==i]
    a=plt.scatter(isub.pc_income,isub.count_pc,alpha=0.5,label=borocods[borocods['int']==i]['name'].values[0],c=colorlist[i-1],**scatter_kwargs)
    #plt.errorbar(isub.pc_income,isub.approxcount_pc,yerr=nintypercentile*isub.error_approxcount/isub.population,xerr=isub.me_pc_income,fmt=None,marker=None,mew=0,ecolor='black',alpha=0.5,label='',**error_kwargs)
#plt.scatter(bsub['pc_income'],bsub['approxcount_pc'],facecolors='none',edgecolors='black',alpha=0.1,label='driver commuter ratio $>$ '+str(round(border,1)))
#plt.scatter(full.loc[outliers|highleverage]['pc_income'],full.loc[outliers|highleverage]['approxcount_pc'],alpha=0.25,label=r'Outliers & high leverage',color='black')
plt.scatter(SIsub.pc_income,SIsub.count_pc,facecolors='none',edgecolors='black',alpha=0.1,label=borocods[borocods['int']==5]['name'].values[0])
#plt.scatter(outliers.pc_income,outliers.count_pc, facecolors='black',edgecolors='black',alpha=0.2, label='Outliers')
#plt.scatter(subset.loc[highleverage]['pc_income'],subset.loc[highleverage]['count_pc'],alpha=0.3,label='High leverage',color='orange')
plt.plot(XX,(XX**fitp.params['pc_income'])*np.exp(fitp.params['const']),color='black',alpha=0.8,lw=2, label='Fit, pseudo-$\mathsf{R}^\mathsf{2}=$ '+str(round(fitp.prsquared,3)) )
#cbar=plt.colorbar(a, ticks=borocods[borocods['int']<5]['int'].tolist())
#cbar.set_label("boro")


plt.yscale('log')
plt.xscale('log')
plt.legend(loc='lower right')
plt.savefig('2013count_dropoffs_vs_income_low_car_poisson.png',dpi=500)
plt.show()


    

In [17]:

    

psub['prediction']=pandas2ri.ri2py(r_model.rx2('fitted.values'))

psub['anscombe_resid']=(3/2)*(psub['count']**(2/3)-psub['prediction']**(2/3))/(psub['prediction']**(1/6))


    

    




/usr/local/lib/python2.7/dist-packages/IPython/kernel/__main__.py:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  if __name__ == '__main__':
/usr/local/lib/python2.7/dist-packages/IPython/kernel/__main__.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  app.launch_new_instance()

In [18]:

    

sm.qqplot(psub['anscombe_resid'],line='s')
plt.ylabel('Anscombe residual quantiles')
plt.savefig('2013quasi-poisson_QQ.png',dpi=500)
plt.show()


    

In [19]:

    

plt.scatter(psub['pc_income'],psub['anscombe_resid'])
plt.xscale('log')
plt.show()