# -*- coding: utf-8 -*- """ X03: 観光統計 4 指標のペアプロット解析 — 散布図行列・PCA・階層クラスタで広島県市町を読み解く ================================================================================ DoBoX が公開する観光関連データセット 4 件: #1567 広島県クルーズ船観光客分析データ (令和6年度) #1282 瀬戸内しまたびライン利用状況 #1280 せとうちモニタークルーズ実施結果 (2022年12月) #1447 インフラツーリズム_施設情報 これらを **市町別 4 次元ベクトル** に統合し、 1) pandas の散布図行列 (scatter_matrix) で全 4×4 ペアを俯瞰 2) 単回帰直線 + R² を各ペアに重ね描き 3) 相関ヒートマップで相関行列を 1 枚に圧縮 4) **PCA** で 4 次元 → 2 次元、PC1×PC2 散布 + ローディング棒 5) **重回帰** (OLS) でインフラツーリズム掲載数を他 3 指標から予測 6) **階層クラスタリング** (Ward 法) で市町を 3〜4 群に分ける の 6 分析を 10 セクション L01 水準で展開する。 データの実態: - #1567 / #1282 / #1280 は DoBoX 限定で resource_download が無く、 ローカル CSV として取得できない (S35/S47/S57 でも同じ制約を確認済み)。 - そこで本記事では **「実データ + カタログメタ + 物理インフラの代理指標」** の 3 層で 4 指標を構築する: X1 クルーズ寄港強度 = sea_route.csv (#1278) の市町別寄港地数 (広島県内のみ) X2 しまたび航路強度 = sea_route.csv のうち港名 or 住所に「島」を含む寄港地数 (#1282 「瀬戸内しまたび」 = 島嶼間航路の代理) X3 モニター実施強度 = sea_route.csv の市町別住所末端ユニーク数 (#1280 「モニタークルーズ」 = 多地点実証の代理) X4 インフラツーリズム = infra_tourism.csv (#1447, 40 施設) を最近傍寄港地に 紐付け、市町別の施設数 すべての X1-X4 は **sea_route.csv (実データ) + infra_tourism.csv (実データ) + dataset_index.csv のカタログメタ** から派生しており、 数値そのものは本物の物理インフラ・施設配置を反映する。 実行: cd "2026 DoBoX 教材" py -X utf8 lessons/X03_tourism_4indicators_pairplot.py 制約: - 並列処理なし、メモリ 50MB 未満 - DoBoX 限定データの取得は試行のみ (resource_download なしを確認後フォールバック) """ from __future__ import annotations import re import math import json import warnings from pathlib import Path import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.patches import Patch from pandas.plotting import scatter_matrix from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LinearRegression from scipy.cluster.hierarchy import linkage, fcluster, dendrogram from scipy.spatial.distance import pdist from _common import ROOT, ASSETS, LESSONS, render_lesson, code, figure, data_recipe plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False warnings.filterwarnings("ignore", category=UserWarning) RNG = np.random.default_rng(42) # === パス定義 ================================================================= IDX_CSV = ROOT / "data" / "dataset_index.csv" ROUTE_CSV = ROOT / "data" / "extras" / "sea_route.csv" INFRA_CSV = ROOT / "data" / "extras" / "infra_tourism.csv" # 中間 CSV (assets 配下 = HTML から直リンク) OUT_INDICATORS = ASSETS / "X03_indicators_matrix.csv" # 市町×4指標 OUT_CORR = ASSETS / "X03_correlation_matrix.csv" # 4×4 相関行列 OUT_REG = ASSETS / "X03_regression_pairs.csv" # 6 ペアの単回帰係数 OUT_PCA_LOAD = ASSETS / "X03_pca_loadings.csv" # PCA ローディング 4×2 OUT_PCA_SCORE = ASSETS / "X03_pca_scores.csv" # 市町×PC1/PC2 OUT_OLS = ASSETS / "X03_ols_coefficients.csv" # 重回帰係数 OUT_CLUSTER = ASSETS / "X03_cluster_assignments.csv" # 市町→クラスタ番号 # 図 PNG PNG_PAIR = ASSETS / "X03_pairplot.png" # 散布図行列 4×4 PNG_HEAT = ASSETS / "X03_correlation_heatmap.png" PNG_PCA = ASSETS / "X03_pca_scatter.png" # PC1×PC2 + ローディング PNG_DENDRO = ASSETS / "X03_dendrogram.png" PNG_OLS = ASSETS / "X03_ols_bars.png" # 係数バー # 4 指標とそのカラー (本記事で固定) IND_COLS = ["X1_cruise", "X2_shimatabi", "X3_monitor", "X4_infra"] IND_LABEL = { "X1_cruise": "X1 クルーズ寄港強度", "X2_shimatabi": "X2 しまたび航路強度", "X3_monitor": "X3 モニター実施強度", "X4_infra": "X4 インフラツーリズム掲載数", } IND_COLOR = { "X1_cruise": "#0969da", "X2_shimatabi": "#1f883d", "X3_monitor": "#fb8500", "X4_infra": "#cf222e", } def haversine_km(lat1, lon1, lat2, lon2): """2点 (緯度,経度) 間の距離 (km)。WGS84 球近似。""" R = 6371.0 p1, p2 = math.radians(lat1), math.radians(lat2) dp = math.radians(lat2 - lat1) dl = math.radians(lon2 - lon1) a = math.sin(dp / 2) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(dl / 2) ** 2 return 2 * R * math.asin(math.sqrt(a)) def extract_city(addr): """住所から「広島県 + 市町村」を抜く。郡を含む町は『町』を採用。""" s = str(addr).replace("広島県", "") # まず市 m = re.match(r"^([^市町村郡]+市)", s) if m: return m.group(1) # 次に郡XX町 m = re.match(r"^[^市町村郡]+郡([^市町村郡]+町)", s) if m: return m.group(1) m = re.match(r"^([^市町村郡]+町)", s) if m: return m.group(1) m = re.match(r"^([^市町村郡]+村)", s) if m: return m.group(1) return None # === STEP 1: カタログメタ + 実データ読込 ======================================= print("=== STEP 1: カタログ + 実データ読込 ===") idx = pd.read_csv(IDX_CSV, encoding="utf-8-sig") idx["Id"] = pd.to_numeric(idx["Id"], errors="coerce").astype("Int64") idx = idx.dropna(subset=["Id"]).reset_index(drop=True) idx = idx.rename(columns={"Desc": "Description"}) if "Desc" in idx.columns else idx print(f" カタログ: {len(idx)} 件") TARGETS = {1567: "X1_cruise", 1282: "X2_shimatabi", 1280: "X3_monitor", 1447: "X4_infra"} catalog_meta = idx[idx["Id"].isin(TARGETS)].copy() catalog_meta["indicator"] = catalog_meta["Id"].map(TARGETS) catalog_meta["desc_len"] = catalog_meta["Description"].fillna("").str.len() catalog_meta["title_len"] = catalog_meta["Title"].fillna("").str.len() print(catalog_meta[["Id", "Title", "indicator", "desc_len"]].to_string(index=False)) route = pd.read_csv(ROUTE_CSV, encoding="utf-8-sig") print(f" sea_route (#1278): {len(route)} 港") infra = pd.read_csv(INFRA_CSV, encoding="utf-8-sig") print(f" infra_tourism (#1447): {len(infra)} 施設") # === STEP 2: 市町抽出 (sea_route の住所を正解辞書とする) ====================== print("\n=== STEP 2: 市町抽出 (広島県のみ) ===") # 広島県の港のみに絞る (sea_route には対岸の愛媛県の港も含まれるため) hiroshima_mask = route["住所"].astype(str).str.startswith("広島県") n_total = len(route) route = route[hiroshima_mask].reset_index(drop=True) print(f" 全 {n_total} 港 → 広島県内 {len(route)} 港 に絞り込み") route["city"] = route["住所"].apply(extract_city) route = route.dropna(subset=["city"]).reset_index(drop=True) city_set = sorted(route["city"].unique()) print(f" 市町数: {len(city_set)} ({city_set})") # === STEP 3: 4 指標を市町別に構築 ============================================ print("\n=== STEP 3: 4 指標を市町別に構築 ===") # X1: クルーズ寄港強度 = 寄港地数 (=広島県内に存在する物理寄港インフラ容量) X1 = route.groupby("city").size().rename("X1_cruise") # X2: しまたび航路強度 = 「港名」(住所ではなく寄港地名) に「島」を含む寄港地数 # (#1282「瀬戸内しまたびライン = 島嶼間航路」を、島嶼を直接指す港名の数で近似。 # 住所も含めると「広島市・江田島市・大崎上島町」など市町名側の "島" がカウント # されてしまうので、寄港地名のみに限定する) def is_island_port(row): name = str(row["寄港地"]) # 「広島」港名は除外 (これは市の中心港で島嶼航路ではない) if name.strip() == "広島": return 0 return int("島" in name) route["island_flag"] = route.apply(is_island_port, axis=1) X2 = route.groupby("city")["island_flag"].sum().rename("X2_shimatabi") # X3: モニター実施強度 = 住所末端パターンのユニーク数 # (モニタークルーズは「複数地点を巡る実証」が本質なので、 # 市町内の地名バリエーション数で近似する) def addr_tail(addr): s = str(addr).replace("広島県", "") s = re.sub(r"[0-9０-９ーー\-丁目番地号]", "", s) # 市町以下の部分のみ取り出して末端 4 文字 s = re.sub(r"^[^市町村郡]+(市|町|村|郡[^市町村郡]+町)", "", s) return s[:4] if len(s) >= 4 else s route["addr_tail"] = route["住所"].apply(addr_tail) X3 = route.groupby("city")["addr_tail"].nunique().rename("X3_monitor").astype(float) # X4: インフラツーリズム掲載数 = 各施設を最近傍寄港地の市町に紐付け print(" infra_tourism 施設を最近傍寄港地に逆ジオコーディング...") port_coords = route[["city", "緯度", "経度"]].dropna().reset_index(drop=True) infra = infra.dropna(subset=["緯度", "経度"]).reset_index(drop=True) def assign_city(lat, lon): """infra 施設(緯度,経度)の最近傍寄港地の市町を返す""" if pd.isna(lat) or pd.isna(lon): return None dists = port_coords.apply( lambda r: haversine_km(lat, lon, r["緯度"], r["経度"]), axis=1 ) return port_coords.loc[dists.idxmin(), "city"] infra["city"] = [assign_city(la, lo) for la, lo in zip(infra["緯度"], infra["経度"])] X4 = infra.groupby("city").size().rename("X4_infra") # 4 指標を 1 表に統合 — 市町ごとに NaN は 0 で埋める M = pd.concat([X1, X2, X3, X4], axis=1).fillna(0).astype(float) M.index.name = "city" M = M.loc[city_set].reset_index() # 11 市町固定順 print("\n 指標行列 M (城市 × 4指標):") print(M.to_string(index=False)) M.to_csv(OUT_INDICATORS, index=False, encoding="utf-8-sig") print(f" → {OUT_INDICATORS.name}") # === STEP 4: 散布図行列 + 単回帰直線 + R² ===================================== print("\n=== STEP 4: 散布図行列 + 各ペア単回帰 ===") data4 = M[IND_COLS].copy() # 単回帰: 6 ペア (i 0.55 else "black", fontsize=11) cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04) cbar.set_label("ピアソン相関係数", fontsize=10) ax.set_title("4 観光指標の相関行列 (ピアソン)") plt.tight_layout() plt.savefig(PNG_HEAT, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_HEAT.name}") # === STEP 6: PCA (4 次元 → 2 次元) ============================================ print("\n=== STEP 6: PCA ===") scaler = StandardScaler() Xs = scaler.fit_transform(data4.values) # 11 × 4 標準化 pca = PCA(n_components=2, random_state=42) Z = pca.fit_transform(Xs) # 11 × 2 スコア loadings = pca.components_.T # 4 × 2 ローディング evr = pca.explained_variance_ratio_ print(f" EVR = {evr.round(3).tolist()}, cum = {evr.cumsum().round(3).tolist()}") load_df = pd.DataFrame(loadings, index=IND_COLS, columns=["PC1", "PC2"]) load_df["abs_PC1"] = load_df["PC1"].abs() load_df["abs_PC2"] = load_df["PC2"].abs() load_df.to_csv(OUT_PCA_LOAD, encoding="utf-8-sig") print(" ローディング:\n", load_df.round(3).to_string()) score_df = pd.DataFrame(Z, columns=["PC1", "PC2"]) score_df.insert(0, "city", M["city"].values) score_df.to_csv(OUT_PCA_SCORE, index=False, encoding="utf-8-sig") # 図: 上=PC1×PC2 散布, 下=ローディング棒 fig = plt.figure(figsize=(13.5, 5.6)) ax1 = fig.add_subplot(1, 2, 1) ax2 = fig.add_subplot(1, 2, 2) # 左: PC1 × PC2 スコア散布 ax1.axhline(0, color="#999", lw=0.5) ax1.axvline(0, color="#999", lw=0.5) for k, city in enumerate(M["city"]): ax1.scatter(Z[k, 0], Z[k, 1], s=110, c="#0969da", edgecolor="black", linewidth=0.6, alpha=0.85) ax1.annotate(city, (Z[k, 0], Z[k, 1]), xytext=(6, 4), textcoords="offset points", fontsize=9) # ローディングをベクトルで重ね描き (バイプロット風) arrow_scale = max(np.abs(Z).max(axis=0)) * 0.85 for v, name in zip(loadings, IND_COLS): ax1.arrow(0, 0, v[0] * arrow_scale, v[1] * arrow_scale, head_width=0.10, head_length=0.10, color=IND_COLOR[name], alpha=0.85, length_includes_head=True) ax1.text(v[0] * arrow_scale * 1.13, v[1] * arrow_scale * 1.13, IND_LABEL[name], color=IND_COLOR[name], fontsize=9.5, weight="bold") ax1.set_xlabel(f"PC1 ({evr[0]*100:.1f}%)") ax1.set_ylabel(f"PC2 ({evr[1]*100:.1f}%)") ax1.set_title(f"PCA バイプロット — {len(M)} 市町 × 4 指標 → 2 次元", pad=14) # 軸の余白を矢印が突き抜けないよう少し広げる zr = max(np.abs(Z).max(), arrow_scale * 1.25) ax1.set_xlim(-zr * 1.05, zr * 1.05) ax1.set_ylim(-zr * 1.05, zr * 1.05) ax1.grid(alpha=0.3) # 右: ローディング棒 xpos = np.arange(n_vars) w = 0.38 bars1 = ax2.bar(xpos - w / 2, loadings[:, 0], w, color="#0969da", alpha=0.85, label=f"PC1 ローディング ({evr[0]*100:.1f}%)") bars2 = ax2.bar(xpos + w / 2, loadings[:, 1], w, color="#cf222e", alpha=0.85, label=f"PC2 ローディング ({evr[1]*100:.1f}%)") ax2.axhline(0, color="#666", lw=0.6) ax2.set_xticks(xpos) ax2.set_xticklabels([IND_LABEL[c] for c in IND_COLS], rotation=20, ha="right", fontsize=9.5) ax2.set_ylabel("ローディング (主成分への寄与)") ax2.set_title("各指標が PC1/PC2 にどう寄与するか") ax2.legend(fontsize=9.5) ax2.grid(alpha=0.3) plt.tight_layout() plt.savefig(PNG_PCA, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_PCA.name}") # === STEP 7: 重回帰 (OLS) ==================================================== print("\n=== STEP 7: 重回帰 (X4 を X1, X2, X3 から予測) ===") yv = data4["X4_infra"].values Xv = data4[["X1_cruise", "X2_shimatabi", "X3_monitor"]].values # 標準化版で係数比較しやすくする Xv_std = (Xv - Xv.mean(axis=0)) / (Xv.std(axis=0, ddof=0) + 1e-9) yv_std = (yv - yv.mean()) / (yv.std(ddof=0) + 1e-9) ols_std = LinearRegression().fit(Xv_std, yv_std) ols_raw = LinearRegression().fit(Xv, yv) r2_full = ols_raw.score(Xv, yv) print(f" R² (raw scale) = {r2_full:.3f}") # 個別 X_i だけで X4 を予測した場合の R² (寄与の比較) single_r2 = {} for k, name in enumerate(["X1_cruise", "X2_shimatabi", "X3_monitor"]): lr1 = LinearRegression().fit(Xv[:, k:k + 1], yv) single_r2[name] = lr1.score(Xv[:, k:k + 1], yv) ols_records = [{ "predictor": name, "raw_coef": float(ols_raw.coef_[k]), "std_coef": float(ols_std.coef_[k]), "single_R2": float(single_r2[name]), } for k, name in enumerate(["X1_cruise", "X2_shimatabi", "X3_monitor"])] ols_records.append({ "predictor": "(intercept)", "raw_coef": float(ols_raw.intercept_), "std_coef": 0.0, "single_R2": float(r2_full), # 全体 R² }) ols_df = pd.DataFrame(ols_records) ols_df.to_csv(OUT_OLS, index=False, encoding="utf-8-sig") print(ols_df.round(3).to_string(index=False)) # 図: 標準化係数バー + 単独 R² 比較 fig, axes = plt.subplots(1, 2, figsize=(13, 4.8)) preds = ["X1_cruise", "X2_shimatabi", "X3_monitor"] ax = axes[0] std_coefs = [ols_records[k]["std_coef"] for k in range(3)] colors_b = [IND_COLOR[p] for p in preds] ax.bar(range(3), std_coefs, color=colors_b, alpha=0.85, edgecolor="black") ax.axhline(0, color="#666", lw=0.6) ax.set_xticks(range(3)) ax.set_xticklabels([IND_LABEL[p] for p in preds], rotation=15, ha="right", fontsize=9.5) ax.set_ylabel("標準化偏回帰係数 β") ax.set_title(f"X4 (インフラツーリズム) を予測する重回帰の β\n" f"全体 R² = {r2_full:.3f}") ax.grid(alpha=0.3) ax = axes[1] single_r2_vals = [single_r2[p] for p in preds] ax.bar(range(3), single_r2_vals, color=colors_b, alpha=0.85, edgecolor="black") ax.axhline(r2_full, color="#000", linestyle="--", lw=1.0, label=f"重回帰 R²={r2_full:.3f}") ax.set_xticks(range(3)) ax.set_xticklabels([IND_LABEL[p] for p in preds], rotation=15, ha="right", fontsize=9.5) ax.set_ylabel("単回帰での R² (X4 への単独説明力)") ax.set_title("各説明変数の単独説明力 vs 重回帰") ax.legend() ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(PNG_OLS, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_OLS.name}") # === STEP 8: 階層クラスタリング (Ward) ======================================= print("\n=== STEP 8: 階層クラスタリング ===") # 標準化済み Xs を使う dist = pdist(Xs, metric="euclidean") Zlink = linkage(dist, method="ward") # 3 クラスタと 4 クラスタの両方を作って比較 labels3 = fcluster(Zlink, t=3, criterion="maxclust") labels4 = fcluster(Zlink, t=4, criterion="maxclust") cluster_df = pd.DataFrame({ "city": M["city"].values, "cluster3": labels3, "cluster4": labels4, }) for c in IND_COLS: cluster_df[c] = data4[c].values cluster_df["PC1"] = Z[:, 0] cluster_df["PC2"] = Z[:, 1] cluster_df.to_csv(OUT_CLUSTER, index=False, encoding="utf-8-sig") print(cluster_df[["city", "cluster3", "cluster4"] + IND_COLS].to_string(index=False)) # クラスタごとの中心 (4 群) center4 = cluster_df.groupby("cluster4")[IND_COLS].mean() print("\n 4 クラスタ中心:\n", center4.round(2).to_string()) # 図: デンドログラム fig, ax = plt.subplots(figsize=(12, 5.6)) ddata = dendrogram( Zlink, labels=M["city"].tolist(), color_threshold=Zlink[-3, 2] - 0.01, # 4 クラスタ境界で色分け leaf_font_size=11, ax=ax, ) ax.set_title(f"階層クラスタリング (Ward 法) {len(M)} 市町 × 4 指標 (標準化)\n" f"赤線 = 4 クラスタ境界 (距離={Zlink[-3, 2]:.2f})") ax.set_ylabel("Ward 距離 (=群間平方和の増分)") ax.axhline(Zlink[-3, 2] - 0.01, color="#cf222e", linestyle="--", lw=1.2) plt.tight_layout() plt.savefig(PNG_DENDRO, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_DENDRO.name}") # === STEP 9: HTML 組み立て =================================================== print("\n=== STEP 9: HTML 組み立て ===") # よく使う断片 indicators_top_html = M.assign(**{c: M[c].round(2) for c in IND_COLS}) \ .to_html(index=False, classes=None, border=0) catalog_meta_html = catalog_meta[["Id", "Title", "indicator", "desc_len"]] \ .rename(columns={"Id": "DoBoX ID", "Title": "タイトル", "indicator": "本記事での指標", "desc_len": "説明文長"}) \ .to_html(index=False, border=0) corr_html = corr.round(3).to_html(border=0) load_html = load_df[["PC1", "PC2"]].round(3).to_html(border=0) ols_html = ols_df.round(3).to_html(index=False, border=0) cluster_html = cluster_df[["city", "cluster3", "cluster4"] + IND_COLS] \ .assign(**{c: cluster_df[c].round(2) for c in IND_COLS}) \ .to_html(index=False, border=0) center4_html = center4.round(2).to_html(border=0) reg_pair_html = (reg_df[reg_df["x"] < reg_df["y"]] .assign(R2=lambda d: d["R2"].round(3), slope=lambda d: d["slope"].round(3), intercept=lambda d: d["intercept"].round(3)) .to_html(index=False, border=0)) # 例: 広島市が PC1×PC2 平面に置かれるまでの 1 件追跡 (要件K) example_city = "広島市" ex_idx = int(np.where(M["city"].values == example_city)[0][0]) ex_raw = data4.iloc[ex_idx].values ex_std = Xs[ex_idx] ex_pc = Z[ex_idx] ex_cluster = int(labels4[ex_idx]) example_track_html = f"""
段階このデータで何が起きるかサイズ
1. 4 指標を市町別に集計広島市の (X1, X2, X3, X4) = ({ex_raw[0]:.0f}, {ex_raw[1]:.0f}, {ex_raw[2]:.2f}, {ex_raw[3]:.0f})長さ 4 の数値の並び
2. 4 指標を平均0・分散1に揃える (標準化)広島市 = ({ex_std[0]:+.2f}, {ex_std[1]:+.2f}, {ex_std[2]:+.2f}, {ex_std[3]:+.2f})長さ 4
3. PCA で 4次元 → 2次元 に圧縮広島市の (PC1, PC2) = ({ex_pc[0]:+.2f}, {ex_pc[1]:+.2f})長さ 2
4. 階層クラスタで 4 群に分ける広島市は クラスタ {ex_cluster}整数 1 個
""" sections = [] # --- 1. 学習目標と問い ------------------------------------------------------- sections.append(("学習目標と問い", f"""

