Cross-Sectional Model Comparison: neighbayes vs spreg¶

import geopandas as gpd
import libpysal
import numpy as np
import pandas as pd
import spreg
from IPython.display import display
from spreg.sputils import _sp_effects, spmultiplier

from neighbayes.models import SAR, SDEM, SDM, SEM, SLX

columbus = gpd.read_file(libpysal.examples.get_path("columbus.shp"))
columbus.plot("HOVAL", scheme="quantiles")

w = libpysal.graph.Graph.build_contiguity(columbus).transform("r")

y = columbus.HOVAL.values
X = columbus[["INC", "CRIME"]]

spreg_results = {
    "SLX": spreg.OLS(y, X, w=w, slx_lags=1, name_y="HOVAL", name_x=["INC", "CRIME"]),
    "SAR": spreg.ML_Lag(y, X, w=w, name_y="HOVAL", name_x=["INC", "CRIME"]),
    "SEM": spreg.ML_Error(y, X, w=w, name_y="HOVAL", name_x=["INC", "CRIME"]),
    "SDM": spreg.ML_Lag(
        y, X, w=w, slx_lags=1, w_lags=2, name_y="HOVAL", name_x=["INC", "CRIME"]
    ),
    "SDEM": spreg.GM_Error(
        y, X, w=w, slx_lags=1, name_y="HOVAL", name_x=["INC", "CRIME"]
    ),
}

spec = "HOVAL ~  1 + INC + CRIME"
sample_kwargs = dict(draws=1000, tune=1000, chains=4, random_seed=42, progressbar=False)
bayes_models = {
    "SLX": SLX(spec, data=columbus, W=w),
    "SAR": SAR(spec, data=columbus, W=w),
    "SEM": SEM(spec, data=columbus, W=w),
    "SDM": SDM(spec, data=columbus, W=w),
    "SDEM": SDEM(spec, data=columbus, W=w),
}
bayes_idata = {name: model.fit(**sample_kwargs) for name, model in bayes_models.items()}

bayes_models["SAR"].summary()

from neighbayes.diagnostics import spatial_mcmc_diagnostic

spatial_mcmc_diagnostic(bayes_models["SAR"], emit_warnings=False).to_frame()

print(spreg_results["SAR"].summary)

def extract_spreg_params(model):
    names = list(getattr(model, "name_x", []))
    betas = np.asarray(model.betas).flatten()
    if len(betas) == len(names) + 1 and type(model).__name__ == "GM_Lag":
        names = names + ["rho"]
    if len(names) != len(betas):
        names = [f"param_{i}" for i in range(len(betas))]
    return pd.Series(betas, index=names, dtype=float)


def extract_neighbayes_params(model_name, model, idata):
    out = {}
    beta_mean = idata.posterior["beta"].mean(("chain", "draw")).values
    beta_names = (
        model._beta_names()
        if model_name in {"SLX", "SDM", "SDEM"}
        else list(model._feature_names)
    )
    for name, val in zip(beta_names, beta_mean):
        out[name] = float(val)
    if model_name in {"SAR", "SDM"}:
        out["rho"] = float(idata.posterior["rho"].mean())
    if model_name in {"SEM", "SDEM"}:
        out["lambda"] = float(idata.posterior["lam"].mean())
    harmonized = {}
    for k, v in out.items():
        key = k.replace("W*", "W_")
        if key.lower() == "intercept":
            key = "CONSTANT"
        harmonized[key] = v
    return pd.Series(harmonized, dtype=float)


def bayes_effects_df(model):
    eff = model.spatial_effects()
    # spatial_effects() returns a DataFrame indexed by feature names
    return pd.DataFrame(
        {
            "feature": eff.index.tolist(),
            "direct_neighbayes": eff["direct"].values,
            "indirect_neighbayes": eff["indirect"].values,
            "total_neighbayes": eff["total"].values,
        }
    )


