6. Empirical analysis¶

# --- Read empirical data --------------------
data_path = os.path.join(base_path, "EpistasisClass", "Sabela_et_al_2025_N10.csv")

df = pd.read_csv(data_path, index_col=False)

L_raw = ep.Landscape.from_dataframe(df)
display(L_raw)

# Correct Batch effects
L, post = ep.correct_batch_effect(
    L_raw, return_posteriors=True,
    seed=12345678, chains=3, iter_warmup=1000, iter_sampling=1000, link="log",
)
print("Cleaned data")
display(L)
#display(L_batchy)

import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator

# Choose a specific binary state
x = 123  # modify as needed

vals_raw = L_raw.values[x, :]
vals_adj = L.values[x, :]

colors = ["#ffc300", "#118ab2"]  # with batch / without batch

fig, axes = plt.subplots(1, 2, figsize=(7, 3.6), sharey=True)

# ---------------------------------------------------
# LEFT PANEL — WITH BATCH
# ---------------------------------------------------
ax = axes[0]

# Panel label A
ax.text(
    -0.15, 1.05, "A",
    transform=ax.transAxes,
    fontsize=16,
    fontweight="bold",
    va="top",
    ha="left"
)

# Scatter points
ax.scatter(
    range(1, L_raw.R + 1),
    vals_raw,
    marker="^",
    s=140,
    facecolor=colors[0],
    edgecolor="black",
    linewidth=2,
    zorder=3
)

# Horizontal mean line
mean_raw = vals_raw.mean()
ax.axhline(mean_raw, color="gray", linestyle="--", linewidth=2, zorder=2)
ax.text(
    4.5, mean_raw + 0.005,
    r"$\hat{F}$",
    fontsize=18,
    va="bottom",
    ha="left",
    color="gray"
)

ax.set_ylabel(r"Observed $\tilde{F}(\mathbf{\mathcal{x}})$", fontsize=14)

# Replicate labels rotated 45º
ax.set_xticks(range(1, L_raw.R + 1))
ax.set_xticklabels([f"$\\tilde{{F}}_{{{i}}}$" for i in range(1, L.R + 1)],rotation=0,fontsize=15)

# Ticks every 0.2 on y-axis
ax.yaxis.set_major_locator(MultipleLocator(0.1))

# y-limits
ax.set_ylim(0.1, 0.4)

# Label (BATCH)
ax.text(
    0.08, 0.98,
    "BATCH",
    transform=ax.transAxes,
    va="top", ha="left",
    fontsize=12,
    fontweight="bold"
)

ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)


# ---------------------------------------------------
# RIGHT PANEL — NO BATCH
# ---------------------------------------------------
ax = axes[1]

# Panel label B
ax.text(
    -0.15, 1.05, "B",
    transform=ax.transAxes,
    fontsize=16,
    fontweight="bold",
    va="top",
    ha="left"
)

ax.scatter(
    range(1, L.R + 1),
    vals_adj,
    marker="^",
    s=140,
    facecolor=colors[1],
    edgecolor="black",
    linewidth=2,
    zorder=3
)

mean_adj = vals_adj.mean()
ax.axhline(mean_adj, color="gray", linestyle="--", linewidth=2, zorder=2)
ax.text(
    4.5, mean_adj + 0.005,
    r"$\hat{F}$",
    fontsize=18,
    va="bottom",
    ha="left",
    color="gray"
)

ax.set_xticks(range(1, L.R + 1))
ax.set_xticklabels([f"$\\tilde{{F}}_{{{i}}}$" for i in range(1, L.R + 1)],rotation=0,fontsize=15)

ax.yaxis.set_major_locator(MultipleLocator(0.1))
ax.set_ylim(0.08, 0.35)

ax.text(
    0.04, 0.98,
    "NO BATCH",
    transform=ax.transAxes,
    va="top", ha="left",
    fontsize=12,
    fontweight="bold"
)

ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)


# --------------------------------------------------

fig.suptitle(f"{L.states[x]}")

# ---------------------------------------------------
plt.tight_layout()
plt.show()

colors = ["#ffc300", "#118ab2"]  # colors for (with batch, without batch)

