from pathlib import Path

import numpy as np
import pandas as pd
import pyreadstat
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
from IPython.display import display

sav_path = Path.cwd() / "data_combined_20.08.2026.sav"
columns = ["imp3", "E31", "E32", "E33", "cntry", "age", "gender", "degree", "PDI1", "IDV"]
df, meta = pyreadstat.read_sav(
    str(sav_path), usecols=columns, apply_value_formats=False, user_missing=True
)

valid = (
    df[["imp3", "E31", "E32", "E33"]].isin([1, 2, 3, 4, 5]).all(axis=1)
    & df["age"].between(16, 100)
    & df["cntry"].isin(range(1, 10))
    & df["gender"].isin([1, 2, 3])
    & df["degree"].isin([1, 2, 3, 4, 5])
)
analysis = df.loc[valid].copy()
# Pool Hong Kong with China in the analysis; retain original codes in df.
analysis["cntry"] = analysis["cntry"].replace({7: 4})
country_labels = meta.variable_value_labels["cntry"].copy()
country_labels[4] = "China (incl. Hong Kong)"
country_labels.pop(7, None)
analysis["honesty_importance"] = analysis["imp3"]
analysis["C_minus_B"] = analysis["E33"] - analysis["E32"]
analysis["region"] = np.where(analysis["cntry"].isin([1, 2, 9]), "Western", "Eastern")
for variable in ["cntry", "gender", "degree"]:
    analysis[variable] = analysis[variable].astype(int).astype("category")

print(f"Common regression sample: {len(analysis):,} of {len(df):,} rows "
      f"({len(df) - len(analysis):,} excluded by validity / complete-case rules).")

Common regression sample: 2,734 of 7,194 rows (4,460 excluded by validity / complete-case rules).

# Shared presentation; each regression is fitted in its own cell below.
def show_regression(result, title, sample=None, predictors=None, country_lines=False):
    if sample is None:
        sample = analysis
    if predictors is None:
        predictors = ["honesty_importance"]

    ci = result.conf_int()
    table = pd.DataFrame({
        "beta": result.params,
        "HC3 SE": result.bse,
        "p-value": result.pvalues,
        "95% CI lower": ci[0],
        "95% CI upper": ci[1],
    })
    table["significance"] = pd.cut(
        table["p-value"], bins=[-np.inf, .001, .01, .05, .10, np.inf],
        labels=["***", "**", "*", ".", ""], right=False
    )
    table = table.rename(index={
        "PDI1": "Power distance (PDI1)", "IDV": "Individualism (IDV)"
    })
    print(title)
    print(f"N = {int(result.nobs):,} | R² = {result.rsquared:.3f} | "
          f"Adjusted R² = {result.rsquared_adj:.3f} | Standard errors: {result.cov_type}")
    display(table.style.format({
        "beta": "{:.4f}", "HC3 SE": "{:.4f}", "p-value": "{:.3g}",
        "95% CI lower": "{:.4f}", "95% CI upper": "{:.4f}"
    }))

    for predictor in predictors:
        if predictor == "honesty_importance":
            levels = np.arange(1, 6)
            ticks = levels
            xlabel = "Honesty importance (1 = least, 5 = very important)"
            xlim = (.7, 5.3)
            jitter_width = .12
        else:
            # Limit the fitted line to the observed score range.
            levels = np.linspace(sample[predictor].min(), sample[predictor].max(), 40)
            ticks = [1, 2, 3, 4]
            xlabel = "Rejection of dishonesty (ten-item mean; 4 = strongest rejection of dishonesty)"
            xlim = (.8, 4.2)
            jitter_width = .025

        rng = np.random.default_rng(42)
        plot_x = sample[predictor].to_numpy() + rng.uniform(-jitter_width, jitter_width, len(sample))
        plot_y = sample["C_minus_B"].to_numpy() + rng.uniform(-.10, .10, len(sample))
        fig, ax = plt.subplots(figsize=(8, 5))
        for region, color in [("Western", "blue"), ("Eastern", "red")]:
            mask = sample["region"].eq(region).to_numpy()
            ax.scatter(plot_x[mask], plot_y[mask], color=color, s=13, alpha=.18,
                       edgecolors="none", label=region)
        mean_predictions = [
            result.predict(sample.assign(**{predictor: level})).mean()
            for level in levels
        ]
        ax.plot(levels, mean_predictions, color="black", linewidth=2.5,
                label="Model mean prediction", zorder=3)
        if country_lines and predictor == "honest_behavior":
            line_styles = ["-", "--", "-.", ":", (0, (5, 1)), (0, (3, 1, 1, 1))]
            for i, country in enumerate(sample["cntry"].cat.categories):
                country_sample = sample.loc[sample["cntry"].eq(country)]
                if country_sample.empty:
                    continue
                country_levels = np.linspace(
                    country_sample[predictor].min(), country_sample[predictor].max(), 30
                )
                predictions = [
                    result.predict(country_sample.assign(**{predictor: level})).mean()
                    for level in country_levels
                ]
                color = "blue" if country in [1, 2, 9] else "red"
                ax.plot(country_levels, predictions, color=color,
                        linestyle=line_styles[i % len(line_styles)], linewidth=1.8,
                        label=country_labels[country])
        ax.axhline(0, color="gray", linestyle="--", linewidth=1)
        ax.set(xlabel=xlabel,
               ylabel="C_minus_B (reporter − student permitting copying)",
               title=title, xticks=ticks, yticks=range(-4, 5),
               xlim=xlim, ylim=(-4.3, 4.3))
        if country_lines and predictor == "honest_behavior":
            ax.legend(bbox_to_anchor=(1.02, 1), loc="upper left", fontsize=8)
        else:
            ax.legend()
        fig.tight_layout()
        plt.show()
        plt.close(fig)

