<jupyter_start><jupyter_text>Basic data description<jupyter_code>VALUES_REPLACES = {
    "Zle kontrolovaná (sTK nad 160 mmHg)": True,
    "Dobre kontrolovaná": False,
    ">50%": False,
    "40-49%": True,
    "Neprítomná": 0.0,
    "I.\tstupňa": 1.0,
    "II.\tstupňa": 2.0,
    "III.\tstupňa": 3.0,
    "None": False,
}
RENAME_COLS = {"Počet": "i. Ischemická CMP/TIA"}

cont_cols = [
    "Hmotnosť",
    "Výška",
    "Vek",
    "BMI",
    "Tlak krvi systolický",
    "Tlak krvi diastolický",
    "Kreatinin (µmol/l)",
    "Klírens kreatinínu",
    "Trvanie AFib (mesiace)",
    "Počet",  # "i. Ischemická CMP/TIA"
    "i. STE – Infarkt myokardu",
    "ii. NSTE – Infarkt myokardu",
    "b. Diastolická dysfunkcia",
    "c. End-diastolický priemer ľavej komory (mm)",
    "d. Priemer ľavej predsiene v PLAX (parasternal short axis) (mm)",
    "D–dimér [ug/ml]",
    "Fibrinogén [g/l]",
    "CRP [mg/l]",
    "NT-proBNP  [ng/l]",
    "Hs-troponín [ng/l]",
    "Apelin[ng/ml]",
]

bool_cols = [
    "Pohlavie",
    "Diabetes mellitus.1",
    "Fajčenie (> 5 cigariet denne).1",
    "Cievna mozgová príhoda/TIA/systémový tromboembolizmus",
    "Betablokátory.1",
    "Supraventrikulárne extrasystoly",
    "Flutter predsiení",
    "Iná supraventrikulárna arytmia",
    "Komorová tachykardia / komorová fibrilácia",
    "a.\tArteriálna hypertenzia",
    "Arteriálna hypertenzia",
    "b.\tCievna mozgová príhoda/systémový tromboembolizmus",
    "ii. Hemoragická CMP",
    "iii. Systémový tromboembolizmus",
    "c. Pľúcna embólia",
    "d. Hlboká žilová trombóza",
    "e. Periférne artériové ochorenie/aortálny plát",
    "g. Známky srdcového zlyhávania",
    "ii. Ejekčná frakcia ľavej komory",
    "iii. Diastolická dysfunkcia",
    "iv. Hypertrofia ľavej komory",
    "v. Dilatácia ľavej predsiene",
    "i. Stabilná koronárna choroba ",
    "j. Nestabilná angina pectoris",
    "k. Infarkt myokardu ",
    "l. Chronická obštrukčná choroba pľúc",
    "m. Obštrukčné spánkové apnoe",
    "a. Elektrická kardioverzia",
    "b. Farmakologická kardioverzia",
    "c. ACE inhibítor/Blokátor angiotenzínových receptorov (ARB)",
    "d. Spironolakton/Eplerenon",
    "e. Beta blokátor",
    "f. Digoxin",
    "g. Amiodaron",
    "h. Dronedaron",
    "i. Inhibítor protónovej pumpy",
    "j. Antidepresíva / Antipsychotiká",
    "<NAME>",
    "<NAME>",
    "<NAME>",
    "<NAME>",
    "<NAME>",
    "p.\tD<NAME>",
    "<NAME>",
    "<NAME>",
    "<NAME>",
    "<NAME>",
    "f. Interventrikulárne septum/zadná stena > 10mm",
    "Kardiomyopatia",
    "Tyreotoxikóza",
    "Renálna insuficiencia K/DIGO G4-5 (Dialýza/St.p. transplantácii/CrCl < 0,5ml/s)",
    "Ochorenie pečene (Cirhóza/Transaminázy > 3x HLN/Bilirubin > 2x HLN)",
    "Abuzus alkoholu (≥8 „drinkov“ za týždeň)",
    "Tehotenstvo",
    "Kardiochirurgický výkon menej než 3 mesiace pred zaradením do štúdie",
    "Apelin[ng/ml]",
]


def mean_std_str(x):
    return f"{np.nanmean(x):.2f} $\pm$ {np.nanstd(x):.2f}"


def perc_str(x):
    sum_ = np.nansum(x)
    return f"{sum_} ({(sum_ / len(x)) * 100.0:.1f}%)"


## continuous
df_stats = df[cont_cols + ["Typ"]].rename(columns=RENAME_COLS).replace(VALUES_REPLACES)

