"""Supplementary analysis code for: Quantum technology adoption in Japanese industry: prevalence, application interests, timelines, and barriers Python 3.11. Required packages: numpy, pandas, scipy, matplotlib, scikit-learn, openpyxl. The script expects the confidential merged sampling-frame/survey file to be available locally as Merged_all.xlsx. Firm-level data are not distributed with this supplementary code. """ # Configuration DATA_PATH = "Merged_all.xlsx" # merged survey + Corporate Data Bank file (full frame) N_SENT = 3890 # questionnaires mailed SURVEY_YEAR = 2025 RANDOM_SEED = 42 FIG_DIR = "figures_paper2" # --- Adoption stage (Q2); analyses use numeric codes parsed from raw labels --- STAGE_COL = "Q2" STAGE_LABELS_EN = { 6: "No interest", 5: "Interest only", 4: "PoC planned", 3: "PoC ongoing", 2: "Pilot/limited", 1: "Full deployment", } STAGE_ORDER_ASC = [6, 5, 4, 3, 2, 1] # least -> most advanced GROUP_ORDER = ["No interest", "Interest only", "Active (PoC planned+)"] # --- Resources (Q3) --- BUDGET_COL, STAFF_COL = "Q3_1", "Q3_2" BUDGET_LABELS_EN = {1: "No budget", 2: "< JPY 5M", 3: "JPY 5-20M", 4: "JPY 20-50M", 5: "JPY 50-100M", 6: ">= JPY 100M"} STAFF_LABELS_EN = {1: "0", 2: "1-5", 3: "6-10", 4: "11-20", 5: ">= 21"} # --- Q4_1 interest domains --- DOMAIN_LABELS = { "Q4_1_1": "Optimization", "Q4_1_2": "Finance applications", "Q4_1_3": "Supply chain", "Q4_1_4": "Quantum machine learning", "Q4_1_5": "Chemistry / simulation", "Q4_1_6": "CAE", "Q4_1_7": "Quantum computing hardware", "Q4_1_8": "Control / readout systems", "Q4_1_9": "Chip fabrication / packaging", "Q4_1_10": "Development environments", "Q4_1_11": "Theory development", "Q4_1_12": "Technology scouting / roadmapping", "Q4_1_13": "Execution / operations environments", "Q4_1_14": "Quantum communication / networks", "Q4_1_15": "Quantum-safe security", "Q4_1_16": "Quantum sensing", "Q4_1_17": "Other", } # --- Q4_2 roadmap milestones --- MILESTONE_LABELS = { "Q4_2_1": "Information gathering / preparation", "Q4_2_2": "Proof of concept", "Q4_2_3": "Limited commercialization", "Q4_2_4": "Full commercialization", } # The fielded questionnaire used code 5 for "No plans"; code 9 is also accepted for compatibility. ROADMAP_LABELS = {1: "By 2030", 2: "By 2035", 3: "By 2040", 4: "Undecided", 5: "No plans", 9: "No plans"} ROADMAP_ORDER = ["By 2030", "By 2035", "By 2040", "Undecided", "No plans"] # --- Q7 barriers --- BARRIER_LABELS = { "Q7_1": "Unclear return on investment", "Q7_2": "Internal quantum-talent shortage", "Q7_3": "Lack of standards / specifications", "Q7_4": "Implementation horizon too distant", "Q7_5": "Unclear use cases", "Q7_6": "No collaboration partners", "Q7_7": "Cost too high", "Q7_8": "Information-security requirements", "Q7_9": "Insufficient basis for decisions", } # --- Q6 facilitators --- FACILITATOR_LABELS = { 1: "Expert demonstration of business potential", 2: "Government financial support", 3: "Support for internal talent development", 4: "Technical / development-environment support", 5: "Regulation and legal frameworks", 6: "Other", } # --- Q1 respondent characteristics --- ROLE_COL, INVOLVE_COL, DEPT_COL = "Q1_2", "Q1_3", "Q1_1" ROLE_LABELS = {1: "Executive", 2: "Division head", 3: "Manager", 4: "Staff", 5: "Other"} INVOLVE_LABELS = {1: "Works directly with QT", 2: "Involved in QT projects", 3: "Information gathering only", 4: "Not involved"} # --- Corporate Data Bank columns (same conventions as Analysis_Revised.ipynb) --- MARKET_COL = "上場市場名" INDUSTRY_COL = "業種名" ADDRESS_COL = "住所" CAPITAL_KEYWORD = "資本金" EMPLOYEE_KEYWORDS = ["従業員", "社員数", "employees"] FOUNDING_KEYWORDS = ["設立", "創立"] # --- Firm-size classes (consolidated employees) --- EMP_BINS = [0, 300, 1000, 5000, float("inf")] EMP_LABELS = ["< 300", "300-999", "1,000-4,999", ">= 5,000"] # Imports and helper functions import os import re import warnings import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy import stats warnings.filterwarnings("ignore", category=FutureWarning) np.random.seed(RANDOM_SEED) os.makedirs(FIG_DIR, exist_ok=True) pd.set_option("display.width", 160) pd.set_option("display.max_columns", 60) pd.set_option("display.max_rows", 250) ZEN2HAN = str.maketrans("０１２３４５６７８９", "0123456789") def code_of(x): # Numeric code from raw values like "1．該当", "２．非該当", 3, "3.0". # Returns NaN for missing values and "回答なし". if pd.isna(x): return np.nan if isinstance(x, (int, float)) and not isinstance(x, bool): return float(x) s = str(x).strip().translate(ZEN2HAN) if s == "" or "回答なし" in s: return np.nan m = re.match(r"^\s*(\d+)", s) return float(m.group(1)) if m else np.nan def codes(series): return series.map(code_of) def wilson_ci(k, n, alpha=0.05): if n is None or n <= 0: return (np.nan, np.nan) z = stats.norm.ppf(1 - alpha / 2) p = k / n denom = 1 + z ** 2 / n centre = (p + z ** 2 / (2 * n)) / denom half = z * np.sqrt(p * (1 - p) / n + z ** 2 / (4 * n ** 2)) / denom return (max(0.0, centre - half), min(1.0, centre + half)) def kish_ess(w): w = np.asarray(w, dtype=float) return w.sum() ** 2 / (w ** 2).sum() if len(w) and (w ** 2).sum() > 0 else 0.0 def weighted_prop_ci(mask, w): # Weighted proportion with a Wilson CI on the Kish effective sample size. mask = np.asarray(mask, dtype=float) w = np.asarray(w, dtype=float) if w.sum() <= 0: return (np.nan, np.nan, np.nan) p = float((mask * w).sum() / w.sum()) n_eff = kish_ess(w) lo, hi = wilson_ci(p * n_eff, n_eff) return (p, lo, hi) def holm_adjust(pvals): # Holm step-down adjusted p-values. p = np.asarray(pvals, dtype=float) order = np.argsort(p) m = len(p) adj = np.empty(m) running = 0.0 for rank, idx in enumerate(order): running = max(running, (m - rank) * p[idx]) adj[idx] = min(1.0, running) return adj def epsilon_squared(h, n, k): # Effect size for Kruskal-Wallis. return (h - k + 1) / (n - k) if n > k else np.nan def fmt_pct(p, lo, hi): if pd.isna(p): return "—" return f"{p * 100:.1f}% [{lo * 100:.1f}-{hi * 100:.1f}]" def prop_row(k, n): lo, hi = wilson_ci(k, n) return fmt_pct(k / n if n else np.nan, lo, hi) def save_fig(fig, name): for ext in ("pdf", "png"): fig.savefig(os.path.join(FIG_DIR, f"{name}.{ext}"), bbox_inches="tight", dpi=300) print(f"[saved] {FIG_DIR}/{name}.pdf /.png") data = pd.read_excel(DATA_PATH) print(f"Data shape: {data.shape}") survey_cols = [c for c in data.columns if re.match(r"^Q\d", str(c))] if "is_respondent" not in data.columns: data["is_respondent"] = data[survey_cols].notna().any(axis=1) print(f"NOTE: no 'is_respondent' column — derived from {len(survey_cols)} " "survey columns (respondent = any survey item non-missing): " f"{int(data['is_respondent'].sum())} respondents identified.") data["is_respondent"] = data["is_respondent"].astype(bool) # Exclude returned questionnaires containing respondent information only and no substantive Q2+ responses. analysis_cols = [c for c in survey_cols if not re.match(r"^Q1(?:_|$)", str(c))] substantive = pd.DataFrame({ c: data[c].notna() & ~data[c].astype(str).str.contains("\u56de\u7b54\u306a\u3057") for c in analysis_cols }) blank_return = data["is_respondent"] & (substantive.sum(axis=1) == 0) data.loc[blank_return, "is_respondent"] = False if int(data["is_respondent"].sum()) != 154: raise ValueError( f"Expected 154 analyzable respondents after excluding blank returns; " f"found {int(data['is_respondent'].sum())}." ) HAS_FRAME = bool((~data["is_respondent"]).sum() > 0) # --- Firm attributes (same logic as Analysis_Revised.ipynb) --- def normalize_market_segment(s): if pd.isna(s): return "Unknown" x = re.sub(r"[\s　]+", "", str(s)) x = (x.replace("東京Ｐ", "東京(P)").replace("東京Ｓ", "東京(S)").replace("東京Ｇ", "東京(G)") .replace("東京（Ｐ）", "東京(P)").replace("東京（Ｓ）", "東京(S)").replace("東京（Ｇ）", "東京(G)")) if "東京(P)" in x: return "TSE Prime" if "東京(S)" in x: return "TSE Standard" if "東京(G)" in x: return "TSE Growth" if any(t in x for t in ["名古屋", "福岡", "札幌"]): return "Other domestic" return "Other/Unknown" INDUSTRY_GROUPS = { "情報・通信": "IT/Services", "サービス業": "IT/Services", "電気・ガス業": "IT/Services", "電気機器": "Electronics/Precision", "精密機器": "Electronics/Precision", "その他製品": "Electronics/Precision", "化学": "Chemicals/Pharma", "医薬品": "Chemicals/Pharma", "ゴム製品": "Chemicals/Pharma", "機械": "Manufacturing", "輸送用機器": "Manufacturing", "鉄鋼": "Manufacturing", "非鉄金属": "Manufacturing", "金属製品": "Manufacturing", "繊維製品": "Manufacturing", "パルプ・紙": "Manufacturing", "食料品": "Manufacturing", "窯業": "Manufacturing", "石油・石炭製品": "Manufacturing", "銀行業": "Finance", "証券、商品先物取引業": "Finance", "保険業": "Finance", "その他金融業": "Finance", } data["market_segment"] = (data[MARKET_COL].apply(normalize_market_segment) if MARKET_COL in data.columns else "Unknown") data["industry_group"] = (data[INDUSTRY_COL].map(INDUSTRY_GROUPS).fillna("Other") if INDUSTRY_COL in data.columns else "Other") def find_column(keywords): for c in data.columns: if any(k in str(c) for k in keywords): return c return None cap_col = find_column([CAPITAL_KEYWORD]) emp_col = find_column(EMPLOYEE_KEYWORDS) fnd_col = find_column(FOUNDING_KEYWORDS) print(f"Capital column: {cap_col}\nEmployee column: {emp_col}\nFounding column: {fnd_col}") data["capital"] = pd.to_numeric(data[cap_col], errors="coerce") if cap_col else np.nan data["employees"] = pd.to_numeric(data[emp_col], errors="coerce") if emp_col else np.nan if fnd_col: founding_year = data[fnd_col].astype(str).str.extract(r"(\d{4})")[0].astype(float) data["company_age"] = SURVEY_YEAR - founding_year else: data["company_age"] = np.nan data["emp_class"] = pd.cut(data["employees"], bins=EMP_BINS, labels=EMP_LABELS, right=False) # --- Headquarters prefecture and region --- PREFS = [ "北海道", "青森県", "岩手県", "宮城県", "秋田県", "山形県", "福島県", "茨城県", "栃木県", "群馬県", "埼玉県", "千葉県", "東京都", "神奈川県", "新潟県", "富山県", "石川県", "福井県", "山梨県", "長野県", "岐阜県", "静岡県", "愛知県", "三重県", "滋賀県", "京都府", "大阪府", "兵庫県", "奈良県", "和歌山県", "鳥取県", "島根県", "岡山県", "広島県", "山口県", "徳島県", "香川県", "愛媛県", "高知県", "福岡県", "佐賀県", "長崎県", "熊本県", "大分県", "宮崎県", "鹿児島県", "沖縄県", ] REGION_OF = {} for p in PREFS: if p == "北海道": REGION_OF[p] = "Hokkaido/Tohoku" elif p in ["青森県", "岩手県", "宮城県", "秋田県", "山形県", "福島県"]: REGION_OF[p] = "Hokkaido/Tohoku" elif p in ["茨城県", "栃木県", "群馬県", "埼玉県", "千葉県", "神奈川県"]: REGION_OF[p] = "Kanto (excl. Tokyo)" elif p == "東京都": REGION_OF[p] = "Tokyo" elif p in ["新潟県", "富山県", "石川県", "福井県", "山梨県", "長野県", "岐阜県", "静岡県", "愛知県"]: REGION_OF[p] = "Chubu" elif p in ["三重県", "滋賀県", "京都府", "大阪府", "兵庫県", "奈良県", "和歌山県"]: REGION_OF[p] = "Kinki" elif p in ["鳥取県", "島根県", "岡山県", "広島県", "山口県", "徳島県", "香川県", "愛媛県", "高知県"]: REGION_OF[p] = "Chugoku/Shikoku" else: REGION_OF[p] = "Kyushu/Okinawa" # Fallbacks for addresses that omit the prefecture: Tokyo's 23 wards and the # major designated cities. Ward names are only matched at the start of the string. TOKYO_WARDS = ["千代田区", "中央区", "港区", "新宿区", "文京区", "台東区", "墨田区", "江東区", "品川区", "目黒区", "大田区", "世田谷区", "渋谷区", "中野区", "杉並区", "豊島区", "北区", "荒川区", "板橋区", "練馬区", "足立区", "葛飾区", "江戸川区"] CITY_TO_PREF = { # designated cities "札幌市": "北海道", "仙台市": "宮城県", "さいたま市": "埼玉県", "千葉市": "千葉県", "横浜市": "神奈川県", "川崎市": "神奈川県", "相模原市": "神奈川県", "新潟市": "新潟県", "静岡市": "静岡県", "浜松市": "静岡県", "名古屋市": "愛知県", "京都市": "京都府", "大阪市": "大阪府", "堺市": "大阪府", "神戸市": "兵庫県", "岡山市": "岡山県", "広島市": "広島県", "北九州市": "福岡県", "福岡市": "福岡県", "熊本市": "熊本県", # remaining prefectural capitals (addresses in the data often omit the prefecture) "青森市": "青森県", "盛岡市": "岩手県", "秋田市": "秋田県", "山形市": "山形県", "福島市": "福島県", "水戸市": "茨城県", "宇都宮市": "栃木県", "前橋市": "群馬県", "富山市": "富山県", "金沢市": "石川県", "福井市": "福井県", "甲府市": "山梨県", "長野市": "長野県", "岐阜市": "岐阜県", "津市": "三重県", "大津市": "滋賀県", "奈良市": "奈良県", "和歌山市": "和歌山県", "鳥取市": "鳥取県", "松江市": "島根県", "山口市": "山口県", "徳島市": "徳島県", "高松市": "香川県", "松山市": "愛媛県", "高知市": "高知県", "佐賀市": "佐賀県", "長崎市": "長崎県", "大分市": "大分県", "宮崎市": "宮崎県", "鹿児島市": "鹿児島県", "那覇市": "沖縄県", } def prefecture_of(addr): if pd.isna(addr): return np.nan s = str(addr).strip() for p in PREFS: if s.startswith(p) or p in s[:8]: return p for city, pref in CITY_TO_PREF.items(): if s.startswith(city): return pref for ward in TOKYO_WARDS: if s.startswith(ward): return "東京都" return np.nan if ADDRESS_COL in data.columns: data["prefecture"] = data[ADDRESS_COL].map(prefecture_of) data["region"] = data["prefecture"].map(REGION_OF) n_addr_missing = int(data[ADDRESS_COL].isna().sum()) n_unparsed = int((data[ADDRESS_COL].notna() & data["prefecture"].isna()).sum()) print(f"Address parsing: {n_addr_missing} missing addresses, " f"{n_unparsed} non-missing but unparsable.") if n_unparsed: bad = data.loc[data[ADDRESS_COL].notna() & data["prefecture"].isna(), ADDRESS_COL].astype(str).str[:6] print("Most common unparsable address prefixes (first 6 chars):") print(bad.value_counts().head(10).to_string()) else: print(f"NOTE: address column '{ADDRESS_COL}' not found — regional analyses skipped.") data["prefecture"] = np.nan data["region"] = np.nan # --- Respondent subset and adoption-stage coding --- resp = data[data["is_respondent"]].copy() resp["stage_code"] = codes(resp[STAGE_COL]) resp["stage_label"] = resp["stage_code"].map(STAGE_LABELS_EN) def stage_group_of(c): if pd.isna(c): return np.nan if c == 6: return "No interest" if c == 5: return "Interest only" return "Active (PoC planned+)" resp["stage_group"] = resp["stage_code"].map(stage_group_of) resp["engaged"] = np.where(resp["stage_code"].isna(), np.nan, (resp["stage_code"] <= 5).astype(float)) resp["active"] = np.where(resp["stage_code"].isna(), np.nan, (resp["stage_code"] <= 4).astype(float)) resp["advanced"] = np.where(resp["stage_code"].isna(), np.nan, (resp["stage_code"] <= 3).astype(float)) # --- STROBE-style flow --- n_frame = len(data) n_resp = len(resp) n_stage = int(resp["stage_code"].notna().sum()) print("\n=== Participant flow ===") print(f"Questionnaires mailed: {N_SENT}") print(f"Rows in merged frame: {n_frame}") print(f"Responding firms: {n_resp} " f"({n_resp / N_SENT * 100:.1f}% crude response rate)") print(f" with non-missing adoption stage: {n_stage}") print(f"Frame includes non-respondents: {HAS_FRAME}") print("\n" + "=" * 62) if n_resp == EXPECTED_RESPONDENTS: print(f"CONSISTENCY CHECK PASSED: {n_resp} respondents " "(matches the companion paper).") else: print(f"*** CONSISTENCY CHECK FAILED: {n_resp} respondents, " f"expected {EXPECTED_RESPONDENTS}. ***") print("*** Outputs below are NOT final — stop and report this run. ***") print("=" * 62) print(f"=== Table 1 — responding firms (n = {len(resp)}) ===") for col, label in [("industry_group", "Industry group"), ("market_segment", "Market segment"), ("emp_class", "Employees (consolidated)"), ("region", "Headquarters region")]: vc = resp[col].value_counts(dropna=False) tab = pd.DataFrame({ "n": vc, "%": (vc / len(resp) * 100).round(1), }) print(f"\n--- {label} ---") print(tab.to_string()) if INDUSTRY_COL in resp.columns: print("\n--- Detailed industry (top 12) ---") vc = resp[INDUSTRY_COL].value_counts().head(12) print(pd.DataFrame({"n": vc, "%": (vc / len(resp) * 100).round(1)}).to_string()) print("\n--- Continuous firm attributes (median [IQR]) ---") for col, label in [("employees", "Employees"), ("capital", "Paid-in capital"), ("company_age", "Company age (years)")]: v = resp[col].dropna() if len(v): print(f"{label:<22}: {v.median():,.0f} " f"[{v.quantile(.25):,.0f}-{v.quantile(.75):,.0f}] (n = {len(v)})") print("\n=== Informant characteristics ===") for col, labels, name in [(ROLE_COL, ROLE_LABELS, "Position (Q1_2)"), (INVOLVE_COL, INVOLVE_LABELS, "Involvement with QT (Q1_3)")]: if col in resp.columns: c = codes(resp[col]).map(labels) vc = c.value_counts(dropna=False) print(f"\n--- {name} ---") print(pd.DataFrame({"n": vc, "%": (vc / len(resp) * 100).round(1)}).to_string()) else: print(f"NOTE: column {col} not found.") resp["ipw"] = 1.0 ipw_raw = None if not HAS_FRAME: print("Frame contains respondents only — IPW skipped; all weights set to 1.") else: from sklearn.linear_model import LogisticRegression from sklearn.metrics import roc_auc_score CAT_PREDICTORS = ["market_segment", "industry_group"] CONT_PREDICTORS = ["employees", "capital", "company_age"] def build_design_matrix(df): X = pd.get_dummies(df[CAT_PREDICTORS], drop_first=True) for c in CONT_PREDICTORS: v = pd.to_numeric(df[c], errors="coerce") X[c] = np.log1p(v.fillna(v.median()).clip(lower=0)) return X.astype(float) X_all = build_design_matrix(data) y_all = data["is_respondent"].astype(int).values prop_model = LogisticRegression(max_iter=5000, C=1e6) prop_model.fit(X_all, y_all) phat = prop_model.predict_proba(X_all)[:, 1] auc = roc_auc_score(y_all, phat) print(f"Response-propensity model: AUC = {auc:.3f} " f"(n = {len(y_all)}, respondents = {y_all.sum()})") phat_s = pd.Series(phat, index=data.index) p_resp = phat_s.loc[resp.index].clip(lower=1e-4) ipw_raw = 1.0 / p_resp # Horvitz-Thompson weights (totals) cap = ipw_raw.quantile(0.99) w = ipw_raw.clip(upper=cap) resp["ipw"] = w / w.mean() # trimmed + normalized (proportions) print("\nPropensity among respondents: " f"median {p_resp.median():.4f}, range {p_resp.min():.4f}-{p_resp.max():.4f}") print(f"Raw weights: median {ipw_raw.median():.1f}, max {ipw_raw.max():.1f}; " f"trimmed at 99th pct = {cap:.1f}") print(f"Kish effective sample size: {kish_ess(resp['ipw']):.1f} " f"(nominal n = {len(resp)})") st = resp.dropna(subset=["stage_code"]).copy() n = len(st) print(f"Respondents with non-missing adoption stage: n = {n}\n") rows = [] for c in STAGE_ORDER_ASC: k = int((st["stage_code"] == c).sum()) lo, hi = wilson_ci(k, n) row = {"Stage": STAGE_LABELS_EN[c], "n": k, "Unweighted % [95% CI]": fmt_pct(k / n, lo, hi)} if HAS_FRAME: pw, lw, hw = weighted_prop_ci((st["stage_code"] == c).values, st["ipw"].values) row["IPW % [95% CI]"] = fmt_pct(pw, lw, hw) rows.append(row) print("=== Adoption-stage distribution ===") print(pd.DataFrame(rows).to_string(index=False)) print("\n=== Binary outcomes ===") for col, label in [("engaged", "Engaged (interest or beyond)"), ("active", "Active (PoC planned or beyond)"), ("advanced", "Advanced (PoC ongoing or beyond)")]: k = int(st[col].sum()) lo, hi = wilson_ci(k, n) line = f"{label:<35}: {k}/{n} {fmt_pct(k / n, lo, hi)}" if HAS_FRAME: pw, lw, hw = weighted_prop_ci(st[col].values, st["ipw"].values) line += f" | IPW: {fmt_pct(pw, lw, hw)}" print(line) def subgroup_prevalence(group_col, title): sub = st.dropna(subset=[group_col]) if st[group_col].isna().any() else st if sub[group_col].nunique() < 2: print(f"\n--- {title}: fewer than 2 groups, skipped ---") return rows = [] for g, grp in sub.groupby(group_col, observed=True): m = len(grp) rows.append({ title: str(g), "n": m, "Engaged % [CI]": prop_row(int(grp["engaged"].sum()), m), "Active % [CI]": prop_row(int(grp["active"].sum()), m), }) print(f"\n--- {title} ---") print(pd.DataFrame(rows).sort_values("n", ascending=False).to_string(index=False)) ct = pd.crosstab(sub[group_col], sub["active"]) if ct.shape[0] >= 2 and ct.shape[1] == 2: chi2, p, dof, exp = stats.chi2_contingency(ct) note = " (caution: expected counts < 5)" if (exp < 5).any() else "" print(f"Chi-square for 'Active' across groups: " f"chi2(df={dof}) = {chi2:.2f}, p = {p:.4f}{note}") for col, title in [("industry_group", "Industry group"), ("emp_class", "Employees (consolidated)"), ("market_segment", "Market segment"), ("region", "Headquarters region")]: subgroup_prevalence(col, title) # --- Figure 1: stage distribution by industry group (stacked proportions) --- pivot = (pd.crosstab(st["industry_group"], st["stage_code"].astype(int), normalize="index") .reindex(columns=STAGE_ORDER_ASC, fill_value=0)) counts = st["industry_group"].value_counts() pivot = pivot.loc[counts.index] # order by group size fig, ax = plt.subplots(figsize=(9, 0.6 * len(pivot) + 2)) left = np.zeros(len(pivot)) cmap = plt.get_cmap("viridis") for i, c in enumerate(STAGE_ORDER_ASC): vals = pivot[c].values ax.barh(pivot.index, vals, left=left, color=cmap(i / 5), label=STAGE_LABELS_EN[c], edgecolor="white") left += vals ax.set_xlabel("Proportion of responding firms") ax.set_xlim(0, 1) ax.invert_yaxis() yticklabels = [f"{g} (n={counts[g]})" for g in pivot.index] ax.set_yticks(range(len(pivot))) ax.set_yticklabels(yticklabels) ax.legend(bbox_to_anchor=(1.02, 1), loc="upper left", frameon=False) ax.set_title("Adoption stage by industry group") save_fig(fig, "fig1_stage_by_industry") plt.show() dom_cols = [c for c in DOMAIN_LABELS if c in resp.columns] missing_cols = [c for c in DOMAIN_LABELS if c not in resp.columns] if missing_cols: print(f"NOTE: missing Q4_1 columns: {missing_cols}") if not dom_cols: print("No Q4_1 columns found — section skipped.") else: dom = pd.DataFrame({c: (codes(resp[c]) == 1) for c in dom_cols}, index=resp.index) nonmiss = pd.DataFrame({c: codes(resp[c]).notna() for c in dom_cols}, index=resp.index) eng_idx = resp.index[resp["engaged"] == 1] rows = [] n_all, n_eng = len(resp), len(eng_idx) for c in dom_cols: k_all = int(dom[c].sum()) k_eng = int(dom.loc[eng_idx, c].sum()) rows.append({ "Domain": DOMAIN_LABELS[c], "n selected": k_all, "% of all resp. [CI]": prop_row(k_all, n_all), "% of engaged [CI]": prop_row(k_eng, n_eng), "n non-missing": int(nonmiss[c].sum()), }) tab = pd.DataFrame(rows).sort_values("n selected", ascending=False) print(f"=== Interest domains (all respondents n = {n_all}; engaged n = {n_eng}) ===") print(tab.to_string(index=False)) # Number of domains per firm, by stage group resp["n_domains"] = dom.sum(axis=1) print("\n--- Number of selected domains per firm, by stage group ---") desc = (resp.dropna(subset=["stage_group"]) .groupby("stage_group")["n_domains"] .agg(["count", "mean", "median", "max"]) .reindex(GROUP_ORDER).round(2)) print(desc.to_string()) # Figure 2: domain prevalence among engaged firms plot_tab = tab[tab["Domain"] != "Other"].copy() ks = [int(dom.loc[eng_idx, c].sum()) for c in dom_cols if DOMAIN_LABELS[c] != "Other"] labels = [DOMAIN_LABELS[c] for c in dom_cols if DOMAIN_LABELS[c] != "Other"] order = np.argsort(ks) fig, ax = plt.subplots(figsize=(8, 0.4 * len(ks) + 1.5)) ps = [k / n_eng * 100 for k in ks] errs = np.array([[k / n_eng * 100 - wilson_ci(k, n_eng)[0] * 100, wilson_ci(k, n_eng)[1] * 100 - k / n_eng * 100] for k in ks]) ax.barh(np.arange(len(ks)), np.array(ps)[order], xerr=errs[order].T, color="#3b6ea5", capsize=2) ax.set_yticks(np.arange(len(ks))) ax.set_yticklabels(np.array(labels)[order]) ax.set_xlabel(f"% of engaged firms (n = {n_eng}), Wilson 95% CI") ax.set_title("Quantum technology interest domains") save_fig(fig, "fig2_interest_domains") plt.show() # Industry x domain heatmap (groups with n >= 10) big_groups = [g for g, m in resp["industry_group"].value_counts().items() if m >= 10] H = pd.DataFrame(index=big_groups, columns=[DOMAIN_LABELS[c] for c in dom_cols if DOMAIN_LABELS[c] != "Other"], dtype=float) for g in big_groups: idx = resp.index[resp["industry_group"] == g] for c in dom_cols: if DOMAIN_LABELS[c] == "Other": continue H.loc[g, DOMAIN_LABELS[c]] = dom.loc[idx, c].mean() * 100 fig, ax = plt.subplots(figsize=(12, 0.55 * len(big_groups) + 2.5)) im = ax.imshow(H.values.astype(float), aspect="auto", cmap="YlGnBu") ax.set_xticks(range(H.shape[1])) ax.set_xticklabels(H.columns, rotation=45, ha="right") ax.set_yticks(range(H.shape[0])) ax.set_yticklabels([f"{g} (n={resp['industry_group'].value_counts()[g]})" for g in H.index]) for i in range(H.shape[0]): for j in range(H.shape[1]): ax.text(j, i, f"{H.values[i, j]:.0f}", ha="center", va="center", fontsize=7, color="white" if H.values[i, j] > 50 else "black") fig.colorbar(im, ax=ax, label="% of firms in group") ax.set_title("Interest domains by industry group") save_fig(fig, "fig3_domain_by_industry") plt.show() # Co-occurrence (Jaccard) among domains selected by >= 10 firms top_cols = [c for c in dom_cols if dom[c].sum() >= 10 and DOMAIN_LABELS[c] != "Other"] if len(top_cols) >= 3: J = pd.DataFrame(index=[DOMAIN_LABELS[c] for c in top_cols], columns=[DOMAIN_LABELS[c] for c in top_cols], dtype=float) for a in top_cols: for b in top_cols: inter = (dom[a] & dom[b]).sum() union = (dom[a] | dom[b]).sum() J.loc[DOMAIN_LABELS[a], DOMAIN_LABELS[b]] = ( inter / union if union else np.nan) print("\n--- Domain co-occurrence (Jaccard), domains selected by >= 10 firms ---") print(J.round(2).to_string()) ms_cols = [c for c in MILESTONE_LABELS if c in resp.columns] if not ms_cols: print("No Q4_2 columns found — section skipped.") else: eng = resp[resp["engaged"] == 1] print(f"Denominator: engaged firms (n = {len(eng)}); " "percentages among those answering each item.\n") CATS = ROADMAP_ORDER + ["Unmapped"] dist_rows, cum_rows = [], [] props = {} for c in ms_cols: v = codes(eng[c]) answered = v.dropna() m = len(answered) counts = {lab: 0 for lab in ROADMAP_ORDER} for code, lab in ROADMAP_LABELS.items(): counts[lab] += int((answered == code).sum()) counts["Unmapped"] = m - sum(counts.values()) if counts["Unmapped"]: unk = answered[~answered.isin(list(ROADMAP_LABELS))] raw_vals = eng.loc[unk.index, c].astype(str).value_counts() print(f"WARNING: {MILESTONE_LABELS[c]}: {counts['Unmapped']} answers " f"with codes outside ROADMAP_LABELS {sorted(set(unk))}.") print(" Raw values (add them to ROADMAP_LABELS in the config cell):") print(" " + raw_vals.head(8).to_string().replace("\n", "\n ")) row = {"Milestone": MILESTONE_LABELS[c], "n answered": m} row.update({lab: f"{counts[lab]} ({counts[lab] / m * 100:.0f}%)" if m else "—" for lab in CATS}) dist_rows.append(row) props[c] = ([counts[lab] / m for lab in CATS] if m else [0] * len(CATS)) if m: c30 = counts["By 2030"] c35 = c30 + counts["By 2035"] c40 = c35 + counts["By 2040"] cum_rows.append({ "Milestone": MILESTONE_LABELS[c], "Cum. % by 2030": f"{c30 / m * 100:.0f}%", "Cum. % by 2035": f"{c35 / m * 100:.0f}%", "Cum. % by 2040": f"{c40 / m * 100:.0f}%", "Undecided": f"{counts['Undecided'] / m * 100:.0f}%", "No plans": f"{counts['No plans'] / m * 100:.0f}%", }) print("=== Milestone timing distribution ===") print(pd.DataFrame(dist_rows).to_string(index=False)) print("\n=== Cumulative expectation (among answering firms) ===") print(pd.DataFrame(cum_rows).to_string(index=False)) fig, ax = plt.subplots(figsize=(9, 3.5)) cmap = plt.get_cmap("plasma") colors = [cmap(i / (len(CATS) - 1)) for i in range(len(CATS) - 1)] + ["#bbbbbb"] ylabels = [MILESTONE_LABELS[c] for c in ms_cols] left = np.zeros(len(ms_cols)) for j, lab in enumerate(CATS): vals = np.array([props[c][j] for c in ms_cols]) if vals.sum() == 0: continue # keep empty categories (e.g. Unmapped) out of the legend ax.barh(ylabels, vals, left=left, color=colors[j], label=lab, edgecolor="white") left += vals ax.set_xlim(0, 1) ax.invert_yaxis() ax.set_xlabel("Proportion of answering engaged firms") ax.legend(bbox_to_anchor=(1.02, 1), loc="upper left", frameon=False) ax.set_title("Expected timing of adoption milestones") save_fig(fig, "fig4_roadmap") plt.show() bar_cols = [c for c in BARRIER_LABELS if c in resp.columns] if not bar_cols: print("No Q7 columns found — section skipped.") else: B = pd.DataFrame({c: codes(resp[c]) for c in bar_cols}, index=resp.index) rows = [] for c in bar_cols: v = B[c].dropna() m = len(v) k = int((v >= 4).sum()) rows.append({"Barrier": BARRIER_LABELS[c], "n": m, "Mean": round(v.mean(), 2), "SD": round(v.std(ddof=1), 2), "Median": v.median(), "% rating 4-5 [CI]": prop_row(k, m)}) tab = pd.DataFrame(rows).sort_values("Mean", ascending=False) print("=== Barrier ratings (all responding firms) ===") print(tab.to_string(index=False)) # Figure 5: % rating 4-5 with CI ks = [int((B[c].dropna() >= 4).sum()) for c in bar_cols] ms = [int(B[c].notna().sum()) for c in bar_cols] ps = [k / m * 100 if m else np.nan for k, m in zip(ks, ms)] errs = np.array([[p - wilson_ci(k, m)[0] * 100, wilson_ci(k, m)[1] * 100 - p] for p, k, m in zip(ps, ks, ms)]) order = np.argsort(ps) labels = np.array([BARRIER_LABELS[c] for c in bar_cols]) fig, ax = plt.subplots(figsize=(8, 0.45 * len(bar_cols) + 1.5)) ax.barh(np.arange(len(bar_cols)), np.array(ps)[order], xerr=errs[order].T, color="#a53b3b", capsize=2) ax.set_yticks(np.arange(len(bar_cols))) ax.set_yticklabels(labels[order]) ax.set_xlabel("% rating the item a barrier (4-5), Wilson 95% CI") ax.set_title("Perceived barriers to quantum technology adoption") save_fig(fig, "fig5_barriers") plt.show() # Stage-group comparison — uses only groups with n >= 5 per barrier # (in the real data, barrier items were answered only by engaged firms, # so the "No interest" group is typically empty). g = resp["stage_group"] group_sizes = {grp: int(B.loc[g == grp].notna().any(axis=1).sum()) for grp in GROUP_ORDER} print(f"\nFirms answering any barrier item, by stage group: {group_sizes}") kw_rows, pvals, ok = [], [], [] for c in bar_cols: present = [(grp, B.loc[g == grp, c].dropna()) for grp in GROUP_ORDER] present = [(grp, s) for grp, s in present if len(s) >= 5] if len(present) < 2: print(f"NOTE: {BARRIER_LABELS[c]}: fewer than 2 groups with n >= 5 — skipped.") continue samples = [s for _, s in present] k = len(samples) h, p = stats.kruskal(*samples) ntot = sum(len(s) for s in samples) row = {"Barrier": BARRIER_LABELS[c], "groups": k, "H": round(h, 2), "p": p, "eps^2": round(epsilon_squared(h, ntot, k), 3)} for grp in GROUP_ORDER: s = B.loc[g == grp, c].dropna() row[grp] = round(s.mean(), 2) if len(s) >= 5 else np.nan kw_rows.append(row) pvals.append(p) ok.append((c, present)) if kw_rows: kw = pd.DataFrame(kw_rows) kw["p (Holm)"] = holm_adjust(pvals).round(4) kw["p"] = kw["p"].round(4) print("\n=== Barrier ratings by stage group (means; Kruskal-Wallis on " "groups with n >= 5; 2-group case equals Mann-Whitney) ===") print(kw.to_string(index=False)) sig = [(c, present) for (c, present), padj in zip(ok, holm_adjust(pvals)) if padj < 0.05] for c, present in sig: if len(present) < 3: continue # omnibus already pairwise print(f"\n--- Pairwise Mann-Whitney (Holm within barrier): {BARRIER_LABELS[c]} ---") pairs = [(a, b) for i, a in enumerate(present) for b in present[i + 1:]] praw = [] for (ga, sa), (gb, sb) in pairs: _, p = stats.mannwhitneyu(sa, sb, alternative="two-sided") praw.append(p) padj = holm_adjust(praw) for ((ga, _), (gb, _)), p, q in zip(pairs, praw, padj): print(f" {ga} vs {gb}: p = {p:.4f}, Holm p = {q:.4f}") q6_cols = [c for c in resp.columns if re.match(r"^Q6", str(c)) and not re.search(r"\u81ea\u7531|free|text|\u305d\u306e\u4ed6", str(c), re.IGNORECASE)] eng_q6 = resp[resp["engaged"] == 1].copy() if not q6_cols: print("No coded Q6 column found; facilitator analysis skipped.") elif len(q6_cols) == 1: col = q6_cols[0] facilitator = codes(eng_q6[col]).map(FACILITATOR_LABELS) facilitator = facilitator.dropna() n_q6 = len(facilitator) vc = facilitator.value_counts().reindex(FACILITATOR_LABELS.values(), fill_value=0) rows = [] for label, k in vc.items(): rows.append({"Facilitator": label, "n": int(k), "% [CI]": prop_row(int(k), n_q6)}) print(f"=== Facilitators among engaged firms answering Q6 (n = {n_q6}) ===") print(pd.DataFrame(rows).sort_values("n", ascending=False).to_string(index=False)) else: n_q6 = len(eng_q6) rows = [] for i, col in enumerate(sorted(q6_cols), start=1): ind = codes(eng_q6[col]) == 1 k = int(ind.sum()) rows.append({"Facilitator": FACILITATOR_LABELS.get(i, col), "n": k, "% [CI]": prop_row(k, n_q6)}) print(f"=== Facilitators among engaged firms (n = {n_q6}) ===") print(pd.DataFrame(rows).sort_values("n", ascending=False).to_string(index=False)) for col, labels, name in [(BUDGET_COL, BUDGET_LABELS_EN, "Quantum budget (Q3_1)"), (STAFF_COL, STAFF_LABELS_EN, "Dedicated personnel (Q3_2)")]: if col not in resp.columns: print(f"NOTE: column {col} not found — skipped.") continue v = codes(resp[col]) lab = v.map(labels) m = int(v.notna().sum()) vc = lab.value_counts().reindex(labels.values(), fill_value=0) print(f"\n=== {name} (n answered = {m}; answered predominantly by engaged " "firms — denominator is answering firms) ===") print(pd.DataFrame({"n": vc, "%": (vc / m * 100).round(1)}).to_string()) sub = resp.dropna(subset=["stage_group"]).copy() sub["_lab"] = codes(sub[col]).map(labels) ct = pd.crosstab(sub["stage_group"], sub["_lab"]).reindex(GROUP_ORDER) ct = ct.reindex(columns=[x for x in labels.values() if x in ct.columns]) pct = ct.div(ct.sum(axis=1), axis=0).mul(100).round(1) print("\n--- Counts by stage group ---") print(ct.to_string()) print("--- Row % by stage group ---") print(pct.to_string()) any_col = "any_budget" if col == BUDGET_COL else "any_staff" resp[any_col] = np.where(v.isna(), np.nan, (v >= 2).astype(float)) vv = resp[any_col].dropna() k = int(vv.sum()) line = f"\nAny {('budget' if col == BUDGET_COL else 'dedicated staff')}: " f"{k}/{len(vv)} {prop_row(k, len(vv))}" if HAS_FRAME: w = resp.loc[vv.index, "ipw"].values pw, lw, hw = weighted_prop_ci(vv.values, w) line += f" | IPW: {fmt_pct(pw, lw, hw)}" print(line) # Sensitivity: stage-known no-interest firms that skipped the item are # assumed to have no budget / no dedicated staff (denominator = all # stage-known respondents with observed or imputable values). base = resp.dropna(subset=["stage_code"]) v2 = codes(base[col]) v2 = v2.where(~((base["stage_code"] == 6) & v2.isna()), 1.0) vv2 = (v2[v2.notna()] >= 2).astype(float) k2 = int(vv2.sum()) print(f" Sensitivity (no-interest non-answerers assumed '{labels[1]}'): " f"{k2}/{len(vv2)} {prop_row(k2, len(vv2))}") if not HAS_FRAME or ipw_raw is None: print("Frame-level weights unavailable — population estimates skipped.") else: targets = [("engaged", "Engaged (interest or beyond)"), ("active", "Active (PoC planned or beyond)"), ("advanced", "Advanced (PoC ongoing or beyond)"), ("any_budget", "Any quantum budget"), ("any_staff", "Any dedicated quantum staff")] print(f"Frame size: {len(data)} firms; mailed: {N_SENT}\n") rows = [] for col, label in targets: if col not in resp.columns: continue v = resp[col].dropna() w = ipw_raw.loc[v.index] est = float((w * v).sum()) rows.append({"Outcome": label, "Respondents with outcome": int(v.sum()), "HT estimated firms (of frame)": round(est), "HT estimated %": f"{est / len(data) * 100:.1f}%"}) print(pd.DataFrame(rows).to_string(index=False)) print("\nNOTE: firms with missing outcome data are excluded from each estimate " "(conservative); weights are untrimmed Horvitz-Thompson weights.") print("=" * 64) print("SUMMARY — Descriptive epidemiology of quantum technology adoption") print("=" * 64) print(f"Frame: {len(data)} firms (mailed {N_SENT}); respondents: {len(resp)}; " f"with adoption stage: {int(resp['stage_code'].notna().sum())}") stx = resp.dropna(subset=["stage_code"]) for col, label in [("engaged", "Engaged"), ("active", "Active"), ("advanced", "Advanced")]: k = int(stx[col].sum()) print(f"{label:<10}: {k}/{len(stx)} {prop_row(k, len(stx))}") if "n_domains" in resp.columns: print(f"Median interest domains per firm: {resp['n_domains'].median():.0f}") bar_cols_s = [c for c in BARRIER_LABELS if c in resp.columns] if bar_cols_s: means = {BARRIER_LABELS[c]: codes(resp[c]).mean() for c in bar_cols_s} top = max(means, key=means.get) print(f"Highest-rated barrier: {top} (mean {means[top]:.2f})") print(f"\nFigures saved under: {FIG_DIR}/") print("Done.")