model_1 = smf.ols(
    "C_minus_B ~ honesty_importance", data=analysis, missing="raise"
).fit(cov_type="HC3")
show_regression(model_1, "Model 1: honesty importance")

Model 1: honesty importance
N = 2,734 | R² = 0.000 | Adjusted R² = -0.000 | Standard errors: HC3

model_2 = smf.ols(
    "C_minus_B ~ honesty_importance + C(cntry)", data=analysis, missing="raise"
).fit(cov_type="HC3")
show_regression(model_2, "Model 2: honesty importance + country fixed effects")

Model 2: honesty importance + country fixed effects
N = 2,734 | R² = 0.211 | Adjusted R² = 0.209 | Standard errors: HC3

model_3 = smf.ols(
    "C_minus_B ~ honesty_importance + C(cntry) + age + C(gender) + C(degree)",
    data=analysis, missing="raise"
).fit(cov_type="HC3")
show_regression(model_3, "Model 3: country fixed effects + demographic controls")

Model 3: country fixed effects + demographic controls
N = 2,734 | R² = 0.218 | Adjusted R² = 0.214 | Standard errors: HC3

e5_items = [f"E5{i}" for i in range(1, 11)]
e5_raw, e5_meta = pyreadstat.read_sav(
    str(sav_path), usecols=e5_items, apply_value_formats=False, user_missing=True
)
# Row indices refer to the same source file loaded above.
e5_sample = e5_raw.loc[analysis.index].copy()
item_names = {
    "E51": "Buying stolen goods",
    "E52": "Keeping money found in the street",
    "E53": "Having an affair while married",
    "E54": "Not reporting damage to a parked vehicle",
    "E55": "Fare evasion",
    "E56": "Tax cheating",
    "E57": "Accepting a bribe",
    "E58": "Fabricating a job application",
    "E59": "Promising a child a reward that will not be given",
    "E510": "Lying about a dangerous wolf to keep a child out of the forest",
}
response_labels = {
    1: "1 Always justified", 2: "2 Sometimes justified",
    3: "3 Rarely justified", 4: "4 Never justified", 5: "5 Don't know"
}

response_counts = pd.DataFrame({
    label: e5_sample.eq(value).sum() for value, label in response_labels.items()
})
response_counts["Missing"] = e5_sample.isna().sum()
response_counts["Invalid code"] = (
    e5_sample.notna() & ~e5_sample.isin([1, 2, 3, 4, 5])
).sum()
response_percent = 100 * response_counts / len(e5_sample)
response_table = response_counts.copy().astype(str)
for column in response_counts:
    response_table[column] = [
        f"{count:,} ({percent:.1f}%)"
        for count, percent in zip(response_counts[column], response_percent[column])
    ]
