import%20marimo%0A%0A__generated_with%20%3D%20%220.23.9%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Practical%3A%20Does%20Industrial%20Activity%20Affect%20Local%20Air%20Pollution%3F%0A%0A%20%20%20%20%3E%20%E2%9A%A0%EF%B8%8F%20**Teaching%20example%2C%20not%20a%20real%20causal%20study.**%20This%20notebook%20is%20a%20**simplified%2C%20illustrative**%20walk-through%20of%20spatial%20and%20causal%20methods.%20The%20numbers%20are%20**not**%20validated%20findings%20about%20livestock%20and%20ammonia%20and%20should%20not%20be%20cited%3B%20read%20every%20estimate%20as%20a%20demonstration%20of%20how%20a%20method%20behaves%20(and%20where%20it%20breaks).%0A%0A%20%20%20%20%23%23%20From%20spatial%20correlation%20to%20causal%20reasoning%2C%20in%20four%20steps%0A%0A%20%20%20%20%7C%20Step%20%7C%20Question%20%7C%0A%20%20%20%20%7C------%7C----------%7C%0A%20%20%20%20%7C%20**1**%20%7C%20Are%20livestock%20farms%20and%20NH%E2%82%83%20spatially%20associated%3F%20(OLS)%20%7C%0A%20%20%20%20%7C%20**2**%20%7C%20Can%20spatial%20models%20fix%20the%20residual%20clustering%3F%20(OLS%20vs%20SAR)%20%7C%0A%20%20%20%20%7C%20**3**%20%7C%20Can%20farm%20gains%20identify%20a%20causal%20effect%3F%20(Difference-in-Differences)%20%7C%0A%20%20%20%20%7C%20**4**%20%7C%20What%20happens%20when%20pollution%20spills%20across%20cells%3F%20(spatial%20model%20%2B%20DiD)%20%7C%0A%0A%20%20%20%20%3E%20**A%20map%20shows%20you%20*where*.%20A%20regression%20shows%20you%20*what%20correlates*.%0A%20%20%20%20%3E%20A%20design%20tells%20you%20*what%20would%20need%20to%20be%20true*%20for%20a%20causal%20claim.**%0A%0A%20%20%20%20Everything%20works%20at%20the%20**1%20%C3%97%201%20km%20grid%20cell**%20level%2C%20the%20same%0A%20%20%20%20resolution%20as%20the%20official%20RIVM%20pollution%20maps.%0A%0A%20%20%20%20**Data%3A**%20RIVM%20NH%E2%82%83%20concentration%20grids%20(2018%E2%80%932024)%2C%20Bureau%20van%20Dijk%20Orbis%0A%20%20%20%20firm%20data%2C%20and%20CBS%20neighbourhood%20statistics.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Setup%3A%20load%20libraries%20and%20a%20few%20small%20helpers%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20import%20contextlib%2C%20io%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20quietly%20silence%20spreg's%20chatty%20output%0A%20%20%20%20import%20warnings%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20hide%20non-essential%20warnings%0A%20%20%20%20import%20numpy%20as%20np%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20arrays%20%2B%20fast%20maths%0A%20%20%20%20import%20pandas%20as%20pd%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20data%20tables%20(our%20grid%20is%20one%20big%20DataFrame)%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%20%20%20%20%20%20%20%20%20%20%20%20%23%20plotting%0A%20%20%20%20import%20statsmodels.formula.api%20as%20smf%20%20%20%20%20%20%23%20regression%20with%20R-style%20formulas%20(OLS%2C%20DiD)%0A%20%20%20%20from%20libpysal.weights%20import%20DistanceBand%20%20%23%20build%20the%20spatial%20weights%20matrix%20W%0A%20%20%20%20from%20esda.moran%20import%20Moran%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20Moran's%20I%3A%20spatial%20clustering%20in%20the%20residuals%3F%0A%20%20%20%20from%20spreg%20import%20ML_Lag%2C%20ML_Error%20%20%20%20%20%20%20%20%20%23%20spatial%20regressions%3A%20SAR%2FSDM%20%3D%20ML_Lag%2C%20SEM%20%3D%20ML_Error%0A%20%20%20%20from%20pathlib%20import%20Path%0A%0A%20%20%20%20warnings.filterwarnings(%22ignore%22)%0A%20%20%20%20pd.set_option(%22display.float_format%22%2C%20lambda%20x%3A%20f%22%7Bx%3A.3f%7D%22)%0A%20%20%20%20plt.rcParams%5B%22figure.dpi%22%5D%20%3D%20120%0A%20%20%20%20plt.rcParams%5B%22figure.facecolor%22%5D%20%3D%20%22white%22%0A%0A%20%20%20%20def%20style_map_axes(*axes)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Make%20scatter%20'maps'%20look%20like%20maps%3A%20equal%20aspect%2C%20no%20frame%2Fticks.%22%22%22%0A%20%20%20%20%20%20%20%20for%20_ax%20in%20axes%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.set_aspect(%22equal%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.set_axis_off()%0A%0A%20%20%20%20def%20fit_spatial(model_cls%2C%20y%2C%20X%2C%20w%2C%20name_y%2C%20name_x%2C%20**kwargs)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Fit%20a%20spreg%20spatial%20model%20quietly%3B%20return%20(model%2C%20residuals%2C%20Moran's%20I).%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20SAR%20-%3E%20ML_Lag(...)%20%20%20%20%20%20%20%20%20%20%20%20%20lag%20of%20the%20OUTCOME%20%20%20%20%20%20%20%20%20%20%20%20(rho%20*%20W%20y)%0A%20%20%20%20%20%20%20%20%20%20%20%20SDM%20-%3E%20ML_Lag(...%2C%20slx_lags%3D1)%20also%20lag%20of%20X%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20(%2B%20theta%20*%20W%20X)%0A%20%20%20%20%20%20%20%20%20%20%20%20SEM%20-%3E%20ML_Error(...)%20%20%20%20%20%20%20%20%20%20%20spatially%20correlated%20errors%20%20%20(lambda)%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20with%20contextlib.redirect_stdout(io.StringIO())%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20m%20%3D%20model_cls(y%2C%20X%2C%20w%3Dw%2C%20name_y%3Dname_y%2C%20name_x%3Dname_x%2C%20method%3D%22lu%22%2C%20**kwargs)%0A%20%20%20%20%20%20%20%20resid%20%3D%20(m.e_filtered%20if%20model_cls%20is%20ML_Error%20else%20m.u).flatten()%0A%20%20%20%20%20%20%20%20return%20m%2C%20resid%2C%20Moran(resid%2C%20w)%0A%0A%20%20%20%20def%20spatial_impacts(model)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Tidy%20Direct%20%2F%20Indirect%20%2F%20Total%20effects%20table%20read%20off%20the%20model%20summary.%22%22%22%0A%20%20%20%20%20%20%20%20_rows%2C%20_capture%20%3D%20%5B%5D%2C%20False%0A%20%20%20%20%20%20%20%20for%20_line%20in%20model.summary.splitlines()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20%22IMPACTS%22%20in%20_line%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_capture%20%3D%20True%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20continue%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20_capture%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_parts%20%3D%20_line.split()%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20len(_parts)%20%3D%3D%204%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_rows.append(%5B_parts%5B0%5D%2C%20float(_parts%5B1%5D)%2C%20float(_parts%5B2%5D)%2C%20float(_parts%5B3%5D)%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20except%20ValueError%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pass%0A%20%20%20%20%20%20%20%20return%20pd.DataFrame(_rows%2C%20columns%3D%5B%22variable%22%2C%20%22direct%22%2C%20%22indirect%22%2C%20%22total%22%5D).set_index(%22variable%22)%0A%0A%20%20%20%20DATA%20%3D%20Path(%22..%2Fdata%2Ffinal%22)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20DATA%2C%0A%20%20%20%20%20%20%20%20DistanceBand%2C%0A%20%20%20%20%20%20%20%20ML_Error%2C%0A%20%20%20%20%20%20%20%20ML_Lag%2C%0A%20%20%20%20%20%20%20%20Moran%2C%0A%20%20%20%20%20%20%20%20fit_spatial%2C%0A%20%20%20%20%20%20%20%20np%2C%0A%20%20%20%20%20%20%20%20pd%2C%0A%20%20%20%20%20%20%20%20plt%2C%0A%20%20%20%20%20%20%20%20smf%2C%0A%20%20%20%20%20%20%20%20spatial_impacts%2C%0A%20%20%20%20%20%20%20%20style_map_axes%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%E2%94%80%20SECTOR%20SELECTOR%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%23%20Pick%20which%20Orbis%20sector%20to%20study%20against%20NH%E2%82%83.%20The%20WHOLE%20analysis%20below%0A%20%20%20%20%23%20re-runs%20reactively%20when%20you%20change%20this.%20Default%20is%20livestock.%20%0A%20%20%20%20%23%20The%20legends%20will%20keep%20saying%20%22farms%22%20even%20if%20you%20choose%20the%20sector%0A%20%20%20%20SECTORS%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%22Livestock%20(NACE%20014x)%22%3A%20%22n_livestock%22%2C%0A%20%20%20%20%20%20%20%20%22Agriculture%2C%20Horticulture%20%26%20Livestock%22%3A%20%22n_agriculture_horticulture_lives%22%2C%0A%20%20%20%20%20%20%20%20%22Food%20%26%20Tobacco%20Manufacturing%22%3A%20%22n_food_tobacco_manufacturing%22%2C%0A%20%20%20%20%20%20%20%20%22Chemicals%2C%20Petroleum%2C%20Rubber%20%26%20Plastic%22%3A%20%22n_chemicals_petroleum_rubber_pla%22%2C%0A%20%20%20%20%20%20%20%20%22Business%20Services%22%3A%20%22n_business_services%22%2C%0A%20%20%20%20%20%20%20%20%22Wholesale%22%3A%20%22n_wholesale%22%2C%0A%20%20%20%20%20%20%20%20%22Construction%22%3A%20%22n_construction%22%2C%0A%20%20%20%20%20%20%20%20%22Transport%2C%20Freight%20%26%20Storage%22%3A%20%22n_transport_freight_storage%22%2C%0A%20%20%20%20%20%20%20%20%22Retail%22%3A%20%22n_retail%22%2C%0A%20%20%20%20%7D%0A%20%20%20%20SECTOR%20%3D%20%22Livestock%20(NACE%20014x)%22%0A%20%20%20%20return%20SECTOR%2C%20SECTORS%0A%0A%0A%40app.cell%0Adef%20_(DATA%2C%20SECTOR%2C%20SECTORS%2C%20np%2C%20pd)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Load%20the%20grid%20and%20build%20generic%20n_sector_YYYY%20columns%20for%20the%20pick%20%E2%94%80%E2%94%80%0A%20%20%20%20_prefix%20%3D%20SECTORS%5BSECTOR%5D%0A%20%20%20%20sector_label%20%3D%20%22livestock%20farms%22%20if%20SECTOR.startswith(%22Livestock%22)%20else%20SECTOR.lower()%0A%0A%20%20%20%20%23%20Create%20a%20column%20n_sector_YEAR%20to%20use%20in%20next%20steps%0A%20%20%20%20grid%20%3D%20pd.read_csv(DATA%20%2F%20%22workshop_grid_1km.csv%22)%0A%20%20%20%20for%20_yr%20in%20range(2018%2C%202025)%3A%0A%20%20%20%20%20%20%20%20_src%20%3D%20f%22%7B_prefix%7D_%7B_yr%7D%22%0A%20%20%20%20%20%20%20%20grid%5Bf%22n_sector_%7B_yr%7D%22%5D%20%3D%20grid%5B_src%5D%20if%20_src%20in%20grid.columns%20else%200%0A%0A%20%20%20%20grid%5B%22log_pop_density%22%5D%20%3D%20np.log1p(grid%5B%22population_density%22%5D.clip(lower%3D0))%0A%20%20%20%20grid%5B%22rurality%22%5D%20%3D%20grid%5B%22urbanity_index%22%5D.round().fillna(0).astype(int)%0A%20%20%20%20grid%5B%22sector_net_change%22%5D%20%3D%20grid%5B%22n_sector_2024%22%5D%20-%20grid%5B%22n_sector_2020%22%5D%0A%0A%20%20%20%20_n_sec%20%3D%20int(grid%5B%22n_sector_2024%22%5D.sum())%0A%20%20%20%20_n_cells%20%3D%20int((grid%5B%22n_sector_2024%22%5D%20%3E%200).sum())%0A%20%20%20%20_n_gain%20%3D%20int((grid%5B%22sector_net_change%22%5D%20%3E%201).sum())%0A%20%20%20%20print(f%22Sector%3A%20%7BSECTOR%7D%22)%0A%20%20%20%20print(f%22Grid%3A%20%7Blen(grid)%3A%2C%7D%20cells%20(%7Bgrid.shape%5B1%5D%7D%20columns)%22)%0A%20%20%20%20print(f%22%20%20%3E%20%7B_n_sec%3A%2C%7D%20%7Bsector_label%7D%20active%20in%202024%20across%20%7B_n_cells%3A%2C%7D%20cells%22)%0A%20%20%20%20print(f%22%20%20%3E%20%7B_n_gain%3A%2C%7D%20cells%20gained%20%3E%201%20%7Bsector_label%7D%20between%202020%20and%202024%22)%0A%20%20%20%20return%20grid%2C%20sector_label%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Part%200.%20Meet%20the%20data%0A%0A%20%20%20%20Before%20any%20analysis%2C%20get%20a%20feel%20for%20the%20grid.%20Each%20**row%20is%20one%0A%20%20%20%201%20%C3%97%201%20km%20cell**%3B%20the%20columns%20hold%20the%20NH%E2%82%83%20pollution%20series%2C%20firm%20counts%0A%20%20%20%20per%20sector%20and%20year%2C%20and%20a%20few%20neighbourhood%20covariates.%0A%0A%20%20%20%20Look%20at%3A%0A%0A%20%20%20%20-%20the%20**shape**%3A%20how%20many%20cells%2C%20how%20many%20columns%0A%20%20%20%20-%20the%20**NH%E2%82%83%20series**%20%60nh3_2018%60%20%E2%80%A6%20%60nh3_2024%60%0A%20%20%20%20-%20the%20**sector%20counts**%20%60n_sector_2018%60%20%E2%80%A6%20%60n_sector_2024%60%0A%20%20%20%20-%20the%20**covariates**%20%60population_density%60%2C%20%60urbanity_index%60%0A%20%20%20%20-%20the%20helper%20columns%20built%20in%20setup%3A%20%60n_sector_2024%60%2C%0A%20%20%20%20%20%20%60sector_net_change%60%2C%20%60log_pop_density%60%2C%20%60rurality%60%0A%0A%20%20%20%20%3E%20Keep%20one%20question%20in%20mind%3A%20at%20the%20**1%20km**%20scale%2C%20what%20does%20a%0A%20%20%20%20%3E%20%22sector%20cell%22%20actually%20mean%2C%20and%20would%20that%20change%20at%20250%20m%2C%20or%20at%0A%20%20%20%20%3E%20the%20municipality%20level%3F%20(That%20is%20the%20MAUP%20from%20the%20lecture.)%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(grid%2C%20mo%2C%20sector_label)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Part%200.%20A%20quick%20look%20at%20the%20grid%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_peek%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%22x_rd%22%2C%20%22y_rd%22%2C%0A%20%20%20%20%20%20%20%20%22nh3_2020%22%2C%20%22nh3_2024%22%2C%0A%20%20%20%20%20%20%20%20%22n_sector_2020%22%2C%20%22n_sector_2024%22%2C%20%22sector_net_change%22%2C%0A%20%20%20%20%20%20%20%20%22population_density%22%2C%20%22urbanity_index%22%2C%0A%20%20%20%20%5D%0A%20%20%20%20print(f%22Grid%20shape%3A%20%7Bgrid.shape%5B0%5D%3A%2C%7D%20cells%20x%20%7Bgrid.shape%5B1%5D%7D%20columns%22)%0A%20%20%20%20print(f%22Cells%20with%20%3E%3D%201%20%7Bsector_label%7D%20in%202024%3A%20%7B(grid%5B'n_sector_2024'%5D%20%3E%200).sum()%3A%2C%7D%22)%0A%20%20%20%20print(f%22Cells%20that%20gained%20%3E%201%20%7Bsector_label%7D%20(2020%20-%3E%202024)%3A%20%7B(grid%5B'sector_net_change'%5D%20%3E%201).sum()%3A%2C%7D%22)%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22**A%20few%20representative%20columns%20(first%20rows)%3A**%22)%2C%0A%20%20%20%20%20%20%20%20grid%5B_peek%5D.head()%2C%0A%20%20%20%20%20%20%20%20mo.md(%22**Distribution%20of%20the%20outcome%2C%20treatment%20intensity%2C%20and%20a%20covariate%3A**%22)%2C%0A%20%20%20%20%20%20%20%20grid%5B%5B%22nh3_2020%22%2C%20%22nh3_2024%22%2C%20%22n_sector_2024%22%2C%20%22sector_net_change%22%2C%20%22population_density%22%5D%5D.describe().round(2)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Which%20sectors%20locate%20in%20high-NH%E2%82%83%20areas%3F%0A%0A%20%20%20%20Before%20any%20regression%2C%20a%20descriptive%20question%3A%20**which%20BvD%20sectors%0A%20%20%20%20tend%20to%20be%20in%20places%20with%20more%20NH%E2%82%83%3F**%0A%0A%20%20%20%20The%20bar%20chart%20shows%20the%20mean%20NH%E2%82%83%20concentration%20at%20firm%20locations%20by%0A%20%20%20%20sector.%20**Livestock%20(NACE%20014x)**%20comes%20from%20a%20separate%20Orbis%20extract%0A%20%20%20%20without%20an%20employee-size%20threshold%20(the%20main%20firm%20data%20only%20has%20firms%0A%20%20%20%20with%20%E2%89%A5%2010%20employees).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(grid%2C%20np%2C%20pd%2C%20plt)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Which%20sectors%20sit%20in%20high-NH3%20places%3F%20(a%20first%20descriptive%20look)%20%E2%94%80%E2%94%80%0A%20%20%20%20_rows%20%3D%20%5B%5D%0A%20%20%20%20for%20_sector%20in%20%5B%0A%20%20%20%20%20%20%20%20%22livestock%22%2C%0A%20%20%20%20%20%20%20%20%22agriculture_horticulture_lives%22%2C%0A%20%20%20%20%20%20%20%20%22food_tobacco_manufacturing%22%2C%0A%20%20%20%20%20%20%20%20%22chemicals_petroleum_rubber_pla%22%2C%0A%20%20%20%20%20%20%20%20%22business_services%22%2C%0A%20%20%20%20%20%20%20%20%22wholesale%22%2C%0A%20%20%20%20%20%20%20%20%22construction%22%2C%0A%20%20%20%20%20%20%20%20%22transport_freight_storage%22%2C%0A%20%20%20%20%20%20%20%20%22retail%22%0A%20%20%20%20%20%20%20%20%5D%3A%0A%20%20%20%20%20%20%20%20_s_cells%20%3D%20grid%5Bgrid%5Bf%22n_%7B_sector%7D_2024%22%5D%20%3E%200%5D%0A%20%20%20%20%20%20%20%20_rows.append(%7B%22sector%22%3A%20_sector.replace(%22_%22%2C%22%20%22).title()%2C%20%22mean_nh3%22%3A%20_s_cells%5B%22nh3_2024%22%5D.mean()%7D)%0A%0A%20%20%20%20_sector_nh3%20%3D%20pd.DataFrame(_rows).sort_values(%22mean_nh3%22)%0A%20%20%20%20_overall_mean%20%3D%20grid%5B%22nh3_2024%22%5D.mean()%0A%0A%20%20%20%20_is_agri%20%3D%20_sector_nh3%5B%22sector%22%5D.str.contains(%22Livestock%7CAgri%22)%0A%20%20%20%20_bar_colors%20%3D%20np.where(_is_agri%2C%20%22%23d95f0e%22%2C%20%22%233b528b%22)%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(10%2C%206))%0A%20%20%20%20_ax.barh(_sector_nh3%5B%22sector%22%5D%2C%20_sector_nh3%5B%22mean_nh3%22%5D%2C%20color%3D_bar_colors)%0A%20%20%20%20_ax.axvline(_overall_mean%2C%20color%3D%22grey%22%2C%20ls%3D%22--%22%2C%20lw%3D1%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Overall%20grid%20mean%20(%7B_overall_mean%3A.1f%7D%20ug%2Fm3)%22)%0A%20%20%20%20_ax.set_xlabel(%22Mean%20NH3%20at%20firm%20locations%20(ug%2Fm3%2C%202024)%22)%0A%20%20%20%20_ax.set_title(%22Which%20sectors%20locate%20in%20high-NH3%20areas%3F%22)%0A%20%20%20%20_ax.legend(loc%3D%22lower%20right%22)%0A%20%20%20%20_fig.tight_layout()%0A%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Maps%3A%20NH%E2%82%83%20and%20sector%20firm%20locations%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(SECTOR%2C%20grid%2C%20plt%2C%20style_map_axes)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Maps%3A%20where%20is%20NH3%20high%2C%20and%20where%20are%20the%20farms%3F%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(14%2C%205))%0A%0A%20%20%20%20_sc%20%3D%20_axes%5B0%5D.scatter(grid%5B%22x_rd%22%5D%2C%20grid%5B%22y_rd%22%5D%2C%20c%3Dgrid%5B%22nh3_2024%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cmap%3D%22YlOrRd%22%2C%20s%3D1%2C%20vmin%3D0%2C%20vmax%3D20)%0A%20%20%20%20_axes%5B0%5D.set_title(%22NH3%20concentration%202024%20(ug%2Fm3)%22)%0A%20%20%20%20_fig.colorbar(_sc%2C%20ax%3D_axes%5B0%5D%2C%20shrink%3D0.7%2C%20label%3D%22NH3%20(ug%2Fm3)%22)%0A%0A%20%20%20%20_axes%5B1%5D.scatter(grid%5B%22x_rd%22%5D%2C%20grid%5B%22y_rd%22%5D%2C%20c%3D%22lightgrey%22%2C%20s%3D0.5%2C%20label%3D%22no%20farms%22)%0A%20%20%20%20_has_farms%20%3D%20grid%5Bgrid%5B%22n_sector_2024%22%5D%20%3E%200%5D%0A%20%20%20%20_sc2%20%3D%20_axes%5B1%5D.scatter(_has_farms%5B%22x_rd%22%5D%2C%20_has_farms%5B%22y_rd%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20c%3D_has_farms%5B%22n_sector_2024%22%5D%2C%20cmap%3D%22Greens%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20s%3D4%2C%20vmin%3D0%2C%20vmax%3D5%2C%20label%3D%22has%20farms%22)%0A%20%20%20%20_axes%5B1%5D.set_title(f%22%7BSECTOR%7D%20per%20cell%20(2024)%22)%0A%20%20%20%20_fig.colorbar(_sc2%2C%20ax%3D_axes%5B1%5D%2C%20shrink%3D0.7%2C%20label%3D%22farms%20in%20cell%22)%0A%20%20%20%20_axes%5B1%5D.legend(loc%3D%22upper%20left%22%2C%20markerscale%3D5%2C%20fontsize%3D9%2C%20frameon%3DTrue)%0A%0A%20%20%20%20style_map_axes(*_axes)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Step%201.%20Are%20sector%20firms%20and%20NH%E2%82%83%20spatially%20associated%3F%20(OLS)%0A%0A%20%20%20%20We%20ask%3A%20*Do%20grid%20cells%20with%20more%20firms%20also%20have%20higher%20NH%E2%82%83%3F*%0A%20%20%20%20This%20is%20a%20**cross-sectional%20association**%2C%20not%20a%20causal%20claim.%0A%0A%20%20%20%20**The%20five%20questions%20for%20this%20step**%0A%0A%20%20%20%20%7C%20Question%20%7C%20Answer%20for%20Step%201%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20**treatment**%20%7C%20number%20of%20firms%20in%20a%20cell%20(2024)%20%7C%0A%20%20%20%20%7C%20**estimand**%20%7C%20the%20*association*%20of%20firm%20count%20with%20NH%E2%82%83%20(descriptive)%20%7C%0A%20%20%20%20%7C%20**comparison**%20%7C%20cells%20with%20more%20vs%20fewer%20firms%2C%20same%20year%20%7C%0A%20%20%20%20%7C%20**credible%20if**%20%7C%20observed%20controls%20capture%20all%20relevant%20differences%20between%20cells%20%7C%0A%20%20%20%20%7C%20**fails%20when**%20%7C%20unmeasured%20place%20characteristics%20impact%20both%20treatment%20and%20outcome%20%7C%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(grid%2C%20np%2C%20plt%2C%20sector_label%2C%20smf%2C%20style_map_axes)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%201a.%20Naive%20OLS%2C%20and%20the%20same%20model%20with%20controls%20(HC3%20robust%20SEs)%20%E2%94%80%E2%94%80%0A%20%20%20%20reg_data%20%3D%20grid.dropna(subset%3D%5B%22population_density%22%2C%20%22urbanity_index%22%5D).copy()%0A%20%20%20%20model_naive%20%3D%20smf.ols(%22nh3_2024%20~%20n_sector_2024%22%2C%20data%3Dreg_data).fit(cov_type%3D%22HC3%22)%20%0A%20%20%20%20%23%20Regression%20with%20controls%0A%20%20%20%20model_adj%20%3D%20smf.ols(%22nh3_2024%20~%20n_sector_2024%20%2B%20log_pop_density%20%2B%20C(rurality)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20data%3Dreg_data).fit(cov_type%3D%22HC3%22)%0A%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(14%2C%205)%2C%20gridspec_kw%3D%7B%22width_ratios%22%3A%20%5B1%2C%201.2%5D%7D)%0A%0A%20%20%20%20_jitter%20%3D%20np.random.default_rng(42).normal(0%2C%200.15%2C%20size%3Dlen(reg_data))%0A%20%20%20%20_axes%5B0%5D.scatter(reg_data%5B%22n_sector_2024%22%5D%20%2B%20_jitter%2C%20reg_data%5B%22nh3_2024%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20alpha%3D0.04%2C%20s%3D4%2C%20color%3D%22%233b528b%22)%0A%20%20%20%20_x_line%20%3D%20np.linspace(0%2C%20reg_data%5B%22n_sector_2024%22%5D.max()%2C%20100)%0A%20%20%20%20_y_line%20%3D%20model_naive.params%5B%22Intercept%22%5D%20%2B%20model_naive.params%5B%22n_sector_2024%22%5D%20*%20_x_line%0A%20%20%20%20_axes%5B0%5D.plot(_x_line%2C%20_y_line%2C%20color%3D%22%23d95f0e%22%2C%20lw%3D2.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Naive%20slope%20%3D%20%7Bmodel_naive.params%5B'n_sector_2024'%5D%3A.2f%7D%20(p%20%3C%200.001)%22)%0A%20%20%20%20_axes%5B0%5D.set_xlabel(f%22Number%20of%20%7Bsector_label%7D%20in%20cell%22)%0A%20%20%20%20_axes%5B0%5D.set_ylabel(%22NH3%20(ug%2Fm3%2C%202024)%22)%0A%20%20%20%20_axes%5B0%5D.set_title(%22Naive%20association%22)%0A%20%20%20%20_axes%5B0%5D.legend(fontsize%3D10)%0A%0A%20%20%20%20reg_data%5B%22residual%22%5D%20%3D%20model_adj.resid%0A%20%20%20%20_sc%20%3D%20_axes%5B1%5D.scatter(reg_data%5B%22x_rd%22%5D%2C%20reg_data%5B%22y_rd%22%5D%2C%20c%3Dreg_data%5B%22residual%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cmap%3D%22coolwarm%22%2C%20s%3D1%2C%20vmin%3D-5%2C%20vmax%3D5)%0A%20%20%20%20_axes%5B1%5D.set_title(%22OLS%20residuals%20(with%20controls)%2C%20still%20clustered%22)%0A%20%20%20%20_fig.colorbar(_sc%2C%20ax%3D_axes%5B1%5D%2C%20shrink%3D0.7%2C%20label%3D%22residual%20(ug%2Fm3)%22)%0A%20%20%20%20style_map_axes(_axes%5B1%5D)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%20model_adj%2C%20model_naive%2C%20reg_data%0A%0A%0A%40app.cell%0Adef%20_(model_adj%2C%20model_naive%2C%20pd)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%201b.%20OLS%20coefficients%20(sector%20AND%20the%20controls)%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_coef_table%20%3D%20pd.DataFrame(%7B%0A%20%20%20%20%20%20%20%20%22naive%22%3A%20model_naive.params%2C%0A%20%20%20%20%20%20%20%20%22%2B%20controls%22%3A%20model_adj.params%2C%0A%20%20%20%20%7D)%0A%20%20%20%20_coef_table%5B%22p%20(%2Bcontrols)%22%5D%20%3D%20model_adj.pvalues%0A%0A%20%20%20%20_b_naive%20%3D%20model_naive.params%5B%22n_sector_2024%22%5D%0A%20%20%20%20_b_adj%20%3D%20model_adj.params%5B%22n_sector_2024%22%5D%0A%20%20%20%20print(f%22Sector%20coefficient%3A%20%7B_b_naive%3A.3f%7D%20(naive)%20%20-%3E%20%20%7B_b_adj%3A.3f%7D%20(with%20controls).%22)%0A%20%20%20%20print(%22Adding%20controls%20shrinks%20it.%20But%20these%20are%20example%20of%20BAD%20controls%20for%20causal%20inference!%22)%0A%20%20%20%20print(%22But%20the%20residual%20map%20is%20still%20spatially%20clustered%3A%20OLS%20SEs%20are%20too%22)%0A%20%20%20%20print(%22optimistic%2C%20and%20unmeasured%20spatial%20causes%20may%20remain.%22)%0A%20%20%20%20_coef_table.round(3)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%3E%20%23%23%23%23%20%F0%9F%94%A7%20Your%20turn%20(Step%201)%0A%20%20%20%20%3E%0A%20%20%20%20%3E%201.%20**Switch%20the%20sector.**%20In%20*marimo*%2C%20pick%20another%20sector%20from%20the%20dropdown%20at%0A%20%20%20%20%3E%20%20%20%20the%20top%20(e.g.%20Chemicals%2C%20Transport)%3B%20in%20*Jupyter*%2C%20change%0A%20%20%20%20%3E%20%20%20%20%60value%3D%22Livestock%20(NACE%20014x)%22%60%20in%20the%20dropdown%20cell%20to%20another%20option%20and%0A%20%20%20%20%3E%20%20%20%20re-run.%20Does%20the%20na%C3%AFve%20slope%20stay%20positive%20for%20every%20sector%3F%0A%20%20%20%20%3E%0A%20%20%20%20%3E%20*Five-question%20check%3A*%20which%20of%20the%20five%20questions%20did%20dropping%20a%20control%20touch%2C%0A%20%20%20%20%3E%20the%20**estimand**%2C%20the%20**credibility**%2C%20or%20the%20**failure%20mode**%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Step%202.%20Can%20spatial%20models%20fix%20the%20residual%20clustering%3F%0A%0A%20%20%20%20The%20OLS%20residuals%20are%20clustered%20in%20space.%20Three%20spatial%20models%20describe%20that%0A%20%20%20%20clustering%20in%20different%20ways%3A%20choose%20based%20on%20*where%20you%20think%20the%20spillover%20lives*%3A%0A%0A%20%20%20%20%7C%20Model%20%7C%20Extra%20term%20%7C%20Use%20it%20when%E2%80%A6%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20**SAR**%20(spatial%20lag)%20%7C%20%24y%20%3D%20%5Crho%5C%2CWy%20%2B%20%5Cbeta%20X%20%2B%20%5Cvarepsilon%24%20%7C%20the%20**outcome%20itself**%20spreads%20(a%20gas%20like%20NH%E2%82%83)%2C%20**our%20choice%20here**%20%7C%0A%20%20%20%20%7C%20**SEM**%20(spatial%20error)%20%7C%20%24y%20%3D%20%5Cbeta%20X%20%2B%20u%2C%5C%3B%5C%3B%20u%20%3D%20%5Clambda%20Wu%20%2B%20%5Cvarepsilon%24%20%7C%20you%20mainly%20want%20**honest%20SEs**%3B%20spatial%20structure%20is%20a%20nuisance%20%7C%0A%20%20%20%20%7C%20**SDM**%20(spatial%20Durbin)%20%7C%20%24y%20%3D%20%5Crho%5C%2CWy%20%2B%20%5Cbeta%20X%20%2B%20%5Ctheta%5C%2CWX%20%2B%20%5Cvarepsilon%24%20%7C%20your%20**neighbours'%20X**%20enter%20your%20equation%20directly%20%7C%0A%0A%20%20%20%20SAR%20is%20the%20special%20case%20of%20SDM%20with%20%24%5Ctheta%20%3D%200%24.%20**SEM%20is%20*not*%20just%20a%20coefficient%0A%20%20%20%20restriction%20of%20SDM**%3A%20the%20dependence%20sits%20in%20the%20**error%20term**%20(%24%5Clambda%20Wu%24).%0A%20%20%20%20We%20**adopt%20SAR**%20here%3A%20NH%E2%82%83%20is%20a%20gas%2C%20so%20the%20**outcome**%20drifts%20between%0A%20%20%20%20cells%20(%24%5Crho%20Wy%24).%0A%0A%20%20%20%20**What%20is%20W%3F**%20The%20spatial%20weights%20matrix%20says%20*who%20is%20a%20neighbour%20of%20whom%2C%20and%20how%0A%20%20%20%20strongly*.%20Here%20%60DistanceBand(threshold%3D3000%2C%20alpha%3D-3)%60%20%3D%20neighbours%20within%20**3%20km**%2C%0A%20%20%20%20weighted%20by%20**1%2Fdistance%C2%B3**%3B%20%60transform%3D%22R%22%60%20**row-standardises**%20it%20so%20each%20cell's%0A%20%20%20%20weights%20sum%20to%201%2C%20making%20%24Wy%24%20a%20neighbour%20**average**.%20W%20is%20a%20*choice%20you%20make*%2C%20not%0A%20%20%20%20data%2C%20and%20a%20different%20W%20can%20change%20the%20story.%0A%0A%20%20%20%20We%20fit%20on%20a%20small%20**contiguous%20eastern%20subset**%3A%20fitting%20spatial%20models%20is%20computationally%20heav.%0A%0A%20%20%20%20**The%20five%20questions%20for%20this%20step**%0A%0A%20%20%20%20%7C%20Question%20%7C%20Answer%20for%20Step%202%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20**treatment**%20%7C%20same%20firm%20count%20as%20Step%201%20%7C%0A%20%20%20%20%7C%20**estimand**%20%7C%20the%20same%20association%2C%20now%20with%20spatially%20honest%20SEs%20%7C%0A%20%20%20%20%7C%20**comparison**%20%7C%20the%20same%20cross-section%2C%20modelling%20spatial%20dependence%20%7C%0A%20%20%20%20%7C%20**credible%20if**%20%7C%20the%20spatial%20model%20is%20the%20right%20description%20of%20the%20dependence%2C%20no%20unmeasured%20confounders%20%7C%0A%20%20%20%20%7C%20**fails%20when**%20%7C%20dependence%20%E2%89%A0%20confounding%3B%20unmeasured%20spatial%20causes%20remain%20%7C%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20DistanceBand%2C%0A%20%20%20%20ML_Error%2C%0A%20%20%20%20ML_Lag%2C%0A%20%20%20%20Moran%2C%0A%20%20%20%20fit_spatial%2C%0A%20%20%20%20np%2C%0A%20%20%20%20pd%2C%0A%20%20%20%20reg_data%2C%0A%20%20%20%20smf%2C%0A)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%202a.%20Fit%20OLS%2C%20SAR%2C%20SEM%2C%20SDM%20on%20a%20contiguous%20subset%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20sample%20%3D%20reg_data.loc%5B(reg_data%5B%22grid_row%22%5D%20%3E%20200)%20%26%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20(reg_data%5B%22grid_col%22%5D%20%3E%20150)%5D.copy().reset_index(drop%3DTrue)%0A%20%20%20%20print(f%22Geographic%20subset%3A%20%7Blen(sample)%3A%2C%7D%20cells%22)%0A%0A%20%20%20%20_w%20%3D%20DistanceBand(sample%5B%5B%22x_rd%22%2C%20%22y_rd%22%5D%5D.values%2C%20threshold%3D3000%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20binary%3DFalse%2C%20alpha%3D-3%2C%20silence_warnings%3DTrue)%0A%20%20%20%20_w.transform%20%3D%20%22R%22%0A%0A%20%20%20%20_ols_s%20%3D%20smf.ols(%22nh3_2024%20~%20n_sector_2024%20%2B%20log_pop_density%20%2B%20C(rurality)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20data%3Dsample).fit(cov_type%3D%22HC3%22)%20%23heteroskedasticity-robust%20standard%20errors.%0A%20%20%20%20sample%5B%22ols_resid%22%5D%20%3D%20_ols_s.resid%0A%20%20%20%20_moran_ols%20%3D%20Moran(sample%5B%22ols_resid%22%5D.values%2C%20_w)%0A%0A%20%20%20%20_rur_dum%20%3D%20pd.get_dummies(sample%5B%22rurality%22%5D.astype(str)%2C%20prefix%3D%22rur%22%2C%20drop_first%3DTrue)%0A%20%20%20%20_X%20%3D%20pd.concat(%5Bsample%5B%5B%22n_sector_2024%22%2C%20%22log_pop_density%22%5D%5D%2C%20_rur_dum%5D%2C%20axis%3D1).astype(float)%0A%20%20%20%20_y%20%3D%20sample%5B%5B%22nh3_2024%22%5D%5D.values%0A%20%20%20%20_xn%20%3D%20list(_X.columns)%0A%0A%20%20%20%20%23%20We%20use%20SAR%3B%20SEM%20and%20SDM%20are%20shown%20for%20the%20menu.%0A%20%20%20%20m_sar%2C%20_sar_resid%2C%20_mr_lag%20%3D%20fit_spatial(ML_Lag%2C%20_y%2C%20_X.values%2C%20_w%2C%20%22nh3_2024%22%2C%20_xn%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20spat_impacts%3D%22full%22)%0A%20%20%20%20sample%5B%22sar_resid%22%5D%20%3D%20_sar_resid%0A%20%20%20%20_m_err%2C%20_%2C%20_mr_err%20%3D%20fit_spatial(ML_Error%2C%20_y%2C%20_X.values%2C%20_w%2C%20%22nh3_2024%22%2C%20_xn)%0A%20%20%20%20_m_sdm%2C%20_%2C%20_mr_sdm%20%3D%20fit_spatial(ML_Lag%2C%20_y%2C%20_X.values%2C%20_w%2C%20%22nh3_2024%22%2C%20_xn%2C%20slx_lags%3D1)%0A%0A%20%20%20%20def%20_liv_beta(model)%3A%0A%20%20%20%20%20%20%20%20return%20float(np.ravel(model.betas)%5B1%5D)%0A%0A%20%20%20%20_comparison%20%3D%20pd.DataFrame(%5B%0A%20%20%20%20%20%20%20%20%7B%22Model%22%3A%20%22OLS%22%2C%20%22sector%20coef%22%3A%20_ols_s.params%5B%22n_sector_2024%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%22resid%20Moran%20I%22%3A%20_moran_ols.I%2C%20%22spatial%20param%22%3A%20%22n%2Fa%22%7D%2C%0A%20%20%20%20%20%20%20%20%7B%22Model%22%3A%20%22SAR%22%2C%20%22sector%20coef%22%3A%20_liv_beta(m_sar)%2C%0A%20%20%20%20%20%20%20%20%20%22resid%20Moran%20I%22%3A%20_mr_lag.I%2C%20%22spatial%20param%22%3A%20f%22rho%20%3D%20%7Bfloat(m_sar.rho)%3A.2f%7D%22%7D%2C%0A%20%20%20%20%20%20%20%20%7B%22Model%22%3A%20%22SEM%22%2C%20%22sector%20coef%22%3A%20_liv_beta(_m_err)%2C%0A%20%20%20%20%20%20%20%20%20%22resid%20Moran%20I%22%3A%20_mr_err.I%2C%20%22spatial%20param%22%3A%20f%22lam%20%3D%20%7Bfloat(_m_err.lam)%3A.2f%7D%22%7D%2C%0A%20%20%20%20%20%20%20%20%7B%22Model%22%3A%20%22SDM%22%2C%20%22sector%20coef%22%3A%20_liv_beta(_m_sdm)%2C%0A%20%20%20%20%20%20%20%20%20%22resid%20Moran%20I%22%3A%20_mr_sdm.I%2C%20%22spatial%20param%22%3A%20f%22rho%20%3D%20%7Bfloat(_m_sdm.rho)%3A.2f%7D%22%7D%2C%0A%20%20%20%20%5D).set_index(%22Model%22)%0A%0A%20%20%20%20print(%22Moran's%20I%20drops%20toward%200%20(clustering%20absorbed).%20The%20sector%20coef%20stays%20clearly%22)%0A%20%20%20%20print(%22positive%20-%3E%20a%20spatially%20honest%20CORRELATION%2C%20still%20not%20a%20cause%3A%20the%20models%22)%0A%20%20%20%20print(%22cannot%20separate%20'farms%20raise%20NH3'%20from%20'farms%20and%20NH3%20share%20spatial%20causes'.%22)%0A%20%20%20%20_comparison.round(3)%0A%20%20%20%20return%20m_sar%2C%20sample%0A%0A%0A%40app.cell%0Adef%20_(plt%2C%20sample%2C%20style_map_axes)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%202b.%20Residual%20maps%3A%20OLS%20vs%20SAR%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(14%2C%206))%0A%20%20%20%20for%20_ax%2C%20_col%2C%20_title%20in%20zip(%0A%20%20%20%20%20%20%20%20_axes%2C%20%5B%22ols_resid%22%2C%20%22sar_resid%22%5D%2C%0A%20%20%20%20%20%20%20%20%5B%22OLS%20residuals%20(strong%20spatial%20clustering)%22%2C%20%22SAR%20residuals%20(spatial%20dependence%20absorbed)%22%5D%2C%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20_sc%20%3D%20_ax.scatter(sample%5B%22x_rd%22%5D%2C%20sample%5B%22y_rd%22%5D%2C%20c%3Dsample%5B_col%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cmap%3D%22coolwarm%22%2C%20s%3D5%2C%20vmin%3D-5%2C%20vmax%3D5)%0A%20%20%20%20%20%20%20%20_ax.set_title(_title%2C%20fontsize%3D12)%0A%20%20%20%20style_map_axes(*_axes)%0A%20%20%20%20_fig.colorbar(_sc%2C%20ax%3D_axes%2C%20shrink%3D0.7%2C%20label%3D%22Residual%20(ug%2Fm3)%22)%0A%20%20%20%20_fig.suptitle(%22Spatial%20models%20absorb%20the%20clustering%20in%20OLS%20residuals%22%2C%20fontsize%3D13%2C%20y%3D1.02)%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(m_sar%2C%20spatial_impacts)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%202c.%20SAR%20effects%3A%20direct%2C%20indirect%2C%20total%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_impacts%20%3D%20spatial_impacts(m_sar)%0A%20%20%20%20_rho%20%3D%20float(m_sar.rho)%0A%20%20%20%20print(f%22Outcome%20feedback%20rho%20%3D%20%7B_rho%3A.2f%7D%20%20-%3E%20%20spatial%20multiplier%201%2F(1-rho)%20%3D%20%7B1%2F(1-_rho)%3A.2f%7D%22)%0A%0A%20%20%20%20_liv%20%3D%20_impacts.loc%5B%22n_sector_2024%22%5D%0A%20%20%20%20print(f%22%5CnSector%20effect%20on%20NH3%20(ug%2Fm3)%3A%20direct%20%7B_liv%5B'direct'%5D%3A.2f%7D%2C%20%22%0A%20%20%20%20%20%20%20%20%20%20f%22indirect%20%7B_liv%5B'indirect'%5D%3A.2f%7D%2C%20total%20%7B_liv%5B'total'%5D%3A.2f%7D.%22)%0A%20%20%20%20print(%22%20%20direct%20%20%20%3D%20effect%20on%20a%20cell's%20OWN%20NH3%22)%0A%20%20%20%20print(%22%20%20indirect%20%3D%20spillover%20to%20OTHER%20cells%20through%20the%20outcome%20feedback%20(rho%20*%20W%20y)%22)%0A%20%20%20%20print(%22%20%20total%20%20%20%20%3D%20direct%20%2B%20indirect%22)%0A%20%20%20%20print(%22NOTE%3A%20this%20is%20a%20DIFFERENT%20ESTIMAND%20from%20Step%204.%20Here%20the%20total%20is%20a%20cross-sectional%22)%0A%20%20%20%20print(%22LEVEL%20association%20(confounded%20by%20time-invariant%20place%20traits)%3B%20Step%204%20first-differences%22)%0A%20%20%20%20print(%22those%20away%2C%20so%20its%20total%20is%20a%20change-on-change%20effect.%20The%20two%20totals%20are%20not%20comparable.%22)%0A%20%20%20%20_impacts.head(1).round(3)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%3E%20%23%23%23%23%20%F0%9F%94%A7%20Your%20turn%20(Step%202)%0A%20%20%20%20%3E%0A%20%20%20%20%3E%20**W%20is%20a%20choice%2C%20not%20data.**%20In%20cell%20*2a*%2C%20change%20the%20weights%3A%0A%20%20%20%20%3E%20%60threshold%3D3000%60%20%E2%86%92%20%605000%60%2C%20and%2For%20%60alpha%3D-3%60%20%E2%86%92%20%60-1%60%20(gentler%20distance%20decay).%0A%20%20%20%20%3E%20Re-run%20and%20watch%20two%20things%20in%20the%20comparison%20table%3A%0A%20%20%20%20%3E%0A%20%20%20%20%3E%201.%20the%20spatial%20parameter%20**%CF%81**%20(does%20a%20wider%20%2F%20gentler%20W%20raise%20it%3F)%3B%0A%20%20%20%20%3E%202.%20the%20**residual%20Moran's%20I**%20(does%20it%20still%20collapse%20toward%200%3F).%0A%20%20%20%20%3E%0A%20%20%20%20%3E%20*Five-question%20check%3A*%20a%20different%20W%20gives%20a%20different%20%CF%81%20and%20a%20different%20%22total%22.%0A%20%20%20%20%3E%20Which%20of%20the%20five%20questions%20does%20that%20uncertainty%20live%20in%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Step%203.%20Can%20farm%20gains%20identify%20a%20causal%20effect%3F%20(DiD)%0A%0A%20%20%20%20Cross-sectional%20analysis%20can't%20separate%20causation%20from%20correlation.%20We%20need%20a%0A%20%20%20%20**shock**%20that%20changes%20exposure%20over%20time.%20Between%202020%20and%202024%20some%20cells%20gained%0A%20%20%20%20**more%20than%20one**%20firm%2C%20and%20if%20firms%20drive%20NH%E2%82%83%20those%20cells'%20NH%E2%82%83%20should%20*rise%20relative%0A%20%20%20%20to*%20cells%20whose%20firm%20count%20stayed%20flat%20or%20fell.%0A%0A%20%20%20%20The%20**DiD%20estimate**%20is%20the%20*extra*%20change%20in%20treated%20cells%20beyond%20the%20control%20trend%3A%0A%0A%20%20%20%20%24%24%5Ctext%7BDiD%7D%20%3D%20(%5Cbar%7BY%7D%5E%7B%5C%2Cpost%7D_%7Btreated%7D%20-%20%5Cbar%7BY%7D%5E%7B%5C%2Cpre%7D_%7Btreated%7D)%20-%20(%5Cbar%7BY%7D%5E%7B%5C%2Cpost%7D_%7Bcontrol%7D%20-%20%5Cbar%7BY%7D%5E%7B%5C%2Cpre%7D_%7Bcontrol%7D)%24%24%0A%0A%20%20%20%20%7C%20Group%20%7C%20Definition%20%7C%0A%20%20%20%20%7C-------%7C-----------%7C%0A%20%20%20%20%7C%20**Treated**%20%7C%20has%20firms%20in%202024%2C%20**gained%20%3E%201**%20since%202020%20%7C%0A%20%20%20%20%7C%20**Control**%20%7C%20has%20firms%20in%202024%2C%20**no%20gain**%20(net%20change%20%E2%89%A4%200)%20%7C%0A%0A%20%20%20%20Note%20that%20%22**no%20gain**%22%20lumps%20together%20stable%20cells%20**and%20cells%20that%20lost%20farms**%3B%20if%0A%20%20%20%20losing%20farms%20lowers%20NH%E2%82%83%2C%20the%20control%20trend%20tilts%20down%20and%20the%20DiD%20is%20inflated%20a%20little.%0A%0A%20%20%20%20**Key%20assumption%2C%20parallel%20trends%3A**%20absent%20the%20gains%2C%20treated%20and%20control%20cells%0A%20%20%20%20would%20have%20followed%20the%20**same**%20NH%E2%82%83%20trend.%20We%20check%20it%20with%20a%20**pre-trend%20test**%0A%20%20%20%20(2018%E2%80%932020)%3B%20a%20significant%20pre-trend%20means%20the%20estimate%20is%20biased%2C%20a%20real%0A%20%20%20%20limitation%2C%20not%20a%20bug.%0A%0A%20%20%20%20**A%20note%20on%20power.**%20With%20only%20~130%20treated%20cells%20the%20DiD%20is%20**not%20significant**%0A%20%20%20%20(p%20%E2%89%88%200.07).%0A%0A%20%20%20%20**The%20five%20questions%20for%20this%20step**%0A%0A%20%20%20%20%7C%20Question%20%7C%20Answer%20for%20Step%203%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20**treatment**%20%7C%20net%20farm%20gain%20in%20a%20cell%2C%202020%20%E2%86%92%202024%20%7C%0A%20%20%20%20%7C%20**estimand**%20%7C%20Effect%20of%20gaining%201%2B%20firms%20on%20cell%20NH%E2%82%83%20%7C%0A%20%20%20%20%7C%20**comparison**%20%7C%20gain%20vs%20no-gain%20%C3%97%20pre%2Fpost%20%7C%0A%20%20%20%20%7C%20**credible%20if**%20%7C%20parallel%20trends%20hold%20%7C%0A%20%20%20%20%7C%20**fails%20when**%20%7C%20NH%E2%82%83%20drifts%20to%20controls%20(interference%20%E2%86%92%20Step%204)%20%7C%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(plt)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20How%20Difference-in-Differences%20works%20(schematic%2C%20not%20real%20data)%20%E2%94%80%E2%94%80%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(9%2C%205))%0A%20%20%20%20_pre%2C%20_post%20%3D%200%2C%201%0A%20%20%20%20_c0%2C%20_c1%20%3D%209.1%2C%207.6%20%20%20%20%20%20%20%20%20%20%23%20control%3A%20pre%20-%3E%20post%0A%20%20%20%20_t0%2C%20_t1%20%3D%2011.1%2C%2010.1%20%20%20%20%20%20%20%20%23%20treated%3A%20pre%20-%3E%20post%20(observed)%0A%20%20%20%20_cf%20%3D%20_t0%20%2B%20(_c1%20-%20_c0)%20%20%20%20%20%20%23%20treated%20COUNTERFACTUAL%0A%0A%20%20%20%20_ax.plot(%5B_pre%2C%20_post%5D%2C%20%5B_c0%2C%20_c1%5D%2C%20%22o-%22%2C%20color%3D%22%232b8cbe%22%2C%20lw%3D2.5%2C%20ms%3D10%2C%20label%3D%22Control%22)%0A%20%20%20%20_ax.plot(%5B_pre%2C%20_post%5D%2C%20%5B_t0%2C%20_t1%5D%2C%20%22s-%22%2C%20color%3D%22%23d95f0e%22%2C%20lw%3D2.5%2C%20ms%3D10%2C%20label%3D%22Treated%20(observed)%22)%0A%20%20%20%20_ax.plot(%5B_pre%2C%20_post%5D%2C%20%5B_t0%2C%20_cf%5D%2C%20%22s--%22%2C%20color%3D%22grey%22%2C%20lw%3D2%2C%20ms%3D10%2C%20mfc%3D%22white%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22Treated%20counterfactual%5Cn(control%20trend%20applied%20to%20treated)%22)%0A%20%20%20%20_ax.annotate(%22%22%2C%20xy%3D(_post%2C%20_t1)%2C%20xytext%3D(_post%2C%20_cf)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20arrowprops%3Ddict(arrowstyle%3D%22%3C-%3E%22%2C%20color%3D%22black%22%2C%20lw%3D2))%0A%20%20%20%20_ax.text(_post%20%2B%200.03%2C%20(_t1%20%2B%20_cf)%20%2F%202%2C%20f%22DiD%20%3D%20%7B_t1%20-%20_cf%3A%2B.1f%7D%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20va%3D%22center%22%2C%20fontsize%3D13%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20_ax.set_xticks(%5B_pre%2C%20_post%5D)%0A%20%20%20%20_ax.set_xticklabels(%5B%222020%20(pre)%22%2C%20%222024%20(post)%22%5D)%0A%20%20%20%20_ax.set_ylabel(%22Mean%20NH3%20(ug%2Fm3)%22)%0A%20%20%20%20_ax.set_title(%22DiD%20%3D%20the%20extra%20change%20in%20treated%20cells%2C%20beyond%20the%20control%20trend%22)%0A%20%20%20%20_ax.set_xlim(-0.15%2C%201.45)%0A%20%20%20%20_ax.legend(loc%3D%22lower%20left%22%2C%20fontsize%3D10)%0A%20%20%20%20_ax.spines%5B%5B%22top%22%2C%20%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(grid%2C%20sector_label)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%203a.%20Define%20treated%20and%20control%20cells%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%23%20Keep%20only%20cells%20that%20HAVE%20firms%20in%202024%2C%20then%20split%20by%20the%202020-%3E2024%20shock%3A%0A%20%20%20%20%23%20%20%20%20%20treated%20%3D%20gained%20MORE%20than%201%20%20%20(sector_net_change%20%3E%201)%0A%20%20%20%20%23%20%20%20%20%20control%20%3D%20no%20gain%20%20%20%20%20%20%20%20%20%20%20%20%20%20(sector_net_change%20%3C%3D%200)%0A%20%20%20%20%23%20Cells%20with%20a%20net%20change%20of%20exactly%20%2B1%20are%20dropped%20as%20an%20ambiguous%20middle%20group.%0A%20%20%20%20did_cells%20%3D%20grid%5Bgrid%5B%22n_sector_2024%22%5D%20%3E%200%5D.copy()%0A%20%20%20%20did_cells%5B%22treated%22%5D%20%3D%20(did_cells%5B%22sector_net_change%22%5D%20%3E%201).astype(int)%0A%20%20%20%20did_cells%20%3D%20did_cells%5B(did_cells%5B%22sector_net_change%22%5D%20%3C%3D%200)%20%7C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20(did_cells%5B%22sector_net_change%22%5D%20%3E%201)%5D.copy()%0A%0A%20%20%20%20_n_t%20%3D%20int(did_cells%5B%22treated%22%5D.sum())%0A%20%20%20%20_n_c%20%3D%20int((did_cells%5B%22treated%22%5D%20%3D%3D%200).sum())%0A%20%20%20%20print(f%22Treated%20(gained%20%3E%201%20%7Bsector_label%7D%2C%202020-2024)%3A%20%7B_n_t%3A%2C%7D%20cells%22)%0A%20%20%20%20print(f%22Control%20(no%20gain%2C%20net%20change%20%3C%3D%200)%3A%20%20%20%20%20%20%20%20%20%20%20%20%20%7B_n_c%3A%2C%7D%20cells%22)%0A%20%20%20%20return%20(did_cells%2C)%0A%0A%0A%40app.cell%0Adef%20_(SECTOR%2C%20did_cells%2C%20plt)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%203b.%20Parallel-trends%20check%20%2B%20build%20the%20long%20panel%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_years%20%3D%20list(range(2018%2C%202025))%0A%20%20%20%20_nh3_cols%20%3D%20%5Bf%22nh3_%7B_y%7D%22%20for%20_y%20in%20_years%5D%0A%20%20%20%20_means%20%3D%20did_cells.groupby(%22treated%22)%5B_nh3_cols%5D.mean()%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(9%2C%205))%0A%20%20%20%20_ax.plot(_years%2C%20_means.loc%5B0%5D%2C%20%22o-%22%2C%20color%3D%22%232b8cbe%22%2C%20lw%3D2%2C%20ms%3D7%2C%20label%3D%22Control%20(no%20gain)%22)%0A%20%20%20%20_ax.plot(_years%2C%20_means.loc%5B1%5D%2C%20%22s-%22%2C%20color%3D%22%23d95f0e%22%2C%20lw%3D2%2C%20ms%3D7%2C%20label%3D%22Treated%20(gained%20%3E%201)%22)%0A%20%20%20%20_ax.axvline(2020.5%2C%20color%3D%22black%22%2C%20ls%3D%22--%22%2C%20lw%3D1.5%2C%20label%3D%22Treatment%20window%22)%0A%20%20%20%20_ax.axvspan(2020.5%2C%202024.5%2C%20alpha%3D0.07%2C%20color%3D%22grey%22)%0A%20%20%20%20_ax.set_xlabel(%22Year%22)%0A%20%20%20%20_ax.set_ylabel(%22Mean%20NH3%20(ug%2Fm3)%22)%0A%20%20%20%20_ax.set_title(f%22Parallel%20trends%20check%3A%20%7BSECTOR%7D%22)%0A%20%20%20%20_ax.legend(fontsize%3D11)%0A%20%20%20%20_fig.tight_layout()%0A%0A%20%20%20%20panel%20%3D%20did_cells.melt(id_vars%3D%5B%22grid_row%22%2C%20%22grid_col%22%2C%20%22treated%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20value_vars%3D_nh3_cols%2C%20var_name%3D%22year%22%2C%20value_name%3D%22nh3%22)%0A%20%20%20%20panel%5B%22year%22%5D%20%3D%20panel%5B%22year%22%5D.str.replace(%22nh3_%22%2C%20%22%22).astype(int)%0A%20%20%20%20panel%5B%22cell_id%22%5D%20%3D%20panel%5B%22grid_row%22%5D.astype(str)%20%2B%20%22_%22%20%2B%20panel%5B%22grid_col%22%5D.astype(str)%0A%20%20%20%20_fig%0A%20%20%20%20return%20(panel%2C)%0A%0A%0A%40app.cell%0Adef%20_(panel%2C%20pd%2C%20smf)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%203c.%20The%20DiD%20regression%20(%2B%20pre-trend%20test)%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_did_data%20%3D%20panel%5Bpanel%5B%22year%22%5D.isin(%5B2020%2C%202024%5D)%5D.copy()%0A%20%20%20%20_did_data%5B%22post%22%5D%20%3D%20(_did_data%5B%22year%22%5D%20%3D%3D%202024).astype(int)%0A%20%20%20%20_model_did%20%3D%20smf.ols(%22nh3%20~%20treated%20*%20post%22%2C%20data%3D_did_data).fit(%0A%20%20%20%20%20%20%20%20cov_type%3D%22cluster%22%2C%20cov_kwds%3D%7B%22groups%22%3A%20_did_data%5B%22cell_id%22%5D%7D)%0A%20%20%20%20_did_est%20%3D%20_model_did.params%5B%22treated%3Apost%22%5D%0A%20%20%20%20_did_p%20%3D%20_model_did.pvalues%5B%22treated%3Apost%22%5D%0A%0A%20%20%20%20%23%20Test%20parallel%20slopes%20assumption%0A%20%20%20%20_pre%20%3D%20panel%5Bpanel%5B%22year%22%5D%20%3C%3D%202020%5D.copy()%0A%20%20%20%20_pre%5B%22year_idx%22%5D%20%3D%20_pre%5B%22year%22%5D%20-%202018%0A%20%20%20%20_m_pre%20%3D%20smf.ols(%22nh3%20~%20treated%20*%20year_idx%22%2C%20data%3D_pre).fit(%0A%20%20%20%20%20%20%20%20cov_type%3D%22cluster%22%2C%20cov_kwds%3D%7B%22groups%22%3A%20_pre%5B%22cell_id%22%5D%7D)%0A%0A%20%20%20%20_did_table%20%3D%20pd.DataFrame(%5B%0A%20%20%20%20%20%20%20%20%7B%22Test%22%3A%20%22DiD%20(treated%20x%20post)%22%2C%20%22Estimate%22%3A%20_did_est%2C%20%22p-value%22%3A%20_did_p%7D%2C%0A%20%20%20%20%20%20%20%20%7B%22Test%22%3A%20%22Pre-trend%20(treated%20x%20year)%22%2C%20%22Estimate%22%3A%20_m_pre.params%5B%22treated%3Ayear_idx%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%22p-value%22%3A%20_m_pre.pvalues%5B%22treated%3Ayear_idx%22%5D%7D%2C%0A%20%20%20%20%5D).set_index(%22Test%22)%0A%0A%20%20%20%20if%20_m_pre.pvalues%5B%22treated%3Ayear_idx%22%5D%20%3C%200.05%3A%0A%20%20%20%20%20%20%20%20print(%22WARNING%3A%20significant%20pre-trend%20-%3E%20parallel%20trends%20is%20doubtful%20here.%22)%0A%20%20%20%20%20%20%20%20print(%22Treated%20cells%20were%20already%20on%20a%20different%20path%20before%20the%20gains%2C%22)%0A%20%20%20%20%20%20%20%20print(%22so%20the%20DiD%20estimate%20is%20likely%20biased%20(a%20real%20limitation%2C%20not%20a%20bug).%22)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20print(%22Pre-trend%20not%20significant%3A%20parallel%20trends%20is%20at%20least%20plausible.%22)%0A%20%20%20%20print(%22%5CnStandard%2Frobust%20DiD%20alternative%3A%20pyfixest%20feols(...%20%7C%20cell_id%20%2B%20year).%22)%0A%20%20%20%20_did_table.round(3)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%3E%20%23%23%23%23%20%F0%9F%94%A7%20Your%20turn%20(Step%203)%0A%20%20%20%20%3E%0A%20%20%20%20%3E%20**Redefine%20the%20treatment.**%20In%20cell%20*3a*%2C%20change%20the%20cutoff%20from%20%60%3E%201%60%20to%20%60%3E%200%60%0A%20%20%20%20%3E%20(treat%20*any*%20gain%20as%20treatment).%20The%20%60%3E%201%60%20appears%20twice%2C%20in%20the%20%60treated%20%3D%20...%60%0A%20%20%20%20%3E%20line%20and%20in%20the%20filter%20line%20just%20below%2C%20change%20**both**.%20Re-run%20and%20compare%3A%20the%0A%20%20%20%20%3E%20treated-cell%20count%2C%20the%20DiD%20estimate%2C%20and%20its%20**p-value**.%20Change%20also%20the%20controls.%0A%20%20%20%20%3E%0A%20%20%20%20%3E%20A%20looser%20definition%20adds%20more%20(smaller-gain)%20treated%20cells.%20Does%20the%20estimate%20get%0A%20%20%20%20%3E%20bigger%20or%20smaller%3F%20More%20or%20less%20significant%3F%0A%20%20%20%20%3E%0A%20%20%20%20%3E%20P-HACKING%20WARNING%3A%20Do%20not%20do%20this%20in%20your%20own%20research!%20Pre-register%20your%20analysis.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Step%204.%20What%20happens%20when%20pollution%20spills%20across%20cells%3F%0A%0A%20%20%20%20Standard%20DiD%20assumes%20controls%20are%20**not**%20affected%20by%20treatment%20(no%20interference).%0A%20%20%20%20But%20NH%E2%82%83%20is%20a%20gas%2C%20it%20**drifts**.%20A%20gain%20in%20one%20cell%20can%20raise%20NH%E2%82%83%20in%20its%0A%20%20%20%20neighbours%2C%20**contaminating%20the%20controls**%20and%20biasing%20the%20DiD%20toward%20zero.%0A%0A%20%20%20%20**First%2C%20why%20is%20a%20first%20difference%20a%20DiD%3F**%20With%20only%20two%20periods%2C%20take%20each%20cell's%0A%20%20%20%20change%20%24%5CDelta%5Ctext%7BNH%7D_3%20%3D%20%5Ctext%7BNH%7D_3%5E%7B2024%7D%20-%20%5Ctext%7BNH%7D_3%5E%7B2020%7D%24.%20Differencing%0A%20%20%20%20cancels%20everything%20time-invariant%20about%20a%20cell%20(its%20location%2C%20soil%2C%20baseline%20level%3A%0A%20%20%20%20the%20cell%20fixed%20effect).%20Regress%20that%20change%20on%20the%20treated%20dummy%3A%0A%0A%20%20%20%20%24%24%5CDelta%5Ctext%7BNH%7D_3%7B%7D_%7B%2Ci%7D%20%3D%20%5Calpha%20%2B%20%5Ctau%5C%2CT_i%20%2B%20%5Cvarepsilon_i%24%24%0A%0A%20%20%20%20The%20intercept%20%24%5Calpha%24%20is%20the%20**average%20change%20in%20control%20cells**%2C%20the%20common%20trend%2C%0A%20%20%20%20i.e.%20our%20counterfactual%20for%20what%20treated%20cells%20would%20have%20done%20anyway%20(%60post%60).%20%24%5Ctau%24%20is%20the%0A%20%20%20%20**extra**%20change%20in%20treated%20cells%3A%20this%20is%20*exactly*%20the%20difference-in-differences%0A%20%20%20%20estimate%20(algebraically%20identical%20to%20the%20%60treated%3Apost%60%20term%20in%20Step%203).%0A%20%20%20%20The%20%60treated%60%20term%20dissapears%20with%20the%20difference.%0A%20%20%20%20The%20spatial-lag%20(SAR)%20model%20then%20**adds%20one%20spatial%20term%20on%20top%20of%20that%20DiD**%3A%0A%20%20%20%20the%20outcome%20feedback%20%24%5Crho%5C%2CW%5CDelta%20Y%24.%0A%0A%20%20%20%20**Option%201.%20Model%20the%20spillover**%20(this%20notebook).%20Fit%20a%20**spatial-lag%20(SAR)**%20model%20on%20the%0A%20%20%20%20**first%20difference**%20(%24%5CDelta%5Ctext%7BNH%7D_3%20%3D%202024%20-%202020%24)%3A%0A%0A%20%20%20%20%24%24%5CDelta%20Y_i%20%3D%20%5Calpha%20%2B%20%5Crho%5C%2C(W%5CDelta%20Y)_i%20%2B%20%5Ctau%5C%2CT_i%20%2B%20%5Cvarepsilon_i%24%24%0A%0A%20%20%20%20-%20%24%5Ctau%24%3A%20my%20**own**%20farm%20gain%3B%20a%20neighbour's%20gain%20reaches%20me%20**only%20through%20the%20outcome**%2C%20since%20their%20farm%20raises%20*their*%20%24%5CDelta%5Ctext%7BNH%7D_3%24%2C%20which%20drifts%20to%20me%20via%20%24%5Crho%5C%2C(W%5CDelta%20Y)%24%0A%20%20%20%20-%20the%20%24%5Crho%5C%2C(W%5CDelta%20Y)%24%20term%20feeds%20each%20cell's%20change%20back%20through%20its%20neighbours%2C%20so%20the%20**total**%20effect%20carries%20the%20spatial%20multiplier%20%241%2F(1-%5Crho)%24.%0A%20%20%20%20-%20*Cost%3A*%20the%20answer%20is%20only%20as%20good%20as%20%24W%24%2C%20and%20a%20high%20%24%5Crho%24%20**amplifies**%20everything.%0A%20%20%20%20-%20*Identifying%20assumption%3A*%20SAR%20loads%20**all**%20the%20spatial%20structure%20of%20%24%5CDelta%5Ctext%7BNH%7D_3%24%20onto%20outcome%20feedback%20%24%5Crho%5C%2C(W%5CDelta%20Y)%24.%20If%20part%20of%20it%20is%20instead%20**common%20spatial%20shocks**%20(weather%2C%20regional%20trends%20that%20hit%20neighbours%20alike)%2C%20%24%5Crho%24%20absorbs%20those%20too%20and%20the%20**total**%20is%20overstated.%0A%0A%20%20%20%20**Option%202.%20Redesign%20with%20distance%20rings**%20(the%20lecture's%20alternative).%20Keep%20the%20treated%0A%20%20%20%20cells%2C%20**drop%20%22controls%22%20near%20a%20treated%20cell**%20(partly%20exposed)%2C%20and%20compare%20only%20to%0A%20%20%20%20**far**%20controls.%20The%20assumption%3A%20spillovers%20vanish%20beyond%20the%20buffer.%0A%0A%20%20%20%20**The%20five%20questions%20for%20this%20step**%0A%0A%20%20%20%20%7C%20Question%20%7C%20Answer%20for%20Step%204%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20**treatment**%20%7C%20own%20gain%3B%20the%20spillover%20enters%20only%20through%20the%20outcome%20%24%5Crho%5C%2C(W%5CDelta%20Y)%24%20%7C%0A%20%20%20%20%7C%20**estimand**%20%7C%20direct%20%2B%20spillover%20%3D%20**total**%20effect%20%7C%0A%20%20%20%20%7C%20**comparison**%20%7C%20treated%20vs%20control%20change%2C%20with%20outcome%20feedback%20across%20cells%20%7C%0A%20%20%20%20%7C%20**credible%20if**%20%7C%20%24W%24%20and%20the%20spillover%20range%20are%20right%20%7C%0A%20%20%20%20%7C%20**fails%20when**%20%7C%20the%20wrong%20%24W%24%20inflates%20the%20%241%2F(1-%5Crho)%24%20multiplier%20%7C%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20DistanceBand%2C%0A%20%20%20%20ML_Lag%2C%0A%20%20%20%20did_cells%2C%0A%20%20%20%20fit_spatial%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20pd%2C%0A%20%20%20%20smf%2C%0A%20%20%20%20spatial_impacts%2C%0A)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%204a.%20Spillover-aware%20DiD%3A%20a%20first-difference%20spatial%20lag%20(SAR)%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20fd%20%3D%20did_cells.copy()%0A%20%20%20%20fd%5B%22delta_nh3%22%5D%20%3D%20fd%5B%22nh3_2024%22%5D%20-%20fd%5B%22nh3_2020%22%5D%0A%20%20%20%20fd%20%3D%20fd.dropna(subset%3D%5B%22delta_nh3%22%5D).copy()%0A%20%20%20%20fd%20%3D%20fd.loc%5Bfd%5B%22n_sector_2024%22%5D%20%3E%200%5D%0A%20%20%20%20fd%20%3D%20fd.loc%5B(fd%5B%22grid_row%22%5D%20%3E%2050)%20%26%20(fd%5B%22grid_col%22%5D%20%3E%20100)%5D.copy().reset_index(drop%3DTrue)%20%23%20reduce%20the%20grid%0A%20%20%20%20print(f%22Geographic%20subset%20for%20the%20spillover%20model%3A%20%7Blen(fd)%3A%2C%7D%20cells%22)%0A%0A%20%20%20%20m_fd_ols%20%3D%20smf.ols(%22delta_nh3%20~%20treated%22%2C%20data%3Dfd).fit(cov_type%3D%22HC3%22)%0A%20%20%20%20print(f%22Plain%20DiD%20(own%20cell%20only)%3A%20treated%20%3D%20%7Bm_fd_ols.params%5B'treated'%5D%3A%2B.2f%7D%20ug%2Fm3%20%22%0A%20%20%20%20%20%20%20%20%20%20f%22(p%20%3D%20%7Bm_fd_ols.pvalues%5B'treated'%5D%3A.3f%7D)%22)%0A%0A%20%20%20%20_w_fd%20%3D%20DistanceBand(fd%5B%5B%22x_rd%22%2C%20%22y_rd%22%5D%5D.values%2C%20threshold%3D3000%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20binary%3DFalse%2C%20alpha%3D-3%2C%20silence_warnings%3DTrue)%0A%20%20%20%20_w_fd.transform%20%3D%20%22R%22%0A%20%20%20%20X_fd%20%3D%20fd%5B%5B%22treated%22%5D%5D.astype(float)%0A%20%20%20%20y_fd%20%3D%20fd%5B%5B%22delta_nh3%22%5D%5D.values%0A%0A%20%20%20%20_m_sar_fd%2C%20_%2C%20_%20%3D%20fit_spatial(ML_Lag%2C%20y_fd%2C%20X_fd.values%2C%20_w_fd%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22delta_nh3%22%2C%20list(X_fd.columns)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20spat_impacts%3D%22full%22)%0A%0A%20%20%20%20def%20_stars(p)%3A%0A%20%20%20%20%20%20%20%20return%20%22***%22%20if%20p%20%3C%200.01%20else%20%22**%22%20if%20p%20%3C%200.05%20else%20%22*%22%20if%20p%20%3C%200.10%20else%20%22%22%0A%0A%20%20%20%20_nice%20%3D%20%7B%22CONSTANT%22%3A%20%22intercept%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%22treated%22%3A%20%22treated%20%20(own%20gain%2C%20tau)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%22W_delta_nh3%22%3A%20%22W*delta_NH3%20%20(rho%2C%20outcome%20feedback)%22%7D%0A%20%20%20%20_reg_table%20%3D%20pd.DataFrame(%7B%0A%20%20%20%20%20%20%20%20%22term%22%3A%20%5B_nice.get(n%2C%20n)%20for%20n%20in%20_m_sar_fd.name_x%5D%2C%0A%20%20%20%20%20%20%20%20%22coef%22%3A%20_m_sar_fd.betas.flatten()%2C%0A%20%20%20%20%20%20%20%20%22std%20err%22%3A%20np.ravel(_m_sar_fd.std_err)%2C%0A%20%20%20%20%20%20%20%20%22p-value%22%3A%20%5Bfloat(z%5B1%5D)%20for%20z%20in%20_m_sar_fd.z_stat%5D%2C%0A%20%20%20%20%7D)%0A%20%20%20%20_reg_table%5B%22sig%22%5D%20%3D%20_reg_table%5B%22p-value%22%5D.apply(_stars)%0A%0A%20%20%20%20_rho_fd%20%3D%20float(_m_sar_fd.rho)%0A%20%20%20%20_names%20%3D%20list(_m_sar_fd.name_x)%0A%20%20%20%20_tau%20%3D%20float(_m_sar_fd.betas.flatten()%5B_names.index(%22treated%22)%5D)%0A%20%20%20%20_impacts_treated%20%3D%20spatial_impacts(_m_sar_fd).loc%5B%5B%22treated%22%5D%5D%0A%20%20%20%20_direct%20%3D%20float(_impacts_treated%5B%22direct%22%5D.iloc%5B0%5D)%0A%20%20%20%20_total%20%3D%20float(_impacts_treated%5B%22total%22%5D.iloc%5B0%5D)%0A%0A%20%20%20%20print(f%22tau%20(own%20gain)%20%3D%20%7B_tau%3A.2f%7D%2C%20rho%20(outcome%20feedback)%20%3D%20%7B_rho_fd%3A.2f%7D%22)%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22**Spatial%20lag%20%2F%20SAR%20(first%20difference)%2C%20raw%20coefficients**%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22(sig%3A%20***%20p%3C0.01%20%20**%20p%3C0.05%20%20*%20p%3C0.10)%3A%22)%2C%0A%20%20%20%20%20%20%20%20_reg_table.round(3)%2C%0A%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22The%20treatment%20story%20is%20in%20two%20coefficients%3A%20**%CF%84%20%3D%20%7B_tau%3A.2f%7D**%20(my%20*own*%20gain)%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22and%20**%CF%81%20%3D%20%7B_rho_fd%3A.2f%7D**%20(outcome%20feedback%2C%20multiplier%201%2F(1%E2%88%92%CF%81)%20%3D%20%7B1%2F(1-_rho_fd)%3A.2f%7D).%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%5Cn%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22Raw%20coefficients%20are%20**not**%20the%20marginal%20effect%3A%20with%20feedback%20%CF%81%2C%20a%20gain%20ripples%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22out%20to%20neighbours%20and%20loops%20back.%20The%20**impact%20decomposition**%20converts%20(%CF%84%2C%20%CF%81)%20into%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22the%20effect%20of%20the%20treatment%20(the%20control's%20row%20is%20omitted)%3A%22%0A%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20_impacts_treated.round(3)%2C%0A%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22-%20**direct%20%3D%20%7B_direct%3A.2f%7D**%3A%20effect%20on%20a%20cell's%20*own*%20%CE%94NH%E2%82%83%2C%20the%20own%20gain%20%CF%84%20**plus**%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22the%20feedback%20that%20loops%20back%20to%20the%20same%20cell%2C%20so%20it%20is%20a%20little%20larger%20than%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22%CF%84%20%3D%20%7B_tau%3A.2f%7D%20(that%20is%20why%20it%20does%20**not**%20match%20the%20coefficient%20table).%20It%20is%2C%20however%2C%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22**below**%20the%20plain%20own-cell%20DiD%20(%7Bm_fd_ols.params%5B'treated'%5D%3A%2B.2f%7D)%3A%20adding%20%CF%81(W%CE%94Y)%20reassigns%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22part%20of%20what%20the%20plain%20regression%20loaded%20onto%20%60treated%60%20to%20the%20spatial-lag%20term.%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22-%20**indirect**%3A%20the%20**spillover**%2C%20my%20own%20gain%20recirculated%20to%20neighbours%20(and%20back)%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22through%20the%20outcome%20feedback%20%CF%81.%20With%20SAR%20this%20is%20the%20*only*%20spillover%20channel.%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22-%20**total%20%3D%20direct%20%2B%20indirect%20%3D%20%7B_total%3A.2f%7D**%3A%20the%20change%20if%20treatment%20rose%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22everywhere%20at%20once%2C%20modestly%20above%20the%20~0.5%20own-cell%20DiD.%22%0A%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20(fd%2C)%0A%0A%0A%40app.cell%0Adef%20_(SECTOR%2C%20fd%2C%20plt%2C%20style_map_axes)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%204b.%20Maps%3A%20the%20DiD%20design%2C%20and%20where%20delta_NH3%20actually%20moved%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fd_plot%20%3D%20fd.dropna(subset%3D%5B%22nh3_2024%22%2C%20%22nh3_2020%22%5D).copy()%0A%20%20%20%20_fd_plot%5B%22delta_nh3%22%5D%20%3D%20_fd_plot%5B%22nh3_2024%22%5D%20-%20_fd_plot%5B%22nh3_2020%22%5D%0A%0A%20%20%20%20_pad%20%3D%203000%0A%20%20%20%20_xlim%20%3D%20(_fd_plot%5B%22x_rd%22%5D.min()%20-%20_pad%2C%20_fd_plot%5B%22x_rd%22%5D.max()%20%2B%20_pad)%0A%20%20%20%20_ylim%20%3D%20(_fd_plot%5B%22y_rd%22%5D.min()%20-%20_pad%2C%20_fd_plot%5B%22y_rd%22%5D.max()%20%2B%20_pad)%0A%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(14%2C%206))%0A%20%20%20%20_ctrl%20%3D%20_fd_plot%5B_fd_plot%5B%22treated%22%5D%20%3D%3D%200%5D%0A%20%20%20%20_trt%20%3D%20_fd_plot%5B_fd_plot%5B%22treated%22%5D%20%3D%3D%201%5D%0A%20%20%20%20_axes%5B0%5D.scatter(_ctrl%5B%22x_rd%22%5D%2C%20_ctrl%5B%22y_rd%22%5D%2C%20c%3D%22%232b8cbe%22%2C%20s%3D4%2C%20alpha%3D0.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Control%20(%7Blen(_ctrl)%3A%2C%7D)%22)%0A%20%20%20%20_axes%5B0%5D.scatter(_trt%5B%22x_rd%22%5D%2C%20_trt%5B%22y_rd%22%5D%2C%20c%3D%22%23d95f0e%22%2C%20s%3D12%2C%20alpha%3D0.9%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Treated%20(%7Blen(_trt)%3A%2C%7D)%22)%0A%20%20%20%20_axes%5B0%5D.set_title(f%22DiD%20design%3A%20%7BSECTOR%7D%22)%0A%20%20%20%20_axes%5B0%5D.legend(markerscale%3D2%2C%20fontsize%3D11)%0A%0A%20%20%20%20_sc%20%3D%20_axes%5B1%5D.scatter(_fd_plot%5B%22x_rd%22%5D%2C%20_fd_plot%5B%22y_rd%22%5D%2C%20c%3D_fd_plot%5B%22delta_nh3%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cmap%3D%22RdBu_r%22%2C%20s%3D4%2C%20vmin%3D-3%2C%20vmax%3D3)%0A%20%20%20%20_axes%5B1%5D.set_title(%22Change%20in%20NH3%2C%202024%20-%202020%22)%0A%20%20%20%20_fig.colorbar(_sc%2C%20ax%3D_axes%5B1%5D%2C%20shrink%3D0.7%2C%20label%3D%22delta_NH3%20(ug%2Fm3)%22)%0A%20%20%20%20for%20_ax%20in%20_axes%3A%0A%20%20%20%20%20%20%20%20_ax.set_xlim(*_xlim)%0A%20%20%20%20%20%20%20%20_ax.set_ylim(*_ylim)%0A%20%20%20%20style_map_axes(*_axes)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%3E%20%23%23%23%23%20%F0%9F%94%A7%20Your%20turn%20(Step%204)%0A%20%20%20%20%3E%0A%20%20%20%20%3E%20**Make%20all%20the%20models%20comparable%20by%20running%20them%20on%0A%20%20%20%20%3E%0A%20%20%20%20%60_fa%20%3D%20_fa.loc%5B(_fa%5B%22grid_row%22%5D%20%3E%20200)%20%26%20(_fa%5B%22grid_col%22%5D%20%3E%20150)%5D.copy().reset_index(drop%3DTrue)%60%0A%20%20%20%20%3E%0A%20%20%20%20%3E%201.%20How%20do%20the%20models%20compare%3F%0A%20%20%20%20%3E%0A%20%20%20%20%3E%20The%20*Aside*%20below%20does%20exactly%20this%20across%20five%20W.%20Does%20the%20choice%20of%20W%20impact%20the%20results%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Wrap-up%3A%20assumptions%2C%20how%20they%20fail%20here%2C%20and%20how%20you'd%20fix%20them%0A%0A%20%20%20%20%7C%20Assumption%20%7C%20How%20it%20fails%20in%20*this*%20analysis%20%7C%20Hypothetical%20fix%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20**Exchangeability**%20%7C%20sector%20cells%20differ%20in%20soil%2C%20climate%2C%20land%20use%20and%20history%2C%20much%20unobserved%20and%20spatially%20clustered%20%7C%20a%20timing%20shock%20(Step%203%20DiD)%2C%20a%20regulation%20boundary%2C%20or%20a%20wind%20instrument%20%7C%0A%20%20%20%20%7C%20**Consistency**%20%7C%20%22**firm%20count**%22%20is%20a%20crude%20treatment%3A%20it%20ignores%20herd%20size%2C%20emission%20technology%2C%20manure%20handling%20and%20exact%20timing%20%7C%20measure%20animal%20numbers%20%2F%20N-emissions%20%2F%20permits%2C%20not%20firm%20counts%20%7C%0A%20%20%20%20%7C%20**Positivity**%20%7C%20treated%20(gain)%20cells%20sit%20in%20a%20special%20part%20of%20the%20country%3B%20few%20truly%20comparable%20controls%20live%20there%20%7C%20restrict%20to%20a%20common-support%20region%3B%20match%20on%20covariates%20%7C%0A%20%20%20%20%7C%20**No%20interference**%20%7C%20NH%E2%82%83%20drifts%2C%20so%20%22control%22%20cells%20are%20partly%20treated%20%7C%20model%20it%20(Step%204%20SAR)%20or%20use%20**far-ring**%20controls%20%7C%0A%0A%20%20%20%20%3E%20**A%20map%20shows%20you%20*where*.%20A%20regression%20shows%20you%20*what%20correlates*.%20A%20design%20tells%20you%20*what%20would%20need%20to%20be%20true*%20for%20a%20causal%20claim.**%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Aside%20(a%20demonstration%2C%20not%20an%20exercise)%3A%20the%20SAR%20%22total%22%20rides%20on%20%CF%81%0A%0A%20%20%20%20We%20take%20the%20same%20first-difference%20spatial-lag%20(SAR)%20model%20from%20Step%204%20and%20swap%20**only%20the%0A%20%20%20%20weights%20matrix%20W**%20across%20five%20choices%3A%20Queen%20contiguity%20(1st%20and%202nd%20order)%20and%20the%203%20km%0A%20%20%20%20distance%20band%20at%20three%20decays%20(%241%2Fd%24%2C%20%241%2Fd%5E%7B2%7D%24%2C%20%241%2Fd%5E%7B3%7D%24).%20For%20speed%20we%20fit%20it%20on%20the%0A%20%20%20%20**smaller%20Step-2%20window**%20(the%20dense%20eastern%20livestock%20belt)%2C%20not%20the%20full%20Step-4%20subset.%0A%0A%20%20%20%20%3E%20%E2%9A%A0%EF%B8%8F%20The%20livestock%E2%86%92NH%E2%82%83%20effect%20is%20**spatially%20concentrated**%2C%20so%20on%20this%20small%20eastern%0A%20%20%20%20%3E%20window%20the%20local%20effect%20is%20**much%20larger**%20than%20the%20Step-4%20headline%20(**direct%20%E2%89%881.3%20here%0A%20%20%20%20%3E%20vs%20%E2%89%880.3**%20on%20the%20broader%20Step-4%20subset).%20That%20is%20a%20*different%20estimand%20on%20a%20different%0A%20%20%20%20%3E%20region*%2C%20**not**%20an%20error%3A%20watch%20the%20**pattern**%2C%20not%20the%20magnitude.%0A%0A%20%20%20%20-%20**direct**%20%3D%20a%20cell's%20**own**%20treatment%20on%20its%20**own**%20%24%5CDelta%5Ctext%7BNH%7D_3%24%3B%0A%20%20%20%20-%20**total**%20%3D%20the%20change%20if%20treatment%20rose%20**everywhere%20at%20once**%3B%0A%20%20%20%20-%20**indirect**%20%3D%20total%20%E2%88%92%20direct%20%3D%20the%20spillover%2C%20carried%20by%20the%20outcome%20feedback%20%CF%81.%0A%0A%20%20%20%20The%20**direct**%20effect%20barely%20budges%20across%20all%20five%20W.%20The%20**total**%20**rides%20on%20%CF%81**%20(fixed%0A%20%20%20%20by%20your%20choice%20of%20W)%3A%20same%20data%2C%20wildly%20different%20totals.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(DistanceBand%2C%20ML_Lag%2C%20did_cells%2C%20fit_spatial%2C%20pd%2C%20smf%2C%20spatial_impacts)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Aside%3A%20re-fit%20the%20Step-4%20SAR%20under%20five%20different%20W%20(only%20W%20changes)%20%E2%94%80%E2%94%80%0A%20%20%20%20%23%20Run%20it%20on%20the%20SMALLER%20Step-2%20window%20(row%3E200%20%26%20col%3E150)%20so%20the%20five%20fits%20are%0A%20%20%20%20%23%20fast%20(~seconds%20vs%20minutes).%20NOTE%3A%20this%20dense%20eastern%20window%20has%20a%20much%20LARGER%0A%20%20%20%20%23%20local%20effect%20than%20the%20Step-4%20headline%20(the%20livestock-%3ENH3%20effect%20is%20spatially%0A%20%20%20%20%23%20concentrated)%3B%20the%20takeaway%20here%20is%20the%20PATTERN%2C%20not%20the%20magnitude.%0A%20%20%20%20from%20libpysal.weights%20import%20W%20%20%20%23%20generic%20weights%20from%20a%20%7Bcell%3A%20%5Bneighbours%5D%7D%20dict%0A%20%20%20%20_fa%20%3D%20did_cells.copy()%0A%20%20%20%20_fa%5B%22delta_nh3%22%5D%20%3D%20_fa%5B%22nh3_2024%22%5D%20-%20_fa%5B%22nh3_2020%22%5D%0A%20%20%20%20_fa%20%3D%20_fa.dropna(subset%3D%5B%22delta_nh3%22%5D).copy()%0A%20%20%20%20_fa%20%3D%20_fa.loc%5B_fa%5B%22n_sector_2024%22%5D%20%3E%200%5D%0A%0A%20%20%20%20_fa%20%3D%20_fa.loc%5B(_fa%5B%22grid_row%22%5D%20%3E%20200)%20%26%20(_fa%5B%22grid_col%22%5D%20%3E%20150)%5D.copy().reset_index(drop%3DTrue)%0A%20%20%20%20_X%20%3D%20_fa%5B%5B%22treated%22%5D%5D.astype(float)%0A%20%20%20%20_y%20%3D%20_fa%5B%5B%22delta_nh3%22%5D%5D.values%0A%20%20%20%20print(f%22Aside%20subset%20(Step-2%20window)%3A%20%7Blen(_fa)%3A%2C%7D%20cells%2C%20%7Bint(_fa%5B'treated'%5D.sum())%7D%20treated%22)%0A%0A%20%20%20%20def%20_queen_from_grid(df%2C%20order%3D1)%3A%0A%20%20%20%20%20%20%20%20_pos%20%3D%20%7B(r%2C%20c)%3A%20i%20for%20i%2C%20(r%2C%20c)%20in%20enumerate(zip(df.grid_row%2C%20df.grid_col))%7D%0A%20%20%20%20%20%20%20%20_offs%20%3D%20%5B(dr%2C%20dc)%20for%20dr%20in%20range(-order%2C%20order%20%2B%201)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20dc%20in%20range(-order%2C%20order%20%2B%201)%20if%20(dr%2C%20dc)%20!%3D%20(0%2C%200)%5D%0A%20%20%20%20%20%20%20%20_neigh%20%3D%20%7Bi%3A%20%5B_pos%5B(r%20%2B%20dr%2C%20c%20%2B%20dc)%5D%20for%20dr%2C%20dc%20in%20_offs%20if%20(r%20%2B%20dr%2C%20c%20%2B%20dc)%20in%20_pos%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20(r%2C%20c)%2C%20i%20in%20_pos.items()%7D%0A%20%20%20%20%20%20%20%20return%20W(_neigh%2C%20id_order%3Dlist(range(len(df)))%2C%20silence_warnings%3DTrue)%0A%0A%20%20%20%20def%20_dist_band(alpha)%3A%0A%20%20%20%20%20%20%20%20return%20DistanceBand(_fa%5B%5B%22x_rd%22%2C%20%22y_rd%22%5D%5D.values%2C%20threshold%3D3000%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20binary%3DFalse%2C%20alpha%3Dalpha%2C%20silence_warnings%3DTrue)%0A%0A%20%20%20%20_W_family%20%3D%20%7B%22Queen-1%22%3A%20_queen_from_grid(_fa%2C%201)%2C%20%22Queen-2%22%3A%20_queen_from_grid(_fa%2C%202)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Dist%20a%3D-1%22%3A%20_dist_band(-1)%2C%20%22Dist%20a%3D-2%22%3A%20_dist_band(-2)%2C%20%22Dist%20a%3D-3%22%3A%20_dist_band(-3)%7D%0A%20%20%20%20for%20_w%20in%20_W_family.values()%3A%0A%20%20%20%20%20%20%20%20_w.transform%20%3D%20%22R%22%0A%0A%20%20%20%20_rows%20%3D%20%5B%5D%0A%20%20%20%20for%20_name%2C%20_w%20in%20_W_family.items()%3A%0A%20%20%20%20%20%20%20%20_m%2C%20_%2C%20_%20%3D%20fit_spatial(ML_Lag%2C%20_y%2C%20_X.values%2C%20_w%2C%20%22delta_nh3%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20list(_X.columns)%2C%20spat_impacts%3D%22full%22)%0A%20%20%20%20%20%20%20%20_eff%2C%20_rho%20%3D%20spatial_impacts(_m).loc%5B%22treated%22%5D%2C%20float(_m.rho)%0A%20%20%20%20%20%20%20%20_rows.append(%7B%22W%22%3A%20_name%2C%20%22avg%20neigh%22%3A%20round(_w.mean_neighbors%2C%201)%2C%20%22rho%22%3A%20_rho%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%221%2F(1-rho)%22%3A%201%20%2F%20(1%20-%20_rho)%2C%20%22direct%22%3A%20_eff%5B%22direct%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22indirect%22%3A%20_eff%5B%22indirect%22%5D%2C%20%22total%22%3A%20_eff%5B%22total%22%5D%7D)%0A%0A%20%20%20%20wsens%20%3D%20pd.DataFrame(_rows).set_index(%22W%22)%0A%20%20%20%20own%20%3D%20smf.ols(%22delta_nh3%20~%20treated%22%2C%20data%3D_fa).fit().params%5B%22treated%22%5D%0A%20%20%20%20print(f%22Own-cell%20DiD%20on%20this%20window%20(the%20anchor)%3A%20%7Bown%3A%2B.2f%7D%20ug%2Fm3%22)%0A%20%20%20%20print(f%22SAR%20DIRECT%20stays%20~%7Bwsens%5B'direct'%5D.min()%3A.2f%7D-%7Bwsens%5B'direct'%5D.max()%3A.2f%7D%20%20%22%0A%20%20%20%20%20%20%20%20%20%20f%22(robust%20across%20W%2C%20sits%20next%20to%20the%20own-cell%20DiD)%22)%0A%20%20%20%20print(f%22SAR%20TOTAL%20swings%20%7Bwsens%5B'total'%5D.min()%3A.2f%7D-%7Bwsens%5B'total'%5D.max()%3A.2f%7D%20%20%22%0A%20%20%20%20%20%20%20%20%20%20f%22(rides%201%2F(1-rho)%20-%3E%20a%20property%20of%20W%2C%20not%20of%20the%20data)%22)%0A%20%20%20%20wsens.round(2)%0A%20%20%20%20return%20own%2C%20wsens%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20own%2C%20plt%2C%20wsens)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Aside%3A%20see%20it%3A%20direct%20is%20flat%2C%20total%20tracks%20the%20spatial%20multiplier%20%E2%94%80%E2%94%80%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(9%2C%205))%0A%20%20%20%20_xpos%20%3D%20np.arange(len(wsens))%0A%20%20%20%20_ax.bar(_xpos%2C%20wsens%5B%22total%22%5D%2C%20width%3D0.55%2C%20color%3D%22%23c6dbef%22%2C%20edgecolor%3D%22%233182bd%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22SAR%20total%20(own%20%2B%20spillover%20via%20rho)%22)%0A%20%20%20%20_ax.bar(_xpos%2C%20wsens%5B%22direct%22%5D%2C%20width%3D0.55%2C%20color%3D%22%233182bd%22%2C%20label%3D%22SAR%20direct%20(own%20cell)%22)%0A%20%20%20%20_ax.axhline(own%2C%20color%3D%22%23d95f0e%22%2C%20lw%3D2%2C%20ls%3D%22--%22%2C%20label%3Df%22plain%20own-cell%20DiD%20%3D%20%7Bown%3A.2f%7D%22)%0A%20%20%20%20for%20_i%2C%20(_%2C%20_r)%20in%20enumerate(wsens.iterrows())%3A%0A%20%20%20%20%20%20%20%20_ax.text(_i%2C%20_r%5B%22total%22%5D%20%2B%200.08%2C%20f%22rho%3D%7B_r%5B'rho'%5D%3A.2f%7D%22%2C%20ha%3D%22center%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20va%3D%22bottom%22%2C%20fontsize%3D9%2C%20color%3D%22%233182bd%22)%0A%20%20%20%20_ax.set_xticks(_xpos)%0A%20%20%20%20_ax.set_xticklabels(wsens.index)%0A%20%20%20%20_ax.set_ylabel(%22effect%20on%20delta_NH3%20per%20treated%20unit%20(ug%2Fm3)%22)%0A%20%20%20%20_ax.set_title(%22Same%20SAR%2C%20five%20W%3A%20the%20own-cell%20(direct)%20effect%20is%20stable%3B%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22the%20'total'%20rides%20the%201%2F(1-rho)%20multiplier%22)%0A%20%20%20%20_ax.legend(loc%3D%22upper%20left%22%2C%20fontsize%3D8.5)%0A%20%20%20%20_ax.margins(y%3D0.15)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
a81a2a599859a72ad801613f69ed3a10