IIabm/imputation_robustness_check.ipynb

Missing-data Imputation Robustness Check

This notebook audits where regression-based imputation is applied and evaluates robustness from two angles:

1. Firm size distribution: all firms (with imputed values) vs firms with original observed revenue only.
2. Network structure robustness: simulated firm network from all firms vs from observed-only firms.

In [1]:

import numpy as np
import pandas as pd
import networkx as nx
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme(style='whitegrid')

In [2]:

# Data loading
firm_amended = pd.read_csv('Firm_amended.csv')
firm_original = pd.read_csv('Firm.csv')

# The notebook `AmendFirm_20230216.ipynb` uses OLS:
# Revenue_Log ~ Num_Employ_Log
# and fills missing Revenue_Log by predicted values from Num_Employ_Log.

observed_mask = firm_original['Revenue'].notna()
firm_observed_only = firm_amended.loc[observed_mask].copy()

print(f'All firms (with imputation): {len(firm_amended)}')
print(f'Observed-only firms: {len(firm_observed_only)}')
print(f'Imputed-count firms: {len(firm_amended) - len(firm_observed_only)}')

All firms (with imputation): 171
Observed-only firms: 111
Imputed-count firms: 60

In [3]:

# Size distribution comparison: Revenue_Log
size_all = firm_amended['Revenue_Log'].dropna()
size_obs = firm_observed_only['Revenue_Log'].dropna()

quantiles = [0.1, 0.25, 0.5, 0.75, 0.9]
summary = pd.DataFrame({
    'all_firms': size_all.quantile(quantiles),
    'observed_only': size_obs.quantile(quantiles)
})
summary.index = [f'Q{int(q*100)}' for q in quantiles]
summary['delta_all_minus_observed'] = summary['all_firms'] - summary['observed_only']
print(summary)

mw = stats.mannwhitneyu(size_all, size_obs, alternative='two-sided')
ks = stats.ks_2samp(size_all, size_obs)
print('\nMann-Whitney U:', mw)
print('Kolmogorov-Smirnov:', ks)

     all_firms  observed_only  delta_all_minus_observed
Q10  18.266132      19.364587                 -1.098455
Q25  19.126847      20.995426                 -1.868579
Q50  21.273558      22.654014                 -1.380456
Q75  24.023677      25.818061                 -1.794384
Q90  26.562162      26.958298                 -0.396136

Mann-Whitney U: MannwhitneyuResult(statistic=6881.5, pvalue=9.665367898714309e-05)
Kolmogorov-Smirnov: KstestResult(statistic=0.23502449818239293, pvalue=0.0009254559883309539, statistic_location=20.22851328801836, statistic_sign=1)

In [4]:

# ECDF plot for distribution-shape comparison
plt.figure(figsize=(8, 5))
sns.ecdfplot(size_all, label='All firms (with imputation)')
sns.ecdfplot(size_obs, label='Observed-only firms')
plt.xlabel('Revenue_Log')
plt.ylabel('ECDF')
plt.title('Revenue_Log Distribution Comparison')
plt.legend()
plt.tight_layout()
plt.show()

Direct Comparison Using count_with_gfirm.csv

This section extracts full experimental firm networks from g_firm in the large count_with_gfirm.csv file and compares them against non-imputed networks generated with the same logic in model.py.

To control runtime and memory, it uses the first N unique s_id networks.

In [5]:

# Prepare network-comparison dependencies
# 1) lightweight check for count_with_gfirm.csv
# 2) build_firm_graph() matching model.py logic for observed-only simulation

import os
import json

if os.path.exists('count_with_gfirm.csv'):
    probe = pd.read_csv('count_with_gfirm.csv', usecols=['s_id', 'g_firm'], nrows=1)
    print('count_with_gfirm.csv exists; probe rows =', len(probe))
    if len(probe) > 0 and pd.notna(probe.loc[0, 'g_firm']):
        G_probe = nx.adjacency_graph(json.loads(probe.loc[0, 'g_firm']))
        print('One experimental graph loaded: nodes =', G_probe.number_of_nodes(), 'edges =', G_probe.number_of_edges())
    else:
        print('g_firm is empty in probe row; verify extraction in analysis_firm_risk_component.ipynb.')
else:
    print('count_with_gfirm.csv not found. Please extract Sample.g_firm first.')