# compute agg
stats = df_stats.groupby("Typ").agg(mean_std_str).T
# p-values
for col in df_stats:
    if (col == "Typ") or (col == "Trvanie AFib (mesiace)"):
        continue
    if pg.normality(df_stats[col]).iloc[0, 2]:
        pval = pg.ttest(
            df_stats[df_stats["Typ"] == "Bez fibrilácie"][col],
            df_stats[df_stats["Typ"] == "S fibriláciou"][col],
        ).iloc[0, 3]
        stats.loc[col, "p-value"] = pval
        stats.loc[col, "test"] = "t-test"
    else:
        pval = pg.mwu(
            df_stats[df_stats["Typ"] == "Bez fibrilácie"][col],
            df_stats[df_stats["Typ"] == "S fibriláciou"][col],
        ).iloc[0, 2]
        stats.loc[col, "p-value"] = pval
        stats.loc[col, "test"] = "MWU"


## boolean
df_stats = df[bool_cols + ["Typ"]].replace(VALUES_REPLACES)
df_stats = df_stats.replace({"Muž": True, "Žena": False})
# compute agg
bool_stats = df_stats.groupby("Typ").agg(perc_str).T
# p-values
for col in df_stats.columns:
    if col == "Typ":
        continue
    _, _, stats_ = pg.chi2_independence(x=col, y="Typ", data=df_stats)
    pval = stats_["pval"].min()
    bool_stats.loc[col, "p-value"] = pval
    bool_stats.loc[col, "test"] = "chi2"

stats = pd.concat([stats, bool_stats], axis=0)
if SAVE_RESULTS:
    stats.round(3).to_csv(os.path.join(RESULTS_FOLDER, "basic_stats.csv"))
stats.round(3)
df_compare = df[["Pohlavie", "Vek", "Apelin[ng/ml]", "Typ"]]
df_compare["Pohlavie"] = df_compare["Pohlavie"].apply(lambda x: x[0])
df_compare = df_compare.rename(columns={"Apelin[ng/ml]": "Apelin [ng/ml]"})
healthy_controls["Typ"] = "control"
df_compare = pd.concat([df_compare, healthy_controls], axis=0)

cont_cols = ["Vek", "Apelin [ng/ml]"]
bool_cols = ["Pohlavie"]

stats2 = df_compare[cont_cols + ["Typ"]].groupby("Typ").agg(mean_std_str).T
for col in cont_cols:
    pval = pg.anova(data=df_compare, dv=col, between="Typ").iloc[0, 4]
    stats2.loc[col, "p-value"] = pval
    stats2.loc[col, "test"] = "ANOVA"
# stats2.round(4)

_, _, chi_stats = pg.chi2_independence(data=df_compare, x="Pohlavie", y="Typ")
stats2_bool = (
    df_compare[bool_cols + ["Typ"]]
    .replace({"M": 1, "Ž": 0})
    .groupby("Typ")
    .agg(perc_str)
    .T
)
stats2_bool.loc["Pohlavie", "p-value"] = chi_stats["pval"].min()
stats2_bool.loc["Pohlavie", "test"] = "chi2"
stats2 = pd.concat([stats2, stats2_bool], axis=0)
if SAVE_RESULTS:
    stats2.round(4).to_csv(os.path.join(RESULTS_FOLDER, "basic_stats_w_control.csv"))
stats2.round(4)<jupyter_output><empty_output><jupyter_text>Apelin differences between groups* basic stats* normality test (parametric vs. non-parametric ANOVA)* equal variances test (Welch vs. classic ANOVA)* ANOVA* post-hoc* boxplot with plotted significance<jupyter_code>df_apelin = df_compare[["Apelin [ng/ml]", "Typ"]]
df_apelin.groupby("Typ").describe()
# normality - Shapiro-Wilk
df_apelin.loc[df_apelin["Typ"] == "S fibriláciou", "idx"] = np.arange(
    len(df_apelin[df_apelin["Typ"] == "S fibriláciou"])
)
df_apelin.loc[df_apelin["Typ"] == "Bez fibrilácie", "idx"] = np.arange(
    len(df_apelin[df_apelin["Typ"] == "Bez fibrilácie"])
)
df_apelin.loc[df_apelin["Typ"] == "control", "idx"] = np.arange(
    len(df_apelin[df_apelin["Typ"] == "control"])
)
df_apelin = df_apelin.pivot(values="Apelin [ng/ml]", columns="Typ", index="idx")
df_apelin = df_apelin.rename(
    columns={"S fibriláciou": "AFib", "Bez fibrilácie": "no AFib"}
)
pg.normality(df_apelin)
# normality Q-Q
_, axs = plt.subplots(1, 3, sharex=True, sharey=True)
for i, col in enumerate(df_apelin):
    pg.qqplot(df_apelin[col], "norm", ax=axs[i])
    sns.despine(trim=True)
    axs[i].set_title(col)
if SAVE_RESULTS:
    plt.savefig(
        os.path.join(RESULTS_FOLDER, "apelin_qqplot" + PLOTS_EXT),
        dpi=DPI,
        transparent=True,
        bbox_inches="tight",
    )