response_table.insert(0, "Question", pd.Series(item_names))
response_table.index.name = "Item"
print(f"Response table denominator: {len(e5_sample):,} regression-sample respondents per row.")
with pd.option_context("display.max_columns", None, "display.max_colwidth", None):
    display(response_table)

# Raw Cronbach's alpha: k/(k-1) * (1 - sum(item variances)/variance(total)).
e5_valid = e5_sample.where(e5_sample.isin([1, 2, 3, 4]))
e5_complete = e5_valid.dropna()
k = len(e5_items)
if len(e5_complete) < 2:
    raise ValueError("Too few complete respondents to estimate Cronbach's alpha.")
item_variances = e5_complete.var(ddof=1)
total_variance = e5_complete.sum(axis=1).var(ddof=1)
if total_variance <= 0:
    raise ValueError("Cronbach's alpha is undefined: total scores have no variance.")
cronbach_alpha = k / (k - 1) * (1 - item_variances.sum() / total_variance)
alpha_summary = pd.DataFrame({
    "Items": [k],
    "Starting sample N": [len(e5_sample)],
    "Complete ten-item N": [len(e5_complete)],
    "Excluded for missing / don't know / invalid": [len(e5_sample) - len(e5_complete)],
    "Cronbach's alpha": [cronbach_alpha],
})
display(alpha_summary.style.format({"Cronbach's alpha": "{:.3f}"}))
print(f"Cronbach's alpha = {cronbach_alpha:.3f}; N = {len(e5_complete):,} complete respondents.")

Response table denominator: 2,734 regression-sample respondents per row.

Cronbach's alpha = 0.719; N = 2,390 complete respondents.

analysis_4 = analysis.loc[e5_complete.index].copy()
analysis_4["honest_behavior"] = e5_complete.mean(axis=1)
print(f"Model 4 sample: {len(analysis_4):,}; "
      f"{len(analysis) - len(analysis_4):,} fewer respondents than Model 3.")

model_4 = smf.ols(
    "C_minus_B ~ honesty_importance + honest_behavior + C(cntry) + age + C(gender) + C(degree)",
    data=analysis_4, missing="raise"
).fit(cov_type="HC3")
show_regression(
    model_4, "Model 4: both honesty measures + country and demographic controls",
    sample=analysis_4, predictors=["honesty_importance", "honest_behavior"]
)

Model 4 sample: 2,390; 344 fewer respondents than Model 3.
Model 4: both honesty measures + country and demographic controls
N = 2,390 | R² = 0.249 | Adjusted R² = 0.244 | Standard errors: HC3

# Retain the existing index units; require both additional controls for Model 5 only.
analysis_5 = analysis_4.loc[
    np.isfinite(analysis_4[["PDI1", "IDV"]]).all(axis=1)
].copy()
for variable in ["cntry", "gender", "degree"]:
    analysis_5[variable] = analysis_5[variable].cat.remove_unused_categories()

print(f"Model 5 sample: {len(analysis_5):,}; "
      f"{len(analysis_4) - len(analysis_5):,} fewer respondents than Model 4.")
absent_countries = sorted(set(analysis_4["cntry"]) - set(analysis_5["cntry"]))
print("Countries lost relative to Model 4:",
      [country_labels[code] for code in absent_countries] or "None")

model_5 = smf.ols(
    "C_minus_B ~ honesty_importance + honest_behavior + C(cntry) "
    "+ age + C(gender) + C(degree) + PDI1 + IDV",
    data=analysis_5, missing="raise"
).fit(cov_type="HC3")

show_regression(
    model_5, "Model 5: add power distance and individualism",
    sample=analysis_5, predictors=["honesty_importance", "honest_behavior"]
)

Model 5 sample: 2,217; 173 fewer respondents than Model 4.
Countries lost relative to Model 4: ['India']
Model 5: add power distance and individualism
N = 2,217 | R² = 0.273 | Adjusted R² = 0.267 | Standard errors: HC3

import re

comparison_models = {
    "Model 1": model_1, "Model 2": model_2,
    "Model 3": model_3, "Model 4": model_4,
    "Model 5": model_5,
}


