Clustered standard errors differ between `DuckRegression` and `statsmodels`

[simolasalla]

When running OLS with clustered data, the standard errors from DuckRegression and statsmodels differ by a non-trivial amount.

Expected behavior:

• Point estimates match (or differ by a tiny amount)
• Standard errors match (or differ by a tiny amount)

Actual behavior:

• Point estimates match (or differ by a tiny amount)
• Standard error differ by a non-trivial amount. In my run with example dataset, statsmodels reports SE ≈ 0.032 for the slope and duckreg reports SE ≈ 0.0251.

What I tried:

• Reproduced on various real data (one example dataset).
• Reproduced on synthetic data (script below).
• Varying n_bootstraps

Replication with simulated data:

# synthetic_repro.py
import numpy as np
import pandas as pd
import duckdb
from duckreg.estimators import DuckRegression
import statsmodels.formula.api as smf

np.random.seed(123)
n_clusters = 50
cluster_size = 10
N = n_clusters * cluster_size

cluster_ids = np.repeat(np.arange(n_clusters), cluster_size)
class_size = np.random.normal(30, 5, size=N)
# cluster random intercept
u = np.repeat(np.random.normal(0, 2, size=n_clusters), cluster_size)
eps = np.random.normal(0, 5, size=N)
# outcome
id_score = 50 + 0.5 * class_size + u + eps

data = pd.DataFrame({
    "id_score": id_score,
    "class_size": class_size,
    "class_id": cluster_ids
})

# statsmodels
model = smf.ols("id_score ~ class_size", data=data).fit()
ms = model.get_robustcov_results(cov_type='cluster', groups=data['class_id'])
print("statsmodels cluster SE:\n", ms.summary())

# duckreg
conn = duckdb.connect("database.db")
conn.execute("DROP TABLE IF EXISTS test_data")
conn.execute("CREATE TABLE test_data AS SELECT * FROM data")
m = DuckRegression(
    db_name="database.db",
    table_name="test_data",
    formula="id_score ~ class_size",
    cluster_col="class_id",
    n_bootstraps=200,
    seed=21,
)
m.fit()
print("duckreg summary:")
display(m.summary())


pd.DataFrame([[ms.params[1], ms.bse[1]], [m.summary()['point_estimate'][1], m.summary()['standard_error'][1]]], columns=['Point estimate', 'Standard error'], index=['Statsmodels', 'DuckRegression'])

[apoorvalal]

Hi @simolasalla , thanks for the bugreport! This turned out to be a real insidious bug that arose out of my misunderstanding of some sql operators: turns out IN (...) implicitly dedupes an array that is being sampled, so the duplicate indices required for a cluster bootstrap were getting dropped, so we were actually doing some (ad-hoc) subsampling. I've pushed a patch in #20. Here is the crucial change - I now use a CTE to compute cluster multiplicities, then apply them as weights.

If you're working on something time sensitive, you're welcome to install from the branch [uv pip install git...@patch/clustboot ]. I will expand the PR to cover the panel data estimators as well and merge it after.

---

Details - Here be dragons

The desired behaviour is

resampled_rows = [3, 7, 7, 2, 44, ...]  # cluster 7 appears twice!

But the SQL query uses:

WHERE cluster_id IN (SELECT unnest([3, 7, 7, 2, 44, ...]))

IN (3, 7, 7, 2, 44) becomes IN (3, 7, 2, 44).

Expected:
50 cluster draws with replacement → ~31.6 unique clusters, but some appear 2-3+ times
Actual: Only ~31.6 unique clusters, each appearing exactly once
Result: Effective sample size is ~63% of what it should be
SE underestimation: √0.63 ≈ 0.794
This matches your observed SE ratio.

I expanded your minimal example (thanks for that, btw) to a larger script that varies (n_clusters, cluster_size, n_boot) and first replicated that your discrepancy magnitude was fairly stable, which tipped me off to a potential scaling issue.

I then worked through the bootstrap query with an LLM and diagnosed the issue. Rerunning the diagnostic gives much more reassuring results.

---

"""
diagnostic script
Tests SE discrepancy with varying configurations based on salla's example
"""

# %%
import numpy as np
import pandas as pd
import duckdb
from duckreg.estimators import DuckRegression
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import os

# %%


