<jupyter_start><jupyter_text>CFT – ECG Analysis  Setup and Helper Functions<jupyter_code>import json
import re
from pathlib import Path

import pandas as pd
import numpy as np
import pingouin as pg

import matplotlib.pyplot as plt
import seaborn as sns

from fau_colors import cmaps
import biopsykit as bp
from biopsykit.utils.dataframe_handling import multi_xs

from cft_analysis.datasets import CftDatasetProcessed

from IPython.display import Markdown

%load_ext autoreload
%autoreload 2
%matplotlib widget
sns.color_palette(cmaps.fau)
[list((int(c * 255) for c in x)) for x in cmaps.fau]
plt.close("all")

palette = sns.color_palette(cmaps.faculties)
sns.set_theme(context="notebook", style="ticks", palette=palette)

plt.rcParams["figure.figsize"] = (8, 4)
plt.rcParams["pdf.fonttype"] = 42
plt.rcParams["mathtext.default"] = "regular"

pg.options["round"] = 4

palette<jupyter_output><empty_output><jupyter_text>Import Data<jupyter_code># get path to analysis results
base_path = Path("../../data")
results_path = base_path.joinpath("../results")
stats_path = results_path.joinpath("statistics")
tex_path = stats_path.joinpath("tex_tables")
bp.utils.file_handling.mkdirs([results_path, stats_path, tex_path])

paper_path = Path("../paper_path.json")
paper_tex_path = None
if paper_path.exists():
    paper_path = Path(json.load(paper_path.open(encoding="utf-8"))["paper_path"])
    paper_tex_path = paper_path.joinpath("tab")
    bp.utils.file_handling.mkdirs([paper_tex_path])
dataset = CftDatasetProcessed(base_path, exclude_subjects=True)
dataset
hue_order = ["Control", "CFT"]<jupyter_output><empty_output><jupyter_text>Response to the MIST  Descriptive Analysis To show the general efficacy of the MIST we compute the mean and maximum heart rate increase for each subject and MIST phase, respectively.  Afterwards, we compute mean and standard deviation over all subjects:<jupyter_code>hr_ensemble = dataset.heart_rate_ensemble
hr_ensemble.head()
hr_increase = hr_ensemble.groupby("phase").agg(["mean", "max"])
hr_increase = hr_increase.stack().agg(["mean", "std"], axis=1)
hr_increase.round(2)<jupyter_output><empty_output><jupyter_text>Statistical Analysis **Population**: Control condition**Analysis** (per measure):1. *Increase BL-AT*: Check whether Arithmetic Tasks of the MIST causes significant HR(V) responses.    * Procedure: Paired t-tests between subphases *BL* and *AT* for each MIST phase    * Expected Result: Significant HR(V) responses for each MIST phase2. *Increase over MIST*: Check whether the HR response increases over time (i.e., over MIST phases).     * Procedure: Repeated-measures ANOVA for each subphase    * *(only for HR)*: Post-hoc test to check which pairs of MIST phases show significant differences    * Expected Result: Significant main effect *MIST Phase* for each subphase**Findings**: * Each MIST Phase induces stress, indicated by significant differences between $BL_{Loc}$ and $AT$ in each individual MIST Phase* Stress levels increase over time, indicated by significant main effect *MIST Phase* and significant differences between MIST phases, especially between *MIST1* and *MIST3*  Increase BL-AT<jupyter_code>display(Markdown("##### Prepare Data"))

condition = "Control"
subphases = ["BL", "AT"]
measures = ["HR_Norm", "HRV_RMSSD", "HRV_pNN50"]

hr_analysis = dataset.get_subset(condition=condition, subphase=subphases).hr_hrv
hr_analysis = multi_xs(hr_analysis, measures, level="type")

display(hr_analysis.head())


################################################
display(Markdown("##### Statistics"))

steps = [("prep", "normality"), ("test", "pairwise_ttests")]
params = {
    "dv": "data",
    "within": "subphase",
    "subject": "subject",
    "groupby": ["type", "phase"],
    "multicomp": {"levels": "phase"},
}