# Build BOM graph once (for simulated network generation)
bom = pd.read_csv('BomCateNet.csv', index_col=0).fillna(0)
G_bom = nx.from_pandas_adjacency(bom.T, create_using=nx.MultiDiGraph())

def build_firm_graph(firm_df, seed=0, netw_prf_n=2, prf_size=True):
    rng = np.random.default_rng(seed)
    F = firm_df.copy().fillna(0)
    F['Code'] = F['Code'].astype(str)

    attrs = F[['Code', 'Name', 'Type_Region', 'Revenue_Log']].copy()
    products = []
    for _, row in F.loc[:, '1':].iterrows():
        products.append(row[row == 1].index.to_list())
    attrs['Product_Code'] = products
    attrs.set_index('Code', inplace=True)

    G = nx.MultiDiGraph()
    G.add_nodes_from(F['Code'].tolist())
    nx.set_node_attributes(G, {code: attrs.loc[code].to_dict() for code in G.nodes()})

    for node in list(G.nodes()):
        pred_products = []
        for p in G.nodes[node].get('Product_Code', []):
            if p in G_bom:
                pred_products += list(G_bom.predecessors(p))
        pred_products = sorted(set(pred_products))

        for pred_p in pred_products:
            pred_firms = F['Code'][F[pred_p] == 1].to_list() if pred_p in F.columns else []
            n = min(netw_prf_n, len(pred_firms))
            if n <= 0:
                continue

            if prf_size:
                sizes = np.array([max(1e-9, float(G.nodes[f].get('Revenue_Log', 0))) for f in pred_firms], dtype=float)
                probs = sizes / sizes.sum() if sizes.sum() > 0 else None
                chosen = rng.choice(pred_firms, size=n, replace=False, p=probs)
            else:
                chosen = rng.choice(pred_firms, size=n, replace=False)

            G.add_edges_from([(str(f), node, {'Product': pred_p}) for f in chosen])

    for node in list(G.nodes()):
        if G.degree(node) == 0:
            for p in G.nodes[node].get('Product_Code', []):
                if p not in G_bom:
                    continue
                for succ_p in G_bom.successors(p):
                    succ_firms = F['Code'][F[succ_p] == 1].to_list() if succ_p in F.columns else []
                    n = min(netw_prf_n, len(succ_firms))
                    if n <= 0:
                        continue
                    if prf_size:
                        sizes = np.array([max(1e-9, float(G.nodes[f].get('Revenue_Log', 0))) for f in succ_firms], dtype=float)
                        probs = sizes / sizes.sum() if sizes.sum() > 0 else None
                        chosen = rng.choice(succ_firms, size=n, replace=False, p=probs)
                    else:
                        chosen = rng.choice(succ_firms, size=n, replace=False)
                    G.add_edges_from([(node, str(f), {'Product': p}) for f in chosen])

    return G

count_with_gfirm.csv exists; probe rows = 1
One experimental graph loaded: nodes = 171 edges = 818

In [6]:

import json
from scipy import stats
from tqdm.auto import tqdm

N = None  # None means use all available s_id; set an integer to cap for faster testing

# 1) Extract one g_firm per s_id from large file
seen = set()
exp_graph_pairs = []
for chunk in tqdm(
    pd.read_csv('count_with_gfirm.csv', usecols=['s_id', 'g_firm'], chunksize=200000),
    desc='Reading chunks from count_with_gfirm.csv',
    unit='chunk'
):
    chunk = chunk.dropna(subset=['s_id', 'g_firm'])
    for sid, g in zip(chunk['s_id'].tolist(), chunk['g_firm'].tolist()):
        sid = int(sid)
        if sid in seen:
            continue
        seen.add(sid)
        exp_graph_pairs.append((sid, g))
        if N is not None and len(exp_graph_pairs) >= N:
            break
    if N is not None and len(exp_graph_pairs) >= N:
        break

print('Extracted experimental networks:', len(exp_graph_pairs))

