Buckets:
| #!/usr/bin/env python3 | |
| """Pre-specified paired analysis for user_study_v1 (protocol.md §11) on the long-format questionnaire table. | |
| Inputs : --questionnaire questionnaire.csv (data_schema.md §4: participant_id, session_index, condition_code, form, item_id, alert_id, raw_value, ...) | |
| --sessions sessions.csv (data_schema.md §2; needs participant_id, session_index, condition_code and, for ES7, alerts_per_hour) | |
| --participants participants.csv (optional; order_code for the §11.3 order check) | |
| --key condition_key.txt (optional; one line `X=C` or `X=S`; without it results are reported as X minus Y, i.e. BLINDED) | |
| --resamples 2000 --seed 20260922 --out results.csv | |
| --no-prompt-neutral sensitivity analysis NP-neutral (protocol §11.1): EP0 / ES1 / ES2 of a session with | |
| sessions.csv no_prompt_box_ticked = 1 are set to the neutral 4 instead of missing | |
| --selftest run on SYNTHETIC random data generated in memory (clearly labelled; nothing is written unless --out) | |
| --check-scoring unit test of the scoring rules: all favourable poles -> usefulness = satisfying = +2 (and the | |
| Likert constructs 7), all unfavourable -> -2 / 1; exit 1 on failure (also run by --selftest) | |
| Output : one row per endpoint: n pairs, mean per condition, mean and median paired difference (C - S, or X - Y when blinded), | |
| percentile bootstrap 95 % CI over participants, Wilcoxon signed-rank p (two-sided), verdict above / below / inconclusive, | |
| and the per-order means of the difference (§11.3) with a two-sample bootstrap CI of their difference. | |
| Scoring (questionnaire_zh_en.md §4): Likert reverse 8 - x for APP3, TIM2, TRU3. Van der Laan: raw_value is entered with the LEFT pole = +2 | |
| and the RIGHT pole = -2 (standard coding); the mirrored pairs VDL3, VDL6, VDL8 carry the favourable word on the right, so score = -raw for | |
| them and raw for the others -> +2 is always the favourable judgement. Construct = mean of its items if at least half are answered. | |
| Requires numpy, pandas, scipy. | |
| """ | |
| import argparse, sys | |
| import numpy as np | |
| import pandas as pd | |
| from scipy import stats | |
| REVERSE_LIKERT = {"APP3", "TIM2", "TRU3"} | |
| MIRRORED_VDL = {"VDL3", "VDL6", "VDL8"} | |
| CONSTRUCTS = { | |
| "EP0_appropriateness": ["APP1", "APP2", "APP3"], | |
| "ES1_timeliness": ["TIM1", "TIM2"], | |
| "ES2_annoyance": ["ANN1", "ANN2"], | |
| "ES3_trust": ["TRU1", "TRU2", "TRU3"], | |
| "ES4_vdl_usefulness": ["VDL1", "VDL3", "VDL5", "VDL7", "VDL9"], | |
| "ES5_vdl_satisfying": ["VDL2", "VDL4", "VDL6", "VDL8"], | |
| } | |
| PER_ALERT = {"ES8_share_appropriate": "AL1", "ES8_share_timely": "AL2", "ES8_share_wanted": "AL3"} | |
| def score(item, raw): | |
| v = pd.to_numeric(raw, errors="coerce") | |
| if item in REVERSE_LIKERT: | |
| return 8 - v | |
| if item in MIRRORED_VDL: | |
| return -v | |
| return v | |
| # favourable word side per Van der Laan item in questionnaire_zh_en.md (left pole = +2 at entry): L = favourable on the left, R = on the right (mirrored) | |
| VDL_FAVOURABLE_SIDE = {"VDL1": "L", "VDL2": "L", "VDL3": "R", "VDL4": "L", "VDL5": "L", "VDL6": "R", "VDL7": "L", "VDL8": "R", "VDL9": "L"} | |
| def check_scoring(): | |
| """all favourable poles -> usefulness = satisfying = +2 and every Likert construct 7 (annoyance: all 'disagree' = 1); | |
| all unfavourable -> -2 / 1 (annoyance 7). Returns a list of failure strings (empty = pass).""" | |
| fails = [] | |
| for label, fav in (("favourable", True), ("unfavourable", False)): | |
| rows = [] | |
| for name, items in CONSTRUCTS.items(): | |
| for it in items: | |
| if it.startswith("VDL"): | |
| side = VDL_FAVOURABLE_SIDE[it] | |
| raw = (2 if side == "L" else -2) if fav else (-2 if side == "L" else 2) | |
| elif name == "ES2_annoyance": | |
| raw = 1 if fav else 7 # 'annoying' items: favourable = disagree | |
| else: | |
| agree = 7 if fav else 1 | |
| raw = (8 - agree) if it in REVERSE_LIKERT else agree | |
| rows.append({"participant_id": "T", "session_index": 1, "form": "S", "item_id": it, "raw_value": raw}) | |
| cs = construct_scores(pd.DataFrame(rows)).iloc[0] | |
| want = {"EP0_appropriateness": 7 if fav else 1, "ES1_timeliness": 7 if fav else 1, "ES2_annoyance": 1 if fav else 7, "ES3_trust": 7 if fav else 1, | |
| "ES4_vdl_usefulness": 2 if fav else -2, "ES5_vdl_satisfying": 2 if fav else -2} | |
| for k, v in want.items(): | |
| if abs(float(cs[k]) - v) > 1e-9: | |
| fails.append(f"all-{label}: {k} = {cs[k]} (expected {v})") | |
| return fails | |
| def construct_scores(q): | |
| """per (participant, session) construct means from S-form rows""" | |
| s = q[q["form"] == "S"].copy() | |
| s["scored"] = [score(i, r) for i, r in zip(s["item_id"], s["raw_value"])] | |
| out = {} | |
| for name, items in CONSTRUCTS.items(): | |
| sub = s[s["item_id"].isin(items)] | |
| g = sub.groupby(["participant_id", "session_index"])["scored"].agg(["mean", "count"]) | |
| g.loc[g["count"] < len(items) / 2, "mean"] = np.nan # fewer than half answered -> missing | |
| out[name] = g["mean"] | |
| return pd.DataFrame(out) | |
| def per_alert_shares(q): | |
| a = q[(q["form"] == "A") & (q["item_id"].isin(PER_ALERT.values()))].copy() | |
| a["v"] = pd.to_numeric(a["raw_value"], errors="coerce") | |
| out = {} | |
| for name, item in PER_ALERT.items(): | |
| sub = a[a["item_id"] == item] | |
| out[name] = sub.groupby(["participant_id", "session_index"])["v"].apply(lambda x: (x >= 5).mean() if len(x) else np.nan) | |
| return pd.DataFrame(out) | |
| def daily_use(q): | |
| e = q[(q["form"] == "E") & (q["item_id"].isin(["END2a", "END2b"]))].copy() | |
| e["session_index"] = e["item_id"].map({"END2a": 1, "END2b": 2}) | |
| e["v"] = pd.to_numeric(e["raw_value"], errors="coerce") | |
| return e.groupby(["participant_id", "session_index"])["v"].mean().rename("ES6_daily_use").to_frame() | |
| def paired_table(per_session, sessions, code_pos, code_neg): | |
| """per_session: DataFrame indexed by (participant_id, session_index); returns wide table with columns pos, neg, diff""" | |
| m = per_session.join(sessions.set_index(["participant_id", "session_index"])["condition_code"], how="inner").reset_index() | |
| wide = m.pivot_table(index="participant_id", columns="condition_code", values=per_session.columns.tolist(), aggfunc="first") | |
| out = {} | |
| for ep in per_session.columns: | |
| if (ep, code_pos) not in wide.columns or (ep, code_neg) not in wide.columns: | |
| continue | |
| t = pd.DataFrame({"pos": wide[(ep, code_pos)], "neg": wide[(ep, code_neg)]}).dropna() | |
| t["diff"] = t["pos"] - t["neg"] | |
| out[ep] = t | |
| return out | |
| def boot_mean(x, rng, n_res): | |
| x = np.asarray(x, float) | |
| idx = rng.integers(0, len(x), size=(n_res, len(x))) | |
| return x[idx].mean(axis=1) | |
| def analyse(tables, orders, n_res, seed, label_pos, label_neg): | |
| rows = [] | |
| for ep, t in tables.items(): | |
| n = len(t) | |
| if n < 2: | |
| rows.append({"endpoint": ep, "n_pairs": n, "note": "fewer than 2 pairs"}) | |
| continue | |
| rng = np.random.default_rng(seed) | |
| b = boot_mean(t["diff"].values, rng, n_res) | |
| lo, hi = np.percentile(b, [2.5, 97.5]) | |
| d = t["diff"].values | |
| if np.all(d == 0): | |
| p = 1.0 | |
| else: | |
| try: | |
| p = float(stats.wilcoxon(d, zero_method="wilcox", alternative="two-sided", method="auto").pvalue) | |
| except ValueError: | |
| p = float("nan") | |
| verdict = "above" if lo > 0 else ("below" if hi < 0 else "inconclusive") | |
| row = {"endpoint": ep, "n_pairs": n, f"mean_{label_pos}": t["pos"].mean(), f"mean_{label_neg}": t["neg"].mean(), "mean_diff": d.mean(), "median_diff": float(np.median(d)), | |
| "ci95_lo": lo, "ci95_hi": hi, "wilcoxon_p": p, "verdict": verdict, "resamples": n_res, "seed": seed} | |
| if orders is not None: | |
| o = t.join(orders, how="left") | |
| grp = {k: v["diff"].values for k, v in o.groupby("order_code")} | |
| if len(grp) == 2: | |
| (ka, da), (kb, db) = sorted(grp.items()) | |
| rng2 = np.random.default_rng(seed) | |
| bd = boot_mean(da, rng2, n_res) - boot_mean(db, rng2, n_res) | |
| olo, ohi = np.percentile(bd, [2.5, 97.5]) | |
| row.update({f"mean_diff_order_{ka}": da.mean(), f"n_order_{ka}": len(da), f"mean_diff_order_{kb}": db.mean(), f"n_order_{kb}": len(db), | |
| "order_diff_ci95_lo": olo, "order_diff_ci95_hi": ohi, "order_interaction": "CI excludes 0" if (olo > 0 or ohi < 0) else "inconclusive"}) | |
| rows.append(row) | |
| return pd.DataFrame(rows) | |
| def synthetic(seed=1): | |
| """SYNTHETIC data with an arbitrary built-in shift, only to exercise the code path. NOT study data.""" | |
| rng = np.random.default_rng(seed) | |
| pids = [f"P{i:02d}" for i in range(1, 31)] | |
| sessions, q = [], [] | |
| for i, pid in enumerate(pids): | |
| order = ("X", "Y") if i % 2 == 0 else ("Y", "X") | |
| for si, code in zip((1, 2), order): | |
| sessions.append({"participant_id": pid, "session_index": si, "condition_code": code, "alerts_per_hour": rng.uniform(0, 12)}) | |
| shift = 0.6 if code == "X" else 0.0 | |
| for items in CONSTRUCTS.values(): | |
| for it in items: | |
| base = rng.normal(4.5 + shift, 1.3) | |
| if it in REVERSE_LIKERT: | |
| base = 8 - base | |
| if it.startswith("VDL"): | |
| base = np.clip(rng.normal(0.5 + shift / 2, 1.0), -2, 2) | |
| if it in MIRRORED_VDL: | |
| base = -base | |
| q.append({"participant_id": pid, "session_index": si, "form": "S", "item_id": it, "raw_value": int(round(base))}) | |
| continue | |
| q.append({"participant_id": pid, "session_index": si, "form": "S", "item_id": it, "raw_value": int(np.clip(round(base), 1, 7))}) | |
| for k in range(rng.integers(1, 6)): | |
| for it in ("AL1", "AL2", "AL3"): | |
| q.append({"participant_id": pid, "session_index": si, "form": "A", "item_id": it, "alert_id": f"{pid}_s{si}_{k + 1:03d}", | |
| "raw_value": int(np.clip(round(rng.normal(4.5 + shift, 1.5)), 1, 7))}) | |
| for it, si in (("END2a", 1), ("END2b", 2)): | |
| code = order[si - 1] | |
| q.append({"participant_id": pid, "session_index": "", "form": "E", "item_id": it, "raw_value": int(np.clip(round(rng.normal(4.5 + (0.6 if code == "X" else 0), 1.4)), 1, 7))}) | |
| parts = pd.DataFrame({"participant_id": pids, "order_code": ["XY" if i % 2 == 0 else "YX" for i in range(len(pids))]}) | |
| return pd.DataFrame(q), pd.DataFrame(sessions), parts | |
| def main(argv=None): | |
| ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) | |
| ap.add_argument("--questionnaire"); ap.add_argument("--sessions"); ap.add_argument("--participants") | |
| ap.add_argument("--key", help="condition_key.txt (X=C or X=S); omit to stay blinded") | |
| ap.add_argument("--resamples", type=int, default=2000); ap.add_argument("--seed", type=int, default=20260922) | |
| ap.add_argument("--out"); ap.add_argument("--selftest", action="store_true") | |
| ap.add_argument("--check-scoring", action="store_true", help="unit test of the scoring rules (VdL sign convention); exit 1 on failure") | |
| ap.add_argument("--no-prompt-neutral", action="store_true", help="sensitivity analysis NP-neutral (protocol 11.1)") | |
| a = ap.parse_args(argv) | |
| if a.check_scoring or a.selftest: | |
| f = check_scoring() | |
| print("scoring check (all favourable -> VdL +2, Likert 7; all unfavourable -> -2 / 1): " + ("PASS" if not f else "FAIL " + "; ".join(f))) | |
| if f: | |
| return 1 | |
| if a.check_scoring and not a.selftest: | |
| return 0 | |
| if a.selftest: | |
| print("=" * 100 + "\nSYNTHETIC SELF-TEST DATA (random numbers with a built-in shift) -- NOT STUDY RESULTS. Exercises the code path only.\n" + "=" * 100) | |
| q, sessions, parts = synthetic() | |
| else: | |
| if not (a.questionnaire and a.sessions): | |
| ap.error("--questionnaire and --sessions are required (or --selftest)") | |
| q, sessions = pd.read_csv(a.questionnaire, dtype=str), pd.read_csv(a.sessions) | |
| parts = pd.read_csv(a.participants, dtype=str) if a.participants else None | |
| sessions["session_index"] = sessions["session_index"].astype(int) | |
| q["session_index"] = pd.to_numeric(q["session_index"], errors="coerce").astype("Int64") | |
| q_s = q.dropna(subset=["session_index"]).copy(); q_s["session_index"] = q_s["session_index"].astype(int) | |
| # condition codes -> labels | |
| pos, neg, lp, ln = "X", "Y", "X", "Y" | |
| if a.key: | |
| k = open(a.key).read().strip().upper().replace(" ", "") | |
| if k == "X=C": | |
| lp, ln = "C", "S" | |
| elif k == "X=S": | |
| pos, neg, lp, ln = "Y", "X", "C", "S" | |
| else: | |
| ap.error("key must be X=C or X=S") | |
| print(f"UNBLINDED with {a.key}: differences are C - S") | |
| else: | |
| print("BLINDED: differences are X - Y (condition codes); apply --key last") | |
| per_session = construct_scores(q_s) | |
| if a.no_prompt_neutral and "no_prompt_box_ticked" in sessions.columns: | |
| flag = sessions[pd.to_numeric(sessions["no_prompt_box_ticked"], errors="coerce") == 1].set_index(["participant_id", "session_index"]).index | |
| per_session = per_session.reindex(per_session.index.union(flag)) | |
| for ep in ("EP0_appropriateness", "ES1_timeliness", "ES2_annoyance"): | |
| per_session.loc[flag, ep] = per_session.loc[flag, ep].fillna(4.0) | |
| print(f"SENSITIVITY ANALYSIS NP-neutral: {len(flag)} no-prompt sessions scored 4 on EP0/ES1/ES2 (protocol 11.1); not the primary analysis") | |
| per_session = per_session.join(per_alert_shares(q_s), how="outer") | |
| e = q[q["form"] == "E"] | |
| if len(e): | |
| per_session = per_session.join(daily_use(e), how="outer") | |
| if "alerts_per_hour" in sessions.columns: | |
| per_session = per_session.join(sessions.set_index(["participant_id", "session_index"])["alerts_per_hour"].rename("ES7_alerts_per_hour"), how="outer") | |
| tables = paired_table(per_session, sessions, pos, neg) | |
| orders = parts.set_index("participant_id")["order_code"] if parts is not None and "order_code" in parts.columns else None | |
| res = analyse(tables, orders, a.resamples, a.seed, lp, ln) | |
| pd.set_option("display.width", 250); pd.set_option("display.max_columns", 40) | |
| print(res.round(3).to_string(index=False)) | |
| print("\nverdict rule: 'above'/'below' = 95 % participant-bootstrap CI excludes 0; otherwise 'inconclusive' (never 'equivalent'). p-values accompany the CI, never alone.") | |
| if a.out: | |
| res.to_csv(a.out, index=False); print(f"wrote {a.out}") | |
| return 0 | |
| if __name__ == "__main__": | |
| sys.exit(main()) | |
Xet Storage Details
- Size:
- 15.1 kB
- Xet hash:
- b0dab8dd725810d7416f1827df3895ce2ac1b1e5776178b0aaf2e0e664c61f92
·
Xet efficiently stores files, intelligently splitting them into unique chunks and accelerating uploads and downloads. More info.