#Next, we can test the models for violations of the equidispersion assumption of Poisson models#test the hypotehsis of equidispersion (var[mu] = mu) against that of QuasiPoisson (var[mu] = phi * mu)#Results = [phi, tvalue, pvalue]printphi_disp(grav)printphi_disp(prod)printphi_disp(att)printphi_disp(doub)#We can see for all four models there is overdispersion (phi >> 0), #which are statistically significant according the tvalues (large)#and pvalues (essentially zero). It does however decrease as more #constraints are introduced and model fit increases

#As a result we can compare our standard errors and tvalues for a Poisson model to a QuasiPoissonprint'Production-constrained Poisson model standard errors and tvalues'printprod.params[-4:]printprod.std_err[-4:]printprod.tvalues[-4:]#Fit the same model using QuasiPoisson frameworkQuasi=Production(flows,o,d_vars,cost,'exp',Quasi=True)print'Production-constrained QuasiPoisson model standard errors and tvalues'printQuasi.params[-4:]printQuasi.std_err[-4:]printQuasi.tvalues[-4:]#As we can see both models result in the same parameters (first line)#We also see the QuasiPoisson results in larger standard errors (middle line)#Which then results in smaller t-values (bottom line)#We would even consdier rejecting the statistical significant of of the #parameter estimate on destination building square footage because its abolsute # value is less than 1.96 (96% confidence level)

#We can also estimate a local model which subsets the data#For a production constrained model this means each local model#is from one origin to all destinations. Since we get a set of #parameter estimates for each origin, we can then map them.local_prod=prod.local()

#There is a set of local parameter estimates, tvalues, and pvalues for each covariate#And there is a set of local values for each diagnosticlocal_prod.keys()

#Prep geometry for plotting#Read in census tracts for NYCcrs={'datum':'WGS84','proj':'longlat'}tracts=ps.examples.get_path('nyct2010.shp')tracts=gp.read_file(tracts)tracts=tracts.to_crs(crs=crs)#subset manhattan tractsman_tracts=tracts[tracts['BoroCode']=='1'].copy()man_tracts['CT2010S']=man_tracts['CT2010'].astype(int).astype(str)#Get tracts for which there are no onbservationsmt=set(man_tracts.CT2010S.unique())lt=set(np.unique(o))nt=list(mt.difference(lt))no_tracts=pd.DataFrame({'no_tract':nt})no_tracts=man_tracts[man_tracts.CT2010S.isin(nt)].copy()#Join local values to census tractslocal_vals=pd.DataFrame({'betas':local_prod['param3'],'tract':np.unique(o)})local_vals=pd.merge(local_vals,man_tracts[['CT2010S','geometry']],left_on='tract',right_on='CT2010S')local_vals=gp.GeoDataFrame(local_vals)

#Plot local estimates for destination capacity: darker red is larger effect; grey is no datafig=plt.figure(figsize=(12,12))ax=fig.add_subplot(111)local_vals['cap']=local_prod['param3']no_tracts['test']=0no_tracts.plot('test',cmap='copper',ax=ax)local_vals.plot('cap',cmap='Reds',ax=ax)plt.legend()plt.xlim(-74.02,-73.95)plt.ylim(40.7,40.78)

(40.7, 40.78)

#Plot local estimates for # of housing units: darker red is larger effect; grey is no datafig=plt.figure(figsize=(12,12))ax=fig.add_subplot(111)local_vals['house']=local_prod['param2']no_tracts['test']=0no_tracts.plot('test',cmap='copper',ax=ax)local_vals.plot('house',cmap='Reds',ax=ax)plt.legend()plt.xlim(-74.02,-73.95)plt.ylim(40.7,40.78)

(40.7, 40.78)

#Plot local estimates for destination building sq footage: darker red is larger effect; grey is no datafig=plt.figure(figsize=(12,12))ax=fig.add_subplot(111)local_vals['foot']=local_prod['param1']no_tracts['test']=0no_tracts.plot('test',cmap='copper',ax=ax)local_vals.plot('foot',cmap='Reds',ax=ax)plt.legend()plt.xlim(-74.02,-73.95)plt.ylim(40.7,40.78)

(40.7, 40.78)

#Drop NA valueslabels=['start station longitude','start station latitude','end station longitude','end station latitude']bikes=bikes.dropna(subset=labels)

#Origin focused Moran's I of OD pairs as vectors in spacevmo=VecMoran(vecs,wo,permutations=1)vmo.I

0.15734115508365823

#Destination focused Moran's I of OD pairs as vectors in spacevmd=VecMoran(vecs,wd,permutations=1)vmd.I

0.21322002653342809

#No substantial examples to show for spatial interaction weights#Will add them once there is a working SAR Lag spatial interaction#model implementation avaialble#from pysal.weights.spintW import vecW, netW, ODW