In [1]:
import pandas as pd
from sqlalchemy import create_engine
import geopandas as gpd
from dotenv import load_dotenv
load_dotenv()
import os
# Database connection details from zshrc environment variables
# use .env file for local development
db_name = 'colorado_spills'
user = os.getenv('DB_USER')
password = os.getenv('DB_PASSWORD')
host = os.getenv('DB_HOST')
port = os.getenv('DB_PORT')
# Create an engine to connect to the PostgreSQL database
engine = create_engine(f'postgresql+psycopg2://{user}:{password}@{host}:{port}/{db_name}')
# Read in the spills_with_demographics_geog as spills
spills = pd.read_sql_table('spills_with_demographics_geog', engine)
/home/dadams/miniconda3/envs/spatial_env2/lib/python3.10/site-packages/pandas/io/sql.py:1725: SAWarning: Did not recognize type 'geometry' of column 'geometry' self.meta.reflect(bind=self.con, only=[table_name], views=True)
In [2]:
# use longitude and latitude to create a GeoDataFrame from spills data
spills['geometry'] = gpd.points_from_xy(spills['Longitude'], spills['Latitude'])
In [3]:
spills_gdf = gpd.GeoDataFrame(spills, crs='EPSG:4326')
In [4]:
# Ensure date columns are in datetime format
spills_gdf['Date of Discovery'] = pd.to_datetime(spills_gdf['Date of Discovery'])
spills_gdf['Initial Report Date'] = pd.to_datetime(spills_gdf['Initial Report Date'])
In [5]:
# drop 2024 data for date of discovery and initial report date
spills_gdf = spills_gdf[spills_gdf['Date of Discovery'].dt.year < 2024]
spills_gdf = spills_gdf[spills_gdf['Initial Report Date'].dt.year < 2024]
In [6]:
# Separate Historical vs. Recent Spills
historical_spills = spills_gdf[spills_gdf['Spill Type'] == 'Historical']
recent_spills = spills_gdf[spills_gdf['Spill Type'] == 'Recent']
In [7]:
# Compare Reporting Delay
report_delay_summary = spills_gdf.groupby('Spill Type')['Report Delay (Days)'].describe()
print(report_delay_summary)
# Ensure the output directory exists
output_dir = 'analysis/output'
os.makedirs(output_dir, exist_ok=True)
# Save the report delay summary to a markdown file
with open(os.path.join(output_dir, 'report_delay_summary.md'), 'w') as f:
f.write("# Report Delay Summary\n\n")
f.write(report_delay_summary.to_markdown())
count mean std min 25% 50% 75% max Spill Type Historical 5312.0 21.490023 205.615913 0.0 0.0 1.0 2.0 9261.0 Recent 10678.0 3.964506 46.997749 0.0 0.0 1.0 2.0 2232.0
In [8]:
import matplotlib.pyplot as plt
# Compare Spill Trends Over Time
spills_by_year = spills_gdf.groupby(['Report Year', 'Spill Type']).size().unstack()
ax = spills_by_year.plot(kind='bar', figsize=(12,6), title='Spills by Year and Type')
ax.set_xlabel('Year')
ax.set_ylabel('Number of Spills')
ax.legend(title='Spill Type')
# Ensure the output directory exists
output_dir = 'analysis/output'
os.makedirs(output_dir, exist_ok=True)
# Save the plot
plt.savefig(f'{output_dir}/spills_by_year_and_type.png')
# Display the plot
plt.show()
In [9]:
import matplotlib.pyplot as plt
import geopandas as gpd
# Ensure spills_gdf is already a GeoDataFrame with the necessary data
# spills_gdf = gpd.GeoDataFrame(...)
# Create a plot for the spatial distribution of spills
fig, ax = plt.subplots(figsize=(10, 6))
# Plot the GeoDataFrame with optimizations
spills_gdf.plot(
ax=ax,
marker='o',
alpha=0.6, # Adjust transparency for better visibility
column='Spill Type',
legend=True,
legend_kwds={'title': "Spill Type"}, # Correctly specify the legend title
cmap='viridis' # Use a colormap for better visual distinction
)
# Set plot title and labels
ax.set_title('Spatial Distribution of Historical vs. Recent Spills', fontsize=14)
ax.set_xlabel('Longitude', fontsize=12)
ax.set_ylabel('Latitude', fontsize=12)
# Improve layout
plt.tight_layout()
# Display the plot
plt.show()
In [10]:
import matplotlib.pyplot as plt
import contextily as ctx
import numpy as np
# Use your ORIGINAL plot code exactly as it is, just adding the basemap at the end
# This ensures all your data points, colors and categories remain intact
# Starting with your original code:
fig, ax = plt.subplots(figsize=(10, 6))
# Plot the GeoDataFrame with optimizations - YOUR ORIGINAL PLOTTING CODE
spills_gdf.plot(
ax=ax,
marker='o',
alpha=0.6, # Adjust transparency for better visibility
column='Spill Type',
legend=True,
legend_kwds={'title': "Spill Type"}, # Correctly specify the legend title
cmap='viridis' # Use a colormap for better visual distinction
)
# Set plot title and labels
ax.set_title('Spatial Distribution of Historical vs. Recent Spills in Colorado', fontsize=14)
ax.set_xlabel('Longitude', fontsize=12)
ax.set_ylabel('Latitude', fontsize=12)
# Store the current axis limits to restore them after adding basemap
xlim = ax.get_xlim()
ylim = ax.get_ylim()
# Before adding the basemap, convert the axes coordinates to Web Mercator
# This is a crucial step that lets us keep your original points exactly as they are
from pyproj import Transformer
# Create transformer from lat/lon to Web Mercator
transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
# Transform the limits
west, south = transformer.transform(xlim[0], ylim[0])
east, north = transformer.transform(xlim[1], ylim[1])
# Add the basemap - this will be UNDER your existing points
ctx.add_basemap(
ax,
crs=spills_gdf.crs,
source=ctx.providers.OpenStreetMap.Mapnik,
zoom="auto"
)
# Restore original axis limits
ax.set_xlim(xlim)
ax.set_ylim(ylim)
# Improve layout
plt.tight_layout()
# Save the plot to the analysis/output directory
output_dir = 'analysis/output'
os.makedirs(output_dir, exist_ok=True)
plt.savefig(f'{output_dir}/spatial_distribution_spills.png')
# Show the plot
plt.show()
In [11]:
# Add a column to indicate whether the spill occurred before or after 2020
spills_gdf['Period'] = spills_gdf['Report Year'].apply(lambda x: 'Before 2020' if x < 2020 else '2020 and After')
# Group by the new period column and calculate the mean response time
response_time_summary = spills_gdf.groupby('Period')['Report Delay (Days)'].describe()
print(response_time_summary)
# Ensure the output directory exists
output_dir = 'analysis/output'
os.makedirs(output_dir, exist_ok=True)
# Save the response time summary to a markdown file
with open(os.path.join(output_dir, 'response_time_summary.md'), 'w') as f:
f.write("# Response Time Summary\n\n")
f.write(response_time_summary.to_markdown())
count mean std min 25% 50% 75% max Period 2020 and After 7177.0 6.924342 124.965839 0.0 0.0 1.0 1.0 9261.0 Before 2020 8813.0 12.117554 124.705876 0.0 0.0 1.0 2.0 5681.0
In [12]:
# t-test for the difference in response time before and after 2020
from scipy.stats import ttest_ind
before_2020 = spills_gdf[spills_gdf['Period'] == 'Before 2020']['Report Delay (Days)']
after_2020 = spills_gdf[spills_gdf['Period'] == '2020 and After']['Report Delay (Days)']
t_stat, p_value = ttest_ind(before_2020, after_2020, equal_var=False)
print(f'T-Statistic: {t_stat:.4f}')
print(f'P-Value: {p_value:.4f}')
# Check if the difference is statistically significant
if p_value < 0.05:
print('The difference in response time is statistically significant.')
else:
print('The difference in response time is not statistically significant.')
# Create a plot for the response time distribution before and after 2020
fig, ax = plt.subplots(figsize=(10, 6))
# Plot the response time distribution for each period
spills_gdf.boxplot(column='Report Delay (Days)', by='Period', ax=ax)
T-Statistic: 2.6161 P-Value: 0.0089 The difference in response time is statistically significant.
Out[12]:
<Axes: title={'center': 'Report Delay (Days)'}, xlabel='Period'>
In [13]:
%pip install tabulate
Requirement already satisfied: tabulate in /home/dadams/miniconda3/envs/spatial_env2/lib/python3.10/site-packages (0.9.0)
Note: you may need to restart the kernel to use updated packages.
In [14]:
import tabulate
from scipy.stats import ttest_ind
# Ensure the output directory exists
output_dir = 'analysis/output'
os.makedirs(output_dir, exist_ok=True)
# Add a column to indicate whether the spill occurred before or after 2020 for historical and recent spills
historical_spills['Period'] = historical_spills['Report Year'].apply(lambda x: 'Before 2020' if x < 2020 else '2020 and After')
recent_spills['Period'] = recent_spills['Report Year'].apply(lambda x: 'Before 2020' if x < 2020 else '2020 and After')
# Historical Spills Analysis
historical_report_delay_summary = historical_spills['Report Delay (Days)'].describe()
print("Historical Spills Report Delay Summary:")
print(historical_report_delay_summary)
# Recent Spills Analysis
recent_report_delay_summary = recent_spills['Report Delay (Days)'].describe()
print("\nRecent Spills Report Delay Summary:")
print(recent_report_delay_summary)
# Historical Spills Response Time Summary
historical_response_time_summary = historical_spills.groupby('Period')['Report Delay (Days)'].describe()
print("\nHistorical Spills Response Time Summary:")
print(historical_response_time_summary)
# Recent Spills Response Time Summary
recent_response_time_summary = recent_spills.groupby('Period')['Report Delay (Days)'].describe()
print("\nRecent Spills Response Time Summary:")
print(recent_response_time_summary)
# t-test for the difference in response time before and after 2020 for historical spills
historical_before_2020 = historical_spills[historical_spills['Period'] == 'Before 2020']['Report Delay (Days)']
historical_after_2020 = historical_spills[historical_spills['Period'] == '2020 and After']['Report Delay (Days)']
t_stat_historical, p_value_historical = ttest_ind(historical_before_2020, historical_after_2020, equal_var=False)
print(f'\nHistorical Spills T-Statistic: {t_stat_historical:.4f}')
print(f'Historical Spills P-Value: {p_value_historical:.4f}')
# Check if the difference is statistically significant for historical spills
if p_value_historical < 0.05:
print('The difference in response time for historical spills is statistically significant.')
else:
print('The difference in response time for historical spills is not statistically significant.')
# t-test for the difference in response time before and after 2020 for recent spills
recent_before_2020 = recent_spills[recent_spills['Period'] == 'Before 2020']['Report Delay (Days)']
recent_after_2020 = recent_spills[recent_spills['Period'] == '2020 and After']['Report Delay (Days)']
t_stat_recent, p_value_recent = ttest_ind(recent_before_2020, recent_after_2020, equal_var=False)
print(f'\nRecent Spills T-Statistic: {t_stat_recent:.4f}')
print(f'Recent Spills P-Value: {p_value_recent:.4f}')
# Check if the difference is statistically significant for recent spills
if p_value_recent < 0.05:
print('The difference in response time for recent spills is statistically significant.')
else:
print('The difference in response time for recent spills is not statistically significant.')
# Save the summaries and t-test results to a markdown file
with open(os.path.join(output_dir, 'analysis_results.md'), 'w') as f:
f.write("# Analysis Results\n\n")
f.write("## Historical Spills Report Delay Summary:\n")
f.write(historical_report_delay_summary.to_markdown() + "\n\n")
f.write("## Recent Spills Report Delay Summary:\n")
f.write(recent_report_delay_summary.to_markdown() + "\n\n")
f.write("## Historical Spills Response Time Summary:\n")
f.write(historical_response_time_summary.to_markdown() + "\n\n")
f.write("## Recent Spills Response Time Summary:\n")
f.write(recent_response_time_summary.to_markdown() + "\n\n")
f.write("## T-Test Results for Historical Spills\n")
f.write(f"T-Statistic: {t_stat_historical:.4f}\n")
f.write(f"P-Value: {p_value_historical:.4f}\n")
if p_value_historical < 0.05:
f.write('The difference in response time for historical spills is statistically significant.\n')
else:
f.write('The difference in response time for historical spills is not statistically significant.\n')
f.write("\n## T-Test Results for Recent Spills\n")
f.write(f"T-Statistic: {t_stat_recent:.4f}\n")
f.write(f"P-Value: {p_value_recent:.4f}\n")
if p_value_recent < 0.05:
f.write('The difference in response time for recent spills is statistically significant.\n')
else:
f.write('The difference in response time for recent spills is not statistically significant.\n')
Historical Spills Report Delay Summary:
count 5312.000000
mean 21.490023
std 205.615913
min 0.000000
25% 0.000000
50% 1.000000
75% 2.000000
max 9261.000000
Name: Report Delay (Days), dtype: float64
Recent Spills Report Delay Summary:
count 10678.000000
mean 3.964506
std 46.997749
min 0.000000
25% 0.000000
50% 1.000000
75% 2.000000
max 2232.000000
Name: Report Delay (Days), dtype: float64
Historical Spills Response Time Summary:
count mean std min 25% 50% 75% max
Period
2020 and After 2705.0 13.533087 200.352339 0.0 0.0 0.0 1.0 9261.0
Before 2020 2607.0 29.746068 210.659491 0.0 0.0 1.0 2.0 5681.0
Recent Spills Response Time Summary:
count mean std min 25% 50% 75% max
Period
2020 and After 4472.0 2.926878 27.302062 0.0 0.0 1.0 1.0 1329.0
Before 2020 6206.0 4.712214 57.116100 0.0 1.0 1.0 2.0 2232.0
Historical Spills T-Statistic: 2.8723
Historical Spills P-Value: 0.0041
The difference in response time for historical spills is statistically significant.
Recent Spills T-Statistic: 2.1457
Recent Spills P-Value: 0.0319
The difference in response time for recent spills is statistically significant.
/home/dadams/miniconda3/envs/spatial_env2/lib/python3.10/site-packages/geopandas/geodataframe.py:1819: 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: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy super().__setitem__(key, value) /home/dadams/miniconda3/envs/spatial_env2/lib/python3.10/site-packages/geopandas/geodataframe.py:1819: 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: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy super().__setitem__(key, value)
In [15]:
import os
# Ensure the output directory exists
output_dir = 'analysis/output'
os.makedirs(output_dir, exist_ok=True)
# Save the t-test results to a markdown file
with open(os.path.join(output_dir, 't_test_results.md'), 'w') as f:
f.write(f"## T-Test Results for Historical Spills\n")
f.write(f"T-Statistic: {t_stat_historical:.4f}\n")
f.write(f"P-Value: {p_value_historical:.4f}\n")
if p_value_historical < 0.05:
f.write('The difference in response time for historical spills is statistically significant.\n')
else:
f.write('The difference in response time for historical spills is not statistically significant.\n')
f.write(f"\n## T-Test Results for Recent Spills\n")
f.write(f"T-Statistic: {t_stat_recent:.4f}\n")
f.write(f"P-Value: {p_value_recent:.4f}\n")
if p_value_recent < 0.05:
f.write('The difference in response time for recent spills is statistically significant.\n')
else:
f.write('The difference in response time for recent spills is not statistically significant.\n')
# Create and save the plots
fig, ax = plt.subplots(figsize=(10, 6))
historical_spills.boxplot(column='Report Delay (Days)', by='Period', ax=ax)
ax.set_title('Historical Spills Response Time Distribution Before and After 2020', fontsize=14)
ax.set_xlabel('Period', fontsize=12)
ax.set_ylabel('Report Delay (Days)', fontsize=12)
plt.suptitle('')
plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'historical_spills_response_time_distribution.png'))
fig, ax = plt.subplots(figsize=(10, 6))
recent_spills.boxplot(column='Report Delay (Days)', by='Period', ax=ax)
ax.set_title('Recent Spills Response Time Distribution Before and After 2020', fontsize=14)
ax.set_xlabel('Period', fontsize=12)
ax.set_ylabel('Report Delay (Days)', fontsize=12)
plt.suptitle('')
plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'recent_spills_response_time_distribution.png'))
In [16]:
# t-test for the difference in response time before and after 2020 for historical spills
historical_before_2020 = historical_spills[historical_spills['Period'] == 'Before 2020']['Report Delay (Days)']
historical_after_2020 = historical_spills[historical_spills['Period'] == '2020 and After']['Report Delay (Days)']
t_stat, p_value = ttest_ind(historical_before_2020, historical_after_2020, equal_var=False)
print(f'\nHistorical Spills T-Statistic: {t_stat:.4f}')
print(f'Historical Spills P-Value: {p_value:.4f}')
# Check if the difference is statistically significant for historical spills
if p_value < 0.05:
print('The difference in response time for historical spills is statistically significant.')
else:
print('The difference in response time for historical spills is not statistically significant.')
# t-test for the difference in response time before and after 2020 for recent spills
recent_before_2020 = recent_spills[recent_spills['Period'] == 'Before 2020']['Report Delay (Days)']
recent_after_2020 = recent_spills[recent_spills['Period'] == '2020 and After']['Report Delay (Days)']
t_stat, p_value = ttest_ind(recent_before_2020, recent_after_2020, equal_var=False)
print(f'\nRecent Spills T-Statistic: {t_stat:.4f}')
print(f'Recent Spills P-Value: {p_value:.4f}')
# Check if the difference is statistically significant for recent spills
if p_value < 0.05:
print('The difference in response time for recent spills is statistically significant.')
else:
print('The difference in response time for recent spills is not statistically significant.')
# Create a plot for the response time distribution before and after 2020 for historical and recent spills
fig, ax = plt.subplots(figsize=(10, 6))
# Plot the response time distribution for historical spills
historical_spills.boxplot(column='Report Delay (Days)', by='Period', ax=ax, positions=[1, 2])
# Plot the response time distribution for recent spills
recent_spills.boxplot(column='Report Delay (Days)', by='Period', ax=ax, positions=[4, 5])
# Save the plot to the /output directory
plt.savefig('analysis/output/response_time_distribution.png')
Historical Spills T-Statistic: 2.8723 Historical Spills P-Value: 0.0041 The difference in response time for historical spills is statistically significant. Recent Spills T-Statistic: 2.1457 Recent Spills P-Value: 0.0319 The difference in response time for recent spills is statistically significant.
In [ ]:
# Historical only
historical_spills = historical_spills.copy()
historical_spills['Post2020'] = historical_spills['Period'].apply(lambda x: 1 if x == '2020 and After' else 0)
model_hist = smf.ols('Q("Report Delay (Days)") ~ Post2020', data=historical_spills).fit()
print("Historical Spills:\n", model_hist.summary())
# Recent only
recent_spills = recent_spills.copy()
recent_spills['Post2020'] = recent_spills['Period'].apply(lambda x: 1 if x == '2020 and After' else 0)
model_recent = smf.ols('Q("Report Delay (Days)") ~ Post2020', data=recent_spills).fit()
print("\nRecent Spills:\n", model_recent.summary())
Historical Spills:
OLS Regression Results
====================================================================================
Dep. Variable: Q("Report Delay (Days)") R-squared: 0.002
Model: OLS Adj. R-squared: 0.001
Method: Least Squares F-statistic: 8.265
Date: Mon, 31 Mar 2025 Prob (F-statistic): 0.00406
Time: 13:58:43 Log-Likelihood: -35825.
No. Observations: 5312 AIC: 7.165e+04
Df Residuals: 5310 BIC: 7.167e+04
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 29.7461 4.024 7.392 0.000 21.857 37.635
Post2020 -16.2130 5.639 -2.875 0.004 -27.269 -5.157
==============================================================================
Omnibus: 12943.048 Durbin-Watson: 1.944
Prob(Omnibus): 0.000 Jarque-Bera (JB): 189644913.057
Skew: 25.310 Prob(JB): 0.00
Kurtosis: 927.266 Cond. No. 2.64
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Recent Spills:
OLS Regression Results
====================================================================================
Dep. Variable: Q("Report Delay (Days)") R-squared: 0.000
Model: OLS Adj. R-squared: 0.000
Method: Least Squares F-statistic: 3.752
Date: Mon, 31 Mar 2025 Prob (F-statistic): 0.0528
Time: 13:58:43 Log-Likelihood: -56260.
No. Observations: 10678 AIC: 1.125e+05
Df Residuals: 10676 BIC: 1.125e+05
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 4.7122 0.597 7.900 0.000 3.543 5.881
Post2020 -1.7853 0.922 -1.937 0.053 -3.592 0.021
==============================================================================
Omnibus: 28541.449 Durbin-Watson: 1.643
Prob(Omnibus): 0.000 Jarque-Bera (JB): 676851501.928
Skew: 32.390 Prob(JB): 0.00
Kurtosis: 1234.707 Cond. No. 2.47
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [ ]:
# Label both groups, combine
historical_spills['Group'] = 'Historical'
recent_spills['Group'] = 'Recent'
combined = pd.concat([historical_spills, recent_spills], ignore_index=True)
combined['Post2020'] = combined['Period'].apply(lambda x: 1 if x == '2020 and After' else 0)
# Run regression with interactions by group
model = smf.ols('Q("Report Delay (Days)") ~ Post2020 * Group', data=combined).fit()
print(model.summary())
OLS Regression Results
====================================================================================
Dep. Variable: Q("Report Delay (Days)") R-squared: 0.006
Model: OLS Adj. R-squared: 0.006
Method: Least Squares F-statistic: 31.12
Date: Mon, 31 Mar 2025 Prob (F-statistic): 4.78e-20
Time: 13:53:41 Log-Likelihood: -99827.
No. Observations: 15990 AIC: 1.997e+05
Df Residuals: 15986 BIC: 1.997e+05
Df Model: 3
Covariance Type: nonrobust
============================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------------------
Intercept 29.7461 2.438 12.200 0.000 24.967 34.525
Group[T.Recent] -25.0339 2.906 -8.616 0.000 -30.729 -19.339
Post2020 -16.2130 3.417 -4.745 0.000 -22.910 -9.516
Post2020:Group[T.Recent] 14.4276 4.200 3.435 0.001 6.196 22.660
==============================================================================
Omnibus: 45995.174 Durbin-Watson: 1.915
Prob(Omnibus): 0.000 Jarque-Bera (JB): 3525625525.882
Skew: 38.959 Prob(JB): 0.00
Kurtosis: 2302.059 Cond. No. 8.21
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [20]:
import matplotlib.pyplot as plt
# Data
pre_means = [29.75, 4.71]
post_means = [29.75 - 16.21, 4.71 - 1.79]
labels = ['Historical Spills', 'Recent Spills']
x = range(len(labels))
width = 0.35
# Plot
fig, ax = plt.subplots(figsize=(8, 6))
ax.bar([i - width / 2 for i in x], pre_means, width, label='Pre-2020', color='skyblue')
ax.bar([i + width / 2 for i in x], post_means, width, label='Post-2020', color='salmon')
# Labels and aesthetics
ax.set_ylabel('Mean Report Delay (Days)')
ax.set_title('Mean Report Delay Before and After 2020 by Spill Type')
ax.set_xticks(x)
ax.set_xticklabels(labels)
ax.legend()
# Save and show
plt.tight_layout()
plt.savefig('report_delay_barplot.png')
plt.show()
| Group | Pre-2020 Mean | Post-2020 Change (coef) | p-value | Significant (p < 0.05) |
|---|---|---|---|---|
| Historical Spills | 29.75 | -16.21 | 0.004 | Yes |
| Recent Spills | 4.71 | -1.79 | 0.053 | No |
Analysis of Reporting Delays for Historical and Recent Spills¶
We analyzed reporting delays for historical and recent spills before and after 2020 using simple OLS regressions. For historical spills, the average delay decreased significantly by approximately 16.2 days after 2020 (p = 0.004), dropping from a pre-2020 mean of 29.75 days. In contrast, the reduction for recent spills was smaller—about 1.8 days—and only marginally significant (p = 0.053), starting from a much lower pre-2020 average of 4.71 days. These results suggest that the most substantial improvements in timeliness occurred within the historical spill reporting process.
Differences in Reporting Delays by Spill Type¶
In [22]:
import pandas as pd
import statsmodels.formula.api as smf
# Combine the datasets and label the groups
historical_spills['SpillType'] = 0 # historical = 0
recent_spills['SpillType'] = 1 # recent = 1
combined = pd.concat([historical_spills, recent_spills], ignore_index=True)
# Convert 'Period' to binary: 0 = Before 2020, 1 = 2020 and After
combined['Post2020'] = combined['Period'].apply(lambda x: 1 if x == '2020 and After' else 0)
# Interaction term is created implicitly in the regression
model = smf.ols('Q("Report Delay (Days)") ~ Post2020 + SpillType + Post2020:SpillType', data=combined).fit()
print(model.summary())
OLS Regression Results
====================================================================================
Dep. Variable: Q("Report Delay (Days)") R-squared: 0.006
Model: OLS Adj. R-squared: 0.006
Method: Least Squares F-statistic: 31.12
Date: Mon, 31 Mar 2025 Prob (F-statistic): 4.78e-20
Time: 14:02:43 Log-Likelihood: -99827.
No. Observations: 15990 AIC: 1.997e+05
Df Residuals: 15986 BIC: 1.997e+05
Df Model: 3
Covariance Type: nonrobust
======================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------------
Intercept 29.7461 2.438 12.200 0.000 24.967 34.525
Post2020 -16.2130 3.417 -4.745 0.000 -22.910 -9.516
SpillType -25.0339 2.906 -8.616 0.000 -30.729 -19.339
Post2020:SpillType 14.4276 4.200 3.435 0.001 6.196 22.660
==============================================================================
Omnibus: 45995.174 Durbin-Watson: 1.915
Prob(Omnibus): 0.000 Jarque-Bera (JB): 3525625525.882
Skew: 38.959 Prob(JB): 0.00
Kurtosis: 2302.059 Cond. No. 8.21
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
| Term | Coefficient | Std. Error | p-value | Interpretation |
|---|---|---|---|---|
| Intercept | 29.75 | 2.44 | <0.001 | Historical spills, pre-2020 average report delay |
| Post2020 | -16.21 | 3.42 | <0.001 | Historical spills saw a 16.2-day reduction in delay after 2020 |
| Group: Recent Spills | -25.03 | 2.91 | <0.001 | Recent spills had a 25-day shorter delay than historical spills pre-2020 |
| Post2020 × Recent Spills | +14.43 | 4.20 | 0.001 | Recent spills improved 14.4 days less than historical spills post-2020 |
In [24]:
import matplotlib.pyplot as plt
# X-axis positions: 0 = Before 2020, 1 = 2020 and After
x = [0, 1]
# Y-values for each group
historical = [29.75, 29.75 - 16.21] # Before and after for historical
recent = [4.71, 4.71 - 1.79] # Before and after for recent
# Plot
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(x, historical, marker='o', label='Historical Spills', linewidth=2)
ax.plot(x, recent, marker='o', label='Recent Spills', linewidth=2)
# Aesthetics
ax.set_xticks(x)
ax.set_xticklabels(['Before 2020', '2020 and After'])
ax.set_ylabel('Mean Report Delay (Days)')
ax.set_title('Difference-in-Differences: Report Delay by Spill Type')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()
In [25]:
import matplotlib.pyplot as plt
improvements = [-16.21, -1.79]
labels = ['Historical Spills', 'Recent Spills']
fig, ax = plt.subplots(figsize=(7, 5))
ax.bar(labels, improvements, color=['steelblue', 'orange'])
ax.axhline(0, color='black', linewidth=0.8)
ax.set_ylabel('Change in Report Delay (Days)')
ax.set_title('Change in Report Delay After 2020')
plt.tight_layout()
plt.show()
Interpretation of Difference-in-Differences Model
This Difference-in-Differences model estimates how report delays changed after 2020 for two groups—historical and recent spills—and compares the magnitude of change between them. The pre-2020 average delay for historical spills was approximately 29.75 days. After 2020, historical spills saw a statistically significant reduction of about 16.2 days (p < 0.001).
Recent spills, by contrast, already had much shorter average delays (about 25 days less than historical spills before 2020). The DiD interaction term (+14.4 days, p = 0.001) shows that the improvement in timeliness for recent spills was significantly smaller than that for historical spills. In other words, while both groups saw reductions in report delays after 2020, the gains were substantially greater for historical spills.
This suggests that any policy or process changes implemented around 2020 had a disproportionately stronger effect on addressing delays in the historical spill category.