# equal variance test - Levene
pg.homoscedasticity(df_apelin.dropna()).round(3)
# ANOVA
df_apelin_anova_format = df_apelin.melt()
anova = pg.anova(data=df_apelin_anova_format, dv="value", between="Typ", detailed=False)
anova
# compute power of ANOVA
pg.power_anova(eta=anova["np2"], k=df_apelin.shape[1], n=df_apelin.shape[0])
posthoc = pg.pairwise_tukey(
    data=df_apelin_anova_format, dv="value", between="Typ", effsize="cohen"
)
if SAVE_RESULTS:
    posthoc.round(3).to_csv(os.path.join(RESULTS_FOLDER, "apelin_posthoc.csv"))
posthoc.round(3)
sns.boxplot(
    x="Typ",
    y="value",
    data=df_apelin_anova_format,
    palette=sns.color_palette("Dark2"),
)
plt.ylabel("Apelin[ng/ml]")
plt.xlabel("")
sns.despine(trim=True)

# plot significance from post-hoc
positions = {"AFib": 1, "no AFib": 0, "control": 2}
y = df_apelin_anova_format["value"].max() * 1.02
h = df_apelin_anova_format["value"].max() * 0.05
color = plt.rcParams["text.color"]
for _, row in posthoc[posthoc["p-tukey"] < ALPHA][::-1].iterrows():
    x1 = positions[row["A"]]
    x2 = positions[row["B"]]
    plt.gca().plot([x1, x1, x2, x2], [y, y + h, y + h, y], lw=1.5, c=color)
    text = text = f"p={row['p-tukey']:.3f}"
    plt.text(
        (x1 + x2) * 0.5,
        y + h * 1.1,
        text,
        ha="center",
        va="bottom",
        color=color,
        size=15,
    )
    y = y + h * 3

if SAVE_RESULTS:
    plt.savefig(
        os.path.join(RESULTS_FOLDER, "apelin_boxplot" + PLOTS_EXT),
        dpi=DPI,
        transparent=True,
        bbox_inches="tight",
    )<jupyter_output><empty_output><jupyter_text>Apelin correlations with other risk factors* apelin vs. LA size* apelin vs. NT-proBNP* apelin vs. diastolic dysfunction correlation<jupyter_code>def plot_corrs_as_reg(data, x_col, y_col, hue, rename=None, replace=None):
    df = data.copy()
    rename = rename or {}
    df = df.rename(columns=rename)
    replace = replace or {}
    df = df.replace(replace)
    g = sns.jointplot(data=df, x=x_col, y=y_col, hue=hue, palette="Dark2")
    sns.despine(trim=True)
    for col, (_, gr) in zip(sns.color_palette("Dark2", 2), df.groupby(hue)):
        sns.regplot(
            x=x_col,
            y=y_col,
            data=gr,
            scatter=False,
            ax=g.ax_joint,
            truncate=False,
            color=col,
        )
    if SAVE_RESULTS:
        plt.savefig(
            os.path.join(
                RESULTS_FOLDER, f"corrplot_{x_col[:5]}_vs_{y_col[:5]}" + PLOTS_EXT
            ),
            dpi=DPI,
            transparent=True,
            bbox_inches="tight",
        )


def get_corrs(data, x_col, y_col, hue):
    overall = pg.corr(data[x_col], data[y_col], method="spearman")
    overall.index = ["overall"]
    for gr_name, gr in data.groupby(hue):
        temp = pg.corr(gr[x_col], gr[y_col], method="spearman")
        temp.index = [gr_name]
        overall = pd.concat([overall, temp], axis=0)
    return overall
# LA size
plot_corrs_as_reg(
    data=df,
    x_col="Apelin [ng/ml]",
    y_col="LA size [mm]",
    hue="Type",
    rename={
        "d. <NAME> v PLAX (parasternal short axis) (mm)": "LA size [mm]",
        "Typ": "Type",
        "Apelin[ng/ml]": "Apelin [ng/ml]",
    },
    replace={"Bez fibrilácie": "no-AFib", "S fibriláciou": "AFib"},
)
corrs = get_corrs(
    df,
    x_col="Apelin[ng/ml]",
    y_col="d. Priemer ľavej predsiene v PLAX (parasternal short axis) (mm)",
    hue="Typ",
)
if SAVE_RESULTS:
    corrs.round(3).to_csv(os.path.join(RESULTS_FOLDER, "ap_lasize_corr.csv"))
