complete reader study analysis pipeline

.gitignore

@@ -49,6 +49,7 @@ corpus/experiments/rhythm/reader-study/*
 !corpus/experiments/rhythm/reader-study/app.py
 !corpus/experiments/rhythm/reader-study/index.html
 !corpus/experiments/rhythm/reader-study/test_app.py
+!corpus/experiments/rhythm/reader-study/test_analyze.py
 !corpus/experiments/rhythm/reader-study/analyze.py
 !corpus/experiments/rhythm/reader-study/semantic-review.md
 corpus/experiments/rhythm/reader-study/data/

corpus/experiments/rhythm/reader-study/README.md

@@ -22,6 +22,7 @@
 uv run corpus/experiments/rhythm/reader-study/validate_stimuli.py
 python corpus/experiments/rhythm/reader-study/app.py --host 127.0.0.1 --port 8765
 uv run corpus/experiments/rhythm/reader-study/analyze.py --check
+uv run corpus/experiments/rhythm/reader-study/test_analyze.py
 uv run corpus/experiments/rhythm/reader-study/analyze.py \
   corpus/experiments/rhythm/reader-study/data/responses.jsonl
 ```

corpus/experiments/rhythm/reader-study/analyze.py

@@ -4,6 +4,7 @@
 #   "numpy>=2.0",
 #   "pandas>=2.2",
 #   "scipy>=1.13",
+#   "scikit-learn>=1.5",
 #   "statsmodels>=0.14.4",
 #   "sudachipy>=0.6.8",
 #   "sudachidict-core>=20240409",
@@ -24,6 +25,13 @@ from pathlib import Path
 import numpy as np
 import pandas as pd
 import statsmodels.formula.api as smf
+from sklearn.compose import ColumnTransformer
+from sklearn.linear_model import Ridge
+from sklearn.metrics import root_mean_squared_error
+from sklearn.model_selection import GroupKFold
+from sklearn.pipeline import make_pipeline
+from sklearn.preprocessing import OneHotEncoder, StandardScaler
+from statsmodels.genmod.bayes_mixed_glm import BinomialBayesMixedGLM
 
 HERE = Path(__file__).resolve().parent
 ROOT = HERE.parents[3]
@@ -151,10 +159,88 @@ def fit_rating(data: pd.DataFrame, outcome: str) -> dict:
         "ci95_high": estimate + 1.96 * se,
         "p_value": float(model.pvalues[term]),
         "converged": bool(model.converged),
+        "aic": float(model.aic),
     }
 
 
-def write_report(results: list[dict], audit: pd.DataFrame, data: pd.DataFrame, output: Path) -> None:
+def fit_metric(data: pd.DataFrame, outcome: str, metric: str) -> dict:
+    frame = data.copy()
+    frame["metric_z"] = (frame[metric] - frame[metric].mean()) / frame[metric].std()
+    model = smf.mixedlm(
+        f"{outcome} ~ metric_z",
+        frame,
+        groups=np.ones(len(frame)),
+        vc_formula={"participant": "0 + C(participant_id)", "item": "0 + C(item_id)"},
+        re_formula="0",
+    ).fit(reml=False, method=["lbfgs", "powell", "cg"])
+    estimate, se = float(model.params["metric_z"]), float(model.bse["metric_z"])
+    return {
+        "outcome": outcome,
+        "metric": metric,
+        "estimate_per_sd": estimate,
+        "standard_error": se,
+        "ci95_low": estimate - 1.96 * se,
+        "ci95_high": estimate + 1.96 * se,
+        "p_value": float(model.pvalues["metric_z"]),
+        "aic": float(model.aic),
+        "converged": bool(model.converged),
+    }
+
+
+def cross_validated_metrics(data: pd.DataFrame, outcome: str, metrics: tuple[str, ...]) -> list[dict]:
+    participants = data["participant_id"].nunique()
+    folds = min(10, participants)
+    splitter = GroupKFold(n_splits=folds)
+    rows = []
+    for metric in metrics:
+        errors = []
+        features = data[["participant_id", "item_id", metric]]
+        for train, test in splitter.split(features, data[outcome], groups=data["participant_id"]):
+            transform = ColumnTransformer([
+                ("ids", OneHotEncoder(handle_unknown="ignore"), ["participant_id", "item_id"]),
+                ("metric", StandardScaler(), [metric]),
+            ])
+            model = make_pipeline(transform, Ridge(alpha=1.0))
+            model.fit(features.iloc[train], data[outcome].iloc[train])
+            prediction = model.predict(features.iloc[test])
+            errors.append(root_mean_squared_error(data[outcome].iloc[test], prediction))
+        rows.append({"outcome": outcome, "metric": metric, "grouped_cv_folds": folds, "rmse": float(np.mean(errors))})
+    return rows
+
+
+def fit_comprehension(data: pd.DataFrame) -> dict:
+    if data["comprehension_correct"].nunique() < 2:
+        value = int(data["comprehension_correct"].iloc[0])
+        return {
+            "outcome": "comprehension_correct",
+            "model": "not estimable: outcome has no variation",
+            "uniform_minus_varied_log_odds": None,
+            "posterior_sd": None,
+            "credible95_low": None,
+            "credible95_high": None,
+            "observed_value": value,
+        }
+    model = BinomialBayesMixedGLM.from_formula(
+        "comprehension_correct ~ C(condition, Treatment(reference='varied'))",
+        {"participant": "0 + C(participant_id)", "item": "0 + C(item_id)"},
+        data,
+    ).fit_vb()
+    names = model.model.exog_names
+    term = "C(condition, Treatment(reference='varied'))[T.uniform]"
+    index = names.index(term)
+    estimate, sd = float(model.fe_mean[index]), float(model.fe_sd[index])
+    return {
+        "outcome": "comprehension_correct",
+        "model": "Bayesian logistic mixed model; random intercepts for participant and item",
+        "uniform_minus_varied_log_odds": estimate,
+        "posterior_sd": sd,
+        "credible95_low": estimate - 1.96 * sd,
+        "credible95_high": estimate + 1.96 * sd,
+    }
+
+
+def write_report(results: list[dict], metric_results: list[dict], cv_results: list[dict], comprehension: dict,
+                 audit: pd.DataFrame, data: pd.DataFrame, output: Path) -> None:
     total = len(audit)
     included = int(audit["included"].sum()) if total else 0
     lines = [
@@ -165,12 +251,41 @@ def write_report(results: list[dict], audit: pd.DataFrame, data: pd.DataFrame, o
     ]
     for row in results:
         lines.append(f"| {row['outcome']} | {row['uniform_minus_varied']:.3f} | {row['ci95_low']:.3f}–{row['ci95_high']:.3f} | {row['p_value']:.4f} | {row['converged']} |")
-    lines += ["", "主要評価は monotony。正の値は uniform のほうが単調と評定されたことを示す。", ""]
+    lines += ["", "主要評価は monotony。正の値は uniform のほうが単調と評定されたことを示す。", "",
+              "## 理解度", "",
+              (f"uniform − varied のlog odds: {comprehension['uniform_minus_varied_log_odds']:.3f} "
+               f"（95%信用区間 {comprehension['credible95_low']:.3f}–{comprehension['credible95_high']:.3f}）"
+               if comprehension["uniform_minus_varied_log_odds"] is not None
+               else f"モデル推定不能: 理解度が全件 {comprehension['observed_value']} で変動がない。"), "",
+              "## リズム指標", "", "| 評価 | 指標 | 1 SDあたり係数 | p | AIC | CV RMSE |", "|---|---|---:|---:|---:|---:|"]
+    cv_map = {(row["outcome"], row["metric"]): row["rmse"] for row in cv_results}
+    for row in metric_results:
+        lines.append(f"| {row['outcome']} | {row['metric']} | {row['estimate_per_sd']:.3f} | {row['p_value']:.4f} | {row['aic']:.1f} | {cv_map[(row['outcome'], row['metric'])]:.3f} |")
+    lines.append("")
     (output / "results.md").write_text("\n".join(lines), encoding="utf-8")
-    (output / "model-results.json").write_text(json.dumps(results, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
+    payload = {"condition_models": results, "metric_models": metric_results,
+               "cross_validation": cv_results, "comprehension_model": comprehension}
+    (output / "model-results.json").write_text(json.dumps(payload, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
     audit.to_csv(output / "exclusions.csv", index=False)
     data.to_csv(output / "responses-long.csv", index=False)
 
+    primary = next(row for row in results if row["outcome"] == "monotony")
+    best = min((row for row in cv_results if row["outcome"] == "monotony"), key=lambda row: row["rmse"])
+    metric_model = next(row for row in metric_results if row["outcome"] == "monotony" and row["metric"] == best["metric"])
+    supported = primary["uniform_minus_varied"] > 0 and primary["p_value"] < 0.05
+    metric_supported = metric_model["p_value"] < 0.05
+    decision = ["# リズム検出器の実装判断", "",
+                f"- 主要仮説: {'支持' if supported else '不支持'}",
+                f"- 主要対比 uniform − varied: {primary['uniform_minus_varied']:.3f} "
+                f"（95% CI {primary['ci95_low']:.3f}–{primary['ci95_high']:.3f}, p={primary['p_value']:.4f}）",
+                f"- 単調さを最もよく予測した指標: {best['metric']}（CV RMSE={best['rmse']:.3f}, p={metric_model['p_value']:.4f}）", ""]
+    if supported and metric_supported:
+        decision += [f"判断: `{best['metric']}` を単調さの疑いの根拠として残し、他のリズム指標は削除または探索扱いにする。"]
+    else:
+        decision += ["判断: 文長系列指標を単調さ・自然さの価値判断には使わない。既存のリズム警告から価値判断を外す。"]
+    decision += ["", "このファイルは解析結果から自動生成された判断案である。実装変更時に、効果量、区間、逸脱、パイロット所見も確認する。", ""]
+    (output / "implementation-decision.md").write_text("\n".join(decision), encoding="utf-8")
+
 
 def main() -> None:
     parser = argparse.ArgumentParser()
@@ -189,8 +304,13 @@ def main() -> None:
     if audit.empty or int(audit["included"].sum()) < 2:
         raise SystemExit("解析には有効回答が2人以上必要です")
     args.output.mkdir(parents=True, exist_ok=True)
-    results = [fit_rating(data, outcome) for outcome in ("monotony", "naturalness", "readability")]
-    write_report(results, audit, data, args.output)
+    outcomes = ("monotony", "naturalness", "readability")
+    metrics = ("mora_cv", "adjacent_abs_diff", "rmssd", "lag1_autocorrelation")
+    results = [fit_rating(data, outcome) for outcome in outcomes]
+    metric_results = [fit_metric(data, outcome, metric) for outcome in outcomes for metric in metrics]
+    cv_results = [row for outcome in outcomes for row in cross_validated_metrics(data, outcome, metrics)]
+    comprehension = fit_comprehension(data)
+    write_report(results, metric_results, cv_results, comprehension, audit, data, args.output)
     print(f"{args.output / 'results.md'} を作成しました")
 
 

corpus/experiments/rhythm/reader-study/test_analyze.py

@@ -0,0 +1,54 @@
+# /// script
+# requires-python = ">=3.10"
+# dependencies = [
+#   "numpy>=2.0", "pandas>=2.2", "scipy>=1.13", "scikit-learn>=1.5",
+#   "statsmodels>=0.14.4", "sudachipy>=0.6.8", "sudachidict-core>=20240409",
+# ]
+# ///
+import importlib.util
+import sys
+import unittest
+from pathlib import Path
+
+HERE = Path(__file__).resolve().parent
+spec = importlib.util.spec_from_file_location("reader_analysis", HERE / "analyze.py")
+analysis = importlib.util.module_from_spec(spec)
+sys.modules[spec.name] = analysis
+spec.loader.exec_module(analysis)
+
+
+def record(**overrides):
+    answers = [{
+        "item_id": f"item-{index}", "condition": ("uniform", "varied", "control")[index % 3],
+        "monotony": 3 + index % 2, "naturalness": 4 + index % 2, "readability": 5 + index % 2,
+        "comprehension": 0, "elapsed_ms": 15_000,
+    } for index in range(12)]
+    value = {"participant_id": "participant-0001", "attention_check": 4, "answers": answers}
+    value.update(overrides)
+    return value
+
+
+class AnalyzeTest(unittest.TestCase):
+    def setUp(self):
+        self.key = {f"item-{index}": 0 for index in range(12)}
+
+    def test_valid_record_is_included(self):
+        self.assertEqual(analysis.exclusion_reasons(record(), self.key), [])
+
+    def test_preregistered_exclusions(self):
+        value = record(attention_check=3)
+        for answer in value["answers"]:
+            answer.update(monotony=4, naturalness=4, readability=4, comprehension=1, elapsed_ms=5_000)
+        self.assertEqual(set(analysis.exclusion_reasons(value, self.key)), {
+            "attention_check", "comprehension_below_6", "median_time_below_10s",
+            "straightlining_all_ratings",
+        })
+
+    def test_incomplete_record_stops_other_checks(self):
+        value = record()
+        value["answers"].pop()
+        self.assertEqual(analysis.exclusion_reasons(value, self.key), ["incomplete"])
+
+
+if __name__ == "__main__":
+    unittest.main()