fig, axes = plt.subplots(1, 2, figsize=(10, 3.6), sharey=True)

# ---------------------------
# Left: with batch effects
# ---------------------------
ax = axes[0]
vals = L_raw.values

ax.text(
    0.02, 1,            
    "BATCH",        
    transform=ax.transAxes,
    va="top", ha="left",
    fontsize=13,
    fontweight="bold"          
)


bp1 = ax.boxplot(
    [vals[:, r] for r in range(L_raw.R)],
    labels=[f"Replic {r+1}" for r in range(L_raw.R)],
    patch_artist=True,
    flierprops=dict(
    marker=".",
    markersize=6,
    markerfacecolor="black",
    markeredgecolor="none",
    markeredgewidth=0
)
)

ax.set_ylim(0,0.8)

ax.yaxis.set_major_locator(MultipleLocator(0.4))

# Color boxes
for patch in bp1["boxes"]:
    patch.set(facecolor=colors[0], alpha=1)
    patch.set(edgecolor="black")   
    patch.set(linewidth=2)

# Style medians
for med in bp1["medians"]:
    med.set(color="black", linewidth=2)

ax.set_ylabel(r"Observed $\tilde{F}(\mathbf{x})$")
ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha="right")

# Make whiskers and caps thicker
for whisker in bp1["whiskers"]:
    whisker.set(linewidth=2)   
for cap in bp1["caps"]:
    cap.set(linewidth=2)

# ---------------------------
# Right: after correction
# ---------------------------
ax = axes[1]
vals2 = L.values

ax.text(
    0.02, 1,
    "NO BATCH",
    transform=ax.transAxes,
    va="top", ha="left",
    fontsize=13,
    fontweight="bold"
)

bp2 = ax.boxplot(
    [vals2[:, r] for r in range(L.R)],
    labels=[f"Replic {r+1}" for r in range(L.R)],
    patch_artist=True,
    flierprops=dict(
    marker=".",
    markersize=6,
    markerfacecolor="black",
    markeredgecolor="none",
    markeredgewidth=0,
)
)

ax.set_ylim(0,0.8)

ax.yaxis.set_major_locator(MultipleLocator(0.4))

# Color boxes
for patch in bp2["boxes"]:
    patch.set(facecolor=colors[1], alpha=1)
    patch.set(edgecolor="black")  
    patch.set(linewidth=2)

for med in bp2["medians"]:
    med.set(color="black", linewidth=2)

ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha="right")

# ---------------------------
# Remove top/right spines to match your style
# ---------------------------

for whisker in bp2["whiskers"]:
    whisker.set(linewidth=2)
for cap in bp2["caps"]:
    cap.set(linewidth=2)

# Layout
plt.tight_layout()
plt.show()

from math import ceil

# Extract posterior draws for a_r
a = np.asarray(post["a"])   # shape: (draws, R)
n_draws, R = a.shape

# Two rows, dynamic number of columns
n_rows = 2
n_cols = ceil(R / n_rows)

fig, axes = plt.subplots(
    n_rows, n_cols,
    figsize=(3.0 * n_cols, 2.5 * n_rows),
)
axes = np.atleast_2d(axes)

# Plot posterior histograms for each a_r
for r in range(R):
    row = r // n_cols
    col = r % n_cols
    ax = axes[row, col]
    ax.hist(a[:, r], bins=40, density=True,color="#a9d6e5")
 
    mean_r = a[:, r].mean()
    ax.axvline(mean_r, color="blue", linestyle="--", linewidth=1.2,label="Mean")
    ax.set_xlabel(rf"$a_{r+1}$")
    ax.set_ylabel(rf"$P(a_{r+1}|\text{{data}})$")

    if(r==0):
        ax.legend(frameon=False,loc="best",fontsize=10)

# Hide unused axes
for idx in range(R, n_rows * n_cols):
    row = idx // n_cols
    col = idx % n_cols
    axes[row, col].axis("off")

fig.tight_layout()
plt.show()

# Extract posterior draws for b_r
b = np.asarray(post["b"])   # shape: (draws, R)
n_draws, R = b.shape