corrs
# NT-proBNP  [ng/l]
plot_corrs_as_reg(
    data=df,
    x_col="Apelin [ng/ml]",
    y_col="NT-proBNP [ng/l]",
    hue="Type",
    rename={
        "NT-proBNP  [ng/l]": "NT-proBNP [ng/l]",
        "Typ": "Type",
        "Apelin[ng/ml]": "Apelin [ng/ml]",
    },
    replace={"Bez fibrilácie": "no-AFib", "S fibriláciou": "AFib"},
)
corrs = get_corrs(
    df,
    x_col="Apelin[ng/ml]",
    y_col="NT-proBNP  [ng/l]",
    hue="Typ",
)
if SAVE_RESULTS:
    corrs.round(3).to_csv(os.path.join(RESULTS_FOLDER, "ap_ntprobnp_corr.csv"))
corrs
# diastolic dysfunction
df_dys = df.replace(
    {"Neprítomná": 0, "I.\tstupňa": 1, "II.\tstupňa": 2, "III.\tstupňa": 3}
)
plot_corrs_as_reg(
    data=df_dys,
    x_col="Apelin [ng/ml]",
    y_col="Diastolic dysfunction",
    hue="Type",
    rename={
        "b. Diastolická dysfunkcia": "Diastolic dysfunction",
        "Typ": "Type",
        "Apelin[ng/ml]": "Apelin [ng/ml]",
    },
    replace={"Bez fibrilácie": "no-AFib", "S fibriláciou": "AFib"},
)
corrs = get_corrs(
    df_dys,
    x_col="Apelin[ng/ml]",
    y_col="<NAME>",
    hue="Typ",
)
if SAVE_RESULTS:
    corrs.round(3).to_csv(os.path.join(RESULTS_FOLDER, "ap_diadys_corr.csv"))
corrs<jupyter_output><empty_output><jupyter_text>Apelin predicts AF (AF vs non-AF patients)* AUC analysis<jupyter_code># prepare
df_logres = df[["Apelin[ng/ml]", "Typ"]].dropna()
df_logres = df_logres.replace({"S fibriláciou": 0.0, "Bez fibrilácie": 1.0})


def manual_roc_curve(data, true_col, pred_col, thresholds):
    tpr = []
    fpr = []
    for thresh in thresholds:
        tp, fn, fp, tn = confusion_matrix(
            data[true_col], data[pred_col] > thresh
        ).ravel()
        tpr.append(tp / (tp + fn))
        fpr.append(fp / (fp + tn))
    return np.array(tpr), np.array(fpr)


def get_metrics_for_thresh(data, true_col, pred_col, thresh):
    print(f"Threshold: {thresh}")
    cm = confusion_matrix(data[true_col], data[pred_col] > thresh)
    total = sum(sum(cm))
    tp, fn, fp, tn = cm.ravel()

    display(
        pd.DataFrame(
            cm, columns=["predicted FIB", "predicted NO FIB"], index=["FIB", "NO FIB"]
        )
    )
    print("---------")
    accuracy = (tp + tn) / total
    print("Accuracy : ", accuracy)

    sensitivity = tp / (tp + fn)
    print("Sensitivity : ", sensitivity)

    specificity = tn / (tn + fp)
    print("Specificity : ", specificity)
    
    
def get_metrics_for_prediction(true, predicted):
    cm = confusion_matrix(true, predicted)
    total = sum(sum(cm))
    tn, fp, fn, tp = cm.ravel()

    display(
        pd.DataFrame(
            cm, columns=["predicted FIB", "predicted NO FIB"], index=["FIB", "NO FIB"]
        )
    )
    print("---------")
    accuracy = (tp + tn) / total
    print("Accuracy : ", accuracy)

    sensitivity = tp / (tp + fn)
    print("Sensitivity : ", sensitivity)

    specificity = tn / (tn + fp)
    print("Specificity : ", specificity)
    
    
def get_metrics_for_kfold(true, predicted):
    cm = confusion_matrix(true, predicted)
    total = sum(sum(cm))
    tn, fp, fn, tp = cm.ravel()
    accuracy = (tp + tn) / total
    sensitivity = tp / (tp + fn)
    specificity = tn / (tn + fp)
    return accuracy, sensitivity, specificity
    


TRUE_COL = "Typ"
PRED_COL = "Apelin[ng/ml]"
thresholds = np.linspace(df_logres[PRED_COL].min(), df_logres[PRED_COL].max(), 100)
tpr, fpr = manual_roc_curve(
    df_logres,
    true_col=TRUE_COL,
    pred_col=PRED_COL,
    thresholds=thresholds,
)
auc = np.abs(np.trapz(tpr, fpr))
plt.figure(figsize=(7, 7))
plt.plot(fpr, tpr, label=f"AUC: {auc:.4f}")
plt.plot([0, 1], [0, 1], color="#AAAAAA", linestyle="--")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title("Apelin predicts AFib: patients only")
plt.legend()
sns.despine(trim=True)
if SAVE_RESULTS:
    plt.savefig(
        os.path.join(RESULTS_FOLDER, "apelin_roc_patients" + PLOTS_EXT),
        dpi=DPI,
        transparent=True,
        bbox_inches="tight",
    )