stats = bp.stats.StatsPipeline(steps, params)
stats.apply(hr_analysis)

stats.export_statistics(stats_path.joinpath("stats_hrv_response_bl_at.xlsx"))
stats.display_results(prep=False)


################################################
display(Markdown("##### Latex Output"))

index_rename_map = {"HR_Norm": "HR", "HRV_RMSSD": "RMSSD", "HRV_pNN50": "pRR50"}

index_value_order = ["HR_Norm", "HRV_RMSSD", "HRV_pNN50"]

caption = "HR(V) responses of the Control condition to the MIST. \
Paired t-tests were performed between BL and AT subphases for each individual MIST phase.\
"

df_latex = stats.results_to_latex_table(
    "pairwise_ttests",
    unstack_levels="phase",
    index_kws={
        "index_value_order": index_value_order,
        "index_rename_map": index_rename_map,
        "index_level_names_tex": ["Measure"],
    },
    caption=caption,
    label="tab:hrv_response_mist",
)

for path in [tex_path, paper_tex_path]:
    if path is not None:
        path.joinpath("tab_hrv_response_bl_at.tex").open(mode="w+").write(df_latex)
print(df_latex)<jupyter_output><empty_output><jupyter_text>Increase over MIST<jupyter_code>display(Markdown("##### Prepare Data"))

measures = ["HR_Norm", "HRV_RMSSD", "HRV_pNN50"]
subphases = ["BL", "AT", "FB"]
conditions = ["Control"]

hr_analysis = dataset.get_subset(condition=conditions, subphase=subphases).hr_hrv

# Select data
hr_analysis = multi_xs(hr_analysis, measures, level="type")
display(hr_analysis.head())

#################################################################
display(Markdown("##### Statistics"))

steps = [("prep", "normality"), ("prep", "equal_var"), ("test", "rm_anova"), ("posthoc", "pairwise_ttests")]
params = {
    "dv": "data",
    "within": "phase",
    "subject": "subject",
    "groupby": ["subphase", "type"],
    "multicomp": {"levels": False},
}

stats = bp.stats.StatsPipeline(steps, params)
stats.apply(hr_analysis)

stats.export_statistics(stats_path.joinpath("stats_hrv_response_mist_phases.xlsx"))
stats.display_results(prep=True)  # , sig_only="posthoc")

#################################################################
display(Markdown("##### Latex Output"))
index_rename_map = {"HR_Norm": "HR", "HRV_RMSSD": "RMSSD", "HRV_pNN50": "pRR50"}

index_value_order = {"type": ["HR_Norm", "HRV_RMSSD", "HRV_pNN50"], "subphase": ["BL", "AT", "FB"]}

caption = "HR(V) measures of the Control condition during MIST. \
Repeated-measurement ANOVAs were performed separately for BL, AT, and FB subphases. \
Within-variable: MIST phase."

df_latex = stats.results_to_latex_table(
    "rm_anova",
    unstack_levels=["subphase"],
    index_kws={
        "index_value_order": index_value_order,
        "index_rename_map": index_rename_map,
        "index_level_names_tex": ["Measure"],
    },
    caption=caption,
    label="tab:hrv_increase_mist",
    multicolumn_format="c",
    escape=False,
    multirow=True,
)

for path in [tex_path, paper_tex_path]:
    if path is not None:
        paper_tex_path.joinpath("tab_hrv_response_mist_phases.tex").open(mode="w+").write(df_latex)

print(stats.stats_to_latex("pairwise_ttests"))
print(df_latex)<jupyter_output><empty_output><jupyter_text>Response to the CFT<jupyter_code>features = ["onset_latency", "peak_brady_latency", "peak_brady_percent", "mean_brady_percent"]<jupyter_output><empty_output><jupyter_text>Descriptive Analysis<jupyter_code>cft_analysis = dataset.cft_parameter
cft_analysis.head()
cft_analysis = multi_xs(cft_analysis, features, level="type")
cft_analysis_agg = cft_analysis.unstack()["data"]