def run_single_comparison(
    n_clusters=50,
    cluster_size=10,
    n_bootstraps=200,
    icc_strength=2.0,
    eps_strength=5.0,
    seed=123,
):
    """
    Run a single comparison matching the original example structure.

    Parameters:
    -----------
    n_clusters : int
        Number of clusters
    cluster_size : int
        Observations per cluster
    n_bootstraps : int
        Number of bootstrap replications for duckreg
    icc_strength : float
        Standard deviation of cluster random intercept (controls ICC)
    eps_strength : float
        Standard deviation of idiosyncratic error
    seed : int
        Random seed for reproducibility
    """
    np.random.seed(seed)
    N = n_clusters * cluster_size
    cluster_ids = np.repeat(np.arange(n_clusters), cluster_size)
    class_size = np.random.normal(30, 5, size=N)

    # cluster random intercept
    u = np.repeat(np.random.normal(0, icc_strength, size=n_clusters), cluster_size)
    eps = np.random.normal(0, eps_strength, size=N)

    # outcome
    id_score = 50 + 0.5 * class_size + u + eps

    data = pd.DataFrame(
        {"id_score": id_score, "class_size": class_size, "class_id": cluster_ids}
    )

    # statsmodels
    model = smf.ols("id_score ~ class_size", data=data).fit()
    ms = model.get_robustcov_results(cov_type="cluster", groups=data["class_id"])

    # duckreg
    db_name = f"temp_db_{seed}.db"
    if os.path.exists(db_name):
        os.remove(db_name)

    conn = duckdb.connect(db_name)
    conn.execute("DROP TABLE IF EXISTS test_data")
    conn.execute("CREATE TABLE test_data AS SELECT * FROM data")
    conn.close()

    m = DuckRegression(
        db_name=db_name,
        table_name="test_data",
        formula="id_score ~ class_size",
        cluster_col="class_id",
        n_bootstraps=n_bootstraps,
        seed=seed,
    )
    m.fit()

    # Clean up
    if os.path.exists(db_name):
        os.remove(db_name)

    # Extract results for slope coefficient
    sm_coef = ms.params[1]
    sm_se = ms.bse[1]
    dr_coef = m.summary()["point_estimate"][1]
    dr_se = m.summary()["standard_error"][1]

    return {
        "n_clusters": n_clusters,
        "cluster_size": cluster_size,
        "N": N,
        "n_bootstraps": n_bootstraps,
        "statsmodels_coef": sm_coef,
        "statsmodels_se": sm_se,
        "duckreg_coef": dr_coef,
        "duckreg_se": dr_se,
        "se_ratio": dr_se / sm_se,
        "se_diff": dr_se - sm_se,
        "se_pct_diff": 100 * (dr_se - sm_se) / sm_se,
    }


def vary_clusters_and_size():
    """
    Test varying both number of clusters and cluster size.
    """
    print("\n" + "=" * 70)
    print("VARYING CLUSTERS AND CLUSTER SIZE")
    print("=" * 70 + "\n")

    # Configuration grid
    n_clusters_grid = [20, 30, 50, 75, 100, 150]
    cluster_size_grid = [5, 10, 20, 50]

    results = []

    for n_clusters in n_clusters_grid:
        for cluster_size in cluster_size_grid:
            print(f"Testing G={n_clusters}, n={cluster_size}...", end=" ")

            # Run multiple replications to account for bootstrap randomness
            for rep in range(5):
                seed = 123 + rep + n_clusters * 1000 + cluster_size * 100
                try:
                    res = run_single_comparison(
                        n_clusters=n_clusters,
                        cluster_size=cluster_size,
                        n_bootstraps=200,
                        seed=seed,
                    )
                    res["replication"] = rep
                    results.append(res)
                except Exception as e:
                    print(f"ERROR: {e}")
                    continue

            print("Done.")

    df = pd.DataFrame(results)

    # Aggregate across replications
    summary = (
        df.groupby(["n_clusters", "cluster_size"])
        .agg(
            {
                "N": "first",
                "statsmodels_se": "mean",
                "duckreg_se": "mean",
                "se_ratio": ["mean", "std"],
                "se_diff": ["mean", "std"],
                "se_pct_diff": ["mean", "std"],
            }
        )
        .round(4)
    )

    print("\nSUMMARY STATISTICS:")
    print(summary)

    return df, summary


def vary_bootstrap_replications():
    """
    Test how number of bootstrap replications affects the discrepancy.
    """
    print("\n" + "=" * 70)
    print("VARYING BOOTSTRAP REPLICATIONS")
    print("=" * 70 + "\n")

    bootstrap_grid = [50, 100, 200, 500, 1000]

    results = []

    for n_bootstraps in bootstrap_grid:
        print(f"Testing B={n_bootstraps}...", end=" ")

        # Multiple replications
        for rep in range(10):
            seed = 123 + rep + n_bootstraps
            try:
                res = run_single_comparison(
                    n_clusters=50, cluster_size=10, n_bootstraps=n_bootstraps, seed=seed
                )
                res["replication"] = rep
                results.append(res)
            except Exception as e:
                print(f"ERROR: {e}")
                continue

        print("Done.")

    df = pd.DataFrame(results)

    summary = (
        df.groupby("n_bootstraps")
        .agg(
            {
                "statsmodels_se": "mean",
                "duckreg_se": ["mean", "std"],
                "se_ratio": ["mean", "std"],
                "se_diff": ["mean", "std"],
                "se_pct_diff": ["mean", "std"],
            }
        )
        .round(4)
    )

    print("\nSUMMARY STATISTICS:")
    print(summary)

    return df, summary