ideal_threshold_pat = thresholds[np.argmax(tpr - fpr)]
get_metrics_for_thresh(
    df_logres, true_col=TRUE_COL, pred_col=PRED_COL, thresh=ideal_threshold_pat
)<jupyter_output>Threshold: 0.9685454545454545<jupyter_text>log Apelin predicts AF (AF vs non-AF patients)<jupyter_code>df_logres_log = df[["Apelin[ng/ml]", "Typ"]].dropna()
df_logres_log = df_logres_log.replace({"S fibriláciou": 0.0, "Bez fibrilácie": 1.0})
df_logres_log["Apelin[ng/ml]"] = np.log(df_logres_log["Apelin[ng/ml]"])

thresholds = np.linspace(
    df_logres_log[PRED_COL].min(), df_logres_log[PRED_COL].max(), 100
)
tpr, fpr = manual_roc_curve(
    df_logres_log,
    true_col=TRUE_COL,
    pred_col=PRED_COL,
    thresholds=thresholds,
)
auc = np.abs(np.trapz(tpr, fpr))
plt.figure(figsize=(7, 7))
plt.plot(fpr, tpr, label=f"AUC: {auc:.4f}")
plt.plot([0, 1], [0, 1], color="#AAAAAA", linestyle="--")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title("Log(Apelin) predicts AFib: patients only")
plt.legend()
sns.despine(trim=True)
if SAVE_RESULTS:
    plt.savefig(
        os.path.join(RESULTS_FOLDER, "log_apelin_roc_patients" + PLOTS_EXT),
        dpi=DPI,
        transparent=True,
        bbox_inches="tight",
    )
ideal_threshold_log = thresholds[np.argmax(tpr - fpr)]
get_metrics_for_thresh(
    df_logres_log, true_col=TRUE_COL, pred_col=PRED_COL, thresh=ideal_threshold_log
)<jupyter_output>Threshold: -0.032553544452226024<jupyter_text>Apelin predicts AF (including healthy controls)* AUC analysis<jupyter_code>df_logres_all = pd.DataFrame(columns=["apelin", "AF"])
for col in df_apelin:
    af = 0.0 if col == "AFib" else 1.0
    temp = pd.DataFrame({"apelin": df_apelin[col].dropna().values, "AF": af})
    df_logres_all = pd.concat([df_logres_all, temp], axis=0)

TRUE_COL = "AF"
PRED_COL = "apelin"
thresholds = np.linspace(
    df_logres_all[PRED_COL].min(), df_logres_all[PRED_COL].max(), 100
)
tpr, fpr = manual_roc_curve(
    df_logres_all,
    true_col=TRUE_COL,
    pred_col=PRED_COL,
    thresholds=thresholds,
)
auc = np.abs(np.trapz(tpr, fpr))
plt.figure(figsize=(7, 7))
plt.plot(fpr, tpr, label=f"AUC: {auc:.4f}")
plt.plot([0, 1], [0, 1], color="#AAAAAA", linestyle="--")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title("Apelin predicts AFib: all data")
plt.legend()
sns.despine(trim=True)
if SAVE_RESULTS:
    plt.savefig(
        os.path.join(RESULTS_FOLDER, "apelin_roc_all" + PLOTS_EXT),
        dpi=DPI,
        transparent=True,
        bbox_inches="tight",
    )
ideal_threshold_all = thresholds[np.argmax(tpr - fpr)]
get_metrics_for_thresh(
    df_logres_all, true_col=TRUE_COL, pred_col=PRED_COL, thresh=ideal_threshold_all
)<jupyter_output>Threshold: 0.8529393939393939<jupyter_text>Exploratory analysis of no-AFib differences* take two groups: no-AFib with Apelin < `ideal_threshold`, i.e. patients WITHOUT AFib and WITH low Apelin (those would be classified as AFib) by our dummy classifier* second groupd is NO AFib and high Apelin, i.e. correctly classified as no-AFib by our dummy classifier<jupyter_code>df_nofib = df[df["Typ"] == "Bez fibrilácie"]
df_nofib["Apelin_low"] = df_nofib["Apelin[ng/ml]"] < ideal_threshold_pat

bivar_nofib = bivariate_analysis(
    data=df_nofib,
    columns=get_cols_by_type(df_nofib, nans_thresh=20),
    exclude_cols=["Apelin_low"],
    match_col="Apelin_low",
    match_vals=[False, True],
)
# display ones with p-value < 0.5
if SAVE_RESULTS:
    bivar_nofib.round(3).to_csv(
        os.path.join(RESULTS_FOLDER, "bivar_noAFib_low_vs_high_apelin.csv")
    )