# 2) Build experimental metrics from full g_firm
# (computed with graph_metrics_with_closeness below)
def graph_metrics_with_closeness(G):
    Gw = nx.DiGraph()
    Gw.add_nodes_from(G.nodes())
    for u, v in G.edges():
        if Gw.has_edge(u, v):
            Gw[u][v]['weight'] += 1
        else:
            Gw.add_edge(u, v, weight=1)

    indeg = dict(G.in_degree())
    outdeg = dict(G.out_degree())
    pr = nx.pagerank(Gw, weight='weight') if Gw.number_of_edges() > 0 else {n: 0 for n in Gw.nodes()}
    bt = nx.betweenness_centrality(Gw, weight='weight') if Gw.number_of_edges() > 0 else {n: 0 for n in Gw.nodes()}
    cl = nx.closeness_centrality(Gw, distance='weight') if Gw.number_of_edges() > 0 else {n: 0 for n in Gw.nodes()}

    return pd.DataFrame({
        'id_firm': [int(str(n)) for n in G.nodes()],
        'firm_num_supplier_firms': [indeg[n] for n in G.nodes()],
        'firm_num_customer_firms': [outdeg[n] for n in G.nodes()],
        'firm_weighted_betweenness': [bt.get(n, 0) for n in G.nodes()],
        'firm_weighted_closeness': [cl.get(n, 0) for n in G.nodes()],
        'firm_pagerank': [pr.get(n, 0) for n in G.nodes()]
    })

exp_full_parts = []
for sid, g in tqdm(exp_graph_pairs, desc='Building metrics for experimental networks', unit='network'):
    G_exp = nx.adjacency_graph(json.loads(g))
    exp_full_parts.append(graph_metrics_with_closeness(G_exp).assign(s_id=sid))
exp_full = pd.concat(exp_full_parts, ignore_index=True)

# 3) Build non-imputed comparison set with same N
sim_parts = []
for seed in tqdm(range(len(exp_graph_pairs)), desc='Building metrics for observed-only simulated networks', unit='network'):
    G_sim = build_firm_graph(firm_observed_only, seed=seed, netw_prf_n=2, prf_size=True)
    sim_parts.append(graph_metrics_with_closeness(G_sim).assign(s_id=seed + 1))
sim_full = pd.concat(sim_parts, ignore_index=True)

# 4) Restrict both to observed-only firm ids
obs_ids = set(firm_observed_only['Code'].astype(int).tolist())
exp_full = exp_full[exp_full['id_firm'].isin(obs_ids)].copy()
sim_full = sim_full[sim_full['id_firm'].isin(obs_ids)].copy()

metrics_cmp = [
    'firm_num_supplier_firms',
    'firm_num_customer_firms',
    'firm_weighted_betweenness',
    'firm_weighted_closeness',
    'firm_pagerank'
]

rows = []
for m in metrics_cmp:
    e = exp_full[m].to_numpy()
    s = sim_full[m].to_numpy()

    # Significance tests
    mw_res = stats.mannwhitneyu(e, s, alternative='two-sided')
    mw_u = float(mw_res.statistic)
    mw_p = float(mw_res.pvalue)

    ks_res = stats.ks_2samp(e, s)
    ks_stat = float(ks_res.statistic)
    ks_p = float(ks_res.pvalue)

    # Effect sizes
    n_e = len(e)
    n_s = len(s)
    cliffs_delta = (2.0 * mw_u) / (n_e * n_s) - 1.0
    wasserstein = float(stats.wasserstein_distance(e, s))

    # Cohen's d (pooled SD)
    e_mean, s_mean = float(np.mean(e)), float(np.mean(s))
    e_std, s_std = float(np.std(e, ddof=1)), float(np.std(s, ddof=1))
    pooled_sd = np.sqrt(((n_e - 1) * e_std**2 + (n_s - 1) * s_std**2) / (n_e + n_s - 2))
    cohen_d = (e_mean - s_mean) / pooled_sd if pooled_sd > 0 else np.nan

    q = [0.1, 0.25, 0.5, 0.75, 0.9]
    eq = np.quantile(e, q)
    sq = np.quantile(s, q)
    rows.append({
        'metric': m,
        'exp_q10': eq[0], 'exp_q25': eq[1], 'exp_q50': eq[2], 'exp_q75': eq[3], 'exp_q90': eq[4],
        'sim_q10': sq[0], 'sim_q25': sq[1], 'sim_q50': sq[2], 'sim_q75': sq[3], 'sim_q90': sq[4],
        'mw_u': mw_u,
        'mw_pvalue': mw_p,
        'ks_statistic': ks_stat,
        'ks_pvalue': ks_p,
        'cliffs_delta': cliffs_delta,
        'wasserstein_distance': wasserstein,
        'cohen_d': cohen_d
    })

