# -*- coding: utf-8 -*- import argparse import math import re import warnings from collections import Counter, defaultdict from itertools import combinations from pathlib import Path import numpy as np import pandas as pd try: import networkx as nx except Exception: nx = None try: import statsmodels.api as sm except Exception: sm = None try: from sklearn.linear_model import LogisticRegression except Exception: LogisticRegression = None YEAR_START = 2016 YEAR_END = 2025 PERIODS = [ (1, "2016-2017", 2016, 2017), (2, "2018-2019", 2018, 2019), (3, "2020-2021", 2020, 2021), (4, "2022-2023", 2022, 2023), (5, "2024-2025", 2024, 2025), ] DEFAULT_INPUT_FILE = r"C:/Users/ACER/Desktop/无人机专利数据/tech1_network_output/patent_tech_long_raw.csv" DEFAULT_OUTPUT_DIR = r"C:/Users/ACER/Desktop/三项筛选条件控制变量鲁棒性检验_技术邻近性优化_无弱方向差异" MIN_TECHS_PER_PATENT = 2 MAX_TECHS_PER_PATENT = 5 BASE_COUNT_THRESHOLD = 40 BASE_PERIOD_THRESHOLD = 4 BASE_EDGE_THRESHOLD = 3 RAW_EDGE_THRESHOLD = 1 CORE_VARS = [ "gwesp_proxy", "evcent_level", "evcent_absdiff", "tech_size_level", "tech_size_absdiff", "proximity", "history_log", ] CONTROL_VARS = [ "short_lifecycle_exposure", "tie_strength_log", ] VAR_CN = { "gwesp_proxy": "局部闭合项", "evcent_level": "中心性水平项", "evcent_absdiff": "中心性差异项", "tech_size_level": "技术规模水平项", "tech_size_absdiff": "技术规模差异项", "proximity": "技术邻近性项", "history_log": "历史关联项", "short_lifecycle_exposure": "短生命周期暴露控制项", "tie_strength_log": "当前关系强度控制项", } GROUP_SPECS = [ { "scheme": "G0_raw_no_three_filters", "group_cn": "原始未筛选组", "min_total_count": None, "min_active_periods": None, "edge_threshold": 1, "purpose": "不设置技术总出现频次阈值、不设置最少出现时期数阈值，边关系按共现次数不低于1次判定，用作原始边界对照。" }, { "scheme": "G1_core_baseline_40_4_edge3", "group_cn": "核心基准组", "min_total_count": 40, "min_active_periods": 4, "edge_threshold": 3, "purpose": "主模型设定：技术总出现频次不低于40次，活跃时期数不少于4个时期，边关系按共现次数不低于3次判定。" }, { "scheme": "G2_no_frequency_threshold", "group_cn": "频次条件对照组", "min_total_count": None, "min_active_periods": 4, "edge_threshold": 3, "purpose": "仅取消技术总出现频次阈值，保持活跃时期数不少于4个时期和边阈值不低于3次，用于检验频次筛选条件。" }, { "scheme": "G3_no_period_threshold", "group_cn": "时期条件对照组", "min_total_count": 40, "min_active_periods": None, "edge_threshold": 3, "purpose": "仅取消最少出现时期数阈值，保持总出现频次不低于40次和边阈值不低于3次，用于检验跨期持续性筛选条件。" }, { "scheme": "G4_edge1_threshold", "group_cn": "边阈值对照组", "min_total_count": 40, "min_active_periods": 4, "edge_threshold": 1, "purpose": "保持技术总出现频次不低于40次、活跃时期数不少于4个时期，仅将边阈值从共现不低于3次放宽为不低于1次，用于检验边关系二值化阈值。" }, ] COMPARISONS = [ { "comparison": "C1_频次阈值控制检验", "left_scheme": "G1_core_baseline_40_4_edge3", "right_scheme": "G2_no_frequency_threshold", "meaning": "控制活跃时期数和边阈值不变，仅取消技术总出现频次不低于40次这一条件。" }, { "comparison": "C2_时期阈值控制检验", "left_scheme": "G1_core_baseline_40_4_edge3", "right_scheme": "G3_no_period_threshold", "meaning": "控制总出现频次和边阈值不变，仅取消活跃时期数不少于4个时期这一条件。" }, { "comparison": "C3_边阈值控制检验", "left_scheme": "G1_core_baseline_40_4_edge3", "right_scheme": "G4_edge1_threshold", "meaning": "控制技术节点筛选条件不变，仅将边关系判定从共现不低于3次放宽为不低于1次。" }, { "comparison": "C4_原始边界对照", "left_scheme": "G1_core_baseline_40_4_edge3", "right_scheme": "G0_raw_no_three_filters", "meaning": "同时取消三项筛选条件，用于观察原始边界组与主模型设定之间的差异，不作为单一条件敏感性判断。" }, ] def extract_year(x): if pd.isna(x): return None m = re.search(r"(20\d{2})", str(x)) return int(m.group(1)) if m else None def period_id_from_year(year): for pid, _, start, end in PERIODS: if start <= int(year) <= end: return pid return None def norm_cdf(x): if pd.isna(x): return np.nan return 0.5 * (1.0 + math.erf(float(x) / math.sqrt(2.0))) def norm_log_sf(z): if pd.isna(z): return np.nan z = float(z) if z < 26.0: p = 0.5 * math.erfc(z / math.sqrt(2.0)) if p > 0.0: return math.log(p) if z > 0: inv2 = 1.0 / (z * z) correction = math.log1p(-inv2 + 3.0 * inv2 * inv2) return -0.5 * z * z - math.log(z) - 0.5 * math.log(2.0 * math.pi) + correction p = 0.5 * math.erfc(z / math.sqrt(2.0)) return math.log(min(max(p, 1e-300), 1.0)) def norm_log_cdf(z): if pd.isna(z): return np.nan return norm_log_sf(-float(z)) def logp_to_sci_string(logp, digits=3): if pd.isna(logp): return "" logp = float(logp) if logp == -math.inf: return "0" if logp >= 0: return "1" log10p = logp / math.log(10.0) exponent = math.floor(log10p) mantissa = 10.0 ** (log10p - exponent) rounded = round(mantissa, digits - 1) if rounded >= 10.0: rounded = 1.0 exponent += 1 return ("%." + str(digits - 1) + "fe%d") % (rounded, exponent) def tost_p_from_z(z_lower, z_upper): log_p_lower = norm_log_sf(z_lower) log_p_upper = norm_log_cdf(z_upper) log_p = max(log_p_lower, log_p_upper) tiny_log = math.log(np.finfo(float).tiny) p_value = math.exp(log_p) if pd.notna(log_p) and log_p > tiny_log else 0.0 return p_value, log_p, logp_to_sci_string(log_p) def fmt_num(x, digits=3): if pd.isna(x): return "" return ("%%.%df" % digits) % float(x) def fmt_p(x): if pd.isna(x): return "" if float(x) < 0.001: return "<0.001" return "%.3f" % float(x) def normalize_columns(df): cols = list(df.columns) def pick(names): for name in names: if name in cols: return name lower = {str(c).lower(): c for c in cols} for name in names: if str(name).lower() in lower: return lower[str(name).lower()] return None patent_col = pick(["patent_id", "patent", "申请号", "公开（公告）号", "序号", "id"]) year_col = pick(["year", "申请年", "年份"]) date_col = pick(["申请日", "申请日期", "date", "application_date"]) tech_col = pick(["tech_code", "technology", "IPC", "ipc", "主分类号", "分类号"]) return patent_col, year_col, date_col, tech_col def extract_main_groups(ipc_raw): if pd.isna(ipc_raw): return [] s = str(ipc_raw).upper().strip() parts = re.split(r"[;,；，\n\r\t]+", s) out = [] seen = set() for p in parts: p = re.sub(r"\s+", "", p.strip()) if not p: continue m = re.match(r"^([A-HY]\d{2}[A-Z]\d+)", p) if m: code = m.group(1) if code not in seen: seen.add(code) out.append(code) return out def load_patent_tech_long(input_path, apply_patent_size_filter=True): input_path = Path(input_path) if not input_path.exists(): raise FileNotFoundError(str(input_path)) if input_path.suffix.lower() in [".xlsx", ".xls"]: df = pd.read_excel(input_path) else: df = pd.read_csv(input_path, encoding="utf-8-sig") patent_col, year_col, date_col, tech_col = normalize_columns(df) if patent_col is None or tech_col is None: raise ValueError("Cannot identify patent id or technology column. Current columns: %s" % list(df.columns)) work = df.copy() if year_col is None: if date_col is None: raise ValueError("Cannot identify year or application date column. Current columns: %s" % list(df.columns)) work["year"] = work[date_col].apply(extract_year) else: work["year"] = work[year_col].apply(extract_year) work = work[(work["year"] >= YEAR_START) & (work["year"] <= YEAR_END)].copy() work = work.dropna(subset=[patent_col, tech_col, "year"]).copy() work[patent_col] = work[patent_col].astype(str).str.strip() rows = [] if str(tech_col).lower() in ["tech_code", "technology"] or "tech" in str(tech_col).lower(): for _, row in work.iterrows(): tech = str(row[tech_col]).strip().upper() if tech: rows.append((row[patent_col], int(row["year"]), tech)) else: for _, row in work.iterrows(): for tech in extract_main_groups(row[tech_col]): rows.append((row[patent_col], int(row["year"]), tech)) long_df = pd.DataFrame(rows, columns=["patent_id", "year", "tech_code"]) long_df["period_id"] = long_df["year"].apply(period_id_from_year) long_df = long_df.dropna(subset=["period_id"]).copy() long_df["period_id"] = long_df["period_id"].astype(int) long_df = long_df.drop_duplicates(["patent_id", "period_id", "tech_code"]).copy() before_patents = long_df["patent_id"].nunique() if apply_patent_size_filter: patent_size = ( long_df.groupby(["patent_id", "period_id"])["tech_code"] .nunique() .reset_index(name="n_tech") ) keep = patent_size[ (patent_size["n_tech"] >= MIN_TECHS_PER_PATENT) & (patent_size["n_tech"] <= MAX_TECHS_PER_PATENT) ][["patent_id", "period_id"]] long_df = long_df.merge(keep, on=["patent_id", "period_id"], how="inner") after_patents = long_df["patent_id"].nunique() print("Loaded patent-tech table: rows=%s, patents_before_size_filter=%s, patents_after_size_filter=%s" % (len(long_df), before_patents, after_patents)) print("Patent-size cleaning applied:", apply_patent_size_filter) return long_df def tech_stats(long_df): return ( long_df.groupby("tech_code") .agg( total_count=("patent_id", "nunique"), n_years=("year", "nunique"), active_periods=("period_id", "nunique"), ) .reset_index() ) def build_roster(long_df, min_total_count=None, min_active_periods=None): stat = tech_stats(long_df) df = stat.copy() if min_total_count is not None: df = df[df["total_count"] >= float(min_total_count)] if min_active_periods is not None: df = df[df["active_periods"] >= float(min_active_periods)] return sorted(df["tech_code"].astype(str).tolist()), stat def build_weights_by_period(long_df, roster): roster_set = set(roster) weights = {pid: Counter() for pid, _, _, _ in PERIODS} active = {pid: set() for pid, _, _, _ in PERIODS} grouped = long_df.groupby(["period_id", "patent_id"])["tech_code"].apply(lambda x: sorted(set(x) & roster_set)) for (pid, _patent), techs in grouped.items(): pid = int(pid) if not techs: continue active[pid].update(techs) if len(techs) >= 2: for i, j in combinations(techs, 2): weights[pid][tuple(sorted((i, j)))] += 1 return weights, active def build_edges_by_period(weights_by_period, threshold): return { pid: {pair for pair, w in weights_by_period[pid].items() if w >= threshold} for pid, _, _, _ in PERIODS } def build_cumulative_weights(weights_by_period): cumulative = {} running = Counter() for pid, _, _, _ in PERIODS: running.update(weights_by_period[pid]) cumulative[pid] = Counter(running) return cumulative def build_proximity(long_df, roster, shrinkage=5.0): roster_set = set(roster) tech_patent_count = Counter() pair_patent_count = Counter() grouped = long_df.groupby("patent_id")["tech_code"].apply(lambda x: sorted(set(x) & roster_set)) for _patent, techs in grouped.items(): for t in techs: tech_patent_count[t] += 1 if len(techs) >= 2: for i, j in combinations(techs, 2): pair_patent_count[tuple(sorted((i, j)))] += 1 prox = {} for pair, cnt in pair_patent_count.items(): i, j = pair denom = math.sqrt(tech_patent_count[i] * tech_patent_count[j]) + shrinkage prox[pair] = 0.0 if denom <= 0 else cnt / denom return prox def build_second_order_stable_proximity_profiles(weights_by_period, roster, stable_threshold=3): """ Technical proximity is defined as similarity in stable technological environments, rather than direct pair co-occurrence. This avoids making proximity mechanically overlap with edge1 weak-tie outcomes. """ out = {} running = Counter() roster_set = set(roster) for pid, _, _, _ in PERIODS: strong_edges = { pair: w for pair, w in weights_by_period[pid].items() if w >= stable_threshold and pair[0] in roster_set and pair[1] in roster_set } running.update(strong_edges) profiles = {t: {} for t in roster} for (a, b), w in running.items(): if a not in profiles or b not in profiles: continue val = float(w) profiles[a][b] = profiles[a].get(b, 0.0) + val profiles[b][a] = profiles[b].get(a, 0.0) + val norms = {} for t, prof in profiles.items(): norms[t] = math.sqrt(sum(v * v for v in prof.values())) out[pid] = {"profiles": profiles, "norms": norms} return out def get_second_order_stable_proximity(profile_obj, pid, i, j): pack = profile_obj.get(pid, {}) profiles = pack.get("profiles", {}) norms = pack.get("norms", {}) pi = profiles.get(i, {}) pj = profiles.get(j, {}) ni = norms.get(i, 0.0) nj = norms.get(j, 0.0) if ni <= 0.0 or nj <= 0.0: return 0.0 if len(pi) <= len(pj): dot = sum(pi.get(k, 0.0) * pj.get(k, 0.0) for k in pi.keys() if k != i and k != j) else: dot = sum(pi.get(k, 0.0) * pj.get(k, 0.0) for k in pj.keys() if k != i and k != j) denom = ni * nj return float(dot / denom) if denom > 0.0 else 0.0 def compute_centrality_by_period(edges_by_period, active_by_period, roster): out = {} for pid, _, _, _ in PERIODS: values = {t: 0.0 for t in roster} active = active_by_period[pid] edges = edges_by_period[pid] if not active: out[pid] = values continue if nx is None: deg = Counter() for i, j in edges: deg[i] += 1 deg[j] += 1 denom = max(1, len(active) - 1) for t in active: values[t] = deg[t] / denom out[pid] = values continue g = nx.Graph() g.add_nodes_from(active) g.add_edges_from(edges) try: c = nx.eigenvector_centrality(g, max_iter=1000, tol=1e-8) except Exception: c = nx.degree_centrality(g) for t, v in c.items(): values[t] = float(v) out[pid] = values return out def compute_size_by_period(long_df, roster): roster_set = set(roster) out = {pid: {t: 0.0 for t in roster} for pid, _, _, _ in PERIODS} counts = long_df[long_df["tech_code"].isin(roster_set)].groupby(["period_id", "tech_code"]).size() for (pid, tech), cnt in counts.items(): out[int(pid)][tech] = float(cnt) return out def degree_maps_by_period(edges_by_period, roster): out = {} for pid, _, _, _ in PERIODS: deg = {t: 0 for t in roster} for i, j in edges_by_period[pid]: deg[i] = deg.get(i, 0) + 1 deg[j] = deg.get(j, 0) + 1 out[pid] = deg return out def active_periods_from_maps(active_by_period, roster): out = {t: 0 for t in roster} for pid in active_by_period: for t in active_by_period[pid]: out[t] = out.get(t, 0) + 1 return out def network_summary_for_scheme(scheme, group_cn, roster, active_by_period, edges_by_period, edge_threshold): rows = [] n = len(roster) all_pairs = n * (n - 1) / 2 if n >= 2 else 0 for pid, label, _, _ in PERIODS: active_nodes = len(active_by_period[pid]) active_dyads = active_nodes * (active_nodes - 1) / 2 if active_nodes >= 2 else 0 edges = len(edges_by_period[pid]) density = edges / active_dyads if active_dyads else 0.0 structural_zeros = int(all_pairs - active_dyads) rows.append({ "scheme": scheme, "group_cn": group_cn, "period": label, "roster_nodes": n, "edge_threshold": edge_threshold, "active_nodes": active_nodes, "edges": edges, "density": density, "structural_zeros": structural_zeros, }) return pd.DataFrame(rows) def build_panel(scheme, group_cn, roster, long_df, edge_threshold, proximity_shrinkage=5.0, proximity_mode="second_order_stable"): weights, active = build_weights_by_period(long_df, roster) edges = build_edges_by_period(weights, edge_threshold) cumulative = build_cumulative_weights(weights) centrality = compute_centrality_by_period(edges, active, roster) size = compute_size_by_period(long_df, roster) if proximity_mode == "direct": proximity = build_proximity(long_df, roster, proximity_shrinkage) stable_proximity_profiles = None else: proximity = {} stable_proximity_profiles = build_second_order_stable_proximity_profiles( weights, roster, stable_threshold=BASE_EDGE_THRESHOLD, ) degree = degree_maps_by_period(edges, roster) active_periods = active_periods_from_maps(active, roster) rows = [] for pid in range(1, 5): next_pid = pid + 1 active_t = active[pid] active_next = active[next_pid] edges_t = edges[pid] edges_next = edges[next_pid] neighbor = defaultdict(set) for a, b in edges_t: neighbor[a].add(b) neighbor[b].add(a) for i, j in combinations(roster, 2): pair = tuple(sorted((i, j))) if i not in active_t or j not in active_t or i not in active_next or j not in active_next: continue edge_t = pair in edges_t edge_next = pair in edges_next common_neighbors = len(neighbor[i] & neighbor[j]) deg_i = degree[pid].get(i, 0) deg_j = degree[pid].get(j, 0) closure = common_neighbors / math.sqrt(max(deg_i, 1) * max(deg_j, 1)) ci = float(centrality[pid].get(i, 0.0)) cj = float(centrality[pid].get(j, 0.0)) si = math.log1p(float(size[pid].get(i, 0.0))) sj = math.log1p(float(size[pid].get(j, 0.0))) current_w = weights[pid].get(pair, 0.0) history_weight = max(cumulative[pid].get(pair, 0.0) - current_w, 0.0) if edge_t else cumulative[pid].get(pair, 0.0) short_lifecycle_exposure = ( (1.0 - active_periods.get(i, 0) / len(PERIODS)) + (1.0 - active_periods.get(j, 0) / len(PERIODS)) ) / 2.0 base = { "scheme": scheme, "group_cn": group_cn, "transition": "%s_to_%s" % (pid, next_pid), "from_period_id": pid, "to_period_id": next_pid, "tech_i": i, "tech_j": j, "gwesp_proxy": closure, "evcent_level": ci + cj, "evcent_absdiff": abs(ci - cj), "tech_size_level": si + sj, "tech_size_absdiff": abs(si - sj), "proximity": proximity.get(pair, 0.0) if proximity_mode == "direct" else get_second_order_stable_proximity(stable_proximity_profiles, pid, i, j), "history_log": math.log1p(history_weight), "tie_strength_log": math.log1p(current_w) if edge_t else 0.0, "short_lifecycle_exposure": short_lifecycle_exposure, } if not edge_t: row = dict(base) row["process"] = "关系形成" row["y"] = 1 if edge_next else 0 rows.append(row) if edge_t: row = dict(base) row["process"] = "关系维持" row["y"] = 1 if edge_next else 0 rows.append(row) network = network_summary_for_scheme(scheme, group_cn, roster, active, edges, edge_threshold) return network, pd.DataFrame(rows) def build_scale_params(panel_df): params = {} for process, g in panel_df.groupby("process"): params[process] = {} for var in CORE_VARS + CONTROL_VARS: if var not in g.columns: continue x = g[var].astype(float) std = float(x.std(ddof=0)) if len(x) else 1.0 if pd.isna(std) or std < 1e-12: std = 1.0 params[process][var] = (float(x.mean()) if len(x) else 0.0, std) return params def standardize(df, process, params, candidate_vars): out = df.copy() used = [] for var in candidate_vars: if var not in out.columns: continue x = out[var].astype(float) if x.std(ddof=0) < 1e-12: continue if params is not None and process in params and var in params[process]: mean, std = params[process][var] else: mean = float(x.mean()) std = float(x.std(ddof=0)) if pd.isna(std) or std < 1e-12: std = 1.0 out[var] = (x - mean) / std used.append(var) return out, used def remove_collinear_predictors(X, used_vars): X_work = X.copy() used_work = list(used_vars) removed = [] try: rank = np.linalg.matrix_rank(X_work.to_numpy(dtype=float)) except Exception: return X_work, used_work, removed if rank >= X_work.shape[1]: return X_work, used_work, removed drop_order = [v for v in used_work if v in CONTROL_VARS] + [v for v in used_work if v in CORE_VARS] while drop_order: try: rank = np.linalg.matrix_rank(X_work.to_numpy(dtype=float)) except Exception: break if rank >= X_work.shape[1]: break drop_var = drop_order.pop(0) if drop_var not in X_work.columns: continue X_work = X_work.drop(columns=[drop_var]) if drop_var in used_work: used_work.remove(drop_var) removed.append(drop_var) return X_work, used_work, removed def fit_model(panel_df, scheme, group_cn, scale_params=None): if sm is None and LogisticRegression is None: raise ImportError("statsmodels or sklearn is required.") rows = [] candidate_vars = CORE_VARS + CONTROL_VARS for process, g in panel_df.groupby("process"): needed = ["y"] + [v for v in candidate_vars if v in g.columns] g = g.dropna(subset=needed).copy() if len(g) < 30 or g["y"].nunique() < 2: print("Skipped model: scheme=%s process=%s n=%s events=%s" % (scheme, process, len(g), int(g["y"].sum()) if len(g) else 0)) continue work, used = standardize(g, process, scale_params, candidate_vars) X = sm.add_constant(work[used], has_constant="add") if sm is not None else pd.concat( [pd.Series(1.0, index=work.index, name="const"), work[used]], axis=1 ) X, used, removed = remove_collinear_predictors(X, used) if removed: print("Collinearity guard: scheme=%s process=%s removed=%s" % (scheme, process, ",".join(removed))) y = work["y"].astype(int) params = bse = pvalues = pred_prob = None status = "" if sm is not None: try: with warnings.catch_warnings(record=True) as caught: warnings.simplefilter("always", RuntimeWarning) result = sm.GLM(y, X, family=sm.families.Binomial()).fit(maxiter=1000, cov_type="HC1") bad_warning = any( ("divide by zero" in str(w.message)) or ("invalid value" in str(w.message)) or ("overflow" in str(w.message)) for w in caught ) bad_bse = bool(pd.Series(result.bse).replace([np.inf, -np.inf], np.nan).isna().any()) if bad_warning or bad_bse: raise RuntimeError("glm_numerical_warning_or_invalid_bse") params = result.params bse = result.bse pvalues = result.pvalues pred_prob = np.asarray(result.predict(X), dtype=float) status = "glm_binomial_hc1" except Exception as e1: try: ridge = sm.GLM(y, X, family=sm.families.Binomial()).fit_regularized( alpha=1e-4, L1_wt=0.0, maxiter=5000, ) params = ridge.params bse = pd.Series(np.nan, index=params.index) pvalues = pd.Series(np.nan, index=params.index) pred_prob = np.asarray(ridge.predict(X), dtype=float) status = "ridge_glm_binomial_fallback" except Exception: params = bse = pvalues = pred_prob = None status = "glm_and_ridge_failed:%s" % e1 if params is None: if LogisticRegression is None: print("Model failed: scheme=%s process=%s error=%s" % (scheme, process, status)) continue try: x_cols = [c for c in X.columns if c != "const"] X_np = X[x_cols].to_numpy(dtype=float) clf = LogisticRegression( penalty="l2", C=0.5, solver="lbfgs", max_iter=5000, random_state=20260529, ) clf.fit(X_np, y.to_numpy()) pred_prob = clf.predict_proba(X_np)[:, 1] coef_values = np.concatenate(([float(clf.intercept_[0])], clf.coef_[0].astype(float))) coef_index = ["const"] + x_cols params = pd.Series(coef_values, index=coef_index) bse = pd.Series(np.nan, index=coef_index) pvalues = pd.Series(np.nan, index=coef_index) status = "sklearn_logistic_l2_fallback" except Exception as e2: print("Model failed: scheme=%s process=%s error=%s" % (scheme, process, e2)) continue pred_prob = np.clip(np.asarray(pred_prob, dtype=float), 1e-8, 1 - 1e-8) gradient = pred_prob * (1.0 - pred_prob) for var in CORE_VARS + CONTROL_VARS: if var not in used: continue coef = params.get(var, np.nan) se = bse.get(var, np.nan) if hasattr(bse, "get") else np.nan pval = pvalues.get(var, np.nan) if hasattr(pvalues, "get") else np.nan if pd.notna(coef): indiv_me = float(coef) * gradient ame = float(np.mean(indiv_me)) ame_se = float(np.std(indiv_me, ddof=1) / math.sqrt(len(indiv_me))) if len(indiv_me) > 1 else np.nan else: ame = np.nan ame_se = np.nan rows.append({ "scheme": scheme, "group_cn": group_cn, "process": process, "variable": var, "variable_cn": VAR_CN.get(var, var), "term_type": "core" if var in CORE_VARS else "control", "coef": float(coef) if pd.notna(coef) else np.nan, "std_error": float(se) if pd.notna(se) else np.nan, "p_value": float(pval) if pd.notna(pval) else np.nan, "ame": ame, "ame_se": ame_se, "n_obs": int(len(work)), "n_events": int(y.sum()), "event_rate": float(y.mean()), "model_status": status, }) return pd.DataFrame(rows) def compare_models(coef_df, left_scheme, right_scheme, comparison, comparison_meaning, ame_margin=0.10, alpha=0.05): rows = [] core = coef_df[coef_df["term_type"] == "core"].copy() for (process, variable, variable_cn), g in core.groupby(["process", "variable", "variable_cn"]): left = g[g["scheme"] == left_scheme] right = g[g["scheme"] == right_scheme] if left.empty or right.empty: continue l = left.iloc[0] r = right.iloc[0] lcoef = float(l["coef"]) rcoef = float(r["coef"]) lp = float(l["p_value"]) if pd.notna(l["p_value"]) else np.nan rp = float(r["p_value"]) if pd.notna(r["p_value"]) else np.nan lame = float(l["ame"]) if pd.notna(l["ame"]) else np.nan rame = float(r["ame"]) if pd.notna(r["ame"]) else np.nan lame_se = float(l["ame_se"]) if pd.notna(l["ame_se"]) else np.nan rame_se = float(r["ame_se"]) if pd.notna(r["ame_se"]) else np.nan ame_delta = rame - lame if pd.notna(lame) and pd.notna(rame) else np.nan ame_se_delta = math.sqrt(lame_se ** 2 + rame_se ** 2) if pd.notna(lame_se) and pd.notna(rame_se) else np.nan if pd.notna(ame_delta) and pd.notna(ame_se_delta) and ame_se_delta > 1e-12: z_lower = (ame_delta + ame_margin) / ame_se_delta z_upper = (ame_delta - ame_margin) / ame_se_delta ame_tost_p, ame_tost_log_p, ame_tost_p_exact = tost_p_from_z(z_lower, z_upper) ame_ci90_low = ame_delta - 1.6448536269514722 * ame_se_delta ame_ci90_high = ame_delta + 1.6448536269514722 * ame_se_delta else: ame_tost_p = ame_tost_log_p = ame_ci90_low = ame_ci90_high = np.nan ame_tost_p_exact = "" evidence_left = pd.notna(lp) and lp < 0.10 evidence_right = pd.notna(rp) and rp < 0.10 substantive_left = abs(lcoef) >= 0.20 and evidence_left substantive_right = abs(rcoef) >= 0.20 and evidence_right reversal = bool(substantive_left and substantive_right and abs(rcoef - lcoef) > 0.50 and np.sign(lcoef) != np.sign(rcoef)) if pd.notna(ame_tost_p) and ame_tost_p < alpha: ame_judgement = "AME-TOST通过" elif pd.notna(ame_delta) and abs(ame_delta) <= ame_margin: ame_judgement = "AME点估计等效" else: ame_judgement = "AME超界或不可判定" if reversal: direction_judgement = "实质性方向反转" else: direction_judgement = "无实质方向反转" rows.append({ "comparison": comparison, "comparison_meaning": comparison_meaning, "left_scheme": left_scheme, "left_group_cn": l["group_cn"], "right_scheme": right_scheme, "right_group_cn": r["group_cn"], "process": process, "variable": variable, "variable_cn": variable_cn, "left_coef": lcoef, "right_coef": rcoef, "coef_delta_reference": rcoef - lcoef, "left_ame": lame, "right_ame": rame, "ame_delta": ame_delta, "abs_ame_delta": abs(ame_delta) if pd.notna(ame_delta) else np.nan, "ame_ci90_low": ame_ci90_low, "ame_ci90_high": ame_ci90_high, "ame_equivalence_margin": ame_margin, "ame_tost_p_value": ame_tost_p, "ame_tost_log_p_value": ame_tost_log_p, "ame_tost_p_exact": ame_tost_p_exact, "ame_in_equivalence_margin": bool(pd.notna(ame_delta) and abs(ame_delta) <= ame_margin), "substantive_direction_reversal": reversal, "direction_judgement": direction_judgement, "ame_judgement": ame_judgement, }) return pd.DataFrame(rows) def build_comparison_summary(comparison_df): rows = [] for comp, g in comparison_df.groupby("comparison"): rows.append({ "comparison": comp, "comparison_meaning": g["comparison_meaning"].iloc[0], "n_terms": int(len(g)), "n_no_substantive_reversal": int((g["substantive_direction_reversal"] == False).sum()), "n_ame_point_equivalent": int(g["ame_in_equivalence_margin"].sum()), "n_ame_tost_pass": int((g["ame_judgement"] == "AME-TOST通过").sum()), "n_substantive_reversal": int((g["direction_judgement"] == "实质性方向反转").sum()), }) return pd.DataFrame(rows) def print_report(group_design, network_summary, comparison_summary, comparison_df, output_dir): print("\n" + "=" * 120) print("三项筛选条件控制变量鲁棒性检验") print("=" * 120) print("\n一、组别设定") print(group_design[["scheme", "group_cn", "min_total_count", "min_active_periods", "edge_threshold", "purpose"]].to_string(index=False)) print("\n二、自动汇总") print(comparison_summary.to_string(index=False)) print("\n三、详细比较结果") show = comparison_df.copy() for col in [ "left_coef", "right_coef", "coef_delta_reference", "left_ame", "right_ame", "ame_delta", "abs_ame_delta", "ame_ci90_low", "ame_ci90_high", "ame_tost_p_value" ]: if col in show.columns: show[col] = show[col].map(lambda x: fmt_num(x, 3)) print(show.to_string(index=False)) print("\n输出目录：") print(output_dir) print("=" * 120 + "\n") def main(): parser = argparse.ArgumentParser() parser.add_argument("--input", default=DEFAULT_INPUT_FILE) parser.add_argument("--out", default=DEFAULT_OUTPUT_DIR) parser.add_argument("--ame-equivalence-margin", type=float, default=0.10) parser.add_argument("--alpha", type=float, default=0.05) parser.add_argument("--proximity-shrinkage", type=float, default=5.0) parser.add_argument("--proximity-mode", choices=["second_order_stable", "direct"], default="second_order_stable") parser.add_argument("--no-patent-size-cleaning", action="store_true") parser.add_argument("--save-panels", action="store_true") args = parser.parse_args() output_dir = Path(args.out) output_dir.mkdir(parents=True, exist_ok=True) print("Input file:", args.input) print("Output folder:", output_dir) print("Design: original raw group + one-condition-at-a-time control comparisons.") print("Technical proximity mode:", args.proximity_mode) print("Direction judgement: only substantive reversal is reported; weak direction difference is not calculated.") print("Default output is a single folder on Desktop.") long_df = load_patent_tech_long( args.input, apply_patent_size_filter=not args.no_patent_size_cleaning, ) panels = {} coef_tables = [] network_tables = [] sample_rows = [] group_rows = [] for spec in GROUP_SPECS: scheme = spec["scheme"] group_cn = spec["group_cn"] print("\nBuilding:", scheme, group_cn) roster, stat_df = build_roster( long_df, min_total_count=spec["min_total_count"], min_active_periods=spec["min_active_periods"], ) group_rows.append({ "scheme": scheme, "group_cn": group_cn, "min_total_count": "不设阈值" if spec["min_total_count"] is None else spec["min_total_count"], "min_active_periods": "不设阈值" if spec["min_active_periods"] is None else spec["min_active_periods"], "edge_threshold": spec["edge_threshold"], "roster_nodes": len(roster), "proximity_mode": args.proximity_mode, "purpose": spec["purpose"], }) sample_rows.append({ "scheme": scheme, "group_cn": group_cn, "roster_nodes": len(roster), "min_total_count": spec["min_total_count"], "min_active_periods": spec["min_active_periods"], "edge_threshold": spec["edge_threshold"], "proximity_mode": args.proximity_mode, }) network_df, panel_df = build_panel( scheme, group_cn, roster, long_df, edge_threshold=spec["edge_threshold"], proximity_shrinkage=args.proximity_shrinkage, proximity_mode=args.proximity_mode, ) panels[scheme] = panel_df network_tables.append(network_df) pd.DataFrame({"tech_code": roster}).to_csv(output_dir / ("roster_%s.csv" % scheme), index=False, encoding="utf-8-sig") if args.save_panels: panel_df.to_csv(output_dir / ("panel_%s.csv" % scheme), index=False, encoding="utf-8-sig") group_design = pd.DataFrame(group_rows) sample_summary = pd.DataFrame(sample_rows) network_summary = pd.concat(network_tables, ignore_index=True) baseline_panel = panels["G1_core_baseline_40_4_edge3"] scale_params = build_scale_params(baseline_panel) for spec in GROUP_SPECS: scheme = spec["scheme"] group_cn = spec["group_cn"] print("Fitting:", scheme, group_cn) coef_tables.append(fit_model(panels[scheme], scheme, group_cn, scale_params=scale_params)) coef_df = pd.concat(coef_tables, ignore_index=True) comp_tables = [] for comp in COMPARISONS: comp_tables.append( compare_models( coef_df, comp["left_scheme"], comp["right_scheme"], comp["comparison"], comp["meaning"], ame_margin=args.ame_equivalence_margin, alpha=args.alpha, ) ) comparison_df = pd.concat(comp_tables, ignore_index=True) comparison_summary = build_comparison_summary(comparison_df) group_design.to_csv(output_dir / "01_group_design.csv", index=False, encoding="utf-8-sig") sample_summary.to_csv(output_dir / "02_sample_summary.csv", index=False, encoding="utf-8-sig") network_summary.to_csv(output_dir / "03_network_summary.csv", index=False, encoding="utf-8-sig") coef_df.to_csv(output_dir / "04_model_coefficients_proximity_optimized_no_weak_direction.csv", index=False, encoding="utf-8-sig") comparison_df.to_csv(output_dir / "05_control_variable_comparison_detail_no_weak_direction.csv", index=False, encoding="utf-8-sig") comparison_summary.to_csv(output_dir / "06_control_variable_comparison_summary_no_weak_direction.csv", index=False, encoding="utf-8-sig") with pd.ExcelWriter(output_dir / "三项筛选条件控制变量鲁棒性检验结果_技术邻近性优化_无弱方向差异.xlsx", engine="openpyxl") as writer: group_design.to_excel(writer, sheet_name="group_design", index=False) sample_summary.to_excel(writer, sheet_name="sample_summary", index=False) network_summary.to_excel(writer, sheet_name="network_summary", index=False) coef_df.to_excel(writer, sheet_name="coefficients", index=False) comparison_df.to_excel(writer, sheet_name="comparison_detail", index=False) comparison_summary.to_excel(writer, sheet_name="comparison_summary", index=False) print_report(group_design, network_summary, comparison_summary, comparison_df, output_dir) if __name__ == "__main__": main()