bivar_nofib[bivar_nofib["p-value"] < 0.5].round(3)<jupyter_output><empty_output><jupyter_text>Compare only low Apelin group* take only patients with Apelin < ideal_threshold (low apelin group) and compare AFib vs. no-AFib<jupyter_code>df_lowap = df[df["Apelin[ng/ml]"] < ideal_threshold_pat]


bivar_low_apelin = bivariate_analysis(
    data=df_lowap,
    columns=get_cols_by_type(df_lowap, nans_thresh=5),
    exclude_cols=["Typ", "Trvanie AFib (mesiace)"],
)
if SAVE_RESULTS:
    bivar_low_apelin.round(3).to_csv(
        os.path.join(RESULTS_FOLDER, "bivar_lowApelin_AFib_vs_noAFib.csv")
    )
# display ones with p-value < 0.5
bivar_low_apelin[bivar_low_apelin["p-value"] < 0.5].round(3)<jupyter_output><empty_output><jupyter_text>Logistic regression* bivariate analysis* variables with significant (p-value < $\alpha$) difference between AFib and no-AFib groups are added as predictors* logistic regression of all data to find significant predictors* ML-style logistic regression with Stratified K-Fold<jupyter_code># bivariate analysis
df_nans = df.replace(
    {
        "Neprítomná": 0,
        "I.\tstupňa": 1,
        "II.\tstupňa": 2,
        "III.\tstupňa": 3,
    }
)

overall = bivariate_analysis(
    data=df_nans,
    columns=get_cols_by_type(df_nans, nans_thresh=5),
    exclude_cols=["Typ", "Trvanie AFib (mesiace)"],
)
if SAVE_RESULTS:
    overall.round(3).to_csv(os.path.join(RESULTS_FOLDER, "bivar_AFib_vs_noAFib.csv"))
# display ones with p-value < 0.3
overall[overall["p-value"] < 0.3].round(3)<jupyter_output><empty_output><jupyter_text>Logistic regression with risk factors All data with detailed results* "Medical feature selection"<jupyter_code>risk_factors = [
    "Apelin[ng/ml]",
    "BMI",
    "NT-proBNP  [ng/l]",
    "Tlak krvi diastolický",
    "Tlak krvi systolický",
    "g. Známky srdcového zlyhávania",
    "l. Chronická obštrukčná choroba pľúc",
    "d. Priemer ľavej predsiene v PLAX (parasternal short axis) (mm)",
    "Vaskulárne ochorenie (1 bod)",
    "Diabetes mellitus (1 bod)",
    "Tyreotoxikóza",
    "CRP [mg/l]",
    "Kreatinin (µmol/l)",
    "D–dimér [ug/ml]",
    "Vek",
    "Pohlavie",
    "Diastolická dysfunkcia (I-III. stupňa).1",
]

predictors = set(risk_factors + ["Typ"])
predictors = list(
    predictors
    - set(
        [
            "Trvanie AFib (mesiace)",
            "Fructosamine[mmol/l]",
            "RAGE[pg/ml]",
            "AGEs[g/l]",
            "AGEs[rfu]",
        ]
    )
)
df_risk = prepare_df(df, predictors)
print(f"Predictors: {', '.join(df_risk.columns)}")
result = logres_statsmodels(df_risk)
if SAVE_RESULTS:
    statsmodels_summary_to_df(result).round(3).to_csv(
        os.path.join(RESULTS_FOLDER, "logres_medicalFS_summary.csv")
    )
result.summary()
fpr, tpr, thresh = roc_curve(df_risk["Afib"], result.predict(df_risk))
plt.figure(figsize=(7, 7))
plt.plot(
    fpr,
    tpr,
    label=f"AUC: {roc_auc_score(df_risk['Afib'], result.predict(df_risk)):.4f}",
)
plt.plot([0, 1], [0, 1], color="#AAAAAA", linestyle="--")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title("Manual feature selection")
plt.legend()
sns.despine(trim=True)
if SAVE_RESULTS:
    plt.savefig(
        os.path.join(RESULTS_FOLDER, "logres_all_data_medicalFS" + PLOTS_EXT),
        dpi=DPI,
        transparent=True,
        bbox_inches="tight",
    )
ideal_threshold_manual = thresh[np.argmax(tpr - fpr)]
print(f"Threshold: {ideal_threshold_manual}")
get_metrics_for_prediction(
    df_risk["Afib"], result.predict(df_risk) > ideal_threshold_manual
)<jupyter_output>Threshold: 0.6332172114058333<jupyter_text>CV Stratified K-Fold<jupyter_code>df_logres = prepare_df(df, predictors)
print(f"Predictors: {', '.join(df_logres.columns)}")
x_data = df_logres.loc[:, df_logres.columns != "Afib"].values
y_data = df_logres["Afib"].values

