#!/usr/bin/env python3 """Reproducible secondary analysis for the CTW manuscript. The script reads the CC BY 4.0 participant-level dataset deposited by Tanevska, Winkle, and Castellano (2025; https://doi.org/10.5281/zenodo.15648216), reshapes six repeated scenes per participant, and separates between-person trust from within-person trust change. It writes analysis-ready data, tables, figures, and a machine-readable results summary. No simulated or imputed participant data are used. """ from __future__ import annotations import json import math from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from scipy import stats SCRIPT_DIR = Path(__file__).resolve().parent SOURCE_NAME = "prolific_AllVariables_Zenodo.csv" # Works both in this project (script under analysis/) and when the submission # copy is placed at package root. The archived source itself is not redistributed. ROOT = SCRIPT_DIR.parent if (SCRIPT_DIR.parent / "source_data" / SOURCE_NAME).exists() else SCRIPT_DIR SOURCE = ROOT / "source_data" / SOURCE_NAME OUT = ROOT / "analysis" / "outputs" OUT.mkdir(parents=True, exist_ok=True) SEED = 20260826 def zscore(series: pd.Series) -> pd.Series: return (series - series.mean()) / series.std(ddof=1) def cluster_ols(y: np.ndarray, X: pd.DataFrame, clusters: np.ndarray) -> pd.DataFrame: """OLS coefficients with CR1 cluster-robust covariance and cluster t tests.""" x = X.to_numpy(dtype=float) y = np.asarray(y, dtype=float) xtx_inv = np.linalg.pinv(x.T @ x) beta = xtx_inv @ x.T @ y resid = y - x @ beta unique = np.unique(clusters) meat = np.zeros((x.shape[1], x.shape[1])) for g in unique: idx = clusters == g xg = x[idx] eg = resid[idx] score = xg.T @ eg meat += np.outer(score, score) n, k, G = len(y), x.shape[1], len(unique) correction = (G / (G - 1)) * ((n - 1) / (n - k)) vcov = correction * xtx_inv @ meat @ xtx_inv se = np.sqrt(np.diag(vcov)) tval = beta / se pval = 2 * stats.t.sf(np.abs(tval), df=G - 1) crit = stats.t.ppf(0.975, df=G - 1) ssr = np.sum(resid**2) sst = np.sum((y - y.mean()) ** 2) result = pd.DataFrame( { "term": X.columns, "estimate": beta, "cluster_se": se, "t": tval, "df": G - 1, "p": pval, "ci_low": beta - crit * se, "ci_high": beta + crit * se, "n_observations": n, "n_participants": G, "r_squared": 1 - ssr / sst, } ) result.attrs["vcov"] = vcov result.attrs["terms"] = list(X.columns) return result def linear_contrast(model: pd.DataFrame, weights: dict[str, float], label: str) -> dict: terms = model.attrs["terms"] vcov = model.attrs["vcov"] b = model.set_index("term")["estimate"].reindex(terms).to_numpy() c = np.array([weights.get(term, 0.0) for term in terms], dtype=float) estimate = float(c @ b) se = float(np.sqrt(c @ vcov @ c)) df = int(model.df.iloc[0]) tval = estimate / se p = float(2 * stats.t.sf(abs(tval), df=df)) crit = stats.t.ppf(.975, df) return { "contrast": label, "estimate": estimate, "cluster_se": se, "t": tval, "df": df, "p": p, "ci_low": estimate - crit * se, "ci_high": estimate + crit * se, } def holm_adjust(pvals: list[float]) -> list[float]: p = np.asarray(pvals, dtype=float) order = np.argsort(p) adjusted_sorted = np.maximum.accumulate((len(p) - np.arange(len(p))) * p[order]) adjusted_sorted = np.minimum(adjusted_sorted, 1.0) adjusted = np.empty_like(adjusted_sorted) adjusted[order] = adjusted_sorted return adjusted.tolist() def bootstrap_spearman(x: np.ndarray, y: np.ndarray, rng: np.random.Generator, b: int = 5000): rho, p = stats.spearmanr(x, y) n = len(x) vals = [] for _ in range(b): idx = rng.integers(0, n, n) r, _ = stats.spearmanr(x[idx], y[idx]) if np.isfinite(r): vals.append(r) low, high = np.quantile(vals, [0.025, 0.975]) return rho, p, low, high def build_long(df: pd.DataFrame) -> pd.DataFrame: rows = [] id_map = {str(value): f"P{i + 1:03d}" for i, value in enumerate(df["ID"].tolist())} for _, r in df.iterrows(): for code, scenario in [("H", "Highway"), ("S", "Suburbs")]: for scene in (1, 2, 3): rows.append( { "ID": id_map[str(r["ID"])], "scenario": scenario, "scene": scene, "loa": float(r[f"{code}Sc{scene}_LoA"]), "trust": float(r[f"{code}Sc{scene}_Trust"]), "comfort": float(r[f"{code}Sc{scene}_Feel"]), "confidence": float(r[f"{code}Sc{scene}_Conf"]), "high_information": int(r["InformationLevel"] == 1), "scenario_order": "Highway first" if r["ScenarioOrder"] == 1 else "Suburbs first", "first_exposure": int( (r["ScenarioOrder"] == 1 and code == "H") or (r["ScenarioOrder"] == 2 and code == "S") ), "age": float(r["Age"]), "sex": r["Sex"], "gender": r["Gender"], "dbq": float(r["DBQ_AVG"]), "avnars": float(r["AVNARS_AVG"]), "driving_frequency": r["DriveFrequency"], "years_active_driver": r["YearsActiveDriver"], } ) long = pd.DataFrame(rows).sort_values(["ID", "scenario", "scene"]).reset_index(drop=True) for v in ["trust", "comfort", "confidence"]: long[f"{v}_between"] = long.groupby("ID")[v].transform("mean") long[f"{v}_within"] = long[v] - long[f"{v}_between"] person = df.copy() person["ID_public"] = person["ID"].astype(str).map(id_map) person = person.set_index("ID_public") for v, source in [("age_z", "Age"), ("avnars_z", "AVNARS_AVG"), ("dbq_z", "DBQ_AVG")]: mapping = zscore(person[source]).to_dict() long[v] = long["ID"].map(mapping) long["highway"] = (long["scenario"] == "Highway").astype(int) long["scene2"] = (long["scene"] == 2).astype(int) long["scene3"] = (long["scene"] == 3).astype(int) long["trust_within_x_highway"] = long["trust_within"] * long["highway"] long["trust_within_x_high_info"] = long["trust_within"] * long["high_information"] return long def main() -> None: df = pd.read_csv(SOURCE) if df.shape != (206, 153): raise ValueError(f"Unexpected source shape: {df.shape}; expected (206, 153)") key = ["ID", "InformationLevel", "ScenarioOrder", "Age", "DBQ_AVG", "AVNARS_AVG"] if df[key].isna().any().any(): raise ValueError("Missing values found in required participant fields") long = build_long(df) long.to_csv(OUT / "analysis_ready_long.csv", index=False) participants = pd.DataFrame( { "measure": [ "Participants", "Age, mean", "Age, SD", "Age, minimum", "Age, maximum", "Female (sex), n", "Male (sex), n", "High-information interface, n", "Low-information interface, n", "Highway-first order, n", "Suburbs-first order, n", "DBQ average, mean", "DBQ average, SD", "AV-NARS average, mean", "AV-NARS average, SD", ], "value": [ len(df), df.Age.mean(), df.Age.std(ddof=1), df.Age.min(), df.Age.max(), int((df.Sex == "Female").sum()), int((df.Sex == "Male").sum()), int((df.InformationLevel == 1).sum()), int((df.InformationLevel == 2).sum()), int((df.ScenarioOrder == 1).sum()), int((df.ScenarioOrder == 2).sum()), df.DBQ_AVG.mean(), df.DBQ_AVG.std(ddof=1), df.AVNARS_AVG.mean(), df.AVNARS_AVG.std(ddof=1), ], } ) participants.to_csv(OUT / "table_1_participants.csv", index=False) trajectories = ( long.groupby(["scenario", "scene"]) .agg( n=("ID", "size"), loa_mean=("loa", "mean"), loa_sd=("loa", "std"), trust_mean=("trust", "mean"), trust_sd=("trust", "std"), comfort_mean=("comfort", "mean"), comfort_sd=("comfort", "std"), confidence_mean=("confidence", "mean"), confidence_sd=("confidence", "std"), ) .reset_index() ) for v in ["loa", "trust", "comfort", "confidence"]: trajectories[f"{v}_ci_low"] = trajectories[f"{v}_mean"] - stats.t.ppf(.975, trajectories.n - 1) * trajectories[f"{v}_sd"] / np.sqrt(trajectories.n) trajectories[f"{v}_ci_high"] = trajectories[f"{v}_mean"] + stats.t.ppf(.975, trajectories.n - 1) * trajectories[f"{v}_sd"] / np.sqrt(trajectories.n) trajectories.to_csv(OUT / "table_2_scene_trajectories.csv", index=False) base_terms = [ "trust_between", "trust_within", "highway", "scene2", "scene3", "high_information", "first_exposure", "trust_within_x_highway", "trust_within_x_high_info", "age_z", "avnars_z", "dbq_z", ] X1 = long[base_terms].copy() X1.insert(0, "Intercept", 1.0) m1 = cluster_ols(long.loa.to_numpy(), X1, long.ID.to_numpy()) m1.insert(0, "model", "M1: trust and design") extra_terms = ["comfort_between", "comfort_within"] X2 = long[base_terms + extra_terms].copy() X2.insert(0, "Intercept", 1.0) m2 = cluster_ols(long.loa.to_numpy(), X2, long.ID.to_numpy()) m2.insert(0, "model", "M2: plus comfort") contrasts = [] for label, weights in [ ("Within-person trust slope: Suburbs, low information", {"trust_within": 1}), ("Within-person trust slope: Highway, low information", {"trust_within": 1, "trust_within_x_highway": 1}), ("Within-person trust slope: Suburbs, high information", {"trust_within": 1, "trust_within_x_high_info": 1}), ("Within-person trust slope: Highway, high information", {"trust_within": 1, "trust_within_x_highway": 1, "trust_within_x_high_info": 1}), ("Within-person minus between-person trust slope", {"trust_within": 1, "trust_between": -1}), ]: row = linear_contrast(m1, weights, label) row["model"] = "M1: trust and design" contrasts.append(row) pd.DataFrame(contrasts).to_csv(OUT / "table_3b_conditional_slopes.csv", index=False) m1_for_export = m1.copy(); m1_for_export.attrs = {} m2_for_export = m2.copy(); m2_for_export.attrs = {} models = pd.concat([m1_for_export, m2_for_export], ignore_index=True) models.to_csv(OUT / "table_3_cluster_robust_models.csv", index=False) changes = [] rng = np.random.default_rng(SEED) for scenario in ["Highway", "Suburbs"]: wide = long[long.scenario == scenario].pivot(index="ID", columns="scene", values=["loa", "trust"]) for start, end, label in [(1, 2, "unexpected event (1 to 2)"), (2, 3, "resolution (2 to 3)")]: dx = (wide[("trust", end)] - wide[("trust", start)]).to_numpy() dy = (wide[("loa", end)] - wide[("loa", start)]).to_numpy() rho, p, low, high = bootstrap_spearman(dx, dy, rng) changes.append( { "scenario": scenario, "transition": label, "n": len(dx), "mean_trust_change": dx.mean(), "sd_trust_change": dx.std(ddof=1), "mean_loa_change": dy.mean(), "sd_loa_change": dy.std(ddof=1), "spearman_rho_change_alignment": rho, "p_raw": p, "rho_boot_ci_low": low, "rho_boot_ci_high": high, } ) change_df = pd.DataFrame(changes) change_df["p_holm"] = holm_adjust(change_df.p_raw.tolist()) change_df.to_csv(OUT / "table_4_transition_alignment.csv", index=False) # Sensitivity: repeat M1 after standardizing the two focal variables. sensitivity = long.copy() sensitivity["loa_z"] = zscore(sensitivity.loa) sensitivity["trust_between_z"] = zscore(sensitivity.trust_between) sensitivity["trust_within_z"] = sensitivity.trust_within / sensitivity.trust_within.std(ddof=1) sensitivity["trust_within_z_x_highway"] = sensitivity.trust_within_z * sensitivity.highway sensitivity["trust_within_z_x_high_info"] = sensitivity.trust_within_z * sensitivity.high_information sens_terms = [ "trust_between_z", "trust_within_z", "highway", "scene2", "scene3", "high_information", "first_exposure", "trust_within_z_x_highway", "trust_within_z_x_high_info", "age_z", "avnars_z", "dbq_z", ] XS = sensitivity[sens_terms].copy(); XS.insert(0, "Intercept", 1.0) sens = cluster_ols(sensitivity.loa_z.to_numpy(), XS, sensitivity.ID.to_numpy()) sens.insert(0, "model", "Sensitivity: standardized M1") sens.to_csv(OUT / "table_s1_standardized_model.csv", index=False) # Figure 1: scenario trajectories, deliberately using separate panels rather than a dual axis. colors = {"Highway": "#1F4E79", "Suburbs": "#B54A4A"} fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.3), constrained_layout=True) for scenario in ["Highway", "Suburbs"]: part = trajectories[trajectories.scenario == scenario] for ax, metric, label in [(axes[0], "trust", "Trust score (−2 to 4)"), (axes[1], "loa", "Conferred LoA (0 to 3)")]: ax.errorbar(part.scene, part[f"{metric}_mean"], yerr=[part[f"{metric}_mean"] - part[f"{metric}_ci_low"], part[f"{metric}_ci_high"] - part[f"{metric}_mean"]], marker="o", linewidth=1.8, capsize=3, color=colors[scenario], label=scenario) ax.set_xlabel("Scene") ax.set_ylabel(label) ax.set_xticks([1, 2, 3]) ax.axvline(2, color="#777777", linestyle="--", linewidth=.8) ax.grid(axis="y", alpha=.22) axes[0].legend(frameon=False, loc="best") fig.savefig(OUT / "figure_1_scene_trajectories.png", dpi=600, bbox_inches="tight") fig.savefig(OUT / "figure_1_scene_trajectories.pdf", bbox_inches="tight") plt.close(fig) # Figure 2: binned within-person trust deviations and observed LoA means. plot = long.copy() plot["trust_within_bin"] = pd.qcut(plot.trust_within, q=7, duplicates="drop") binned = plot.groupby(["scenario", "trust_within_bin"], observed=True).agg( trust_within_mean=("trust_within", "mean"), loa_mean=("loa", "mean"), n=("ID", "size"), loa_sd=("loa", "std") ).reset_index() binned["loa_se"] = binned.loa_sd / np.sqrt(binned.n) fig, ax = plt.subplots(figsize=(4.6, 3.5), constrained_layout=True) for scenario in ["Highway", "Suburbs"]: part = binned[binned.scenario == scenario] ax.errorbar(part.trust_within_mean, part.loa_mean, yerr=1.96 * part.loa_se, marker="o", linewidth=1.7, capsize=3, color=colors[scenario], label=scenario) ax.axvline(0, color="#777777", linestyle="--", linewidth=.8) ax.set_xlabel("Within-person trust deviation") ax.set_ylabel("Mean conferred LoA") ax.grid(alpha=.22) ax.legend(frameon=False) fig.savefig(OUT / "figure_2_within_person_alignment.png", dpi=600, bbox_inches="tight") fig.savefig(OUT / "figure_2_within_person_alignment.pdf", bbox_inches="tight") plt.close(fig) # Compact machine-readable summary used to populate the manuscript. def term(model: pd.DataFrame, name: str) -> dict: r = model.loc[model.term == name].iloc[0] return {k: float(r[k]) for k in ["estimate", "cluster_se", "p", "ci_low", "ci_high"]} summary = { "source": { "doi": "10.5281/zenodo.15648216", "license": "CC BY 4.0", "participants": int(len(df)), "observations": int(len(long)), }, "sample": { "age_mean": float(df.Age.mean()), "age_sd": float(df.Age.std(ddof=1)), "age_min": int(df.Age.min()), "age_max": int(df.Age.max()), "female_n": int((df.Sex == "Female").sum()), "male_n": int((df.Sex == "Male").sum()), "high_information_n": int((df.InformationLevel == 1).sum()), "low_information_n": int((df.InformationLevel == 2).sum()), }, "model1": { "r_squared": float(m1.r_squared.iloc[0]), "trust_between": term(m1, "trust_between"), "trust_within": term(m1, "trust_within"), "within_x_highway": term(m1, "trust_within_x_highway"), "within_x_high_information": term(m1, "trust_within_x_high_info"), "highway": term(m1, "highway"), "avnars_z": term(m1, "avnars_z"), }, "model2": { "r_squared": float(m2.r_squared.iloc[0]), "trust_between": term(m2, "trust_between"), "trust_within": term(m2, "trust_within"), "within_x_highway": term(m2, "trust_within_x_highway"), "within_x_high_information": term(m2, "trust_within_x_high_info"), }, "conditional_slopes": contrasts, "transitions": change_df.to_dict(orient="records"), } (OUT / "results_summary.json").write_text(json.dumps(summary, indent=2), encoding="utf-8") print(json.dumps(summary, indent=2)) if __name__ == "__main__": main()