cft_analysis_agg = (
    cft_analysis_agg.groupby("phase")
    .agg(
        {
            "onset_latency": ["mean", "std"],
            "peak_brady_latency": ["mean", "std"],
            "peak_brady_percent": ["min", "mean", "std"],
            "mean_brady_percent": ["min", "mean", "std"],
        }
    )
    .T
)
cft_analysis_agg
cft_time_bl_glo = dataset.get_subset(condition="CFT", subphase="RP_CFI").time_above_baseline
cft_time_bl_glo = cft_time_bl_glo.xs("HR", level="type")
cft_time_bl_glo = cft_time_bl_glo.groupby("phase").agg(["mean", "std"]).T
cft_time_bl_glo.round(2)<jupyter_output><empty_output><jupyter_text>Statistical Analysis **Population**: CFT condition  **Analysis** (per measure):1. *Decrease BL-RP_CFI*: Check whether Cold Face Intervention causes significant HR(V) responses.    * Procedure: Paired t-tests between subphases *BL* and *RP_CFI* for each MIST phase    * Expected Result: Significant HR(V) responses during each MIST phase**Findings**: * Each MIST Phase induces stress, indicated by significant differences between $BL_{Loc}$ and $AT$ in each individual MIST Phase* Stress levels increase over time, indicated by significant main effect *MIST Phase* and significant differences between MIST phases, especially between *MIST1* and *MIST3*  Decrease *BL-RP_CFI*<jupyter_code>display(Markdown("##### Prepare Data"))

subphases = ["BL", "RP_CFI"]
conditions = ["CFT"]
measures = ["HR_Norm", "HRV_RMSSD", "HRV_pNN50"]

hr_analysis = dataset.get_subset(condition=conditions, subphase=subphases).hr_hrv
hr_analysis = multi_xs(hr_analysis, measures, level="type")
display(hr_analysis.head())

################################################################
display(Markdown("##### Statistics"))

steps = [("prep", "normality"), ("test", "pairwise_ttests")]
params = {
    "dv": "data",
    "within": "subphase",
    "subject": "subject",
    "groupby": ["type", "phase"],
    "multicomp": {"levels": "phase"},
}

stats = bp.stats.StatsPipeline(steps, params)
stats.apply(hr_analysis)

stats.export_statistics(stats_path.joinpath("stats_hr_response_bl_cfi.xlsx"))
stats.display_results(prep=False)

################################################################
display(Markdown("##### Latex Output"))

index_rename_map = {"HR_Norm": "HR", "HRV_RMSSD": "RMSSD", "HRV_pNN50": "pRR50"}

index_value_order = ["HR_Norm", "HRV_RMSSD", "HRV_pNN50"]

caption = "HR(V) responses of the CFI condition to the Cold Face Test. \
Paired t-tests were performed between BL and RP/CFI subphases for each MIST phase separately."

df_latex = stats.results_to_latex_table(
    "pairwise_ttests",
    unstack_levels="phase",
    index_kws={
        "index_value_order": index_value_order,
        "index_rename_map": index_rename_map,
        "index_level_names_tex": ["Measure"],
    },
    caption=caption,
    label="tab:hrv_response_cft",
)

for path in [tex_path, paper_tex_path]:
    if path is not None:
        path.joinpath("tab_hrv_response_bl_cfi.tex").open(mode="w+").write(df_latex)