def spreg_effects_df(model_name, sp_model):
    if model_name == "SEM":
        # SEM has no spatial multiplier on the systematic component,
        # so direct = coefficient, indirect = 0, total = coefficient.
        # neighbayes excludes the intercept from effects, so we
        # only compare non-constant covariates here.
        sp_params = extract_spreg_params(sp_model)
        features = [f for f in ["INC", "CRIME"] if f in sp_params.index]
        return pd.DataFrame(
            {
                "feature": features,
                "direct_spreg": [float(sp_params[f]) for f in features],
                "indirect_spreg": [0.0 for _ in features],
                "total_spreg": [float(sp_params[f]) for f in features],
            }
        )

    output = sp_model.output.copy()
    output["regime"] = "global"

    def classify_var(var_name):
        if var_name == "CONSTANT":
            return "o"
        if var_name in {"lambda", "W_HOVAL"}:
            return "rho"
        if var_name.startswith("W_"):
            return "wx"
        return "x"

    output["var_type"] = output["var_names"].map(classify_var)
    sp_model.output = output

    vars_x = output.query("var_type == 'x'")
    rho = (
        float(np.squeeze(getattr(sp_model, "rho", 0.0)))
        if hasattr(sp_model, "rho")
        else 0.0
    )
    multipliers = spmultiplier(w, rho)
    slx_lags = getattr(sp_model, "slx_lags", 0) or (
        1 if model_name in {"SLX", "SDEM"} else 0
    )
    slx_vars = getattr(sp_model, "slx_vars", None)
    if slx_vars is None:
        slx_vars = "All"

    total_spreg, direct_spreg, indirect_spreg = _sp_effects(
        sp_model,
        vars_x,
        multipliers,
        slx_lags=slx_lags,
        slx_vars=slx_vars,
    )

    return pd.DataFrame(
        {
            "feature": vars_x["var_names"].tolist(),
            "direct_spreg": np.asarray(direct_spreg).flatten(),
            "indirect_spreg": np.asarray(indirect_spreg).flatten(),
            "total_spreg": np.asarray(total_spreg).flatten(),
        }
    )


coef_rows = []
effect_tables = {}

for model_name in ["SLX", "SAR", "SEM", "SDM", "SDEM"]:
    sp = extract_spreg_params(spreg_results[model_name])
    by = extract_neighbayes_params(
        model_name, bayes_models[model_name], bayes_idata[model_name]
    )

    common = sorted(set(sp.index).intersection(set(by.index)))
    for param in common:
        coef_rows.append(
            {
                "model": model_name,
                "parameter": param,
                "neighbayes_posterior_mean": by[param],
                "spreg_estimate": sp[param],
                "difference": by[param] - sp[param],
            }
        )

    beff = bayes_effects_df(bayes_models[model_name]).copy()
    seff = spreg_effects_df(model_name, spreg_results[model_name]).copy()
    merged = beff.merge(seff, on="feature", how="inner")
    if not merged.empty:
        merged["direct_difference"] = (
            merged["direct_neighbayes"] - merged["direct_spreg"]
        )
        merged["indirect_difference"] = (
            merged["indirect_neighbayes"] - merged["indirect_spreg"]
        )
        merged["total_difference"] = merged["total_neighbayes"] - merged["total_spreg"]
        merged.insert(0, "model", model_name)
    effect_tables[model_name] = merged

comparison = (
    pd.DataFrame(coef_rows).sort_values(["model", "parameter"]).reset_index(drop=True)
)
spatial_effects_comparison = pd.concat(
    [tbl for tbl in effect_tables.values() if not tbl.empty],
    ignore_index=True,
)

display(comparison)

display(spatial_effects_comparison)

from neighbayes import dgp
from neighbayes.models import SAR

# Set random seed for reproducibility
np.random.seed(42)

# True parameters
beta_true = np.array([1.0, -1.0])
rho_true = 0.5
sigma_true = 1.0

gdf = dgp.simulate_sar(
    n=40,
    beta=beta_true,
    rho=rho_true,
    sigma=sigma_true,
    create_gdf=True,
    contiguity="queen",
)

w = libpysal.graph.Graph.build_contiguity(gdf, rook=False).transform("r")

# Fit SAR model
model = SAR("y ~ X_1", W=w, data=gdf)
res = model.fit()

model.summary()

model.spatial_effects()

sim_spreg = spreg.GM_Lag(
    gdf["y"],
    gdf["X_1"],
    w=w,
)

print(sim_spreg.summary)