docs/examples/group_pipeline/3.group_data_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9790169d",
   "metadata": {},
   "source": [
    "# Data Analysis\n",
    "\n",
    "Finally, the fun part! :) `wearablehrv` allows you to easily calculate and plot the most important analyses for validation. It can establish for each device, against the criterion, how valid it is in each of your experimental conditions. To achieve this, this notebook example will guide you through the important statistical analyses that are incorporated in `wearablehrv`:\n",
    "\n",
    "- Mean absolute percentage error (MAPE)\n",
    "- Regression analysis\n",
    "- Intraclass correlation coefficient (ICC)\n",
    "- Bland-Altman analysis\n",
    "\n",
    "Not only this, but it also provides functionality to deal with the non-normality of your data, correct for multiple testing, etc.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5930aa1d",
   "metadata": {},
   "source": [
    "<div style=\"border:1px solid; padding:10px; border-radius:5px; margin:10px 0;\">\n",
    "\n",
    "**Note**: Throughout the example notebooks and also in the code, we used the term \"<u>criterion</u>,\" which refers to the device that the rest of the devices are compared against. This is also referred to as \"reference system,\" \"ground truth,\" and \"gold standard\" in the literature. This is usually an electrocardiography (ECG) device.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ee9fe1c",
   "metadata": {},
   "source": [
    "## Previous Steps\n",
    "\n",
    "If you have not done so, first take a look at the following notebooks:\n",
    "\n",
    "- [Determine the signal quality of your wearables](https://github.com/Aminsinichi/wearable-hrv/blob/master/docs/examples/group_pipeline/2.group_signal_quality.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4929a3d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing Module\n",
    "import wearablehrv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c55c49ac",
   "metadata": {},
   "source": [
    "The code in the following cell has been explained in the previous notebook. Run it, so we can continue with the examples in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "504cffe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "wearablehrv.data.clear_wearablehrv_cache() \n",
    "path = wearablehrv.data.download_data_and_get_path([\"P01.csv\", \"P02.csv\", \"P03.csv\", \"P04.csv\", \"P05.csv\", \"P06.csv\", \"P07.csv\", \"P08.csv\", \"P09.csv\", \"P10.csv\"])\n",
    "conditions = ['sitting', 'arithmetic', 'recovery', 'standing', 'breathing', 'neurotask', 'walking', 'biking'] \n",
    "devices = [\"kyto\", \"heartmath\", \"rhythm\", \"empatica\", \"vu\"] \n",
    "criterion = \"vu\" \n",
    "features = [\"rmssd\", \"hf\",'mean_hr', 'nibi_after_cropping', 'artefact'] \n",
    "data, file_names = wearablehrv.group.import_data (path, conditions, devices, features)\n",
    "data = wearablehrv.group.nan_handling (data, devices, features, conditions) \n",
    "data, features, summary_df, quality_df = wearablehrv.group.signal_quality (data, path, conditions, devices, features, criterion,  file_names, exclude = False, save_as_csv = False, ibi_threshold = 0.20, artefact_threshold = 0.20, manual_missing=False, missing_threshold=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4a813f1",
   "metadata": {},
   "source": [
    "## 1. Mean Absolute Percentage Error\n",
    "\n",
    "The `mape_analysis` function computes the mean absolute percentage error by taking the absolute difference between each device's measurements and those of a criterion device for each participant, dividing by the criterion's measurements to obtain percentage errors, and then averaging these percentages across all participants within each condition and feature for each device. Optionally, if `save_as_csv` is set to `True` and a path is provided, the function will save these results in a CSV file at the specified location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26287bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mape_data = wearablehrv.group.mape_analysis (data, criterion, devices, conditions, features, path=None, alpha=0.95, save_as_csv=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0052a918",
   "metadata": {},
   "source": [
    "For instance, if you now want to see the MAPE of the Empatica for mean heart rate in the biking condition, you can run the following code and retrieve the MAPE and confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9732d182",
   "metadata": {},
   "outputs": [],
   "source": [
    "mape_data[\"empatica\"][\"mean_hr\"][\"biking\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba3324bb",
   "metadata": {},
   "source": [
    "\n",
    "The `mape_plot` function generates grouped bar charts to visualize MAPE values for each device across different conditions for a selected feature. The function allows interactive selection of features for comparison. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58f8542b",
   "metadata": {},
   "outputs": [],
   "source": [
    "wearablehrv.group.mape_plot (mape_data, features, conditions, devices)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99c9a657",
   "metadata": {},
   "source": [
    "Looking at this plot shows you how error appears in different devices and conditions. For instance, you can immediately see how, by increasing movement in the conditions, the error got a lot larger in all devices, yet the magnitude is quite different!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b986dc7",
   "metadata": {},
   "source": [
    "## Check for Normality and Log Transformation (Optional)\n",
    "\n",
    "It might be the case that your data is not normally distributed. It is good to check this before doing any further statistical analysis. There are two functions in place for you:\n",
    "\n",
    "The `check_normality` function assesses the normality of data using the Shapiro-Wilk test. It returns three dictionaries:\n",
    "1. **`normality_results`**: Contains p-values and normality status for each device, feature, and condition.\n",
    "2. **`suggestions`**: Provides transformation recommendations for non-normal data.\n",
    "3. **`transformation_summary`**: Summarizes conditions requiring transformation for each device and feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62bf092e",
   "metadata": {},
   "outputs": [],
   "source": [
    "normality_results, suggestions, transformation_summary = wearablehrv.group.check_normality (data, conditions, devices, features, alpha=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b2f2621",
   "metadata": {},
   "source": [
    "You can then get a sense for each device, for which features you have the most issues with non-normality. For instance, if you check out the \"kyto\" device, you can see that for \"rmssd\" and \"hf,\" there is non-normality in many conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "383ee4a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "transformation_summary[\"kyto\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df004d45",
   "metadata": {},
   "source": [
    "Based on your results, you can decide which features require transformation, and then you can use the following function:\n",
    "\n",
    "The `log_transform_data` function applies a log transformation to specified features ('rmssd', 'hf') in the data for all devices and conditions. \n",
    "\n",
    "**It modifies the `data` dictionary in place**, without returning any value. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f0cf39a",
   "metadata": {},
   "outputs": [],
   "source": [
    "transform_features = [\"rmssd\", \"hf\"]\n",
    "wearablehrv.group.log_transform_data(data, devices, transform_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b293194",
   "metadata": {},
   "source": [
    "## 2. Regression Analysis "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53bc8ea3",
   "metadata": {},
   "source": [
    "An important step in validation is to determine the linear strength between a measure calculated from a given device (e.g., RMSSD from Empatica in sitting condition) against the criterion counterpart (e.g., RMSSD from VU, in sitting condition). The `regression_analysis` function does this for you and returns the slope, intercept, r-value, p-value, and standard error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81256e30",
   "metadata": {},
   "outputs": [],
   "source": [
    "regression_data = wearablehrv.group.regression_analysis (data, criterion, conditions, devices, features, path, save_as_csv=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad42987b",
   "metadata": {},
   "source": [
    "Let's check out and see how Empatica indeed linearly correlates with VU by running the following code. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "834e858d",
   "metadata": {},
   "outputs": [],
   "source": [
    "regression_data[\"empatica\"][\"rmssd\"][\"sitting\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fef3f4c",
   "metadata": {},
   "source": [
    "You will see the output results, and a nonsignificant (0.07) p-value, despite a high r-value. Keep in mind that we have loaded only 10 participants as examples, and only four of them did have Empatica!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "945fe139",
   "metadata": {},
   "source": [
    "### Adjusting P-values with Bonferonni (Optional)\n",
    "\n",
    "Given multiple testing, you can correct the p-values for multiple comparison. The `bonferroni_correction_regression` function applies the Bonferroni correction to p-values in regression analysis results. It returns a dictionary similar to the input but with Bonferroni corrected p-values and a flag indicating significance based on the corrected alpha level.\n",
    "\n",
    "**Note**: Please be aware that the number of tests performed is based on the \"features\" array. Therefore, be mindful of the number of variables you include in this array. If there are non-relevant variables present, redefine your \"features\" list, then rerun the `regression_analysis` function and then correct for Bonferroni. For instance, you can only keep a few relevant features: `features = [\"rmssd\", \"hf\",'mean_hr']`\n",
    "\n",
    "**It modifies the `regression_data` dictionary in place**, and adds `corrected_p_value` and `is_significant` to it.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "740692f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "wearablehrv.group.bonferroni_correction_regression (regression_data, alpha=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "242ed6f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now regression_data dictionary contains two more keys: corrected_p_value and is_significant\n",
    "regression_data[\"empatica\"][\"rmssd\"][\"sitting\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7192dfb",
   "metadata": {},
   "source": [
    "The `regression_plot` function creates scatter plots with regression lines for device and criterion pairs across different conditions and features. It visualizes the correlation and significance (p-value) between devices.\n",
    "\n",
    "**Note:** If `bonferroni_corrected` is set to `True`, the function uses Bonferroni corrected p-values for significance annotations in the plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d57320c",
   "metadata": {},
   "outputs": [],
   "source": [
    "wearablehrv.group.regression_plot(regression_data, data, criterion, conditions, devices, features,\n",
    "                    width=15, height_per_condition=4, \n",
    "                    regression_line_style='-', regression_line_color='black', \n",
    "                    marker_color='gray', font_size=12, \n",
    "                    show_grid=True, background_color=None, \n",
    "                    bonferroni_corrected=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32ea6632",
   "metadata": {},
   "source": [
    "The `heatmap_plot` function generates a correlation heatmap that displays the correlation between each device's data and the criterion device's data for each condition and feature. The heatmap provides a visual representation of how closely the data from each device matches the data from the criterion device for the selected feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "972b2546",
   "metadata": {},
   "outputs": [],
   "source": [
    "wearablehrv.group.heatmap_plot (data, criterion, devices, conditions, features)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5928f459",
   "metadata": {},
   "source": [
    "## 3. Intraclass Correlation Coefficient \n",
    "\n",
    "Intraclass correlation coefficient (ICC) is a crucial statistical tool used in validation studies to measure the consistency or agreement of measurements made by different observers measuring the same entity. This is particularly important in validation because it quantifies how much the measurements taken with different devices, under various conditions, or by different raters, can be reliably compared or interchanged.\n",
    "\n",
    "The `icc_analysis` function computes ICC for different devices, conditions, and features. In this function, ICC is computed using the two-way random-effects model, which considers both the inter-rater and intra-rater variability. Given we use [`pingouin.intraclass_corr`](https://pingouin-stats.org/build/html/generated/pingouin.intraclass_corr.html), you may want to check out the documentation to better understand the outputs such as ICC1, ICC2, ICC3, etc.\n",
    "\n",
    "If `save_as_csv = True`, the output will be saved as a .csv file in your specified path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bafbd56",
   "metadata": {},
   "outputs": [],
   "source": [
    "icc_data = wearablehrv.group.icc_analysis (data, criterion, devices, conditions, features, path, save_as_csv=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "713fd7a6",
   "metadata": {},
   "source": [
    "Let's take a look at the output to see different ICC values for the \"heartmath\" device, for RMSSD in the sitting condition, against the criterion device."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6713325",
   "metadata": {},
   "outputs": [],
   "source": [
    "icc_data[\"heartmath\"][\"rmssd\"][\"sitting\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b79f4c8",
   "metadata": {},
   "source": [
    "### Adjusting P-values with Bonferonni (Optional)\n",
    "\n",
    "Similar to the regression analysis, we can also correct the p-values for multiple testing with Bonferonni. The `bonferroni_correction_icc` function applies the Bonferroni correction to p-values in the ICC analysis results.It returns a dictionary similar to the input, but with Bonferroni corrected p-values. Additionally, it includes a flag indicating significance based on the corrected alpha level.\n",
    "\n",
    "**Note**: Please be aware that the number of tests performed is based on the \"features\" array. Therefore, be mindful of the number of variables you include in this array. If there are non-relevant variables present, redefine your \"features\" list, then rerun the `icc_analysis` function and then correct for Bonferroni. For instance, you can only keep a few relevant features: `features = [\"rmssd\", \"hf\",'mean_hr']`\n",
    "\n",
    "**It modifies the `icc_data` data frame in place**, and adds `corrected_pval` and `is_significant` to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d75d3300",
   "metadata": {},
   "outputs": [],
   "source": [
    "wearablehrv.group.bonferroni_correction_icc (icc_data, alpha = 0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5316b9d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "icc_data[\"heartmath\"][\"rmssd\"][\"sitting\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52488ff1",
   "metadata": {},
   "source": [
    "\n",
    "The `icc_plot` function provides a visual representation of the ICC values in the form of an interactive heatmap. The heatmap displays ICC values for each device (compared to the criterion device) across all conditions for a selected feature. Additionally, the heatmap annotations provide the 95% confidence intervals for each ICC value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc8c337b",
   "metadata": {},
   "outputs": [],
   "source": [
    "wearablehrv.group.icc_plot (icc_data, conditions, devices, features, font_size=9, cmap=\"coolwarm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c16cc8c",
   "metadata": {},
   "source": [
    "## 4. Bland-Altman Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b59b015d",
   "metadata": {},
   "source": [
    "So far, with ICC and regression analysis, we could establish the linear strength between a device and our criterion in a condition. Bland-Altman analysis offers a valuable feature: we can calculate the bias by subtracting a value in a device from the counterpart in the criterion (e.g., RMSSD value in Empatica in biking condition minus RMSSD value in VU in biking condition), and calculate the limits of agreement by taking the mean difference ± 1.96 times the standard deviation of the differences. This tells us a lot about whether there is something systematic in a device going on (e.g., if a device is overestimating values all the time compared to the criterion device).\n",
    "\n",
    "The `blandaltman_analysis` performs Bland-Altman analysis for each device against the criterion device, for all features, and conditions. It calculates the bias, standard deviation (SD), and limits of agreement (LOA) for each combination. If `save_as_csv = True`, the output will be saved as a .csv file at your specified path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adb228be",
   "metadata": {},
   "outputs": [],
   "source": [
    "blandaltman_data = wearablehrv.group.blandaltman_analysis (data, criterion, devices, conditions, features, path, save_as_csv=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0db8971",
   "metadata": {},
   "source": [
    "For instance, you can see empatica on average is overestimating the RMSSD values with 2.74 ms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ae3477d",
   "metadata": {},
   "outputs": [],
   "source": [
    "blandaltman_data[\"empatica\"][\"rmssd\"][\"biking\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6e38b1a",
   "metadata": {},
   "source": [
    "\n",
    "The `blandaltman_plot` creates a plot for each combination of device and condition agianst the criterion, with each plot showing the Bland-Altman analysis for the selected feature in a chosen device. It shows the difference between the measurements from the two devices on the y-axis and the average of the measurements on the x-axis. The LOA are shown as dashed lines on the plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5203d0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "wearablehrv.group.blandaltman_plot (data, criterion, conditions, devices, features, \n",
    "                         width=10, height_per_plot=5, agreement_bound=1.96, confidenceInterval=95, \n",
    "                         percentage=False, mean_diff_color='#FF6347', boundary_color='#20B2AA', pointColour='#8B008B', shade=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90636283",
   "metadata": {},
   "source": [
    "That's it! At this point, you have successfully calculated four major statistical analyses common in validation studies that help you interpret the performance of your devices at the group level, against a criterion. There is one more notebook remaining, which illustrates how to create some nice descriptive plots of your data.\n",
    "\n",
    "## Next Steps\n",
    "\n",
    "You're now ready to move on to the next notebook examples. \n",
    "\n",
    "Continue by consulting: \n",
    "\n",
    "- [Descriptive plots for your group data](https://github.com/Aminsinichi/wearable-hrv/blob/master/docs/examples/group_pipeline/4.group_data_plotting.ipynb)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}