print(df_latex)<jupyter_output><empty_output><jupyter_text>Effect of CFT on Acute Stress Response  Statistical Analysis  Interaction Condition x MIST **Population**: Control vs. CFT condition  **Analysis** (per measure):1. *Interaction Condition x MIST Phase*: Check whether Cold Face Intervention has an effect on the HR(V) responses to the MIST.    * Procedure:        * Mixed-ANOVA to determine interaction effect between *Condition* and *MIST Phase*        * In case of significant interaction effect: Post-hoc test to assess during which MIST Phase the conditions showed HR(V) response differences    * Expected Result: Significant interaction effect between *Condition* and *MIST Phase***Findings**: * Each MIST Phase induces stress, indicated by significant differences between $BL_{Loc}$ and $AT$ in each individual MIST Phase* Stress levels increase over time, indicated by significant main effect *MIST Phase* and significant differences between MIST phases, especially between *MIST1* and *MIST3*<jupyter_code>display(Markdown("##### Prepare Data"))

phases = ["MIST1", "MIST2", "MIST3"]
subphases = ["BL", "RP_CFI"]
hr_measures = ["HR_Norm", "HRV_RMSSD", "HRV_pNN50"]
time_bl_measures = ["HR", "HRV_RMSSD", "HRV_pNN50"]

hr_analysis = dataset.get_subset(phase=phases, subphase=subphases).hr_hrv
hr_analysis = multi_xs(hr_analysis, hr_measures, level="type")
time_bl_analysis = dataset.get_subset(phase=phases, subphase=subphases).time_above_baseline
time_bl_analysis = multi_xs(time_bl_analysis, time_bl_measures, level="type")

hr_analysis = hr_analysis.unstack(["type", "category"]).join(time_bl_analysis.unstack(["type", "category"]))
hr_analysis = hr_analysis.stack(["category", "type"])
display(hr_analysis.head())


#################################################################
display(Markdown("##### Statistics"))

steps = [("prep", "normality"), ("prep", "equal_var"), ("test", "mixed_anova"), ("posthoc", "pairwise_ttests")]
params = {
    "dv": "data",
    "within": "phase",
    "between": "condition",
    "subject": "subject",
    "groupby": ["category", "type", "subphase"],
    "padjust": "bonf",
}

stats = bp.stats.StatsPipeline(steps, params)
stats.apply(hr_analysis)

stats.export_statistics(stats_path.joinpath("stats_effect_cft_interaction.xlsx"))
stats.display_results(prep=False, posthoc=False, sig_only=["test", "posthoc"])

#################################################################
display(Markdown("##### Latex Output"))

index_rename_map = {
    "HR_HR_Norm": "$\Delta HR$",
    "HRV_HRV_RMSSD": "RMSSD",
    "HRV_HRV_pNN50": "pRR50",
    "RP_CFI": "RP/CFI",
    "Time_BL_Glo_HR": "$\hat{t}_{Glo}(HR)$",
    "Time_BL_Glo_HRV_RMSSD": "$\hat{t}_{Glo}(RMSSD)$",
    "Time_BL_Glo_HRV_pNN50": "$\hat{t}_{Glo}(pRR50)$",
}

index_value_order = {
    "type": [
        "HR_HR_Norm",
        "Time_BL_Glo_HR",
        "HRV_HRV_RMSSD",
        "Time_BL_Glo_HRV_RMSSD",
        "HRV_HRV_pNN50",
        "Time_BL_Glo_HRV_pNN50",
    ]
}

#################################################################
display(Markdown("###### Interaction Effect"))

data = stats.results["mixed_anova"]
data = data.xs("BL", level="subphase")
data = data.unstack(-1)
data.index = ["_".join(val) for val in data.index]
data.index.names = ["type"]
data = data.stack(-1)

caption = "Mixed-ANOVA results (interaction effect MIST Phase x Condition) of HR(V) measures during BL subphase."

df_latex = stats.results_to_latex_table(
    "mixed_anova",
    data=data,
    stats_effect_type="Interaction",
    collapse_dof=True,
    index_kws={
        "index_value_order": index_value_order,
        "index_rename_map": index_rename_map,
        "index_level_names_tex": ["Measure"],
    },
    caption=caption,
    label="tab:cft_mist_interaction",
)

for path in [tex_path, paper_tex_path]:
    if path is not None:
        path.joinpath("tab_hrv_response_cft_interaction.tex").open(mode="w+").write(df_latex)