N_FOLDS = 5
N_REPEATS = 20

tpr, fpr, auc = [], [], []
acc, sens, spec = [], [], []
thresholds = np.linspace(0, 1, 100)


def manual_roc_curve(true, pred):
    tpr = []
    fpr = []
    for thresh in thresholds:
        tn, fp, fn, tp = confusion_matrix(true, pred > thresh).ravel()
        tpr.append(tp / (tp + fn))
        fpr.append(fp / (fp + tn))
    return np.array(tpr), np.array(fpr)


# define Repeated Stratified K-Fold
rskf = RepeatedStratifiedKFold(n_splits=N_FOLDS, n_repeats=N_REPEATS)
# define classifier
clf = make_pipeline(
    StandardScaler(),
    LogisticRegression(penalty="elasticnet", solver="saga", l1_ratio=0.5, max_iter=500),
)

for train_index, test_index in rskf.split(x_data, y_data):
    x_train, x_test = x_data[train_index], x_data[test_index]
    y_train, y_test = y_data[train_index], y_data[test_index]
    clf.fit(x_train, y_train)
    probs = clf.predict_proba(x_test)[:, 1]
    tpr_, fpr_ = manual_roc_curve(y_test, probs)
    auc_ = roc_auc_score(y_test, probs)
    ideal_threshold_ = thresholds[np.argmax(tpr_ - fpr_)]
    metrics = get_metrics_for_kfold(y_test, probs > ideal_threshold_)
    tpr.append(tpr_)
    fpr.append(fpr_)
    auc.append(auc_)
    acc.append(metrics[0])
    sens.append(metrics[1])
    spec.append(metrics[2])
plt.figure(figsize=(9, 9))
tpr = np.vstack(tpr)
fpr = np.vstack(fpr)
mean_auc = np.mean(auc)
std_auc = np.std(auc)
plt.plot(
    fpr.mean(axis=0),
    tpr.mean(axis=0),
    label=f"AUC: {mean_auc:.3f} $\pm$ {std_auc:.3f}",
    linewidth=3.5,
)
plt.fill_between(
    fpr.mean(axis=0),
    tpr.mean(axis=0) - tpr.std(axis=0),
    tpr.mean(axis=0) + tpr.std(axis=0),
    alpha=0.3,
)
plt.plot([0, 1], [0, 1], color="#AAAAAA", linestyle="--")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title(
    f"Logistic regression ROC \n Stratified K-Fold; {N_FOLDS} folds with {N_REPEATS} repeats"
)
plt.legend()
sns.despine(trim=True)
if SAVE_RESULTS:
    plt.savefig(
        os.path.join(RESULTS_FOLDER, "logres_KFold_medicalFS" + PLOTS_EXT),
        dpi=DPI,
        transparent=True,
        bbox_inches="tight",
    )
    
print("---------")
print(f"Accuracy: {np.mean(acc):.3f} ± {np.std(acc):.3f}")
print(f"Sensitivity: {np.mean(sens):.3f} ± {np.std(sens):.3f}")
print(f"Specificity: {np.mean(spec):.3f} ± {np.std(spec):.3f}")<jupyter_output>---------
Accuracy: 0.721 ± 0.093
Sensitivity: 0.787 ± 0.243
Specificity: 0.657 ± 0.239<jupyter_text>Logistic regression from bivariate analysis All data with detailed results* Bivariate analysis feature selection<jupyter_code>PREDICTOR_THRESH = 0.1
predictors = set(list(overall[overall["p-value"] < PREDICTOR_THRESH].index) + ["Typ"])
predictors = list(
    predictors
    - set(
        [
            "<NAME> (mesiace)",
            "Fructosamine[mmol/l]",
            "RAGE[pg/ml]",
            "AGEs[g/l]",
            "AGEs[rfu]",
            "i. Inhibítor protónovej pumpy",
            "<NAME>",
        ]
    )
)
df_bivar = prepare_df(df, predictors)
print(f"Predictors: {', '.join(df_bivar.columns)}")
result = logres_statsmodels(df_bivar)
if SAVE_RESULTS:
    statsmodels_summary_to_df(result).round(3).to_csv(
        os.path.join(RESULTS_FOLDER, "logres_bivarFS_summary.csv")
    )
result.summary()
fpr, tpr, thresh = roc_curve(df_bivar["Afib"], result.predict(df_bivar))
plt.figure(figsize=(7, 7))
plt.plot(
    fpr,
    tpr,
    label=f"AUC: {roc_auc_score(df_bivar['Afib'], result.predict(df_bivar)):.4f}",
)
plt.plot([0, 1], [0, 1], color="#AAAAAA", linestyle="--")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title("Bivariate analysis feature selection")
plt.legend()
sns.despine(trim=True)
if SAVE_RESULTS:
    plt.savefig(
        os.path.join(RESULTS_FOLDER, "logres_all_data_bivarFS" + PLOTS_EXT),
        dpi=DPI,
        transparent=True,
        bbox_inches="tight",
    )