def significance_mark(p):
    if p < .001:
        return "***"
    if p < .01:
        return "**"
    if p < .05:
        return "*"
    if p < .10:
        return "·"
    return ""


def readable_term(term):
    names = {
        "Intercept": "Intercept",
        "honesty_importance": "Honesty importance",
        "honest_behavior": "Rejection of dishonesty",
        "age": "Age (years)",
        "PDI1": "Power distance (PDI1)",
        "IDV": "Individualism (IDV)",
    }
    if term in names:
        return names[term]
    match = re.fullmatch(r"C\((cntry|gender|degree)\)\[T\.(\d+)\]", term)
    if match:
        variable, code = match.groups()
        labels = country_labels if variable == "cntry" else meta.variable_value_labels[variable]
        label = labels[int(code)]
        prefix = {"cntry": "", "gender": "Gender: ", "degree": "Education: "}[variable]
        return prefix + label
    return term

# Lead with the research variables, then show all other coefficients.
all_terms = list(dict.fromkeys(
    term for model in comparison_models.values() for term in model.params.index
))
term_order = ["honesty_importance", "honest_behavior", "PDI1", "IDV"] + [
    term for term in all_terms if term not in ["honesty_importance", "honest_behavior", "PDI1", "IDV", "Intercept"]
] + ["Intercept"]
comparison_table = pd.DataFrame(index=term_order, columns=comparison_models, dtype=object)
for name, model in comparison_models.items():
    for term in term_order:
        comparison_table.loc[term, name] = (
            f"{model.params[term]:.3f}{significance_mark(model.pvalues[term])} "
            f"({model.bse[term]:.3f})" if term in model.params else "—"
        )
comparison_table.index = [readable_term(term) for term in term_order]
for name, model in comparison_models.items():
    comparison_table.loc["Observations", name] = f"{int(model.nobs):,}"
    comparison_table.loc["R²", name] = f"{model.rsquared:.3f}"
    comparison_table.loc["Adjusted R²", name] = f"{model.rsquared_adj:.3f}"
comparison_table.index.name = "Predictor / statistic"

display(
    comparison_table.style
    .set_caption("Models 1–5 | beta with significance marks (HC3 standard error)")
    .set_properties(**{"text-align": "center", "padding": "9px 14px"})
    .set_table_styles([
        {"selector": "caption", "props": [("font-size", "17px"), ("font-weight", "bold"), ("text-align", "left")]},
        {"selector": "th", "props": [("text-align", "left"), ("padding", "9px 14px")]},
        {"selector": "thead th", "props": [("background-color", "#edf2f7"), ("color", "#172b4d")]},
    ])
)

from matplotlib.lines import Line2D
import re

figure_dir = Path("presentation_figures")
figure_dir.mkdir(exist_ok=True)
blue, red, ink = "#3569A8", "#C45A53", "#253449"
report_style = {
    "font.family": "DejaVu Sans",
    "font.size": 10,
    "axes.titlesize": 12,
    "axes.titleweight": "semibold",
    "axes.labelsize": 10,
    "axes.labelpad": 8,
    "axes.edgecolor": "#BBC3CC",
    "axes.linewidth": .7,
    "text.color": ink,
    "axes.labelcolor": ink,
    "xtick.color": ink,
    "ytick.color": ink,
    "xtick.labelsize": 9,
    "ytick.labelsize": 9,
    "figure.facecolor": "white",
    "axes.facecolor": "white",
    "pdf.fonttype": 42,
    "ps.fonttype": 42,
}


def save_report_figure(fig, filename):
    # Use the PDF in Overleaf; the PNG is a high-resolution alternative.
    for extension in ["pdf", "png"]:
        path = figure_dir / f"{filename}.{extension}"
        fig.savefig(path, dpi=300, bbox_inches="tight", pad_inches=.08,
                    facecolor="white")