dist_compare = pd.DataFrame(rows)
display(dist_compare)

# 5) Node-level mean consistency and top-node overlap
exp_node_mean = exp_full.groupby('id_firm')[metrics_cmp].mean()
sim_node_mean = sim_full.groupby('id_firm')[metrics_cmp].mean()
common_nodes = sorted(set(exp_node_mean.index) & set(sim_node_mean.index))

node_corr = pd.DataFrame({
    'metric': metrics_cmp,
    'pearson_corr': [
        np.corrcoef(exp_node_mean.loc[common_nodes, m], sim_node_mean.loc[common_nodes, m])[0, 1]
        for m in metrics_cmp
    ]
})
display(node_corr)

for m in ['firm_num_supplier_firms', 'firm_num_customer_firms', 'firm_weighted_betweenness', 'firm_weighted_closeness', 'firm_pagerank']:
    top_exp = set(exp_node_mean.loc[common_nodes, m].sort_values(ascending=False).head(10).index)
    top_sim = set(sim_node_mean.loc[common_nodes, m].sort_values(ascending=False).head(10).index)
    j = len(top_exp & top_sim) / len(top_exp | top_sim)
    print(f'Top10 Jaccard ({m}): {j:.4f}')

c:\Users\ASUS\Documents\Project\IIabm\venv\lib\site-packages\tqdm\auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
Reading chunks from count_with_gfirm.csv: 1chunk [00:07,  7.41s/chunk]

Extracted experimental networks: 12881

Building metrics for experimental networks: 100%|██████████| 12881/12881 [03:17<00:00, 65.11network/s]
Building metrics for observed-only simulated networks: 100%|██████████| 12881/12881 [12:03<00:00, 17.81network/s]

	metric	exp_q10	exp_q25	exp_q50	exp_q75	exp_q90	sim_q10	sim_q25	sim_q50	sim_q75	sim_q90	mw_u	mw_pvalue	ks_statistic	ks_pvalue	cliffs_delta	wasserstein_distance	cohen_d
0	firm_num_supplier_firms	0.000000	0.000000	0.000000	6.000000	16.000000	0.000000	0.000000	0.000000	6.000000	14.000000	1.018947e+12	4.438885e-08	0.030905	0.0	-0.003134	0.580167	0.011394
1	firm_num_customer_firms	1.000000	1.000000	2.000000	3.000000	7.000000	1.000000	1.000000	2.000000	4.000000	11.000000	8.561743e+11	0.000000e+00	0.140991	0.0	-0.162380	1.595933	-0.226000
2	firm_weighted_betweenness	0.000000	0.000000	0.000000	0.000000	0.001357	0.000000	0.000000	0.000000	0.000128	0.002363	9.882519e+11	0.000000e+00	0.063961	0.0	-0.033165	0.000443	-0.145224
3	firm_weighted_closeness	0.000000	0.000000	0.000000	0.029412	0.079619	0.000000	0.000000	0.000000	0.024242	0.111364	1.024913e+12	2.397071e-06	0.034881	0.0	0.002702	0.005074	-0.045781
4	firm_pagerank	0.001826	0.001927	0.002224	0.005013	0.014727	0.002991	0.003069	0.003179	0.005520	0.016487	5.323516e+11	0.000000e+00	0.670247	0.0	-0.479185	0.001240	-0.064583

	metric	pearson_corr
0	firm_num_supplier_firms	0.995273
1	firm_num_customer_firms	0.972093
2	firm_weighted_betweenness	0.983254
3	firm_weighted_closeness	0.991394
4	firm_pagerank	0.994468

Top10 Jaccard (firm_num_supplier_firms): 1.0000
Top10 Jaccard (firm_num_customer_firms): 0.5385
Top10 Jaccard (firm_weighted_betweenness): 0.6667
Top10 Jaccard (firm_weighted_closeness): 1.0000
Top10 Jaccard (firm_pagerank): 1.0000

In [13]:

# 五个网络指标的 ECDF 对比图（多子图，x 轴对数坐标）
import matplotlib
from matplotlib.ticker import LogLocator, FuncFormatter, NullFormatter