def create_visualizations(df_clusters, df_bootstraps):
    """
    Create diagnostic visualizations.
    """
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    fig.suptitle("Clustered SE Diagnostic Results", fontsize=16, fontweight="bold")

    # Plot 1: SE ratio vs number of clusters (for different cluster sizes)
    ax = axes[0, 0]
    summary = (
        df_clusters.groupby(["n_clusters", "cluster_size"])["se_ratio"].mean().unstack()
    )
    for col in summary.columns:
        ax.plot(summary.index, summary[col], marker="o", label=f"n={col}")
    ax.axhline(y=1, color="red", linestyle="--", alpha=0.5, label="Perfect agreement")
    ax.set_xlabel("Number of Clusters (G)")
    ax.set_ylabel("SE Ratio (duckreg/statsmodels)")
    ax.set_title("Effect of Number of Clusters")
    ax.legend()
    ax.grid(True, alpha=0.3)

    # Plot 2: SE ratio vs cluster size (for different numbers of clusters)
    ax = axes[0, 1]
    summary = (
        df_clusters.groupby(["cluster_size", "n_clusters"])["se_ratio"].mean().unstack()
    )
    for col in summary.columns[:4]:  # Plot first 4 to avoid clutter
        ax.plot(summary.index, summary[col], marker="s", label=f"G={col}")
    ax.axhline(y=1, color="red", linestyle="--", alpha=0.5, label="Perfect agreement")
    ax.set_xlabel("Cluster Size (n)")
    ax.set_ylabel("SE Ratio (duckreg/statsmodels)")
    ax.set_title("Effect of Cluster Size")
    ax.legend()
    ax.grid(True, alpha=0.3)

    # Plot 3: Convergence with bootstrap replications
    ax = axes[1, 0]
    summary = df_bootstraps.groupby("n_bootstraps").agg({"se_ratio": ["mean", "std"]})
    ax.errorbar(
        summary.index,
        summary["se_ratio"]["mean"],
        yerr=summary["se_ratio"]["std"],
        marker="o",
        capsize=5,
    )
    ax.axhline(y=1, color="red", linestyle="--", alpha=0.5)
    ax.set_xlabel("Number of Bootstrap Replications (B)")
    ax.set_ylabel("SE Ratio (duckreg/statsmodels)")
    ax.set_title("Convergence with Bootstrap Replications")
    ax.grid(True, alpha=0.3)

    # Plot 4: Percent difference distribution
    ax = axes[1, 1]
    ax.hist(
        df_clusters["se_pct_diff"],
        bins=30,
        alpha=0.7,
        color="steelblue",
        edgecolor="black",
    )
    ax.axvline(x=0, color="red", linestyle="--", linewidth=2, label="Zero difference")
    ax.axvline(
        x=df_clusters["se_pct_diff"].median(),
        color="orange",
        linestyle="--",
        linewidth=2,
        label=f'Median: {df_clusters["se_pct_diff"].median():.1f}%',
    )
    ax.set_xlabel("Percent Difference (%)")
    ax.set_ylabel("Frequency")
    ax.set_title("Distribution of SE Percent Differences")
    ax.legend()
    ax.grid(True, alpha=0.3, axis="y")

    plt.tight_layout()
    return fig


# %%


def main():
    """
    Run the full diagnostic analysis.
    """
    # Run experiments
    df_clusters, summary_clusters = vary_clusters_and_size()
    df_bootstraps, summary_bootstraps = vary_bootstrap_replications()

    # Create visualizations
    print("\nGenerating visualizations...")
    fig = create_visualizations(df_clusters, df_bootstraps)
    fig.savefig(
        "cluster_se_diagnostics.png",
        dpi=300,
        bbox_inches="tight",
    )

    # Save detailed results
    df_clusters.to_csv("results_varying_clusters.csv", index=False)
    df_bootstraps.to_csv("results_varying_bootstraps.csv", index=False)

    summary_clusters.to_csv("summary_varying_clusters.csv")
    summary_bootstraps.to_csv("summary_varying_bootstraps.csv")


# %%
if __name__ == "__main__":
    main()

[apoorvalal]

merged the fix to the main branch, thanks again for the valuable bugreport. Note to self - calculate coverage rates for variance estimators, esp when relying on computational tricks like this.