with plt.rc_context(report_style):
    # Figure 1: raw responses and the Model 5 adjusted prediction.
    fig, ax = plt.subplots(figsize=(7.6, 4.6), layout="constrained")
    rng = np.random.default_rng(42)
    x = analysis_5["honest_behavior"].to_numpy()
    y = analysis_5["C_minus_B"].to_numpy()
    jitter_x = x + rng.uniform(-.025, .025, len(x))
    jitter_y = y + rng.uniform(-.09, .09, len(y))
    order = rng.permutation(len(analysis_5))
    colors = np.where(analysis_5["region"].eq("Western"), blue, red)
    ax.scatter(jitter_x[order], jitter_y[order], c=colors[order],
               s=12, alpha=.25, edgecolors="none", rasterized=True, zorder=2)

    levels = np.linspace(x.min(), x.max(), 100)
    predictions = [
        model_5.predict(analysis_5.assign(honest_behavior=level)).mean()
        for level in levels
    ]
    ax.plot(levels, predictions, color=ink, linewidth=2.2, zorder=4)
    ax.axhline(0, color="#8793A1", linestyle=(0, (4, 3)), linewidth=.9, zorder=1)
    ax.set(
        title="Rejection of dishonesty and reporter evaluation",
        xlabel="Rejection of dishonesty (ten-item average, 1–4)",
        ylabel="Relative reporter evaluation\n(C − B)",
        xlim=(.9, 4.1), ylim=(-4.3, 4.3),
        xticks=[1, 2, 3, 4], yticks=[-4, -2, 0, 2, 4],
    )
    ax.set_title(ax.get_title(), loc="center", pad=12)
    ax.grid(axis="y", color="#E7EBEF", linewidth=.6)
    ax.set_axisbelow(True)
    for side in ["top", "right"]:
        ax.spines[side].set_visible(False)
    # In the Model 5 sample, the Western/Eastern groups are Europe/Asia.
    ax.legend(handles=[
        Line2D([0], [0], marker="o", color=blue, linestyle="none",
               markersize=5, label="Europe"),
        Line2D([0], [0], marker="o", color=red, linestyle="none",
               markersize=5, label="Asia"),
        Line2D([0], [0], color=ink, linewidth=2.2,
               label="Adjusted prediction"),
    ], loc="upper left", fontsize=9, frameon=True, facecolor="white",
       edgecolor="none", framealpha=.95, borderpad=.6)
    save_report_figure(fig, "honesty_associations_model5")
    plt.show()
    plt.close(fig)

    # Figure 2: Model 5 country coefficients and pointwise 95% HC3 intervals.
    country_terms = sorted(
        [term for term in model_5.params.index if term.startswith("C(cntry)")],
        key=lambda term: model_5.params[term],
    )
    fig, ax = plt.subplots(figsize=(7.6, 4.2), layout="constrained")
    row_labels = []
    for row, term in enumerate(country_terms):
        code = int(re.search(r"T\.(\d+)", term).group(1))
        beta = model_5.params[term]
        lower, upper = model_5.conf_int().loc[term]
        row_labels.append(country_labels[code])
        ax.errorbar(
            beta, row, xerr=[[beta - lower], [upper - beta]],
            fmt="o", color=blue if code in [1, 2, 9] else red,
            markersize=5.5, capsize=3.5, linewidth=1.6, zorder=3,
        )
    ax.set_yticks(range(len(row_labels)), row_labels)
    ax.set_ylim(len(row_labels) - .5, -.65)
    ax.set_title("Adjusted country differences", pad=12)
    ax.set_xlabel("Difference in relative reporter evaluation\n(reference: Germany)")
    ax.axvline(0, color="#8793A1", linestyle=(0, (4, 3)), linewidth=.9)
    ax.grid(axis="x", color="#E7EBEF", linewidth=.6)
    ax.set_axisbelow(True)
    ax.margins(x=.10)
    ax.tick_params(axis="y", length=0, pad=8)
    for side in ["top", "right", "left"]:
        ax.spines[side].set_visible(False)
    ax.legend(handles=[
        Line2D([0], [0], marker="o", color=blue, linestyle="none",
               markersize=5, label="Europe"),
        Line2D([0], [0], marker="o", color=red, linestyle="none",
               markersize=5, label="Asia"),
    ], loc="upper right", ncol=2, frameon=False, fontsize=9)
    save_report_figure(fig, "country_differences_model5")
    plt.show()
    plt.close(fig)