ideal_threshold_bivar = thresh[np.argmax(tpr - fpr)]
print(f"Threshold: {ideal_threshold_bivar}")
get_metrics_for_prediction(
    df_bivar["Afib"], result.predict(df_bivar) > ideal_threshold_bivar
)<jupyter_output>Threshold: 0.5721753523810448<jupyter_text>CV with Stratified K-Fold<jupyter_code>df_logres = prepare_df(df, predictors)
print(f"Predictors: {', '.join(df_logres.columns)}")
x_data = df_logres.loc[:, df_logres.columns != "Afib"].values
y_data = df_logres["Afib"].values

N_FOLDS = 5
N_REPEATS = 20

tpr, fpr, auc = [], [], []
acc, sens, spec = [], [], []
thresholds = np.linspace(0, 1, 100)


def manual_roc_curve(true, pred):
    tpr = []
    fpr = []
    for thresh in thresholds:
        tn, fp, fn, tp = confusion_matrix(true, pred > thresh).ravel()
        tpr.append(tp / (tp + fn))
        fpr.append(fp / (fp + tn))
    return np.array(tpr), np.array(fpr)


# define Repeated Stratified K-Fold
rskf = RepeatedStratifiedKFold(n_splits=N_FOLDS, n_repeats=N_REPEATS)
# define classifier
clf = make_pipeline(
    StandardScaler(),
    LogisticRegression(penalty="elasticnet", solver="saga", l1_ratio=0.5, max_iter=500),
)

for train_index, test_index in rskf.split(x_data, y_data):
    x_train, x_test = x_data[train_index], x_data[test_index]
    y_train, y_test = y_data[train_index], y_data[test_index]
    clf.fit(x_train, y_train)
    probs = clf.predict_proba(x_test)[:, 1]
    tpr_, fpr_ = manual_roc_curve(y_test, probs)
    auc_ = roc_auc_score(y_test, probs)
    ideal_threshold_ = thresholds[np.argmax(tpr_ - fpr_)]
    metrics = get_metrics_for_kfold(y_test, probs > ideal_threshold_)
    tpr.append(tpr_)
    fpr.append(fpr_)
    auc.append(auc_)
    acc.append(metrics[0])
    sens.append(metrics[1])
    spec.append(metrics[2])
plt.figure(figsize=(9, 9))
tpr = np.vstack(tpr)
fpr = np.vstack(fpr)
mean_auc = np.mean(auc)
std_auc = np.std(auc)
plt.plot(
    fpr.mean(axis=0),
    tpr.mean(axis=0),
    label=f"AUC: {mean_auc:.3f} $\pm$ {std_auc:.3f}",
    linewidth=3.5,
)
plt.fill_between(
    fpr.mean(axis=0),
    tpr.mean(axis=0) - tpr.std(axis=0),
    tpr.mean(axis=0) + tpr.std(axis=0),
    alpha=0.3,
)
plt.plot([0, 1], [0, 1], color="#AAAAAA", linestyle="--")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title(
    f"Logistic regression ROC \n Stratified K-Fold; {N_FOLDS} folds with {N_REPEATS} repeats"
)
plt.legend()
sns.despine(trim=True)
if SAVE_RESULTS:
    plt.savefig(
        os.path.join(RESULTS_FOLDER, "logres_KFold_bivarFS" + PLOTS_EXT),
        dpi=DPI,
        transparent=True,
        bbox_inches="tight",
    )
    
    
print("---------")
print(f"Accuracy: {np.mean(acc):.3f} ± {np.std(acc):.3f}")
print(f"Sensitivity: {np.mean(sens):.3f} ± {np.std(sens):.3f}")
print(f"Specificity: {np.mean(spec):.3f} ± {np.std(spec):.3f}")<jupyter_output>---------
Accuracy: 0.767 ± 0.098
Sensitivity: 0.851 ± 0.209
Specificity: 0.685 ± 0.250
{'fraction_non_alphanumeric': 0.1653466382889641, 'fraction_numerical': 0.015142528174173964, 'mean_word_length': 3.7873726407215633, 'pattern_counts': {'":': 50, '<': 74, '<?xml version=': 0, '>': 73, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '531748', 'license': 'MIT', 'n_tokens_mistral': 12387, 'n_tokens_neox': 11742, 'n_words': 2084, 'path': 'AFISBIO_paper_final.ipynb', 'repo_name': 'jajcayn/apelin_vs_atrial_fibrillation', 'size': '28815'}