# 设置中文字体（含 fallback），并处理负号/数学文本显示
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = ['SimSun', 'Microsoft YaHei', 'DejaVu Sans']
matplotlib.rcParams['axes.unicode_minus'] = False
matplotlib.rcParams['mathtext.fontset'] = 'dejavusans'
matplotlib.rcParams['mathtext.default'] = 'regular'

metric_name_map = {
    'firm_num_supplier_firms': '入度中心性',
    'firm_num_customer_firms': '出度中心性',
    'firm_weighted_betweenness': '加权介数中心性',
    'firm_weighted_closeness': '加权接近中心性',
    'firm_pagerank': 'PageRank值'
}

plot_metrics = list(metric_name_map.keys())
epsilon = 1e-6

fig, axes = plt.subplots(2, 3, figsize=(12, 8), dpi=300, sharey=True)
axes = axes.flatten()

for i, metric in enumerate(plot_metrics):
    ax = axes[i]

    exp_vals = exp_full[metric].to_numpy() + epsilon
    sim_vals = sim_full[metric].to_numpy() + epsilon

    sns.ecdfplot(
        exp_vals,
        ax=ax,
        label='实验网络（含插值）',
        linewidth=1.8,
        alpha=0.85
    )
    sns.ecdfplot(
        sim_vals,
        ax=ax,
        label='对照网络（仅原始样本）',
        linewidth=1.8,
        alpha=0.85
    )

    ax.set_xscale('log')

    # 自定义对数坐标刻度格式，避免 U+2212 字形问题
    ax.xaxis.set_major_locator(LogLocator(base=10.0))
    ax.xaxis.set_minor_formatter(NullFormatter())
    ax.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: f"1e{int(np.log10(x))}" if x > 0 else ""))

    # 子图边框统一为黑色
    for spine in ax.spines.values():
        spine.set_color('black')
        spine.set_linewidth(1.0)

    ax.set_title(metric_name_map[metric], fontsize=12)
    ax.set_xlabel('指标值（对数坐标）', fontsize=11)
    if i % 3 == 0:
        ax.set_ylabel('经验累积分布函数', fontsize=11)
    else:
        ax.set_ylabel('')
    ax.grid(True, alpha=0.3)

# 每个子图显示各自图例
for i in range(len(plot_metrics)):
    axes[i].legend(loc='best', fontsize=9, frameon=False)

# 隐藏第 6 个空子图
axes[-1].axis('off')

plt.tight_layout()
plt.show()

1) 规模分布层面：插值样本与原始样本存在统计差异，但方向清晰且可解释。
Revenue_Log 在含插值样本中的 Q10/Q25/Q50/Q75 仍低于仅原始观测样本（约下降 0.4-1.9），与补充样本中中小企业占比更高这一设定一致。该差异提示：插值主要改变了样本的边际规模分布，而非引入无方向的随机噪声。

2) 网络指标分布层面：全量样本下“统计显著”与“实质效应”仍需区分。
在直接比较（count_with_gfirm.csv 全部 s_id vs 无插值生成网络）中，5 个指标的 MW/KS 检验均显著，但效应量显示明显分化：

• firm_num_supplier_firms 与 firm_weighted_closeness 的 Cliff's delta 接近 0（约 -0.003 与 0.003），firm_weighted_betweenness 也为小效应（约 -0.033）；
• firm_num_customer_firms 为小到中等效应（Cliff's delta 约 -0.162）；
• firm_pagerank 仍为相对最强差异（Cliff's delta 约 -0.479，KS statistic 约 0.670）。 因此，更稳妥的解读是：多数指标体现为“超大样本下显著但效应有限”，而 pagerank 分布存在较明显偏移。

3) 节点相对排序层面：核心结构关系仍高度稳健。
尽管总体分布存在差异，节点均值相关性仍维持在高水平（Pearson: 0.972-0.995；pagerank=0.994）。关键节点识别方面，Top10 Jaccard 分别为：supplier=1.0000、customer=0.5385、betweenness=0.6667、closeness=1.0000、pagerank=1.0000。说明核心节点集合在主要指标上总体保持稳定，尤其是 supplier、closeness 与 pagerank 的关键节点识别几乎完全一致。