# Two rows, dynamic number of columns
n_rows = 2
n_cols = ceil(R / n_rows)

fig, axes = plt.subplots(
    n_rows, n_cols,
    figsize=(3.0 * n_cols, 2.5 * n_rows),
)
axes = np.atleast_2d(axes)

# Plot posterior histograms for each b_r
for r in range(R):
    row = r // n_cols
    col = r % n_cols
    ax = axes[row, col]

    # Histogram
    ax.hist(b[:, r], bins=40, density=True, color="#ffbf69")

    # Posterior mean line
    mean_r = b[:, r].mean()
    ax.axvline(mean_r, color="red", linestyle="--", linewidth=1.2,label="Mean")

    # Labels
    ax.set_xlabel(rf"$b_{r+1}$")
    ax.set_ylabel(rf"$P(b_{r+1}|\text{{data}})$")

    if(r==0):
        ax.legend(frameon=False,loc="best",fontsize=10)

# Hide unused axes
for idx in range(R, n_rows * n_cols):
    row = idx // n_cols
    col = idx % n_cols
    axes[row, col].axis("off")

fig.tight_layout()
plt.show()

# Extract posterior draws for sigma
sigma = np.asarray(post["sigma"])   # shape: (draws,)

fig, ax = plt.subplots(figsize=(5, 3))

# Histogram
ax.hist(sigma, bins=40, density=True, color="#76c893")

# Vertical posterior mean line
mean_sigma = sigma.mean()
ax.axvline(mean_sigma, color="green", linestyle="--", linewidth=2,label="Mean")

# Axis labels
ax.set_xlabel(r"$\sigma$")
ax.set_ylabel(r"$P(\sigma|\text{data})$")

# Legend
ax.legend(frameon=False,loc="best",fontsize=10)

fig.tight_layout()
plt.show()

df_amp = ep.epistasis_amplitude(
    L,
    B_uncertainty=1000,
    B_null=1000,
    ci_level=0.95,
    multipliers="normal",
    as_dataframe=True,
)

display(df_amp)

fig, axes, data = ep.plot_variance_and_amplitude(
    L,
    B_uncertainty=100,
    B_null=100,
    as_fraction=True,
    show_components=True,
    rng=np.random.default_rng(125),
    return_data=True,
    figsize=(15, 4),
    ci_method_uncertainty = "percentile"  # {"percentile","studentized","bca"}
)

plt.plot(df_amp["Order"],1-df_amp["P-value order"],"o-")
plt.ylabel("Detection score (1 - p-value)")
plt.xlabel("Epistatic order")

delta2 = ep.focal_effect(
    L,
    i=2,
    missing_policy="drop",       # tolerates missing states
    nan_policy="omit",           # ignore NaNs inside alternating sum
    B_uncertainty=1000,         # bootstrap draws for uncertainty bands
    B_null=1000,                # bootstrap draws for null test
    multipliers="rademacher",  # type of wild multipliers
    ci_level=0.95,
    as_dataframe=True,          # return as DataFrame
)

display(delta2)

display(full_df.loc[(full_df["p-null"] < 0.05) & (full_df["Order"] == 1)])

df_sig = ep.filter_significant_interactions(
    full_df,
    min_snr_null=2.0,
    max_p_value_var=None,
    min_prob_sign=0.0,
    max_p_null=0.05,
    # Only order 2 interactions 
    min_order=2,
    max_order=2,
    min_feature_overlap=2,
    feature_col="Loci names",
    p_col_for_filter="p-null",
)

edges_df = ep.epistasis_to_network(
    df_sig,
    order=2,
    agg="mean_abs",              
    feature_names=L.feature_names 
)

display(edges_df.head())

ep.plot_epistasis_network(
    edges_df,
    use_names=True,  
    node_size=2000,
    font_size=10,
)

# One panel per order, with significance threshold and % significant annotated
ep.plot_epistasis_volcano(
    full_df,
    orders=[1,2,3,4,5,6,7,8,9,10],
    mode="by-order",
    alpha=0.05,
    clip=1e-6,
    colormap="custom",
)