print(df_latex)


#################################################################
display(Markdown("###### Main Effect Condition"))

data = stats.results["mixed_anova"]
data = data.xs("RP_CFI", level="subphase")
data = data.unstack(-1)
data.index = ["_".join(val) for val in data.index]
data.index.names = ["type"]
data = data.stack(-1)

caption = "Mixed-ANOVA results (main effect Condition) of HR(V) measures during RP/CFI subphase."

df_latex = stats.results_to_latex_table(
    "mixed_anova",
    data=data,
    stats_effect_type="condition",
    index_kws={
        "index_value_order": index_value_order,
        "index_rename_map": index_rename_map,
        "index_level_names_tex": ["Measure"],
    },
    caption=caption,
    label="tab:cft_mist_main_effect",
)
for path in [tex_path, paper_tex_path]:
    if path is not None:
        path.joinpath("tab_hrv_response_cft_main_effect.tex").open(mode="w+").write(df_latex)
print(df_latex)

#################################################################
display(Markdown("###### Post-Hoc"))

posthoc_results = stats._filter_effect("posthoc", "interaction")
posthoc_results = multi_xs(posthoc_results, ["HR", "Time_BL_Glo"], level="category")
posthoc_results = multi_xs(posthoc_results, ["HR", "HR_Norm"], level="type")
posthoc_results = posthoc_results.xs("BL", level="subphase")
posthoc_results = posthoc_results.set_index("phase", append=True)

print(stats.stats_to_latex(data=posthoc_results))<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code>category = ["HR"]
feature = ["HR_Norm"]
subphase = ["BL"]

hr_plot = dataset.get_subset(subphase=subphase).heart_rate
hr_plot = multi_xs(hr_plot, feature, level="type")

stats_data = stats.results_cat("posthoc").loc[category, feature, subphase]
box_pairs, pvalues = stats.sig_brackets(stats_data, stats_effect_type="interaction", plot_type="multi", x="phase")

fig, ax = plt.subplots()
bp.plotting.feature_boxplot(
    data=hr_plot.reset_index(),
    x="phase",
    y="data",
    hue="condition",
    hue_order=hue_order,
    stats_kwargs={"box_pairs": box_pairs, "pvalues": pvalues},
    legend_loc="upper left",
    ax=ax,
)

ax.set_xlabel("Phase")
ax.set_ylabel("$\Delta$ HR during $BL_{Loc}$ [%]")

fig.tight_layout()
category = ["Time_BL_Glo"]
feature = ["HR"]
subphase = ["BL"]

hr_plot = dataset.get_subset(subphase=subphase).time_above_baseline
hr_plot = multi_xs(hr_plot, feature, level="type")


stats_data = stats.results_cat("posthoc").loc[category, feature, subphase, :]
box_pairs, pvalues = stats.sig_brackets(stats_data, stats_effect_type="interaction", plot_type="multi", x="phase")

fig, ax = plt.subplots()
bp.plotting.feature_boxplot(
    data=hr_plot.reset_index(),
    x="phase",
    y="data",
    hue="condition",
    hue_order=hue_order,
    stats_kwargs={"box_pairs": box_pairs, "pvalues": pvalues},
    legend_loc="upper left",
    ax=ax,
)


ax.set_xlabel("Phase")
ax.set_ylabel("$\hat{t}_{Glo}(HR)$ during BL [%]")

fig.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18053384285797, 'fraction_numerical': 0.0040530368826356316, 'mean_word_length': 6.096138044371405, 'pattern_counts': {'":': 66, '<': 41, '<?xml version=': 0, '>': 41, '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': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '22531', 'license': 'MIT', 'n_tokens_mistral': 6024, 'n_tokens_neox': 5796, 'n_words': 1224, 'path': 'experiments/2022_scientific_reports/notebooks/analysis/ECG_Analysis.ipynb', 'repo_name': 'mad-lab-fau/cft-analysis', 'size': '17271'}