regional_formula = (
    "C_minus_B ~ honesty_importance + honest_behavior + age + C(gender) + C(degree) + PDI1 + IDV"
)
regional_codes = {"Europe": [1, 2], "Asia": [3, 4, 5, 6, 8]}
regional_samples = {}
regional_models = {}

for region, codes in regional_codes.items():
    sample = analysis_5.loc[analysis_5["cntry"].isin(codes)].copy()
    # Avoid empty-category coefficients after splitting the sample.
    for variable in ["gender", "degree"]:
        sample[variable] = sample[variable].cat.remove_unused_categories()
    regional_samples[region] = sample
    regional_models[region] = smf.ols(
        regional_formula, data=sample, missing="raise"
    ).fit(cov_type="HC3")
    country_counts = sample["cntry"].value_counts().sort_index()
    coverage = ", ".join(
        f"{country_labels[code]}: {int(count):,}"
        for code, count in country_counts.items() if count > 0
    )
    references = ", ".join(
        f"{variable} = {meta.variable_value_labels[variable][sample[variable].cat.categories[0]]}"
        for variable in ["gender", "degree"]
    )
    print(f"{region}: N = {len(sample):,} ({coverage})")
    print(f"Reference groups: {references}")

regional_terms = list(dict.fromkeys(
    term for model in regional_models.values() for term in model.params.index
))
regional_order = ["honesty_importance", "honest_behavior", "PDI1", "IDV", "age"] + [
    term for term in regional_terms
    if term not in ["honesty_importance", "honest_behavior", "PDI1", "IDV", "age", "Intercept"]
] + ["Intercept"]
regional_table = pd.DataFrame(index=regional_order, columns=regional_models, dtype=object)
for region, model in regional_models.items():
    for term in regional_order:
        regional_table.loc[term, region] = (
            f"{model.params[term]:.3f}{significance_mark(model.pvalues[term])} "
            f"({model.bse[term]:.3f})" if term in model.params else "—"
        )
regional_table.index = [readable_term(term) for term in regional_order]
regional_table.index.name = "Predictor / statistic"
for region, model in regional_models.items():
    regional_table.loc["Observations", region] = f"{int(model.nobs):,}"
    regional_table.loc["R²", region] = f"{model.rsquared:.3f}"
    regional_table.loc["Adjusted R²", region] = f"{model.rsquared_adj:.3f}"
    regional_table.loc["Country dummies", region] = "No"

display(
    regional_table.style
    .set_caption("Regional Model 5 regressions | beta with significance marks (HC3 standard error)")
    .set_properties(**{"text-align": "center", "padding": "9px 14px"})
    .set_table_styles([
        {"selector": "caption", "props": [("font-size", "17px"), ("font-weight", "bold"), ("text-align", "left")]},
        {"selector": "th", "props": [("text-align", "left"), ("padding", "9px 14px")]},
        {"selector": "thead th", "props": [("background-color", "#edf2f7"), ("color", "#172b4d")]},
    ])
)
print("*** p < .001; ** p < .01; * p < .05; · p < .10. All standard errors are HC3.")

regional_focus_rows = []
for region, model in regional_models.items():
    for term in ["honesty_importance", "honest_behavior"]:
        lower, upper = model.conf_int().loc[term]
        regional_focus_rows.append({
            "Sample": region,
            "Predictor": readable_term(term),
            "Beta": model.params[term],
            "HC3 SE": model.bse[term],
            "p-value": model.pvalues[term],
            "95% CI lower": lower,
            "95% CI upper": upper,
        })
regional_focus = pd.DataFrame(regional_focus_rows).set_index(["Sample", "Predictor"])
display(regional_focus.style.format({
    "Beta": "{:.3f}", "HC3 SE": "{:.3f}", "p-value": "{:.3g}",
    "95% CI lower": "{:.3f}", "95% CI upper": "{:.3f}"
}))

Europe: N = 626 (Germany: 569, Estonia: 57)
Reference groups: gender = Male, degree = High school
Asia: N = 1,591 (Taiwan: 611, China (incl. Hong Kong): 524, Vietnam: 340, Japan: 116)
Reference groups: gender = Male, degree = High school

*** p < .001; ** p < .01; * p < .05; · p < .10. All standard errors are HC3.