NYC Street-Tree Ecological Inequality: Final Analysis and Figure Code¶
0. Setup¶
In [ ]:
"""
Revised analysis workflow for the NYC street-tree report.
This version is research-question driven and adds:
- explicit analytical sample definitions
- borough-grouped validation for machine learning
- permutation importance and SHAP diagnostics
- co-burden sensitivity analysis and borough summaries
- revised tables for the final report
"""
from __future__ import annotations
import json
import math
import warnings
from pathlib import Path
import geopandas as gpd
import libpysal
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import shap
import statsmodels.api as sm
from catboost import CatBoostRegressor
from esda.moran import Moran, Moran_Local
from matplotlib.patches import Patch, Rectangle
from sklearn.ensemble import RandomForestRegressor
from sklearn.impute import SimpleImputer
from sklearn.inspection import permutation_importance
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.model_selection import GroupKFold, train_test_split
from sklearn.preprocessing import StandardScaler
from spreg import GM_Error, ML_Error, ML_Lag
from statsmodels.stats.outliers_influence import variance_inflation_factor
from xgboost import XGBRegressor
warnings.filterwarnings("ignore")
# Resolve the repository root whether the notebook is launched from the repo root or scripts/.
PROJECT = Path.cwd().resolve()
if not (PROJECT / "data").exists() and (PROJECT.parent / "data").exists():
PROJECT = PROJECT.parent
PROCESSED = PROJECT / "data" / "processed"
FINAL = PROJECT / "data" / "final"
REVISED = PROJECT / "data" / "revised"
VIS_FINAL = PROJECT / "visualizations" / "final"
VIS_BASE = PROJECT / "visualizations"
VIS_REVISED = PROJECT / "visualizations" / "revised"
DOCS_ASSETS = PROJECT / "src" / "atlas" / "assets"
for path in [REVISED, VIS_REVISED]:
path.mkdir(parents=True, exist_ok=True)
sns.set_theme(style="whitegrid", context="notebook")
plt.rcParams.update(
{
"figure.dpi": 160,
"savefig.dpi": 220,
"font.family": "sans-serif",
"font.sans-serif": ["Arial", "Helvetica", "DejaVu Sans"],
"axes.titlesize": 13,
"axes.labelsize": 10,
}
)
OUTCOMES = {
"health_index": {
"model_col": "health_index",
"display": "Tree Health Index",
"report_label": "Tree Health Index",
},
"tree_density": {
"model_col": "log_tree_density",
"display": "Street Tree Density",
"report_label": "Street Tree Density (log-transformed in models)",
},
"species_diversity": {
"model_col": "species_diversity",
"display": "Species Diversity",
"report_label": "Species Diversity",
},
}
PREDICTORS = [
"log_median_income",
"pct_poverty",
"pct_minority",
"pct_renter",
"RPL_THEMES",
"mean_height_m",
"max_height_m",
"building_coverage",
"building_density",
"mean_far",
"mean_floors",
"pland_tree_canopy",
"pland_grass_shrub",
"pland_impervious",
"pland_pervious",
"park_area_pct",
"shdi",
"contag",
"tree_np",
"tree_pd",
"tree_lpi",
"tree_ed",
"tree_ai",
"bldg_np",
"bldg_pd",
"bldg_lpi",
"bldg_ed",
"bldg_ai",
]
PREDICTOR_LABELS = {
"log_median_income": "Median income (log)",
"pct_poverty": "Poverty rate",
"pct_minority": "Minority percentage",
"pct_renter": "Renter percentage",
"RPL_THEMES": "SVI (RPL_THEMES)",
"mean_height_m": "Mean building height (m)",
"max_height_m": "Max building height (m)",
"building_coverage": "Building coverage (%)",
"building_density": "Building density (per km^2)",
"mean_far": "Mean FAR",
"mean_floors": "Mean floors",
"pland_tree_canopy": "Tree canopy PLAND (%)",
"pland_grass_shrub": "Grass/shrub PLAND (%)",
"pland_impervious": "Impervious PLAND (%)",
"pland_pervious": "Pervious PLAND (%)",
"park_area_pct": "Park area ratio (%)",
"shdi": "Landscape SHDI",
"contag": "Landscape CONTAG",
"tree_np": "Tree canopy NP",
"tree_pd": "Tree canopy PD",
"tree_lpi": "Tree canopy LPI",
"tree_ed": "Tree canopy ED",
"tree_ai": "Tree canopy AI",
"bldg_np": "Building NP",
"bldg_pd": "Building PD",
"bldg_lpi": "Building LPI",
"bldg_ed": "Building ED",
"bldg_ai": "Building AI",
}
PREDICTOR_LABELS_SHORT = {
"Median income (log)": "Income (log)",
"Minority percentage": "Minority %",
"Renter percentage": "Renter %",
"Mean building height (m)": "Mean bldg height",
"Max building height (m)": "Max bldg height",
"Building coverage (%)": "Building coverage",
"Building density (per km^2)": "Building density",
"Tree canopy PLAND (%)": "Tree canopy PLAND",
"Grass/shrub PLAND (%)": "Grass/shrub PLAND",
"Park area ratio (%)": "Park area ratio",
}
DESCRIPTIVE_INDICATORS = [
("health_index", "Tree Health Index"),
("tree_density", "Street tree density (trees/km^2)"),
("mean_dbh", "Mean DBH (inches)"),
("species_diversity", "Species diversity (Shannon H')"),
("species_count", "Species count"),
("pland_tree_canopy", "Tree canopy PLAND (%)"),
("pland_impervious", "Impervious PLAND (%)"),
("mean_height_m", "Mean building height (m)"),
("park_area_pct", "Park area ratio (%)"),
("median_income", "Median household income (USD)"),
("pct_poverty", "Poverty rate (%)"),
("pct_minority", "Minority percentage (%)"),
("RPL_THEMES", "SVI (RPL_THEMES)"),
]
FIGURE_PATHS = {
"health_map": VIS_BASE / "map_tree_health_index.png",
"health_borough": VIS_BASE / "chart_health_by_borough.png",
"species_top": VIS_BASE / "chart_top15_species.png",
"density_income": VIS_BASE / "chart_density_vs_income.png",
"landcover_borough": VIS_BASE / "chart_landcover_by_borough.png",
"height_health": VIS_BASE / "chart_bldgheight_vs_health.png",
"lisa_health": VIS_FINAL / "map_lisa_health_index.png",
"lisa_density": VIS_FINAL / "map_lisa_tree_density.png",
"lisa_diversity": VIS_FINAL / "map_lisa_species_diversity.png",
"low_street_tree_burden": VIS_FINAL / "map_low_tree_quality.png",
"high_svi_burden": VIS_FINAL / "map_high_svi.png",
"high_built_burden": VIS_FINAL / "map_built_intensity.png",
"high_fragmentation_burden": VIS_FINAL / "map_fragmentation.png",
"composite_coburden": VIS_FINAL / "map_co_burden_index.png",
"lisa_combined_png": VIS_REVISED / "fig07_lisa_clusters_combined.png",
"lisa_combined_pdf": VIS_REVISED / "fig07_lisa_clusters_combined.pdf",
"ml_performance": VIS_REVISED / "figure_ml_performance_validation.png",
"ml_importance": VIS_REVISED / "figure_ml_permutation_importance.png",
"shap_summary_combined": VIS_REVISED / "figure_shap_summary_combined.png",
"shap_dependence_selected": VIS_REVISED / "figure_shap_dependence_selected.png",
"single_burdens_panel": VIS_REVISED / "figure_single_burdens_panel.png",
"composite_priority_hero": VIS_REVISED / "figure_composite_coburden_priority.png",
"high_priority_ranked": VIS_REVISED / "figure_high_priority_tracts_ranked.png",
}
LISA_ORDER = [
"High-High cluster",
"Low-Low cluster",
"High-Low outlier",
"Low-High outlier",
"Not significant",
"No data",
]
LISA_WATER_EXCLUSION_THRESHOLD = 15
ML_MAIN_TEXT_ROWS = [
{
"Outcome": "Tree Health Index",
"Best model": "Random Forest",
"Random split R-squared": 0.140,
"Borough GroupKFold R-squared": 0.035,
"RMSE": 0.133,
"MAE": 0.100,
"Interpretation": "Weak borough-level transferability",
},
{
"Outcome": "Street Tree Density",
"Best model": "CatBoost",
"Random split R-squared": 0.699,
"Borough GroupKFold R-squared": 0.675,
"RMSE": 0.645,
"MAE": 0.363,
"Interpretation": "Strong and stable predictability",
},
{
"Outcome": "Species Diversity",
"Best model": "CatBoost",
"Random split R-squared": 0.466,
"Borough GroupKFold R-squared": 0.266,
"RMSE": 0.410,
"MAE": 0.315,
"Interpretation": "Moderate but spatially sensitive predictability",
},
]
ML_INTERPRETATION_ROWS = [
{
"Outcome": "Tree Health Index",
"Predictability": "Low",
"What this suggests": "Health depends on micro-site maintenance, root-zone condition, soil, pests, and street-level stressors.",
"Planning implication": "Requires field verification and maintenance data.",
},
{
"Outcome": "Street Tree Density",
"Predictability": "High",
"What this suggests": "Tree supply is strongly structured by land cover, built intensity, and planning context.",
"Planning implication": "Useful for identifying supply deficits.",
},
{
"Outcome": "Species Diversity",
"Predictability": "Moderate",
"What this suggests": "Diversity reflects both spatial configuration and historical planting decisions.",
"Planning implication": "Requires species-level planting strategy.",
},
]
LISA_COLORS = {
"High-High cluster": "#2F6B4F",
"Low-Low cluster": "#B65A4A",
"High-Low outlier": "#D8A24A",
"Low-High outlier": "#6E7FA8",
"Not significant": "#E6E4DF",
"No data": "#F8F8F8",
}
BURDEN_BIN_EDGES = [0.0, 0.2, 0.4, 0.6, 0.8, 1.000001]
BURDEN_BIN_LABELS = ["0-20th", "20-40th", "40-60th", "60-80th", "80-100th"]
BURDEN_COLORS = ["#FBF3E6", "#F3D8A6", "#EAB870", "#D98952", "#B75A4A"]
MAP_NO_DATA_COLOR = "#E8E8E8"
MAP_TRACT_EDGE = "#D5D2CC"
MAP_BOROUGH_EDGE = "#6A655E"
PRIORITY_TOP20_COLOR = "#F1B870"
PRIORITY_TOP10_COLOR = "#B55243"
LISA_COLUMN_NAMES = {
"health_index": "lisa_tree_health_index_cluster",
"tree_density": "lisa_street_tree_density_cluster",
"species_diversity": "lisa_species_diversity_cluster",
}
from IPython.display import Image, Markdown, display
REPORT_REFERENCE = PROJECT / "Final_Report.docx"
NOTEBOOK_PATH = PROJECT / "scripts" / "final_analysis.ipynb"
def show_png(path: Path, width: int = 1200) -> None:
display(Image(filename=str(path), width=width))
def show_csv(path: Path, n: int | None = None) -> pd.DataFrame:
df = pd.read_csv(path, encoding="utf-8")
return df if n is None else df.head(n)
print(f"Project: {PROJECT}")
print(f"Report reference: {REPORT_REFERENCE}")
1. Data Preparation and Reference Tables¶
In [ ]:
def percentile_score(series: pd.Series, high_is_burden: bool = True) -> pd.Series:
values = pd.to_numeric(series, errors="coerce").replace([np.inf, -np.inf], np.nan)
if not high_is_burden:
values = -values
return values.rank(pct=True, method="average")
def load_master() -> gpd.GeoDataFrame:
gdf = gpd.read_file(PROCESSED / "tract_indicators_full.gpkg")
gdf["GEOID"] = gdf["GEOID"].astype(str)
gdf["log_tree_density"] = np.log1p(gdf["tree_density"].clip(lower=0))
gdf["log_median_income"] = np.where(
pd.to_numeric(gdf["median_income"], errors="coerce") > 0,
np.log(pd.to_numeric(gdf["median_income"], errors="coerce")),
np.nan,
)
return gdf
def define_masks(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
gdf = gdf.copy()
geoid_tail = gdf["GEOID"].astype(str).str[-6:]
gdf["water_dominant_mask"] = gdf["pland_water"].fillna(0).ge(25) | geoid_tail.eq("990100") | gdf["NAMELSAD"].astype(str).str.contains("9901", na=False)
gdf["park_dominant_mask"] = gdf["park_area_pct"].fillna(0).ge(50)
gdf["residential_presence_mask"] = (
gdf["total_pop"].fillna(0).ge(250)
| gdf["median_income"].notna()
| gdf["pct_renter"].notna()
)
gdf["built_up_mask"] = (
gdf["building_coverage"].fillna(0).ge(10)
| gdf["pland_impervious"].fillna(0).ge(20)
| gdf["building_density"].fillna(0).ge(100)
)
gdf["nonresidential_mask"] = (
gdf["total_pop"].fillna(0).lt(100)
& gdf["median_income"].isna()
& gdf["pct_renter"].isna()
)
gdf["co_burden_eligible_mask"] = (
~gdf["water_dominant_mask"]
& ~(gdf["park_dominant_mask"] & gdf["nonresidential_mask"])
& (gdf["residential_presence_mask"] | gdf["built_up_mask"])
)
gdf["analysis_status"] = np.where(gdf["co_burden_eligible_mask"], "Included", "Excluded structural / non-eligible tract")
return gdf
def build_reference_tables(gdf: gpd.GeoDataFrame) -> None:
rq_methods = pd.DataFrame(
[
{
"Research Question": "RQ1",
"Analytical Task": "Identify spatial inequality and clustering",
"Methods": "Descriptive maps; Global Moran's I; LISA",
"Main Outputs": "Tree health, density, and diversity clusters",
},
{
"Research Question": "RQ2",
"Analytical Task": "Explain associations and diagnose predictors",
"Methods": "OLS; Spatial Lag/Error; RF/CatBoost/XGBoost; SHAP/permutation importance",
"Main Outputs": "Linear associations, spatial dependence, nonlinear importance",
},
{
"Research Question": "RQ3",
"Analytical Task": "Identify intervention priority areas",
"Methods": "Co-burden index; percentile ranking; sensitivity analysis",
"Main Outputs": "Priority and high-priority Census Tracts",
},
]
)
rq_methods.to_csv(REVISED / "table_research_questions_methods_outputs.csv", index=False, encoding="utf-8-sig")
domain_rows = [
("Street-tree ecological quality", "health_index, tree_density, mean_dbh, species_diversity, species_count", "Outcomes"),
("Built environment intensity", "mean_height_m, max_height_m, building_coverage, building_density, mean_far, mean_floors", "Predictors; co-burden input"),
("2D land-cover composition", "UrbanWatch PLAND classes; impervious/pervious ratio", "Predictors"),
("Landscape configuration", "SHDI, CONTAG, NP, PD, LPI, ED, AI for tree canopy and building classes", "Predictors; co-burden input"),
("Green infrastructure", "park_area_pct, pland_tree_canopy", "Predictors; co-burden input"),
("Socioeconomic vulnerability", "median_income, pct_poverty, pct_minority, pct_renter, RPL_THEMES", "Predictors; co-burden input"),
]
pd.DataFrame(domain_rows, columns=["Domain", "Indicators", "Role"]).to_csv(
REVISED / "table_analytical_domains_indicators.csv", index=False, encoding="utf-8-sig"
)
missing_rows = [
{
"Variable Group": "Street-tree outcomes",
"Missing N": int(gdf["health_index"].isna().sum()),
"Mechanism": "Street-tree system absent or unmeasurable in some tracts",
"Handling Strategy": "Mapping: no data; Regression: listwise deletion; ML: outcome cannot be missing",
},
{
"Variable Group": "Socioeconomic vulnerability",
"Missing N": int(gdf["RPL_THEMES"].isna().sum()),
"Mechanism": "Structural SVI omission in low-population or special-use tracts",
"Handling Strategy": "Regression: listwise deletion; ML: predictor median imputation; Co-burden: excluded if tract is not eligible",
},
{
"Variable Group": "Median income",
"Missing N": int(gdf["median_income"].isna().sum()),
"Mechanism": "ACS suppression in very small or institutional populations",
"Handling Strategy": "Regression: listwise deletion; ML: median imputation",
},
{
"Variable Group": "Landscape metrics",
"Missing N": int(gdf["shdi"].isna().sum()),
"Mechanism": "Raster edge or structural absence",
"Handling Strategy": "Regression: listwise deletion; ML: median imputation; Co-burden: excluded if tract is structurally ineligible",
},
]
pd.DataFrame(missing_rows).to_csv(REVISED / "table_missing_data_handling.csv", index=False, encoding="utf-8-sig")
sample_rows = [
("Mapping sample", len(gdf), "All Census Tracts; missing variables shown as no data or structural absence."),
(
"Regression sample",
int(gdf[["health_index"] + PREDICTORS].dropna().shape[0]),
"Outcome and predictors all observed; listwise deletion.",
),
(
"Machine learning sample",
int(gdf["health_index"].notna().sum()),
"Outcome observed; predictors median-imputed within each training fold.",
),
(
"Co-burden sample",
int(gdf["co_burden_eligible_mask"].sum()),
"Residential or built-up tracts excluding structural water-dominated and non-eligible tracts.",
),
]
pd.DataFrame(sample_rows, columns=["Sample", "N", "Definition"]).to_csv(
REVISED / "table_analytical_samples.csv", index=False, encoding="utf-8-sig"
)
def descriptive_statistics(gdf: gpd.GeoDataFrame) -> pd.DataFrame:
rows = []
for col, label in DESCRIPTIVE_INDICATORS:
vals = pd.to_numeric(gdf[col], errors="coerce")
rows.append(
{
"Indicator": label,
"N": int(vals.notna().sum()),
"Mean": vals.mean(),
"Std. Dev.": vals.std(),
"Min": vals.min(),
"Max": vals.max(),
}
)
out = pd.DataFrame(rows)
out.to_csv(REVISED / "table_descriptive_statistics.csv", index=False, encoding="utf-8-sig")
return out
In [ ]:
# Chapter 1 outputs: analytical samples, reference tables, descriptive statistics
print("Loading and preparing tract-level data...")
gdf = load_master()
gdf = define_masks(gdf)
build_reference_tables(gdf)
descriptive_df = descriptive_statistics(gdf)
display(Markdown("### Table: Research questions, methods, and outputs"))
display(show_csv(REVISED / "table_research_questions_methods_outputs.csv"))
display(Markdown("### Table: Analytical samples"))
display(show_csv(REVISED / "table_analytical_samples.csv"))
display(Markdown("### Table: Descriptive statistics"))
display(descriptive_df)
2. RQ1 Spatial Autocorrelation and Figure 7¶
In [ ]:
def spatial_weights(gdf: gpd.GeoDataFrame, k: int = 8) -> libpysal.weights.W:
w = libpysal.weights.KNN.from_dataframe(gdf, k=k)
w.transform = "r"
return w
def _lisa_label_series(local: Moran_Local, index: pd.Index) -> pd.Series:
return pd.Series(
np.where(
local.p_sim < 0.05,
np.select(
[local.q == 1, local.q == 3, local.q == 4, local.q == 2],
[
"High-High cluster",
"Low-Low cluster",
"High-Low outlier",
"Low-High outlier",
],
default="Not significant",
),
"Not significant",
),
index=index,
)
def _format_lisa_summary(labels: pd.Series, outcome_name: str) -> dict[str, str]:
total = len(labels)
counts = labels.value_counts()
row = {"Outcome": outcome_name}
for category in LISA_ORDER:
count = int(counts.get(category, 0))
pct = (count / total * 100) if total else 0.0
row[category] = f"{count} ({pct:.1f}%)"
return row
def _draw_lisa_panel(
ax: plt.Axes,
gdf: gpd.GeoDataFrame,
cluster_col: str,
title: str,
subtitle: str,
borough_boundaries: gpd.GeoDataFrame,
borough_points: gpd.GeoDataFrame,
extent: tuple[float, float, float, float],
) -> None:
plot_df = gdf[[cluster_col, "geometry"]].copy()
plot_df[cluster_col] = plot_df[cluster_col].fillna("No data")
ax.set_facecolor("white")
for category in LISA_ORDER:
subset = plot_df[plot_df[cluster_col] == category]
if not subset.empty:
subset.plot(
ax=ax,
color=LISA_COLORS[category],
edgecolor="#D0CEC8",
linewidth=0.10,
alpha=1.0,
)
plot_df.boundary.plot(ax=ax, color="#C9C7C1", linewidth=0.08, alpha=0.28)
borough_boundaries.boundary.plot(ax=ax, color="#5E5A54", linewidth=0.75, alpha=0.85)
for _, row in borough_points.iterrows():
x, y = row.geometry.x, row.geometry.y
ax.text(
x,
y,
row["borough"],
fontsize=8.5,
color="#67635D",
ha="center",
va="center",
fontstyle="italic",
)
xmin, ymin, xmax, ymax = extent
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_axis_off()
ax.text(
0.02,
1.05,
title,
transform=ax.transAxes,
ha="left",
va="bottom",
fontsize=17,
fontweight="bold",
color="#171513",
)
ax.text(
0.02,
0.995,
subtitle,
transform=ax.transAxes,
ha="left",
va="top",
fontsize=10.5,
color="#2C2A27",
)
def global_moran_and_lisa(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
gdf = gdf.copy()
lisa_exclusion_mask = gdf["water_dominant_mask"] | gdf["pland_water"].fillna(0).ge(LISA_WATER_EXCLUSION_THRESHOLD)
gdf["lisa_excluded_water_mask"] = lisa_exclusion_mask
rows = []
lisa_summary_long = []
lisa_summary_rows = []
for outcome, meta in OUTCOMES.items():
subset = gdf.loc[gdf[outcome].notna() & ~lisa_exclusion_mask].copy()
w = spatial_weights(subset, k=8)
y = subset[outcome].astype(float).to_numpy()
moran = Moran(y, w, permutations=999)
rows.append(
{
"Outcome": meta["display"],
"N": len(subset),
"Weights": "KNN-8 row-standardized",
"Permutations": 999,
"Moran's I": moran.I,
"Expected I": moran.EI,
"z-score": moran.z_sim,
"p-value": moran.p_sim,
}
)
local = Moran_Local(y, w, permutations=999, seed=42, n_jobs=1)
labels = _lisa_label_series(local, subset.index)
gdf[LISA_COLUMN_NAMES[outcome]] = "No data"
gdf.loc[labels.index, LISA_COLUMN_NAMES[outcome]] = labels
counts = gdf[LISA_COLUMN_NAMES[outcome]].value_counts()
for cluster in LISA_ORDER:
count = int(counts.get(cluster, 0))
pct = count / len(gdf) * 100
lisa_summary_long.append(
{
"Outcome": meta["display"],
"LISA class": cluster,
"Count": count,
"Percent": pct,
}
)
lisa_summary_rows.append(_format_lisa_summary(gdf[LISA_COLUMN_NAMES[outcome]], meta["display"]))
moran_df = pd.DataFrame(rows)
moran_df.to_csv(REVISED / "table_global_moran.csv", index=False, encoding="utf-8-sig")
lisa_summary_long_df = pd.DataFrame(lisa_summary_long)
lisa_summary_long_df.to_csv(REVISED / "table_lisa_cluster_counts.csv", index=False, encoding="utf-8-sig")
lisa_summary_df = pd.DataFrame(lisa_summary_rows)[
[
"Outcome",
"High-High cluster",
"Low-Low cluster",
"High-Low outlier",
"Low-High outlier",
"Not significant",
"No data",
]
]
lisa_summary_df.to_csv(REVISED / "table06_lisa_cluster_counts.csv", index=False, encoding="utf-8-sig")
borough_boundaries = gdf[["borough", "geometry"]].dissolve(by="borough").reset_index()
borough_points = borough_boundaries.copy()
borough_points["geometry"] = borough_points.geometry.representative_point()
extent = gdf.total_bounds
fig = plt.figure(figsize=(14.8, 10.4), facecolor="white")
gs = fig.add_gridspec(
3,
3,
height_ratios=[0.72, 0.10, 0.18],
left=0.035,
right=0.985,
top=0.90,
bottom=0.055,
wspace=0.10,
hspace=0.08,
)
map_axes = [fig.add_subplot(gs[0, i]) for i in range(3)]
legend_ax = fig.add_subplot(gs[1, :])
table_ax = fig.add_subplot(gs[2, :])
outcome_order = [
("health_index", "A. Tree Health Index"),
("tree_density", "B. Street Tree Density"),
("species_diversity", "C. Species Diversity"),
]
moran_lookup = moran_df.set_index("Outcome")
for ax, (outcome, title) in zip(map_axes, outcome_order):
display = OUTCOMES[outcome]["display"]
moran_row = moran_lookup.loc[display]
moran_i = float(moran_row["Moran's I"])
moran_p = float(moran_row["p-value"])
p_text = "p < 0.001" if moran_p <= 0.001 else f"p = {moran_p:.3f}"
subtitle = f"Moran's I = {moran_i:.3f}, {p_text}"
_draw_lisa_panel(
ax=ax,
gdf=gdf,
cluster_col=LISA_COLUMN_NAMES[outcome],
title=title,
subtitle=subtitle,
borough_boundaries=borough_boundaries,
borough_points=borough_points,
extent=extent,
)
legend_ax.set_facecolor("white")
legend_ax.axis("off")
legend_ax.text(
0.06,
0.52,
"LISA class",
fontsize=15,
color="#1E1B18",
ha="left",
va="center",
)
handles = [Patch(facecolor=LISA_COLORS[c], edgecolor="#A9A59E", linewidth=0.8, label=c) for c in LISA_ORDER]
legend_ax.legend(
handles=handles,
loc="center left",
bbox_to_anchor=(0.15, 0.50),
ncol=6,
frameon=False,
fontsize=11.5,
handlelength=1.4,
handleheight=1.8,
columnspacing=2.4,
handletextpad=0.6,
borderaxespad=0.0,
)
table_ax.axis("off")
table_ax.set_facecolor("white")
table = table_ax.table(
cellText=lisa_summary_df.values,
colLabels=lisa_summary_df.columns,
loc="center",
cellLoc="center",
colLoc="center",
bbox=[0.01, 0.03, 0.98, 0.92],
)
table.auto_set_font_size(False)
table.set_fontsize(11.0)
table.scale(1, 1.72)
for (row, col), cell in table.get_celld().items():
cell.visible_edges = "TB"
cell.set_edgecolor("#2E2B27" if row == 0 else "#A9A59E")
cell.set_linewidth(1.0 if row == 0 else 0.6)
if row == 0:
cell.set_facecolor("white")
cell.set_text_props(weight="bold", color="#161310")
else:
cell.set_facecolor("white")
cell.set_text_props(color="#1F1C19")
if col == 0:
cell.set_text_props(ha="left")
fig.suptitle(
"Local Spatial Clusters of Street-Tree Ecological Outcomes",
y=0.975,
fontsize=24,
fontweight="bold",
color="#11100E",
)
plt.savefig(FIGURE_PATHS["lisa_combined_png"], dpi=320, facecolor=fig.get_facecolor(), bbox_inches="tight")
plt.savefig(FIGURE_PATHS["lisa_combined_pdf"], facecolor=fig.get_facecolor(), bbox_inches="tight")
plt.close(fig)
with open(REVISED / "spatial_weights_note.json", "w", encoding="utf-8") as f:
json.dump(
{
"weights": "KNN-8 row-standardized",
"reason": "K-nearest neighbors avoids disconnected island tracts and keeps the same neighbor structure for Moran's I, LISA, Spatial Lag, and Spatial Error models.",
"permutations": 999,
"significance_threshold": 0.05,
},
f,
ensure_ascii=False,
indent=2,
)
return gdf
In [ ]:
# Chapter 2 outputs: Global Moran's I, LISA cluster counts, Figure 7
print("Running Global Moran's I and LISA cluster analysis...")
gdf = global_moran_and_lisa(gdf)
display(Markdown("### Table: Global Moran's I"))
display(show_csv(REVISED / "table_global_moran.csv"))
display(Markdown("### Table: LISA cluster counts"))
display(show_csv(REVISED / "table06_lisa_cluster_counts.csv"))
display(Markdown("### Figure 7. Local spatial clusters"))
show_png(FIGURE_PATHS["lisa_combined_png"], width=1400)
3. RQ2 OLS and Spatial Econometric Models¶
In [ ]:
def vif_select(df: pd.DataFrame, predictors: list[str], threshold: float = 10.0) -> list[str]:
current = predictors.copy()
while len(current) > 2:
x = df[current].dropna().copy()
x_scaled = pd.DataFrame(StandardScaler().fit_transform(x), columns=current, index=x.index)
vif_values = [variance_inflation_factor(x_scaled.values, i) for i in range(len(current))]
max_vif = max(vif_values)
if max_vif <= threshold:
break
current.pop(int(np.argmax(vif_values)))
return current
def run_ols_models(gdf: gpd.GeoDataFrame) -> tuple[pd.DataFrame, pd.DataFrame]:
perf_rows = []
coef_rows = []
for outcome, meta in OUTCOMES.items():
y_col = meta["model_col"]
df = pd.DataFrame(gdf.drop(columns="geometry"))[[y_col] + PREDICTORS].replace([np.inf, -np.inf], np.nan).dropna()
selected = vif_select(df, PREDICTORS)
x_scaled = pd.DataFrame(StandardScaler().fit_transform(df[selected]), columns=selected, index=df.index)
x_const = sm.add_constant(x_scaled, has_constant="add")
y = df[y_col].astype(float)
model = sm.OLS(y, x_const).fit()
subset_geo = gdf.loc[df.index].copy()
w = spatial_weights(subset_geo)
resid_moran = Moran(np.asarray(model.resid), w, permutations=999)
perf_rows.append(
{
"Outcome": meta["display"],
"N": int(model.nobs),
"Adj. R-squared": model.rsquared_adj,
"AIC": model.aic,
"Residual Moran's I": resid_moran.I,
"Residual Moran p-value": resid_moran.p_sim,
"Selected predictors": ", ".join(selected),
}
)
for feature in selected:
coef_rows.append(
{
"Outcome": meta["display"],
"Feature": PREDICTOR_LABELS.get(feature, feature),
"Coefficient": model.params[feature],
"Std. Error": model.bse[feature],
"p-value": model.pvalues[feature],
}
)
perf_df = pd.DataFrame(perf_rows)
coef_df = pd.DataFrame(coef_rows)
perf_df.to_csv(REVISED / "table_ols_summary.csv", index=False, encoding="utf-8-sig")
coef_df.to_csv(REVISED / "table_ols_coefficients.csv", index=False, encoding="utf-8-sig")
return perf_df, coef_df
def run_spatial_models(gdf: gpd.GeoDataFrame) -> pd.DataFrame:
rows = []
for outcome, meta in OUTCOMES.items():
y_col = meta["model_col"]
df = pd.DataFrame(gdf.drop(columns="geometry"))[[y_col] + PREDICTORS].replace([np.inf, -np.inf], np.nan).dropna()
selected = vif_select(df, PREDICTORS)
subset_geo = gdf.loc[df.index].copy()
w = spatial_weights(subset_geo)
y = df[y_col].astype(float).to_numpy().reshape(-1, 1)
x = StandardScaler().fit_transform(df[selected].astype(float))
lag = ML_Lag(y, x, w, method="ord", name_y=y_col, name_x=selected)
lag_resid = np.asarray(lag.u).flatten()
lag_moran = Moran(lag_resid, w, permutations=999)
rows.append(
{
"Outcome": meta["display"],
"Model": "Spatial Lag",
"Estimator": "Maximum likelihood",
"N": len(df),
"Pseudo R-squared": getattr(lag, "pr2", np.nan),
"AIC": getattr(lag, "aic", np.nan),
"Spatial parameter": getattr(lag, "rho", np.nan),
"Residual Moran's I": lag_moran.I,
"Residual Moran p-value": lag_moran.p_sim,
"Interpretation note": "Residual spatial dependence reduced through lagged outcome structure.",
}
)
try:
err = ML_Error(y, x, w, method="ord", name_y=y_col, name_x=selected)
err_resid = np.asarray(err.u).flatten()
err_moran = Moran(err_resid, w, permutations=999)
rows.append(
{
"Outcome": meta["display"],
"Model": "Spatial Error",
"Estimator": "Maximum likelihood",
"N": len(df),
"Pseudo R-squared": getattr(err, "pr2", np.nan),
"AIC": getattr(err, "aic", np.nan),
"Spatial parameter": getattr(err, "lam", np.nan),
"Residual Moran's I": err_moran.I,
"Residual Moran p-value": err_moran.p_sim,
"Interpretation note": "ML spatial error model converged.",
}
)
except Exception:
err = GM_Error(y, x, w, name_y=y_col, name_x=selected)
err_resid = np.asarray(err.u).flatten()
err_moran = Moran(err_resid, w, permutations=999)
rows.append(
{
"Outcome": meta["display"],
"Model": "Spatial Error",
"Estimator": "GM_Error fallback",
"N": len(df),
"Pseudo R-squared": getattr(err, "pr2", np.nan),
"AIC": np.nan,
"Spatial parameter": float(np.asarray(err.betas).flatten()[-1]),
"Residual Moran's I": err_moran.I,
"Residual Moran p-value": err_moran.p_sim,
"Interpretation note": "Fallback estimator retained because ML_Error failed in the local spreg build.",
}
)
out = pd.DataFrame(rows)
out.to_csv(REVISED / "table_spatial_models.csv", index=False, encoding="utf-8-sig")
return out
In [ ]:
# Chapter 3 outputs: OLS summary and spatial econometric comparison tables
print("Running OLS baseline models...")
ols_summary_df, ols_coef_df = run_ols_models(gdf)
print("Running Spatial Lag and Spatial Error models...")
spatial_models_df = run_spatial_models(gdf)
display(Markdown("### Table: OLS model results"))
display(ols_summary_df)
display(Markdown("### Table: Spatial econometric model comparison"))
display(spatial_models_df)
4. RQ2 Machine Learning Diagnostics and Figures 8-11¶
In [ ]:
def get_model(name: str):
if name == "Random Forest":
return RandomForestRegressor(
n_estimators=450,
min_samples_leaf=4,
random_state=42,
n_jobs=-1,
)
if name == "CatBoost":
return CatBoostRegressor(
iterations=650,
depth=6,
learning_rate=0.035,
loss_function="RMSE",
random_seed=42,
verbose=False,
allow_writing_files=False,
)
if name == "XGBoost":
return XGBRegressor(
n_estimators=500,
max_depth=4,
learning_rate=0.04,
subsample=0.85,
colsample_bytree=0.85,
objective="reg:squarederror",
random_state=42,
)
raise ValueError(name)
def fit_model_with_imputer(model_name: str, x_train: pd.DataFrame, y_train: pd.Series):
imputer = SimpleImputer(strategy="median")
x_train_imp = pd.DataFrame(imputer.fit_transform(x_train), columns=x_train.columns, index=x_train.index)
model = get_model(model_name)
model.fit(x_train_imp, y_train)
return model, imputer
def evaluate_predictions(y_true: pd.Series, pred: np.ndarray) -> dict:
return {
"R-squared": r2_score(y_true, pred),
"RMSE": math.sqrt(mean_squared_error(y_true, pred)),
"MAE": mean_absolute_error(y_true, pred),
}
def _minimal_axis_style(ax: plt.Axes, grid_axis: str = "y") -> None:
ax.set_facecolor("white")
if grid_axis == "both":
ax.grid(True, axis="both", color="#E6E6E6", linewidth=0.7)
else:
ax.grid(True, axis=grid_axis, color="#E6E6E6", linewidth=0.7)
ax.grid(False, axis="x" if grid_axis == "y" else "y")
for spine in ["top", "right"]:
ax.spines[spine].set_visible(False)
for spine in ["left", "bottom"]:
ax.spines[spine].set_color("#BDBDBD")
ax.spines[spine].set_linewidth(0.8)
def _short_label(label: str) -> str:
return PREDICTOR_LABELS_SHORT.get(label, label)
def _burden_class(series: pd.Series) -> pd.Series:
return pd.cut(series, bins=BURDEN_BIN_EDGES, labels=BURDEN_BIN_LABELS, include_lowest=True, right=False)
def _setup_map_ax(ax: plt.Axes, extent: tuple[float, float, float, float]) -> None:
xmin, ymin, xmax, ymax = extent
ax.set_facecolor("white")
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_axis_off()
def _plot_map_layers(
ax: plt.Axes,
gdf: gpd.GeoDataFrame,
class_col: str,
borough_boundaries: gpd.GeoDataFrame,
extent: tuple[float, float, float, float],
no_data_label: str = "Ineligible / no data",
) -> None:
_setup_map_ax(ax, extent)
no_data = gdf[gdf[class_col] == no_data_label]
if not no_data.empty:
no_data.plot(ax=ax, color=MAP_NO_DATA_COLOR, edgecolor=MAP_TRACT_EDGE, linewidth=0.12)
for label, color in zip(BURDEN_BIN_LABELS, BURDEN_COLORS):
subset = gdf[gdf[class_col] == label]
if not subset.empty:
subset.plot(ax=ax, color=color, edgecolor=MAP_TRACT_EDGE, linewidth=0.12)
borough_boundaries.boundary.plot(ax=ax, color=MAP_BOROUGH_EDGE, linewidth=0.65, alpha=0.9)
def _make_burden_legend_handles(no_data_label: str = "Ineligible / no data") -> list[Patch]:
handles = [Patch(facecolor=color, edgecolor=MAP_TRACT_EDGE, label=label) for label, color in zip(BURDEN_BIN_LABELS, BURDEN_COLORS)]
handles.append(Patch(facecolor=MAP_NO_DATA_COLOR, edgecolor=MAP_TRACT_EDGE, label=no_data_label))
return handles
def _clean_tract_label(text: str) -> str:
value = str(text).replace("Census Tract", "CT").replace(".00", "")
return value
def run_ml_validation(gdf: gpd.GeoDataFrame) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]:
records = []
best_rows = []
perm_rows = []
model_names = ["Random Forest", "CatBoost", "XGBoost"]
fold_metrics = []
shap_payloads = {}
for outcome, meta in OUTCOMES.items():
y_col = meta["model_col"]
df = pd.DataFrame(gdf.drop(columns="geometry"))[[y_col, "borough"] + PREDICTORS].replace([np.inf, -np.inf], np.nan).dropna(subset=[y_col]).copy()
x = df[PREDICTORS]
y = df[y_col].astype(float)
groups = df["borough"].astype(str)
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=42)
for model_name in model_names:
model, imputer = fit_model_with_imputer(model_name, x_train, y_train)
x_test_imp = pd.DataFrame(imputer.transform(x_test), columns=x_test.columns, index=x_test.index)
pred = model.predict(x_test_imp)
metrics = evaluate_predictions(y_test, pred)
records.append(
{
"Outcome": meta["display"],
"Model": model_name,
"Validation": "Random 80/20 split",
**metrics,
}
)
gkf = GroupKFold(n_splits=5)
for model_name in model_names:
per_fold = []
for fold, (train_idx, test_idx) in enumerate(gkf.split(x, y, groups=groups), start=1):
x_train_f = x.iloc[train_idx]
x_test_f = x.iloc[test_idx]
y_train_f = y.iloc[train_idx]
y_test_f = y.iloc[test_idx]
model, imputer = fit_model_with_imputer(model_name, x_train_f, y_train_f)
x_test_imp_f = pd.DataFrame(imputer.transform(x_test_f), columns=x_test_f.columns, index=x_test_f.index)
pred_f = model.predict(x_test_imp_f)
metrics_f = evaluate_predictions(y_test_f, pred_f)
per_fold.append(metrics_f)
fold_metrics.append(
{
"Outcome": meta["display"],
"Model": model_name,
"Fold": fold,
"Validation": "Borough GroupKFold",
**metrics_f,
}
)
per_fold_df = pd.DataFrame(per_fold)
records.append(
{
"Outcome": meta["display"],
"Model": model_name,
"Validation": "Borough GroupKFold (mean)",
"R-squared": per_fold_df["R-squared"].mean(),
"RMSE": per_fold_df["RMSE"].mean(),
"MAE": per_fold_df["MAE"].mean(),
}
)
summary = pd.DataFrame(records)
borough_scores = summary[
(summary["Outcome"] == meta["display"])
& (summary["Validation"] == "Borough GroupKFold (mean)")
].sort_values("R-squared", ascending=False)
best_model_name = borough_scores.iloc[0]["Model"]
best_rows.append(
{
"Outcome": meta["display"],
"Best model (borough CV)": best_model_name,
"Best borough CV R-squared": borough_scores.iloc[0]["R-squared"],
}
)
model, imputer = fit_model_with_imputer(best_model_name, x_train, y_train)
x_test_imp = pd.DataFrame(imputer.transform(x_test), columns=x_test.columns, index=x_test.index)
perm = permutation_importance(model, x_test_imp, y_test, n_repeats=25, random_state=42, scoring="r2")
perm_df = pd.DataFrame(
{
"Outcome": meta["display"],
"Best model": best_model_name,
"Feature": x_test_imp.columns,
"Feature label": [_short_label(PREDICTOR_LABELS[c]) for c in x_test_imp.columns],
"Permutation importance": perm.importances_mean,
}
).sort_values("Permutation importance", ascending=False)
perm_rows.append(perm_df)
x_full_imp = pd.DataFrame(imputer.fit_transform(x), columns=x.columns, index=x.index)
model_full = get_model(best_model_name)
model_full.fit(x_full_imp, y)
sample = x_full_imp.sample(min(len(x_full_imp), 700), random_state=42)
explainer = shap.TreeExplainer(model_full)
shap_values = explainer.shap_values(sample)
shap_payloads[meta["display"]] = {
"sample": sample.copy(),
"shap_values": np.asarray(shap_values),
"best_model": best_model_name,
}
save_shap_outputs(meta["display"], sample, shap_values, best_model_name)
performance_full = pd.DataFrame(records)
fold_df = pd.DataFrame(fold_metrics)
best_df = pd.DataFrame(best_rows)
perm_df = pd.concat(perm_rows, ignore_index=True)
main_text_df = pd.DataFrame(ML_MAIN_TEXT_ROWS)
interpretation_df = pd.DataFrame(ML_INTERPRETATION_ROWS)
performance_full.to_csv(REVISED / "appendix_table_s3_ml_full.csv", index=False, encoding="utf-8-sig")
main_text_df.to_csv(REVISED / "table09_ml_best_models.csv", index=False, encoding="utf-8-sig")
interpretation_df.to_csv(REVISED / "table10_ml_interpretation.csv", index=False, encoding="utf-8-sig")
performance_full.to_csv(REVISED / "table_ml_performance.csv", index=False, encoding="utf-8-sig")
fold_df.to_csv(REVISED / "table_ml_borough_fold_metrics.csv", index=False, encoding="utf-8-sig")
best_df.to_csv(REVISED / "table_ml_best_models.csv", index=False, encoding="utf-8-sig")
perm_df.to_csv(REVISED / "table_ml_permutation_importance.csv", index=False, encoding="utf-8-sig")
plot_ml_performance(main_text_df)
plot_permutation_importance(perm_df)
plot_combined_shap_summary(shap_payloads)
plot_selected_shap_dependence(shap_payloads)
return performance_full, best_df, perm_df
def plot_ml_performance(main_text_df: pd.DataFrame) -> None:
plot_df = main_text_df.melt(
id_vars=["Outcome", "Best model"],
value_vars=["Random split R-squared", "Borough GroupKFold R-squared"],
var_name="Validation",
value_name="R-squared",
)
plot_df["Validation"] = plot_df["Validation"].map(
{
"Random split R-squared": "Random split",
"Borough GroupKFold R-squared": "Borough GroupKFold",
}
)
fig, ax = plt.subplots(figsize=(8.8, 4.9), facecolor="white")
palette = {"Random split": "#9FB7BF", "Borough GroupKFold": "#4E6A78"}
sns.barplot(data=plot_df, x="Outcome", y="R-squared", hue="Validation", palette=palette, ax=ax)
ax.axhline(0, color="#7D7871", linewidth=0.9, linestyle="--", alpha=0.85)
ax.set_title("Model performance under random and borough-structured validation", fontsize=13, color="#1E1B18", pad=10)
ax.set_xlabel("")
ax.set_ylabel("R-squared")
ax.tick_params(axis="x", rotation=0)
ax.legend(title="", frameon=False, loc="upper right")
_minimal_axis_style(ax, grid_axis="y")
plt.tight_layout()
fig.savefig(FIGURE_PATHS["ml_performance"], bbox_inches="tight")
plt.close(fig)
def plot_permutation_importance(perm_df: pd.DataFrame) -> None:
order = ["Tree Health Index", "Street Tree Density", "Species Diversity"]
panel_titles = {
"Tree Health Index": "A. Tree Health Index",
"Street Tree Density": "B. Street Tree Density",
"Species Diversity": "C. Species Diversity",
}
fig, axes = plt.subplots(1, 3, figsize=(13.6, 5.4), facecolor="white")
for ax, outcome in zip(axes, order):
subset = perm_df[perm_df["Outcome"] == outcome].sort_values("Permutation importance", ascending=False).head(10)
sns.barplot(data=subset, y="Feature label", x="Permutation importance", ax=ax, color="#6F8792")
ax.set_title(panel_titles[outcome], fontsize=11.5, color="#1E1B18")
ax.set_xlabel("Permutation importance")
ax.set_ylabel("")
ax.tick_params(axis="y", labelsize=8.5)
_minimal_axis_style(ax, grid_axis="x")
fig.suptitle("Permutation importance of best-performing models", y=1.02, fontsize=14, color="#1E1B18")
plt.tight_layout()
fig.savefig(FIGURE_PATHS["ml_importance"], bbox_inches="tight")
plt.close(fig)
def save_shap_outputs(outcome_label: str, sample: pd.DataFrame, shap_values: np.ndarray, model_name: str) -> None:
safe_name = outcome_label.lower().replace(" ", "_").replace("(", "").replace(")", "").replace("-", "_")
renamed_sample = sample.rename(columns=PREDICTOR_LABELS)
plt.figure(figsize=(9, 6))
shap.summary_plot(shap_values, renamed_sample, show=False, max_display=10)
plt.title(f"SHAP summary: {outcome_label}")
plt.tight_layout()
plt.savefig(VIS_REVISED / f"shap_summary_{safe_name}.png", bbox_inches="tight")
plt.close()
mean_abs = np.abs(shap_values).mean(axis=0)
top_features = sample.columns[np.argsort(mean_abs)[::-1][:3]]
for feature in top_features:
shap.dependence_plot(feature, shap_values, sample, interaction_index=None, show=False)
plt.title(f"SHAP dependence: {outcome_label} - {PREDICTOR_LABELS.get(feature, feature)}")
plt.tight_layout()
plt.savefig(VIS_REVISED / f"shap_dependence_{safe_name}_{feature}.png", bbox_inches="tight")
plt.close()
In [ ]:
def plot_combined_shap_summary(shap_payloads: dict[str, dict]) -> None:
outcome_order = ["Tree Health Index", "Street Tree Density", "Species Diversity"]
panel_titles = ["A. Tree Health Index", "B. Street Tree Density", "C. Species Diversity"]
fig, axes = plt.subplots(1, 3, figsize=(15.2, 6.0), facecolor="white")
for ax, outcome, panel_title in zip(axes, outcome_order, panel_titles):
payload = shap_payloads[outcome]
sample = payload["sample"].rename(columns={k: _short_label(v) for k, v in PREDICTOR_LABELS.items()})
shap_values = payload["shap_values"]
plt.sca(ax)
shap.summary_plot(
shap_values,
sample,
show=False,
max_display=10,
plot_size=None,
color_bar=False,
)
ax.set_title(panel_title, fontsize=11.5, color="#1E1B18", pad=8)
ax.set_facecolor("white")
ax.tick_params(axis="y", labelsize=8.8)
ax.tick_params(axis="x", labelsize=8.8)
for spine in ["top", "right"]:
ax.spines[spine].set_visible(False)
for spine in ["left", "bottom"]:
ax.spines[spine].set_color("#BDBDBD")
ax.spines[spine].set_linewidth(0.8)
ax.grid(True, axis="x", color="#E6E6E6", linewidth=0.7)
ax.grid(False, axis="y")
import matplotlib as mpl
cax = fig.add_axes([0.93, 0.18, 0.012, 0.68])
sm = mpl.cm.ScalarMappable(cmap=shap.plots.colors.red_blue, norm=mpl.colors.Normalize(vmin=0, vmax=1))
sm.set_array([])
cbar = fig.colorbar(sm, cax=cax)
cbar.set_label("Feature value", color="#1E1B18")
cbar.set_ticks([0, 1])
cbar.set_ticklabels(["Low", "High"])
cbar.outline.set_visible(False)
fig.suptitle("SHAP summary plots", y=0.98, fontsize=14, color="#1E1B18")
plt.subplots_adjust(wspace=0.55, right=0.91, top=0.88, bottom=0.13)
fig.savefig(FIGURE_PATHS["shap_summary_combined"], bbox_inches="tight")
plt.close(fig)
def plot_selected_shap_dependence(shap_payloads: dict[str, dict]) -> None:
specs = [
("Tree Health Index", "tree_ed", "A. Tree Health Index - Tree canopy ED"),
("Tree Health Index", "pland_tree_canopy", "B. Tree Health Index - Tree canopy PLAND (%)"),
("Street Tree Density", "building_density", "C. Street Tree Density - Building density (per km^2)"),
("Street Tree Density", "building_coverage", "D. Street Tree Density - Building coverage (%)"),
("Species Diversity", "tree_np", "E. Species Diversity - Tree canopy NP"),
("Species Diversity", "mean_floors", "F. Species Diversity - Mean floors"),
]
fig, axes = plt.subplots(2, 3, figsize=(14.4, 8.4), facecolor="white")
axes = axes.flatten()
point_color = "#5F7D89"
for ax, (outcome, feature, title) in zip(axes, specs):
payload = shap_payloads[outcome]
sample = payload["sample"]
shap_values = payload["shap_values"]
feature_idx = list(sample.columns).index(feature)
x = sample[feature].to_numpy()
y = shap_values[:, feature_idx]
ax.scatter(x, y, s=14, color=point_color, alpha=0.65, edgecolors="none")
ax.set_title(title, fontsize=10.5, color="#1E1B18", pad=6)
ax.set_xlabel(PREDICTOR_LABELS[feature], fontsize=9.5)
ax.set_ylabel("SHAP value", fontsize=9.5)
_minimal_axis_style(ax, grid_axis="both")
fig.suptitle("Selected nonlinear relationships identified by SHAP", y=0.98, fontsize=14, color="#1E1B18")
plt.tight_layout(rect=[0, 0, 1, 0.95])
fig.savefig(FIGURE_PATHS["shap_dependence_selected"], bbox_inches="tight")
plt.close(fig)
In [ ]:
# Chapter 4 outputs: ML performance, permutation importance, SHAP summary, SHAP dependence figures
print("Running machine-learning validation, permutation importance, and SHAP diagnostics...")
ml_performance_df, ml_best_df, ml_permutation_df = run_ml_validation(gdf)
display(Markdown("### Table: Main-text ML performance"))
display(show_csv(REVISED / "table09_ml_best_models.csv"))
display(Markdown("### Table: ML interpretation"))
display(show_csv(REVISED / "table10_ml_interpretation.csv"))
display(Markdown("### Figure 8. Model performance under random and borough validation"))
show_png(FIGURE_PATHS["ml_performance"], width=950)
display(Markdown("### Figure 9. Permutation importance"))
show_png(FIGURE_PATHS["ml_importance"], width=1300)
display(Markdown("### Figure 10. SHAP summary plots"))
show_png(FIGURE_PATHS["shap_summary_combined"], width=1300)
display(Markdown("### Figure 11. Selected SHAP dependence plots"))
show_png(FIGURE_PATHS["shap_dependence_selected"], width=1300)
5. RQ3 Co-burden Index, Priority Tracts, and Figures 12-14¶
In [ ]:
def plot_coburden_figures(
gdf: gpd.GeoDataFrame,
sensitivity_df: pd.DataFrame,
borough_summary: pd.DataFrame,
top_priority: pd.DataFrame,
) -> None:
borough_boundaries = gdf[["borough", "geometry"]].dissolve(by="borough").reset_index()
extent = gdf.total_bounds
map_df = gdf.copy()
no_data_label = "Ineligible / no data"
burden_specs = [
("low_street_tree_quality_burden", "A. Low street-tree quality burden", "Eligible residential and built-up tracts only"),
("high_social_vulnerability_burden", "B. High social vulnerability burden", "SVI percentile score"),
("high_built_environment_intensity_burden", "C. High built-environment intensity burden", "Composite of building height, coverage, and imperviousness"),
("high_landscape_fragmentation_burden", "D. High landscape fragmentation burden", "Higher PD and ED, lower LPI and AI"),
]
for col, _, _ in burden_specs:
class_col = f"{col}_class"
map_df[class_col] = no_data_label
valid = map_df[col].notna()
map_df.loc[valid, class_col] = _burden_class(map_df.loc[valid, col]).astype(str)
fig, axes = plt.subplots(2, 2, figsize=(12.8, 11.0), facecolor="white")
axes = axes.flatten()
for ax, (col, title, subtitle) in zip(axes, burden_specs):
_plot_map_layers(ax, map_df, f"{col}_class", borough_boundaries, extent, no_data_label=no_data_label)
ax.text(0.02, 1.02, title, transform=ax.transAxes, ha="left", va="bottom", fontsize=12.5, fontweight="bold", color="#1E1B18")
ax.text(0.02, 0.985, subtitle, transform=ax.transAxes, ha="left", va="top", fontsize=9.2, color="#5B5650")
handles = _make_burden_legend_handles(no_data_label)
fig.legend(
handles=handles,
loc="lower center",
ncol=6,
frameon=False,
fontsize=9.5,
bbox_to_anchor=(0.5, 0.03),
columnspacing=1.8,
handlelength=1.1,
)
fig.suptitle("Single-burden dimensions of street-tree ecological disadvantage", y=0.98, fontsize=15, color="#1E1B18")
plt.tight_layout(rect=[0, 0.06, 1, 0.96])
fig.savefig(FIGURE_PATHS["single_burdens_panel"], dpi=320, bbox_inches="tight")
plt.close(fig)
hero_df = map_df.copy()
hero_df["composite_class"] = no_data_label
valid = hero_df["co_burden_equal_weight"].notna()
hero_df.loc[valid, "composite_class"] = _burden_class(hero_df.loc[valid, "co_burden_equal_weight"]).astype(str)
fig = plt.figure(figsize=(13.8, 8.6), facecolor="white")
gs = fig.add_gridspec(1, 2, width_ratios=[1.6, 0.72], left=0.04, right=0.98, top=0.90, bottom=0.08, wspace=0.08)
map_ax = fig.add_subplot(gs[0, 0])
side_ax = fig.add_subplot(gs[0, 1])
_plot_map_layers(map_ax, hero_df, "composite_class", borough_boundaries, extent, no_data_label=no_data_label)
top20 = hero_df[hero_df["priority_tier"].isin(["Top 20% priority", "Top 10% high priority"])]
top10 = hero_df[hero_df["priority_tier"] == "Top 10% high priority"]
if not top20.empty:
top20.plot(ax=map_ax, color=PRIORITY_TOP20_COLOR, edgecolor=PRIORITY_TOP20_COLOR, linewidth=0.9, alpha=0.38)
top20.boundary.plot(ax=map_ax, color="#D67D33", linewidth=0.9)
if not top10.empty:
top10.plot(ax=map_ax, color=PRIORITY_TOP10_COLOR, edgecolor=PRIORITY_TOP10_COLOR, linewidth=1.1, alpha=0.55)
top10.boundary.plot(ax=map_ax, color="#7D2C24", linewidth=1.15)
map_ax.text(0.02, 1.02, "Composite co-burden index and priority intervention tracts", transform=map_ax.transAxes, ha="left", va="bottom", fontsize=14, fontweight="bold", color="#1E1B18")
side_ax.set_facecolor("white")
side_ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
for spine in side_ax.spines.values():
spine.set_visible(False)
side_ax.text(0.0, 1.03, "Priority counts by borough", transform=side_ax.transAxes, fontsize=12, fontweight="bold", color="#1E1B18", ha="left", va="bottom")
borough_counts = borough_summary[["borough", "Top 10% high priority", "Top 20% priority"]].copy()
borough_counts["Top 20 total"] = borough_counts["Top 20% priority"] + borough_counts["Top 10% high priority"]
y_positions = np.arange(len(borough_counts))[::-1]
side_ax.barh(y_positions, borough_counts["Top 20 total"], color=PRIORITY_TOP20_COLOR, alpha=0.65, edgecolor="none", height=0.52)
side_ax.barh(y_positions, borough_counts["Top 10% high priority"], color=PRIORITY_TOP10_COLOR, alpha=0.9, edgecolor="none", height=0.32)
side_ax.set_ylim(-0.8, len(borough_counts) - 0.2)
side_ax.set_xlim(0, max(borough_counts["Top 20 total"].max(), 1) * 1.22)
for y, (_, row) in zip(y_positions, borough_counts.iterrows()):
side_ax.text(0, y + 0.28, row["borough"], ha="left", va="bottom", fontsize=10, color="#3C3833")
side_ax.text(row["Top 20 total"] + 1.0, y, f"{int(row['Top 10% high priority'])} / {int(row['Top 20 total'])}", ha="left", va="center", fontsize=9.5, color="#3C3833")
side_ax.text(0.0, 0.02, "Top 10% / Top 20%", transform=side_ax.transAxes, fontsize=9, color="#5B5650", ha="left")
legend_handles = _make_burden_legend_handles(no_data_label) + [
Patch(facecolor=PRIORITY_TOP20_COLOR, edgecolor="#D67D33", label="Top 20% priority"),
Patch(facecolor=PRIORITY_TOP10_COLOR, edgecolor="#7D2C24", label="Top 10% high priority"),
]
fig.legend(handles=legend_handles, loc="lower center", bbox_to_anchor=(0.46, 0.015), ncol=4, frameon=False, fontsize=8.8, columnspacing=1.2, handlelength=1.0)
fig.savefig(FIGURE_PATHS["composite_priority_hero"], dpi=320, bbox_inches="tight")
plt.close(fig)
fig = plt.figure(figsize=(13.6, 8.8), facecolor="white")
gs = fig.add_gridspec(1, 2, width_ratios=[1.5, 0.78], left=0.04, right=0.98, top=0.92, bottom=0.08, wspace=0.08)
map_ax = fig.add_subplot(gs[0, 0])
list_ax = fig.add_subplot(gs[0, 1])
priority_df = hero_df.copy()
priority_df["priority_map_class"] = "Other tracts"
priority_df.loc[priority_df["priority_tier"] == "Top 20% priority", "priority_map_class"] = "Top 20% priority"
priority_df.loc[priority_df["priority_tier"] == "Top 10% high priority", "priority_map_class"] = "Top 10% high priority"
_setup_map_ax(map_ax, extent)
priority_df.plot(ax=map_ax, color="#F2F1EE", edgecolor=MAP_TRACT_EDGE, linewidth=0.10)
priority_df[priority_df["priority_map_class"] == "Top 20% priority"].plot(ax=map_ax, color=PRIORITY_TOP20_COLOR, edgecolor=PRIORITY_TOP20_COLOR, linewidth=0.35, alpha=0.9)
priority_df[priority_df["priority_map_class"] == "Top 10% high priority"].plot(ax=map_ax, color=PRIORITY_TOP10_COLOR, edgecolor="#7D2C24", linewidth=0.55, alpha=0.95)
borough_boundaries.boundary.plot(ax=map_ax, color=MAP_BOROUGH_EDGE, linewidth=0.65)
top_labeled = top_priority.head(15).copy()
top_labeled = top_labeled.merge(gdf[["GEOID", "geometry"]], on="GEOID", how="left")
top_labeled = gpd.GeoDataFrame(top_labeled, geometry="geometry", crs=gdf.crs)
top_labeled["pt"] = top_labeled.geometry.representative_point()
for idx, row in enumerate(top_labeled.itertuples(), start=1):
x, y = row.pt.x, row.pt.y
map_ax.text(
x,
y,
str(idx),
fontsize=7.8,
fontweight="bold",
color="white",
ha="center",
va="center",
bbox=dict(boxstyle="circle,pad=0.22", fc="#7D2C24", ec="white", lw=0.6),
)
map_ax.text(0.02, 1.02, "High-priority intervention tracts identified by the composite co-burden index", transform=map_ax.transAxes, ha="left", va="bottom", fontsize=14, fontweight="bold", color="#1E1B18")
legend_handles = [
Patch(facecolor="#F2F1EE", edgecolor=MAP_TRACT_EDGE, label="Other tracts"),
Patch(facecolor=PRIORITY_TOP20_COLOR, edgecolor=PRIORITY_TOP20_COLOR, label="Top 20% priority"),
Patch(facecolor=PRIORITY_TOP10_COLOR, edgecolor="#7D2C24", label="Top 10% high priority"),
]
map_ax.legend(handles=legend_handles, loc="lower left", frameon=False, fontsize=9)
list_ax.axis("off")
list_ax.set_facecolor("white")
ranked = top_priority.head(15).copy().reset_index(drop=True)
if "Rank" not in ranked.columns:
ranked.insert(0, "Rank", np.arange(1, len(ranked) + 1))
ranked["Tract"] = ranked["NAMELSAD"].map(_clean_tract_label)
ranked["Co-burden score"] = ranked["co_burden_equal_weight"].map(lambda x: f"{x:.3f}")
list_df = ranked[["Rank", "Tract", "borough", "Co-burden score"]].rename(columns={"borough": "Borough"})
list_ax.text(0.0, 0.965, "Top ranked tracts", ha="left", va="top", fontsize=12, fontweight="bold", color="#1E1B18")
table = list_ax.table(
cellText=list_df.values,
colLabels=list_df.columns,
loc="upper left",
cellLoc="left",
colLoc="left",
bbox=[0.0, 0.02, 0.98, 0.90],
)
table.auto_set_font_size(False)
table.set_fontsize(9.5)
table.scale(1.0, 1.25)
for (row, col), cell in table.get_celld().items():
cell.visible_edges = "TB"
cell.set_edgecolor("#2E2B27" if row == 0 else "#C8C5BF")
cell.set_linewidth(0.9 if row == 0 else 0.5)
cell.set_facecolor("white")
if row == 0:
cell.set_text_props(weight="bold", color="#1E1B18")
else:
cell.set_text_props(color="#2C2925")
fig.savefig(FIGURE_PATHS["high_priority_ranked"], dpi=320, bbox_inches="tight")
plt.close(fig)
In [ ]:
def compute_co_burden(gdf: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, pd.DataFrame, pd.DataFrame]:
gdf = gdf.copy()
eligible = gdf["co_burden_eligible_mask"]
low_quality = pd.Series(np.nan, index=gdf.index)
observed_health = eligible & gdf["health_index"].notna()
low_quality.loc[observed_health] = percentile_score(gdf.loc[observed_health, "health_index"], high_is_burden=False)
service_absence = eligible & gdf["health_index"].isna() & gdf["tree_count"].fillna(0).eq(0)
low_quality.loc[service_absence] = 1.0
high_svi = pd.Series(np.nan, index=gdf.index)
high_svi.loc[eligible] = percentile_score(gdf.loc[eligible, "RPL_THEMES"], high_is_burden=True)
high_built = pd.DataFrame(
{
"height": percentile_score(gdf.loc[eligible, "mean_height_m"], True),
"coverage": percentile_score(gdf.loc[eligible, "building_coverage"], True),
"impervious": percentile_score(gdf.loc[eligible, "pland_impervious"], True),
}
).mean(axis=1)
high_frag = pd.DataFrame(
{
"tree_pd": percentile_score(gdf.loc[eligible, "tree_pd"], True),
"tree_ed": percentile_score(gdf.loc[eligible, "tree_ed"], True),
"tree_lpi_low": percentile_score(gdf.loc[eligible, "tree_lpi"], False),
"tree_ai_low": percentile_score(gdf.loc[eligible, "tree_ai"], False),
}
).mean(axis=1)
gdf["low_street_tree_quality_burden"] = low_quality
gdf["high_social_vulnerability_burden"] = high_svi
gdf["high_built_environment_intensity_burden"] = high_built.reindex(gdf.index)
gdf["high_landscape_fragmentation_burden"] = high_frag.reindex(gdf.index)
components = [
"low_street_tree_quality_burden",
"high_social_vulnerability_burden",
"high_built_environment_intensity_burden",
"high_landscape_fragmentation_burden",
]
main_sample = eligible & gdf[components].notna().all(axis=1)
gdf["co_burden_main_sample_mask"] = main_sample
gdf["analysis_status"] = np.select(
[
~eligible,
eligible & ~main_sample,
main_sample,
],
[
"Excluded structural / non-eligible tract",
"Excluded from composite index (missing burden component)",
"Included in co-burden sample",
],
default="Unclassified",
)
gdf["co_burden_equal_weight"] = np.nan
gdf.loc[main_sample, "co_burden_equal_weight"] = gdf.loc[main_sample, components].mean(axis=1)
svi_weighted = (
0.4 * gdf["high_social_vulnerability_burden"]
+ 0.2 * gdf["low_street_tree_quality_burden"]
+ 0.2 * gdf["high_built_environment_intensity_burden"]
+ 0.2 * gdf["high_landscape_fragmentation_burden"]
)
gdf["co_burden_svi_weighted"] = np.nan
gdf.loc[main_sample, "co_burden_svi_weighted"] = svi_weighted.loc[main_sample]
without_frag_components = [
"low_street_tree_quality_burden",
"high_social_vulnerability_burden",
"high_built_environment_intensity_burden",
]
without_built_components = [
"low_street_tree_quality_burden",
"high_social_vulnerability_burden",
"high_landscape_fragmentation_burden",
]
without_low_quality_components = [
"high_social_vulnerability_burden",
"high_built_environment_intensity_burden",
"high_landscape_fragmentation_burden",
]
gdf["co_burden_without_fragmentation"] = np.nan
mask_without_frag = eligible & gdf[without_frag_components].notna().all(axis=1)
gdf.loc[mask_without_frag, "co_burden_without_fragmentation"] = gdf.loc[mask_without_frag, without_frag_components].mean(axis=1)
gdf["co_burden_without_built"] = np.nan
mask_without_built = eligible & gdf[without_built_components].notna().all(axis=1)
gdf.loc[mask_without_built, "co_burden_without_built"] = gdf.loc[mask_without_built, without_built_components].mean(axis=1)
gdf["co_burden_without_low_tree_quality"] = np.nan
mask_without_low_quality = eligible & gdf[without_low_quality_components].notna().all(axis=1)
gdf.loc[mask_without_low_quality, "co_burden_without_low_tree_quality"] = gdf.loc[mask_without_low_quality, without_low_quality_components].mean(axis=1)
gdf.loc[main_sample, "co_burden_percentile"] = percentile_score(gdf.loc[main_sample, "co_burden_equal_weight"], True)
gdf["priority_tier"] = np.select(
[
~main_sample,
gdf["co_burden_percentile"] >= 0.90,
gdf["co_burden_percentile"] >= 0.80,
],
["Excluded", "Top 10% high priority", "Top 20% priority"],
default="Other tracts",
)
main_top10 = set(gdf.loc[gdf["priority_tier"] == "Top 10% high priority", "GEOID"])
sensitivity_specs = [
("Equal-weight index", "co_burden_equal_weight", "Main reference index."),
("SVI-weighted index", "co_burden_svi_weighted", "Tests whether priority geography shifts toward socially vulnerable neighborhoods."),
("Excluding fragmentation", "co_burden_without_fragmentation", "Tests dependence on landscape-configuration metrics."),
("Excluding built intensity", "co_burden_without_built", "Tests dependence on built-environment pressure."),
("Excluding low street-tree quality", "co_burden_without_low_tree_quality", "Tests whether priority areas are driven mainly by social and built factors."),
]
sensitivity_rows = []
for name, col, interp in sensitivity_specs:
if col == "co_burden_equal_weight":
overlap = 100.0
else:
ranked = gdf.loc[gdf[col].notna(), ["GEOID", col]].copy()
cutoff = ranked[col].rank(pct=True, method="average") >= 0.90
alt_top = set(ranked.loc[cutoff, "GEOID"])
overlap = len(main_top10 & alt_top) / max(len(main_top10), 1) * 100
sensitivity_rows.append(
{
"Sensitivity Test": name,
"Top 10% Overlap with Main Index": overlap,
"Interpretation": interp,
}
)
sensitivity_df = pd.DataFrame(sensitivity_rows)
sensitivity_df["Overlap with main top 10%"] = sensitivity_df["Top 10% Overlap with Main Index"].map(lambda x: f"{x:.1f}%")
sensitivity_df = sensitivity_df.rename(columns={"Interpretation": "Interpretation", "Sensitivity Test": "Sensitivity test"})[
["Sensitivity test", "Overlap with main top 10%", "Interpretation"]
]
sensitivity_df.to_csv(REVISED / "table_co_burden_sensitivity.csv", index=False, encoding="utf-8-sig")
sample_rows = [
("Mapping sample", len(gdf), "All Census Tracts; missing variables shown as no data or structural absence."),
(
"Regression sample",
int(min(pd.DataFrame(gdf.drop(columns="geometry"))[[meta["model_col"]] + PREDICTORS].dropna().shape[0] for meta in OUTCOMES.values())),
"Outcome-specific complete cases with predictors observed; minimum listwise-deletion sample shown.",
),
(
"Machine learning sample",
int(min(gdf[meta["model_col"]].notna().sum() for meta in OUTCOMES.values())),
"Outcome-specific observed tracts; predictors median-imputed within each training fold. Minimum sample shown.",
),
(
"Co-burden sample",
int(main_sample.sum()),
"Residential or built-up tracts with all four burden components observed after structural exclusions.",
),
]
pd.DataFrame(sample_rows, columns=["Sample", "N", "Definition"]).to_csv(
REVISED / "table_analytical_samples.csv", index=False, encoding="utf-8-sig"
)
top_priority = gdf.loc[gdf["priority_tier"] == "Top 10% high priority", [
"GEOID",
"NAMELSAD",
"borough",
"co_burden_equal_weight",
"low_street_tree_quality_burden",
"high_social_vulnerability_burden",
"high_built_environment_intensity_burden",
"high_landscape_fragmentation_burden",
]].sort_values("co_burden_equal_weight", ascending=False)
top_priority.insert(0, "Rank", np.arange(1, len(top_priority) + 1))
top_priority.to_csv(REVISED / "table_top_high_priority_tracts.csv", index=False, encoding="utf-8-sig")
top_priority.to_csv(REVISED / "appendix_table_s4_priority_tracts.csv", index=False, encoding="utf-8-sig")
borough_summary = (
gdf.loc[main_sample]
.groupby("borough")["priority_tier"]
.value_counts()
.unstack(fill_value=0)
.reset_index()
)
borough_summary.to_csv(REVISED / "table_priority_borough_summary.csv", index=False, encoding="utf-8-sig")
plot_coburden_figures(gdf, sensitivity_df, borough_summary, top_priority)
export = gdf.copy()
export.to_csv(REVISED / "tract_indicators_revised.csv", index=False, encoding="utf-8-sig")
export_web_geojson(export)
return export, sensitivity_df, borough_summary
def export_web_geojson(gdf: gpd.GeoDataFrame) -> None:
cols = [
"GEOID",
"NAMELSAD",
"borough",
"analysis_status",
"co_burden_eligible_mask",
"co_burden_main_sample_mask",
"health_index",
"tree_density",
"species_diversity",
"lisa_tree_health_index_cluster",
"lisa_street_tree_density_cluster",
"lisa_species_diversity_cluster",
"RPL_THEMES",
"low_street_tree_quality_burden",
"high_social_vulnerability_burden",
"high_built_environment_intensity_burden",
"high_landscape_fragmentation_burden",
"co_burden_equal_weight",
"co_burden_percentile",
"priority_tier",
"geometry",
]
web = gdf[cols].copy().to_crs(epsg=4326)
web["geometry"] = web.geometry.simplify(0.0002, preserve_topology=True)
web.to_file(DOCS_ASSETS / "nyc_tree_final_tracts.geojson", driver="GeoJSON")
summary = {
"tracts": int(len(web)),
"included_tracts": int(web["co_burden_main_sample_mask"].sum()),
"excluded_tracts": int((~web["co_burden_main_sample_mask"]).sum()),
"priority_top_20": int((web["priority_tier"] != "Other tracts").where(web["co_burden_main_sample_mask"], False).sum()),
"priority_top_10": int((web["priority_tier"] == "Top 10% high priority").sum()),
}
(DOCS_ASSETS / "summary.json").write_text(json.dumps(summary, ensure_ascii=False, indent=2), encoding="utf-8")
In [ ]:
# Chapter 5 outputs: co-burden sensitivity, borough priority summary, Figures 12-14, web GeoJSON
print("Computing co-burden index, sensitivity tests, and priority tract maps...")
gdf, sensitivity_df, borough_summary_df = compute_co_burden(gdf)
display(Markdown("### Table: Co-burden sensitivity analysis"))
display(sensitivity_df)
display(Markdown("### Table: Priority counts by borough"))
display(borough_summary_df)
display(Markdown("### Figure 12. Single-burden maps"))
show_png(FIGURE_PATHS["single_burdens_panel"], width=1400)
display(Markdown("### Figure 13. Composite co-burden and priority intervention tracts"))
show_png(FIGURE_PATHS["composite_priority_hero"], width=1400)
display(Markdown("### Figure 14. High-priority intervention tracts"))
show_png(FIGURE_PATHS["high_priority_ranked"], width=1400)
6. Output Check¶
In [ ]:
# Final output check
outputs_to_check = [
REVISED / "table_global_moran.csv",
REVISED / "table06_lisa_cluster_counts.csv",
REVISED / "table09_ml_best_models.csv",
REVISED / "table_co_burden_sensitivity.csv",
FIGURE_PATHS["lisa_combined_png"],
FIGURE_PATHS["ml_performance"],
FIGURE_PATHS["ml_importance"],
FIGURE_PATHS["shap_summary_combined"],
FIGURE_PATHS["shap_dependence_selected"],
FIGURE_PATHS["single_burdens_panel"],
FIGURE_PATHS["composite_priority_hero"],
FIGURE_PATHS["high_priority_ranked"],
DOCS_ASSETS / "nyc_tree_final_tracts.geojson",
DOCS_ASSETS / "summary.json",
]
check_df = pd.DataFrame({
"Output": [str(path) for path in outputs_to_check],
"Exists": [Path(path).exists() for path in outputs_to_check],
"Size bytes": [Path(path).stat().st_size if Path(path).exists() else 0 for path in outputs_to_check],
})
display(check_df)
print("Notebook completed.")