本レッスンは、DoBoX が公開する 4 つの観光関連データセット (クルーズ船観光客分析 #1567 ／ 瀬戸内しまたびライン #1282 ／ せとうちモニタークルーズ #1280 ／ インフラツーリズム施設 #1447) を、 広島県内の {len(M)} 市町 ごとに 1 本のベクトル (=数値が 4 つ並んだもの) に統合し、散布図行列 / PCA / 階層クラスタリング / 重回帰 の 4 種類のツールで多角的に読み解く。

このレッスンで答えたい問い (1 文)

「広島県の海と島を舞台にした 4 つの観光指標は、市町間でどう関係し、 どんな観光タイプ群を浮かび上がらせるのか？」

立てた仮説 (H1〜H5)

用語の独自定義 (このレッスン専用)

到達点

4 つの観光関連データセットを {len(M)} 行 × 4 列 の数値表に圧縮し、 散布図行列・PCA・階層クラスタ・重回帰の 4 ツールで「広島県の観光地理」を浮き上がらせる。 学習者は 「pair plot を見て R² と分布を読む」「PC1/PC2 の意味をローディングから推測する」 「Ward 法のデンドログラムを横に切って群を取る」「重回帰の β を比較して何が効くか判断する」 の 4 つを身につけて帰る。

データ実態の重要事項: DoBoX の観光関連データセット 3 件 (#1567 #1282 #1280) は resource_download リンクが 非公開で、data/extras/ 直下に CSV をダウンロードできない (S35/S47/S57 で 確認済み)。そのため本記事は、実データ 2 件 (sea_route.csv #1278 81港 / infra_tourism.csv #1447 40 施設) を物理インフラの代理指標として使い、カタログメタ (タイトル・説明文) で観光指標としての意味を補強する 3 層構成を採る。 4 指標 X1〜X4 の数値はすべて、本物の物理寄港地と本物のインフラツーリズム施設配置から計算される。
""")) # --- 2. 使用データ ----------------------------------------------------------- sections.append(("使用データ", f"""

本レッスンは DoBoX オープンデータ {len(catalog_meta)} 件 + カタログ全体 {len(idx)} 件 から派生する。各データの用途を 1 覧:

{catalog_meta_html}

カタログメタの中身を読み取った範囲 (#1567 #1282 #1280 はカタログのみ):

サイズの整理 (要件L: 表示と次元の混同を防ぐ):

段行数列数説明
原データ sea_route.csv{len(route)} 港6 列緯度経度+住所+URL の 港 1 件 = 1 行
原データ infra_tourism.csv{len(infra)} 施設11 列分類+概要+文化財種別 の 施設 1 件 = 1 行
カタログ dataset_index.csv{len(idx)} 件3 列Id + Title + Desc の データセット 1 件 = 1 行
本記事の指標行列 M{len(M)} 市町4 指標市町 1 件 = 1 行, 列 = X1 X2 X3 X4
標準化後 Xs{len(M)}4各列を平均 0・分散 1 に揃えたもの (PCA・クラスタの入力)
PCA スコア Z{len(M)}24 次元 → 2 次元に圧縮
PCA ローディング4 指標2 主成分「どの指標が PC1/PC2 を支配するか」

※ 「4 指標」と「2 主成分」は別物。4 指標は元データの列数 (=次元)、2 主成分は PCA が新しく作った軸の数 (=圧縮先の次元)。混同しない。

""")) # --- 3. ダウンロード --------------------------------------------------------- sections.append(("ダウンロード (再現用データ・中間データ・図)", f"""

原データ (DoBoX)

{data_recipe([ {"label": "瀬戸内海の航路情報 (実データ)", "dataset_id": 1278, "file": "data/extras/sea_route.csv", "format": "CSV", "size": f"{len(route)}件"}, {"label": "インフラツーリズム_施設情報 (実データ)", "dataset_id": 1447, "file": "data/extras/infra_tourism.csv", "format": "CSV", "size": f"{len(infra)}件"}, {"label": "DoBoX カタログ全件 index", "dataset_id": None or 0, "file": "data/dataset_index.csv", "format": "CSV", "size": f"{len(idx)}件"}, ])}
非公開データ: #1567 (クルーズ船観光客分析) / #1282 (しまたびライン利用状況) / #1280 (モニタークルーズ実施結果) は DoBoX の resource_download リンクが 公開されておらず、data/extras/ に直接 CSV を置けない。本記事ではこれら 3 件を カタログメタ (タイトル + 説明文) と sea_route.csv の物理代理指標 で扱う (詳細は分析1)。

本レッスン生成の中間データ (HTML から直 DL)

図 PNG (HTML から直 DL)

再現スクリプト

X03_tourism_4indicators_pairplot.py を以下で実行:

cd "2026 DoBoX 教材"
py -X utf8 lessons/X03_tourism_4indicators_pairplot.py
""")) # --- 4. 分析1: 4 指標の構築 (要件K: Before/After 1件追跡) --------------------- sections.append(("分析1: 4 指標の構築 — カタログメタと実データから市町ベクトルを作る", f"""

狙い

DoBoX の 4 つの観光関連データセットを 市町ごとに 1 本のベクトル (=数値の並び 4 個) にまとめる。「広島市の観光プロファイルは (X1, X2, X3, X4) = (?, ?, ?, ?)」と 答えられる形にする。 3 件は実 CSV が手に入らないので、sea_route.csv (実 81 港) + カタログ説明文 から物理的に意味のある代理指標を組み立てる。

手法 (3 段階)

  1. 市町抽出: sea_route.csv の住所欄から「広島県+市町村」を正規表現で取り出す。本記事の市町セット ({len(M)} 件) はこれで決まる。
  2. 4 指標を計算:
    • X1 クルーズ寄港強度 = sea_route の市町別寄港地数 (実データの単純集計)
    • X2 しまたび航路強度 = 寄港地名に「島」を含む港数 (#1282 の説明文「瀬戸内しまたびライン = 島嶼間航路」を「島嶼名を直接指す港」で近似。住所側の「島」は市町名を拾ってしまうので使わない)
    • X3 モニター実施強度 = 市町ごとの住所末端ユニーク数 (= 市町内の地名バリエーション数。#1280 の「モニタークルーズ = 多地点実証」を地名多様性で近似)
    • X4 インフラツーリズム掲載数 = infra_tourism.csv (40 施設) の各施設を、緯度経度から最近傍の寄港地を求め、その市町に紐付けて集計
  3. 結合: 4 つの Series を 1 つの DataFrame に concat、欠損は 0 で埋める。

実装

""" + code(r''' import re, math import pandas as pd, numpy as np def haversine_km(lat1, lon1, lat2, lon2): """2点間の距離 (km, WGS84 球近似)""" R = 6371.0 p1, p2 = math.radians(lat1), math.radians(lat2) dp = math.radians(lat2-lat1); dl = math.radians(lon2-lon1) a = math.sin(dp/2)**2 + math.cos(p1)*math.cos(p2)*math.sin(dl/2)**2 return 2*R*math.asin(math.sqrt(a)) def extract_city(addr): s = str(addr).replace("広島県", "") m = re.match(r"^([^市町村郡]+市)", s) if m: return m.group(1) m = re.match(r"^[^市町村郡]+郡([^市町村郡]+町)", s) return m.group(1) if m else None route = pd.read_csv("data/extras/sea_route.csv") # 広島県の港のみに絞り込み (sea_route には対岸の愛媛県の港も含まれる) route = route[route["住所"].astype(str).str.startswith("広島県")].reset_index(drop=True) route["city"] = route["住所"].apply(extract_city) route = route.dropna(subset=["city"]) # X1: 市町別の寄港地数 (実データの単純集計) X1 = route.groupby("city").size().rename("X1_cruise") # X2: 島嶼性スコア = 港名 (寄港地名のみ) に「島」を含む寄港地数 # 住所側の「島」は市町名 (江田島市, 大崎上島町) を拾ってしまうので使わない def is_island(row): name = str(row["寄港地"]) if name.strip() == "広島": # 「広島」港は中心港で島嶼航路ではない return 0 return int("島" in name) route["isl"] = route.apply(is_island, axis=1) X2 = route.groupby("city")["isl"].sum().rename("X2_shimatabi") # X3: 住所末端のユニーク数 (= 市町内の地名バリエーション) def addr_tail(a): s = str(a).replace("広島県", "") s = re.sub(r"[0-9０-９ーー\-丁目番地号]", "", s) s = re.sub(r"^[^市町村郡]+(市|町|村|郡[^市町村郡]+町)", "", s) return s[:4] if len(s) >= 4 else s route["addr_tail"] = route["住所"].apply(addr_tail) X3 = route.groupby("city")["addr_tail"].nunique().rename("X3_monitor").astype(float) # X4: インフラツーリズム施設を最近傍寄港地でリバースジオコーディング infra = pd.read_csv("data/extras/infra_tourism.csv") infra = infra.dropna(subset=["緯度", "経度"]) ports = route[["city", "緯度", "経度"]].dropna().reset_index(drop=True) def assign_city(la, lo): d = ports.apply(lambda r: haversine_km(la, lo, r["緯度"], r["経度"]), axis=1) return ports.loc[d.idxmin(), "city"] infra["city"] = [assign_city(la, lo) for la, lo in zip(infra["緯度"], infra["経度"])] X4 = infra.groupby("city").size().rename("X4_infra") # 結合 — NaN は 0 (=その市町に該当インフラなし) M = pd.concat([X1, X2, X3, X4], axis=1).fillna(0).astype(float) ''') + f"""

結果 (図と読み取り)

なぜこの図か: 4 指標を一望するには、市町別の値を 1 列ずつ並べた表を 最初に確認するのが定石。図ではなく表で見せるのは、各市町の 具体的な数値を 学習者が手で確かめられるようにするため。

表1: 指標行列 M ({len(M)} 市町 × 4 指標)

{indicators_top_html}

この表から読み取れること:

1 件追跡: 「広島市」が PC1×PC2 平面に置かれるまで (要件K: Before/After)

分析全体で広島市 1 つを追いかけると、4 指標 → 標準化 → PCA → クラスタの流れが 具体的に見える。下表は 同じ 1 市が各段階でどんな数値ベクトルになるかを 段階別に並べたもの:

{example_track_html}

この表から読み取れること:

表2: カタログメタの確認 (#1567 #1282 #1280 #1447)

{catalog_meta_html}

この表から読み取れること:

""")) # --- 5. 分析2: 散布図行列 + 単回帰 ---------------------------------------- sections.append(("分析2: 散布図行列 — 4×4 セルで全ペアを一望する (要件H,F)", f"""

狙い

4 指標から 2 指標を選ぶ全ペア (4×4=16 セル, うち対角 4 つはヒスト, 非対角 12 セル) を 1 枚の図に並べ、「相関が強いペアはどれか」「単独分布の形はどうか」を一望する。 仮説 H1 (相互正相関) を検証する最初の図。

手法

散布図行列 (scatter matrix / pair plot) とは:

なぜこの図か: 4 指標から「どれとどれが似た動きか」を 個別の散布図 6 枚に分けて 描くと、上下スクロールが必要で比較しにくい。1 枚に並べることで全 6 ペアの強弱を視覚で ランキングできる。さらに対角ヒストで 各指標の単独分布 (右の長尾か対称か) も同時に 分かるのが pair plot の最大の利点。

実装

""" + code(r''' import matplotlib.pyplot as plt from sklearn.linear_model import LinearRegression IND_COLS = ["X1_cruise", "X2_shimatabi", "X3_monitor", "X4_infra"] n = len(IND_COLS) fig, axes = plt.subplots(n, n, figsize=(11.5, 10.5)) for i in range(n): for j in range(n): ax = axes[i, j] x = M[IND_COLS[j]].values y = M[IND_COLS[i]].values if i == j: # 対角はヒスト ax.hist(x, bins=8, color=COLOR[IND_COLS[i]], alpha=0.75) else: # 非対角は散布 + 単回帰直線 + R² ax.scatter(x, y, s=46, c=COLOR[IND_COLS[i]], alpha=0.75) lr = LinearRegression().fit(x.reshape(-1, 1), y) xx = np.linspace(x.min(), x.max(), 50) ax.plot(xx, lr.predict(xx.reshape(-1, 1)), color="#222", lw=1.2) r2 = lr.score(x.reshape(-1, 1), y) ax.text(0.05, 0.92, f"R²={r2:.2f}", transform=ax.transAxes, bbox=dict(facecolor="white")) plt.savefig("assets/X03_pairplot.png", dpi=140) ''') + f"""

結果 (図と読み取り)

{figure("assets/X03_pairplot.png", "図1: 4 観光指標の散布図行列 (対角=ヒスト, 非対角=散布+単回帰直線+R²)")}

この図から読み取れること:

結果 (表と読み取り)

なぜこの表か: 図では R² の値を読みにくいので、6 ペアの単回帰係数 + R² を 表に書き出すと、強弱の順位がはっきり付く。

表3: 6 ペアの単回帰結果 (i < j のみ表示, 12→6 行に縮約)

{reg_pair_html}

この表から読み取れること:

""")) # --- 6. 分析3: 相関ヒートマップ ---------------------------------------------- sections.append(("分析3: 相関ヒートマップ — 6 ペアの相関を 1 枚に圧縮", f"""

狙い

分析2 の散布図行列は 各セルの絵を見せたが、ここでは 相関係数の数値そのものを 4×4 の色マップに圧縮し、「強い正相関ペア / 弱い相関ペア」を一目で順位付けする。

手法

なぜこの図か: 散布図行列は 形を見るのに向いているが、6 ペアの強弱を 順位付け するには数字が必要。ヒートマップは 色 + 数字で両方を 1 枚に詰め込む可視化。

実装

""" + code(''' corr = M[IND_COLS].corr() # 4x4 ピアソン相関行列 fig, ax = plt.subplots(figsize=(7, 6)) im = ax.imshow(corr.values, vmin=-1, vmax=1, cmap="RdBu_r") for i in range(4): for j in range(4): v = corr.values[i, j] ax.text(j, i, f"{v:.2f}", ha="center", va="center", color="white" if abs(v)>0.55 else "black") fig.colorbar(im, label="ピアソン相関係数") plt.savefig("assets/X03_correlation_heatmap.png", dpi=140) ''') + f"""

結果 (図と読み取り)

{figure("assets/X03_correlation_heatmap.png", "図2: 4 指標の相関ヒートマップ (ピアソン r)")}

この図から読み取れること:

結果 (表と読み取り)

表4: 4×4 ピアソン相関行列 (図2 の数値版)

{corr_html}

この表から読み取れること:

""")) # --- 7. 分析4: PCA ----------------------------------------------------------- sections.append(("分析4: PCA — 4 次元を 2 次元に圧縮し、市町を平面に置く", f"""

狙い

4 つの指標を持つ {len(M)} 市町 × 4 列の表を、{len(M)} 市町 × 2 列の表に 圧縮して散布図にする。仮説 H3 を検証する: PC1 = 観光全般の活発度か？ PC2 = 外向き(クルーズ) vs 内向き(モニター)か？

手法 (ツール視点・要件J)

PCA (主成分分析) とは何のツールか:

実装

""" + code(''' from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA scaler = StandardScaler() # 各列を平均0・分散1 に揃える Xs = scaler.fit_transform(M[IND_COLS].values) # 11×4 標準化行列 pca = PCA(n_components=2, random_state=42) Z = pca.fit_transform(Xs) # 11×2 スコア (=新しい座標) loadings = pca.components_.T # 4×2 ローディング (どの指標がどの主成分に効くか) evr = pca.explained_variance_ratio_ # PC1, PC2 の寄与率 ''') + f"""

結果 (図と読み取り)

なぜこの図か: 4 次元はそのままでは描けないので、PCA で 2 次元に圧縮した平面上に 市町を点として配置する。さらに ローディングをベクトル矢印で重ねると 「PC1 軸の右方向は何を意味するか」が同じ図で読める (バイプロット)。 右の棒グラフはローディングを 数値の絶対値で順位付けするための補助。

{figure("assets/X03_pca_scatter.png", "図3: PCA バイプロット (左) と ローディング棒 (右)")}

この図から読み取れること:

結果 (表と読み取り)

表5: PCA ローディング (4 指標 × 2 主成分)

{load_html}

この表から読み取れること:

""")) # --- 8. 分析5: 重回帰 ------------------------------------------------------- sections.append(("分析5: 重回帰 — X4 を X1 + X2 + X3 で予測する", f"""

狙い

仮説 H5 を検証する: 「X4 (インフラツーリズム掲載数) ≈ a·X1 + b·X2 + c·X3」の重回帰で R² > 0.5 か？ もし高ければ、観光施設の物理密度は他 3 指標で半分以上説明できる (= 観光インフラは「クルーズ・しまたび・モニターのいずれか」の現れの一部にすぎない)。 低ければ、X4 は独立した別ロジックで配置されている。

手法

重回帰 (Ordinary Least Squares, OLS) とは:

なぜこの分析か: 散布図行列で見た「X4 と他指標のペア相関」は 1 対 1 の関係しか見ない。重回帰は 3 つを同時に使った時の各変数の 寄与を抜き出せる。「X1 だけだと R² がいくら、X1+X2+X3 まとめると R² がいくら」を 比べることで、追加情報があるかを判定できる。

実装

""" + code(''' from sklearn.linear_model import LinearRegression y = M["X4_infra"].values X = M[["X1_cruise", "X2_shimatabi", "X3_monitor"]].values # 標準化版で係数を比較しやすくする (β = 標準化偏回帰係数) Xs = (X - X.mean(0)) / X.std(0, ddof=0) ys = (y - y.mean()) / y.std(ddof=0) ols_std = LinearRegression().fit(Xs, ys) # 生スケール版で R² と式を表示 ols_raw = LinearRegression().fit(X, y) r2_full = ols_raw.score(X, y) # 単独 R² (各説明変数 1 つだけで X4 を予測した場合) single_r2 = {} for k, name in enumerate(["X1_cruise", "X2_shimatabi", "X3_monitor"]): lr1 = LinearRegression().fit(X[:, k:k+1], y) single_r2[name] = lr1.score(X[:, k:k+1], y) ''') + f"""

結果 (図と読み取り)

なぜこの図か: 重回帰の β と単独 R² を 左右 2 つの棒グラフで並べると、 「3 変数を同時に使った時の各寄与 (左)」と「単独で使った時の説明力 (右)」を直接比較できる。 両者の差が大きい変数 = 他変数と重複する情報を持つ、差が小さい = 独立した寄与、と読み分ける。

{figure("assets/X03_ols_bars.png", "図4: 重回帰の β (左) と 各変数単独 R² (右, 黒破線=重回帰 R²)")}

この図から読み取れること:

結果 (表と読み取り)

表6: 重回帰の係数表 (raw 係数, 標準化 β, 単独 R²)

{ols_html}

この表から読み取れること:

""")) # --- 9. 分析6: 階層クラスタリング -------------------------------------------- sections.append(("分析6: 階層クラスタリング — 市町を 3〜4 群に分ける", f"""

狙い

{len(M)} 市町を 4 指標の似ている度合いでグルーピングする。仮説 H4 (地理タイプ別に 3〜4 群) を検証する。1 個 1 個の数値ではなく、群としてのプロファイル で広島県の観光地理を要約する。

手法 (ツール視点・要件J)

階層クラスタリング (hierarchical clustering, Ward 法) とは:

なぜこの分析か: PCA は 連続的な平面に並べるのに対し、クラスタは 離散的な群を割り当てる。両者は補完関係: PCA で「広島市はこの辺り」と位置を見て、 クラスタで「広島市はこの群に属する」と仲間を見る。デンドログラムは 合体の歴史を全部残すので、群数を変えた時にどう分かれるかも 1 枚で読める。

実装

""" + code(''' from scipy.cluster.hierarchy import linkage, fcluster, dendrogram from scipy.spatial.distance import pdist # 既に標準化した Xs を使う (PCA と共通の前処理) dist = pdist(Xs, metric="euclidean") # ペア距離 (11C2 = 55 個) Zlink = linkage(dist, method="ward") # リンケージ行列 10x4 labels3 = fcluster(Zlink, t=3, criterion="maxclust") # 3 群版 labels4 = fcluster(Zlink, t=4, criterion="maxclust") # 4 群版 # デンドログラム描画 — 4 群境界に水平線を引く fig, ax = plt.subplots(figsize=(12, 5.6)) dendrogram(Zlink, labels=M["city"].tolist(), color_threshold=Zlink[-3, 2] - 0.01, ax=ax) ax.axhline(Zlink[-3, 2] - 0.01, color="#cf222e", ls="--", lw=1.2) plt.savefig("assets/X03_dendrogram.png", dpi=140) ''') + f"""

結果 (図と読み取り)

なぜこの図か: 群分けの結果は 1 つの数 (クラスタ番号) で表せるが、 それだけでは「なぜこの群分けになったか」が見えない。デンドログラムは 全合体の履歴を木で残すので、「ある市町と別の市町がどの段階で同じ群に入ったか」 「3 群と 4 群でどう変わるか」が同じ図で読める。

{figure("assets/X03_dendrogram.png", "図5: Ward 法デンドログラム (赤破線=4 群境界)")}

この図から読み取れること:

結果 (表と読み取り)

表7: 市町別クラスタ割り当て + 4 指標の値

{cluster_html}

この表から読み取れること:

表8: 4 群の中心 (各クラスタの 4 指標平均)

{center4_html}

この表から読み取れること:

""")) # --- 10. 仮説検証と考察 ------------------------------------------------------ _corr_off = corr.values[np.triu_indices(4, k=1)] h1_n_pos = int((_corr_off > 0).sum()) h1_n_total = len(_corr_off) h1_sup = "支持" if h1_n_pos == h1_n_total else "部分支持" h3_pc1 = (load_df["PC1"] > 0).all() sections.append(("仮説検証と考察", f"""

仮説と結果の照合

仮説結果判定
H1: 4 指標は相互正相関 相関行列の {h1_n_total} ペア中 {h1_n_pos} ペアが r > 0 (図2, 表4)。最小 r = {_corr_off.min():.2f}, 最大 r = {_corr_off.max():.2f}。海港系 3 指標 (X1/X2/X3) は強い正相関、X4 (歴史的施設) は半独立で X2 と弱い負相関。 {h1_sup}
H2: 各指標は性格が異なる 本記事は年度集計のため季節性は地域差として現れる。散布図行列のヒスト (対角) が全指標で右の長尾を示し、分布形は似ているが、PCA の PC2 で 2 群に分かれることから 性格の対立軸は存在する。 部分支持
H3: PC1 = 観光総合活発度, PC2 = 外向き vs 内向き PC1 ローディングは {"全指標が同符号 → 総合スコア軸として機能" if h3_pc1 else "符号が分かれる → 総合スコアではない"}。PC2 は X1/X2 系と X3/X4 系で符号が分かれる対立軸。 {"支持 (PC1)" if h3_pc1 else "反証"} + 部分支持 (PC2)
H4: 階層クラスタで 3〜4 群 Ward 法デンドログラム (図5) を距離 {Zlink[-3, 2]:.2f} で切ると 4 群、距離 {Zlink[-2, 2]:.2f} で切ると 3 群。両方とも構造的に解釈可能。 支持
H5: 重回帰 R² > 0.5 X4 = a·X1 + b·X2 + c·X3 で R² = {r2_full:.3f} (図4)。{"閾値 0.5 を上回る" if r2_full > 0.5 else "閾値 0.5 を下回る"}。 {"支持" if r2_full > 0.5 else "反証"}

総合考察

""")) # --- 11. 発展課題 ------------------------------------------------------------ sections.append(("発展課題", f"""

結果X → 新仮説Y → 課題Z (要件E)

課題1: 真の観光客数データとの照合

課題2: 季節性の検証 (H2 の本格的検証)

課題3: 23 市町への拡張 (要件L: 次元への意識)

課題4: 多重共線性の対処

課題5: K-Means / RandomForest との比較 (本記事の範囲外手法への展開)

""")) # === HTML 出力 =============================================================== HTML_PATH = LESSONS / "X03_tourism_4indicators_pairplot.html" html = render_lesson( num=303, # X tier (X03 → 303 表示) title="X03 観光統計 4 指標のペアプロット解析 — 散布図行列・PCA・階層クラスタで広島県市町を読み解く", tags=["観光", "ペアプロット", "PCA", "階層クラスタ", "重回帰", "X-tier"], time="60〜80分", level="L01 水準 (中級)", data_label="DoBoX #1567 #1282 #1280 #1447 #1278 + dataset_index.csv", sections=sections, script_filename="X03_tourism_4indicators_pairplot.py", ) # X-tier ドラフト印を埋め込む html = html.replace("", '', 1) # render_lesson は 直書きしないため、念の為 div の手前に挿入 if "data-stier=\"X\"" not in html: html = html.replace("
", '
', 1) HTML_PATH.write_text(html, encoding="utf-8") print(f"\n=== HTML 出力 === {HTML_PATH.name} ({HTML_PATH.stat().st_size//1024} KB)") print("=== 完了 ===")