ELBW_analysis_ventilation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![alt text](./Cerny_logo_1.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of Cerny ventilation recordings\n",
    "\n",
    "Analysis of ventilator data obtained from ELBW (birth weight <1000 g) infants transferred during the first 24 hours of life among recordings `AL000001 - AL001100` and  `AT000001 - AT000818`\n",
    "\n",
    "#### Author: Dr Gusztav Belteki"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Import the required libraries and set options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import IPython\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import pickle\n",
    "\n",
    "from scipy import stats\n",
    "from pandas import Series, DataFrame\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "%matplotlib inline\n",
    "matplotlib.style.use('classic')\n",
    "matplotlib.rcParams['figure.facecolor'] = 'w'\n",
    "\n",
    "pd.set_option('display.max_rows', 500)\n",
    "pd.set_option('display.max_columns', 100)\n",
    "pd.set_option('mode.chained_assignment', None) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Python version: {}\".format(sys.version))\n",
    "print(\"pandas version: {}\".format(pd.__version__))\n",
    "print(\"matplotlib version: {}\".format(matplotlib.__version__))\n",
    "print(\"NumPy version: {}\".format(np.__version__))\n",
    "print(\"SciPy version: {}\".format(sp.__version__))\n",
    "print(\"IPython version: {}\".format(IPython.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. List and set the working directory and the directory to write out data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Topic of the Notebook which will also be the name of the subfolder containing results\n",
    "TOPIC = 'ELBW'\n",
    "\n",
    "# Path to project folder containing clinical data (current weights only) and for export of results\n",
    "PATH = os.path.join(os.sep, 'Users', 'guszti', 'Library', 'Mobile Documents', 'com~apple~CloudDocs', \n",
    "                            'Documents', 'Research', 'Ventilation')\n",
    "\n",
    "# Name of the external hard drive\n",
    "DRIVE = 'GUSZTI'\n",
    "\n",
    "# Data loaded from both directories\n",
    "DIR_READ_1 = os.path.join(os.sep, 'Volumes', DRIVE, 'data_dump', 'fabian')\n",
    "DIR_READ_2 = os.path.join(os.sep, 'Volumes', DRIVE, 'data_dump', 'fabian_new')\n",
    "\n",
    "# Results will be written in this folder\n",
    "DIR_WRITE =  os.path.join(os.sep, PATH, 'ventilation_fabian_combined', 'Analyses', TOPIC)\n",
    "os.makedirs(DIR_WRITE, exist_ok=True)\n",
    "\n",
    "# Images and raw data will be written on an external hard drive\n",
    "DATA_DUMP = os.path.join(os.sep, 'Volumes', DRIVE, 'data_dump', 'fabian_combined', TOPIC)\n",
    "os.makedirs(DATA_DUMP, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DIR_READ_1, DIR_READ_2, DIR_WRITE, DATA_DUMP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Import clinical data, ventilation modes and blood gases of the final dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These data were exported by the `ELBW_analysis_clinical` Notebook\n",
    "\n",
    "# Computationally collected data\n",
    "with open(os.path.join(DATA_DUMP, 'clin_data_comput_ELBW.pickle'), 'rb') as handle:\n",
    "    clin_data_comput_ELBW = pickle.load(handle)\n",
    "\n",
    "with open(os.path.join(DATA_DUMP, 'ventilation_modes_ELBW.pickle'), 'rb') as handle:\n",
    "    vent_modes_ELBW = pickle.load(handle)\n",
    "\n",
    "with open(os.path.join(DATA_DUMP, 'blood_gases_comput.pickle'), 'rb') as handle:\n",
    "    blood_gases_comput = pickle.load(handle)\n",
    "\n",
    "# Manually collected data\n",
    "# Demographic, clinical and outcome data\n",
    "with open(os.path.join(DATA_DUMP, 'clin_data_combined.pickle'), 'rb') as handle:\n",
    "    clinical_data_combined = pickle.load(handle)\n",
    "\n",
    "# Vital parameters for the final dataset\n",
    "with open(os.path.join(DATA_DUMP, 'vital_parameters_manual.pickle'), 'rb') as handle:\n",
    "    vital_parameters = pickle.load(handle)\n",
    "    \n",
    "# Blood gases for the final dataset\n",
    "with open(os.path.join(DATA_DUMP, 'blood_gases_manual.pickle'), 'rb') as handle:\n",
    "    blood_gases_manual = pickle.load(handle)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Import ventilator data from pickle archives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. Import data on ventilated recordings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(DIR_READ_1, 'data_pars_measurements_ventilated_1_1100.pickle'), 'rb') as handle:\n",
    "    data_pars_measurements_1_1100 = pickle.load(handle)\n",
    "with open(os.path.join(DIR_READ_1, 'data_pars_settings_ventilated_1_1100.pickle'), 'rb') as handle:\n",
    "    data_pars_settings_1_1100 = pickle.load(handle)\n",
    "   \n",
    "with open(os.path.join(DIR_READ_2, 'data_pars_measurements_ventilated_new_1_1305.pickle'), 'rb') as handle:\n",
    "    data_pars_measurements_new_1_1305 = pickle.load(handle)\n",
    "with open(os.path.join(DIR_READ_2, 'data_pars_settings_ventilated_new_1_1305.pickle'), 'rb') as handle:\n",
    "    data_pars_settings_new_1_1305 = pickle.load(handle)\n",
    "\n",
    "data_pars_measurements_ventilated = {**data_pars_measurements_1_1100, **data_pars_measurements_new_1_1305}\n",
    "data_pars_settings_ventilated =     {**data_pars_settings_1_1100, **data_pars_settings_new_1_1305}\n",
    "\n",
    "len(data_pars_measurements_ventilated), len(data_pars_settings_ventilated)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. Import data on noninvasive recordings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(DIR_READ_1, 'data_pars_measurements_noninvasive_1_300.pickle'), 'rb') as handle:\n",
    "    data_pars_measurements_1_300 = pickle.load(handle)\n",
    "with open(os.path.join(DIR_READ_1, 'data_pars_settings_noninvasive_1_300.pickle'), 'rb') as handle:\n",
    "    data_pars_settings_1_300 = pickle.load(handle)\n",
    "\n",
    "with open(os.path.join(DIR_READ_1, 'data_pars_measurements_noninvasive_301_600.pickle'), 'rb') as handle:\n",
    "    data_pars_measurements_301_600 = pickle.load(handle)\n",
    "with open(os.path.join(DIR_READ_1, 'data_pars_settings_noninvasive_301_600.pickle'), 'rb') as handle:\n",
    "    data_pars_settings_301_600 = pickle.load(handle)\n",
    "\n",
    "with open(os.path.join(DIR_READ_1, 'data_pars_measurements_noninvasive_601_900.pickle'), 'rb') as handle:\n",
    "    data_pars_measurements_601_900 = pickle.load(handle)\n",
    "with open(os.path.join(DIR_READ_1, 'data_pars_settings_noninvasive_601_900.pickle'), 'rb') as handle:\n",
    "    data_pars_settings_601_900 = pickle.load(handle)\n",
    "\n",
    "with open(os.path.join(DIR_READ_1, 'data_pars_measurements_noninvasive_901_1100.pickle'), 'rb') as handle:\n",
    "    data_pars_measurements_901_1100 = pickle.load(handle)\n",
    "with open(os.path.join(DIR_READ_1, 'data_pars_settings_noninvasive_901_1100.pickle'), 'rb') as handle:\n",
    "    data_pars_settings_901_1100 = pickle.load(handle)\n",
    "   \n",
    "with open(os.path.join(DIR_READ_2, 'data_pars_measurements_noninvasive_new_1_1305.pickle'), 'rb') as handle:\n",
    "    data_pars_measurements_new_1_1305 = pickle.load(handle)\n",
    "with open(os.path.join(DIR_READ_2, 'data_pars_settings_noninvasive_new_1_1305.pickle'), 'rb') as handle:\n",
    "    data_pars_settings_new_1_1305 = pickle.load(handle)\n",
    "\n",
    "data_pars_measurements_noninvasive = {**data_pars_measurements_1_300, **data_pars_measurements_301_600, \n",
    "                                      **data_pars_measurements_601_900, **data_pars_measurements_901_1100,\n",
    "                                      **data_pars_measurements_new_1_1305}\n",
    "data_pars_settings_noninvasive =     {**data_pars_settings_1_300, **data_pars_settings_301_600,\n",
    "                                      **data_pars_settings_601_900, **data_pars_settings_901_1100,\n",
    "                                      **data_pars_settings_new_1_1305}\n",
    "\n",
    "len(data_pars_measurements_noninvasive), len(data_pars_settings_noninvasive)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Limit ventilator data to the included cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_pars_measurements_ventilated = {rec: data_pars_measurements_ventilated[rec] \n",
    "                for rec in data_pars_measurements_ventilated if rec in clin_data_comput_ELBW.index}\n",
    "data_pars_settings_ventilated = {rec: data_pars_settings_ventilated[rec] \n",
    "                for rec in data_pars_settings_ventilated if rec in clin_data_comput_ELBW.index}\n",
    "\n",
    "len(data_pars_measurements_ventilated), len(data_pars_settings_ventilated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_pars_measurements_noninvasive = {rec: data_pars_measurements_noninvasive[rec] \n",
    "                for rec in data_pars_measurements_noninvasive if rec in clin_data_comput_ELBW.index}\n",
    "data_pars_settings_noninvasive = {rec: data_pars_settings_noninvasive[rec] \n",
    "                for rec in data_pars_settings_noninvasive if rec in clin_data_comput_ELBW.index}\n",
    "\n",
    "len(data_pars_measurements_noninvasive), len(data_pars_settings_noninvasive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(data_pars_measurements_ventilated.keys()) & set(data_pars_measurements_noninvasive.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_pars_measurements = {**data_pars_measurements_ventilated, **data_pars_measurements_noninvasive}\n",
    "data_pars_settings = {**data_pars_settings_ventilated, **data_pars_settings_noninvasive}\n",
    "len(data_pars_measurements), len(data_pars_settings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Keep only selected parameters and correct data types when needed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. Keep only parameters and settings relevant for the study"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters_to_choose_from = {'PIP', 'MAP', 'PEEP',  'MVresp', 'Leak', 'RR', 'FiO2', 'Flow_exp', 'MV_kg', \n",
    "                             'VTimand_kg', 'VTemand_kg', 'VTespon_pat_kg', 'VTemand_resp_kg', }\n",
    "\n",
    "settings_to_choose_from = {'Ventilator_mode', 'PIP_set', 'PEEP_set', 'FiO2_set', 'Flow_insp_set', 'Flow_exp_set', \n",
    "                           'Ti_set', 'Te_set', 'RR_set', 'IE_I_set', 'IE_E_set', 'VG_set_kg', 'VG_state', }\n",
    "\n",
    "for recording, dta in data_pars_measurements.items():\n",
    "    parameters_to_keep = sorted(set(dta.columns) & parameters_to_choose_from)\n",
    "    data_pars_measurements[recording] = dta[parameters_to_keep]\n",
    "    \n",
    "for recording, dta in data_pars_settings.items():\n",
    "    settings_to_keep = sorted(set(dta.columns) & settings_to_choose_from)\n",
    "    data_pars_settings[recording] = dta[settings_to_keep]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. Correct data types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All measured parameters are floats\n",
    "for patient in data_pars_measurements:\n",
    "    data_pars_measurements[patient] = data_pars_measurements[patient].astype('float')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Settings are floats or categorical\n",
    "\n",
    "categorical_settings = ['Ventilator_mode',  'VG_state',]\n",
    "\n",
    "float_settings = ['PIP_set', 'PEEP_set', 'FiO2_set', 'Flow_insp_set', 'Flow_exp_set', \n",
    "                  'Ti_set', 'Te_set', 'RR_set', 'IE_I_set', 'IE_E_set', 'VG_set_kg',]\n",
    "\n",
    "for recording in data_pars_settings:\n",
    "    for par in float_settings:\n",
    "        if par in data_pars_settings[recording]:\n",
    "            data_pars_settings[recording][par] = data_pars_settings[recording][par].astype('float')\n",
    "    for par in categorical_settings:\n",
    "        if par in data_pars_settings[recording]:\n",
    "            data_pars_settings[recording][par] = data_pars_settings[recording][par].astype('category')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Combine parameters and settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_all = {}\n",
    "\n",
    "for case in data_pars_measurements.keys():\n",
    "    vent_pars = data_pars_measurements[case]\n",
    "    vent_settings = data_pars_settings[case]\n",
    "    data_all[case] = pd.concat([vent_pars, vent_settings], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Cleanup ventilator data by removing outliers due to artifacts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify and remove abnormal inflations (presumaby ventilation artifacts)\n",
    "\n",
    "data_points_before_cleanup = {}\n",
    "for patient in data_all:\n",
    "    data_points_before_cleanup[patient] = len(data_all[patient])\n",
    "data_points_before_cleanup;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. Inflations with negative inflating pressure\n",
    "These are artifacts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pinfl_negative = {}\n",
    "for patient in data_all:\n",
    "    if 'PIP'in data_all[patient].columns and 'PEEP' in data_all[patient].columns:\n",
    "        Pinfl_negative[patient] = sum(data_all[patient]['PIP'] < data_all[patient]['PEEP'])\n",
    "Pinfl_negative;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many occurrence in the whole dataset\n",
    "sum(n for n in Pinfl_negative.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. Inflations with VTemand_kg, VTemand_resp_kg or VTespon_respo_kg >25 mL/kg\n",
    "\n",
    "Expired VT > 25 mL/kg is likely artifact, representing a flow sensor problem or an open circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VTemand_kg_over25 = {}\n",
    "for patient in data_all:\n",
    "    if 'VTemand_kg' in data_all[patient].columns:\n",
    "        VTemand_kg_over25[patient] = (sum(data_all[patient]['VTemand_kg'] > 25))\n",
    "VTemand_kg_over25;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many occurrence in the whole dataset\n",
    "sum(n for n in VTemand_kg_over25.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VTemand_resp_kg_over25 = {}\n",
    "for patient in data_all:\n",
    "    if 'VTemand_resp_kg' in data_all[patient].columns:\n",
    "        VTemand_resp_kg_over25[patient] = (sum(data_all[patient]['VTemand_resp_kg'] > 25))\n",
    "VTemand_resp_kg_over25;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many occurrence in the whole dataset\n",
    "sum(n for n in VTemand_resp_kg_over25.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VTespon_pat_kg_over25 = {}\n",
    "for patient in data_all:\n",
    "    if 'VTespon_pat_kg' in data_all[patient].columns:\n",
    "        VTespon_pat_kg_over25[patient] = (sum(data_all[patient]['VTespon_pat_kg'] > 25))\n",
    "VTespon_pat_kg_over25;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many occurrence in the whole dataset\n",
    "sum(n for n in VTespon_pat_kg_over25.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Inflations with MV over 1 L/min/kg\n",
    "\n",
    "This will not occur during conventional ventilation, unless due to an artifact (eg circuit is open). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MV_over_1 = {}\n",
    "for patient in data_all:\n",
    "    if 'MV_kg' in data_all[patient].columns:\n",
    "        MV_over_1[patient] = sum(data_all[patient]['MV_kg'] > 1) \n",
    "MV_over_1;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many occurrence in the whole dataset\n",
    "sum(n for n in MV_over_1.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### D. Inflations with RR  > 130/min\n",
    "\n",
    "These values are likely abnormal representing autocycling due to condensed water in the circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RR_over_130 = {}\n",
    "for patient in data_all:\n",
    "    if 'RR' in data_all[patient].columns:\n",
    "        RR_over_130[patient] = sum(data_all[patient]['RR'] > 130)\n",
    "RR_over_130;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many occurrence in the whole dataset\n",
    "sum(n for n in RR_over_130.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### E. Remove these anomalies from the dataset\n",
    "\n",
    "Keep rows with missing data of the parameters because other parameters at the same time point may be present and valid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "for patient in data_all:\n",
    "    #print(patient)\n",
    "    dta = data_all[patient]\n",
    "    if 'PIP' in dta and 'PEEP' in dta:\n",
    "        dta = dta[dta['PIP'].isna() | dta['PEEP'].isna() | (dta['PIP'] >= dta['PEEP'])]\n",
    "    if 'VTemand_kg' in dta:\n",
    "        dta = dta[dta['VTemand_kg'].isna() | (dta['VTemand_kg'] <= 25)]\n",
    "    if 'VTemand_resp_kg' in dta:\n",
    "        dta = dta[dta['VTemand_resp_kg'].isna() | (dta['VTemand_resp_kg'] <= 25)]\n",
    "    if 'VTespon_pat_kg' in dta:\n",
    "        dta = dta[dta['VTespon_pat_kg'].isna() | (dta['VTespon_pat_kg'] <= 25)]\n",
    "    if 'MV_kg' in dta:\n",
    "        dta = dta[dta['MV_kg'].isna() | (dta['MV_kg'] <= 1)]\n",
    "    if 'RR' in dta:\n",
    "        dta = dta[dta['RR'].isna() | (dta['RR'] <= 130)]\n",
    "    data_all[patient] = dta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_points_after_cleanup = {}\n",
    "for patient in data_all:\n",
    "    data_points_after_cleanup[patient] = len(data_all[patient])\n",
    "data_points_after_cleanup;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cleanup_frme = DataFrame([data_points_before_cleanup, data_points_after_cleanup, Pinfl_negative, \n",
    "    VTemand_kg_over25, VTemand_resp_kg_over25, VTespon_pat_kg_over25, MV_over_1, RR_over_130]).T\n",
    "cleanup_frme.columns =   ['data_points_before_cleanup', 'data_points_after_cleanup', 'Pinfl_negative', \n",
    "    'VTemand_kg_over25', 'VTemand_resp_kg_over25', 'VTespon_pat_kg_over25', 'MV_over_1', 'RR_over_130',]\n",
    "cleanup_frme['rows_removed'] = cleanup_frme['data_points_before_cleanup'] - \\\n",
    "                                cleanup_frme['data_points_after_cleanup']\n",
    "cleanup_frme.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cleanup_frme.sum()/ cleanup_frme['data_points_before_cleanup'].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cleanup_frme.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cleanup_frme.to_csv(os.path.join(DIR_WRITE, f'data_cleanup.csv'))\n",
    "\n",
    "writer = pd.ExcelWriter(os.path.join(DIR_WRITE, f'data_cleanup.xlsx'))\n",
    "cleanup_frme.to_excel(writer, 'cleanup_frme')\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Calculate ventilator parameters which are interpreted differently for SIPPV and SIMV "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. Calculate the expired VT of ventilator inflations (VTe_kg) as appropriate for each ventilator mode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for patient, dta in data_all.items():\n",
    "    # Skip recordings which had only nCPAP\n",
    "    if 'VTemand_kg' not in dta.columns:\n",
    "        continue\n",
    "    \n",
    "    # For SIMV / SIMVPSV recordings take VTe_kg from VTemand_resp_kg, for SIPPV from VTemand_kg\n",
    "    if 'SIMV' in dta['Ventilator_mode'].unique() or 'SIMVPSV' in dta['Ventilator_mode'].unique():\n",
    "        #print(patient)\n",
    "        dta['VTe_kg'] = np.where(dta['Ventilator_mode'].isin(['SIMV', 'SIMVPSV']), \n",
    "            dta['VTemand_resp_kg'], dta['VTemand_kg']) \n",
    "    else:\n",
    "        dta['VTe_kg'] = dta['VTemand_kg']\n",
    "    \n",
    "    data_all[patient] = dta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. Calculate the number of mandatory ventilator inflations "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for patient, dta in data_all.items():\n",
    "    # Skip recordings which had only nCPAP\n",
    "    if 'VTemand_kg' not in dta.columns:\n",
    "        continue\n",
    "    \n",
    "    # For SIMV / SIMVPSV recordings take RRmand from RR_set, for SIPPV from RR\n",
    "    if 'SIMV' in dta['Ventilator_mode'].unique() or 'SIMVPSV' in dta['Ventilator_mode'].unique():\n",
    "        if 'SIPPV' in dta['Ventilator_mode'].unique():\n",
    "            dta['RRmand'] = np.where(dta['Ventilator_mode'].isin(['SIMV', 'SIMVPSV']), \n",
    "                dta['RR_set'], dta['RR']) \n",
    "        else:\n",
    "            dta['RRmand'] = dta['RR_set']    \n",
    "    else:\n",
    "        dta['RRmand'] = dta['RR']\n",
    "    \n",
    "    data_all[patient] = dta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Descriptive statistics in the whole dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. Statistics on individual cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "percentiles = [0.05, 0.25, 0.5, 0.75, 0.95]\n",
    "index_names = ['data_points', 'mean', 'SD', 'min', '5pc', '25pc', 'median', '75pc', '95pc', 'max']\n",
    "\n",
    "stats_all_patients = {} \n",
    "\n",
    "for recording in data_all:\n",
    "    stats_all_patients[recording] = round(data_all[recording].describe(percentiles = percentiles), 2)\n",
    "    stats_all_patients[recording].index = index_names\n",
    "    stats_all_patients[recording].dropna(how='all', subset=index_names[1:], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create table with statistics for all cases and all relevant parameters\n",
    "stats_all_patients_combined = pd.concat(stats_all_patients, axis = 1).T\n",
    "\n",
    "stats_all_patients_combined.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats_all_patients_combined"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. Statistics on selected individual parameters and setting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Selected individual parameters and settings\n",
    "\n",
    "pars_selected = ['FiO2_set', 'MAP', 'PIP', 'PEEP', 'VTe_kg', 'RRmand', 'MV_kg', 'Ti_set', 'Te_set', 'Leak']\n",
    "\n",
    "stats_all_parameters = {}\n",
    "\n",
    "for par in pars_selected:\n",
    "    stats_all_parameters[par] = stats_all_patients_combined.swaplevel(0,1).loc[par].sort_values('mean', ascending = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Unstack table to create table with different configuration\n",
    "stats_all_parameters_combined = pd.concat(stats_all_parameters, axis=1)\n",
    "\n",
    "stats_all_parameters_combined.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Group statistics on selected parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars_for_mean = ['PIP', 'MAP', 'PEEP', 'MV_kg', 'VTe_kg', 'RRmand',]\n",
    "\n",
    "group_stats_all_mean = \\\n",
    "    round(stats_all_parameters_combined.swaplevel(axis=1)['mean'].describe()[pars_for_mean].T, 2)\n",
    "group_stats_all_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars_for_median = ['FiO2_set', 'Leak', 'Ti_set', 'Te_set']\n",
    "\n",
    "group_stats_all_median = \\\n",
    "    round(stats_all_parameters_combined.swaplevel(axis=1)['median'].describe()[pars_for_median].T, 2)\n",
    "group_stats_all_median"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Export statistics to a multisheet Excel file and pickle archive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save statistics into Excel file\n",
    "writer = pd.ExcelWriter(os.path.join(DIR_WRITE, 'stats_all_patients_ELBW.xlsx'))\n",
    "for patient in stats_all_patients:\n",
    "    stats_all_patients[patient].to_excel(writer, patient)\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save statistics into Excel file\n",
    "writer = pd.ExcelWriter(os.path.join(DIR_WRITE, 'stats_all_parameters_ELBW.xlsx'))\n",
    "for parameter in stats_all_parameters:\n",
    "    stats_all_parameters[parameter].to_excel(writer, parameter)\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save statistics into Excel file\n",
    "writer = pd.ExcelWriter(os.path.join(DIR_WRITE, 'stats_all_combined_ELBW.xlsx'))\n",
    "stats_all_patients_combined.to_excel(writer, 'patients')\n",
    "stats_all_parameters_combined.to_excel(writer, 'parameters')\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save group statistics into Excel file\n",
    "writer = pd.ExcelWriter(os.path.join(DIR_WRITE, 'group_stats_all_ELBW.xlsx'))\n",
    "group_stats_all_mean.to_excel(writer, 'parametric')\n",
    "group_stats_all_median.to_excel(writer, 'non_parametric')\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9. Combine data from all recordings, and export overall statistics about it as an Excel file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As this includes all data with various ventilator modes, some parameters should be interpreted with caution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_all_combined = pd.concat(data_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_all_combined.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter(os.path.join(DIR_WRITE, f'descriptive_stats_combined_ELBW.xlsx'))\n",
    "cleanup_frme.to_excel(writer, 'stats_all_combined')\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10. Analyse ventilator modes used"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. Ventilator modes\n",
    "\n",
    "Some modes do not occur in this dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vent_modes_ELBW.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Limit vent mode table to modes which did occur\n",
    "vent_modes_ELBW = vent_modes_ELBW[['SIMV', 'SIMVPSV', 'SIPPV', 'NCPAP', 'total', 'VG_on',]]\n",
    "vent_modes_ELBW.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. Assign recordings where >95% is one mode to that mode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predominantly_sippv = set(vent_modes_ELBW[vent_modes_ELBW['SIPPV'] > vent_modes_ELBW['total'] * 0.95].index)\n",
    "len(predominantly_sippv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vent_modes_ELBW.loc[sorted(predominantly_sippv)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predominantly_simv = set(vent_modes_ELBW[vent_modes_ELBW['SIMV'] > vent_modes_ELBW['total'] * 0.95].index)\n",
    "len(predominantly_simv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vent_modes_ELBW.loc[sorted(predominantly_simv)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predominantly_simv_ps = set(vent_modes_ELBW[vent_modes_ELBW['SIMVPSV'] > vent_modes_ELBW['total'] * 0.95].index)\n",
    "len(predominantly_simv_ps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vent_modes_ELBW.loc[sorted(predominantly_simv_ps)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# They are all fully nCPAP as these babies were not intubated\n",
    "predominantly_ncpap = set(vent_modes_ELBW[vent_modes_ELBW['NCPAP'] > vent_modes_ELBW['total'] * 0.95].index)\n",
    "len(predominantly_ncpap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vent_modes_ELBW.loc[sorted(predominantly_ncpap)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_predominant_mode = set(vent_modes_ELBW.index) - (predominantly_simv | predominantly_simv_ps |\n",
    "                          predominantly_sippv |  predominantly_ncpap)\n",
    "vent_modes_ELBW.loc[sorted(no_predominant_mode)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "both_simv_and_sippv = no_predominant_mode.copy()\n",
    "# Exclude recordings which have SIMV and SIMVPSV\n",
    "both_simv_and_sippv.discard('AL000152')\n",
    "vent_modes_ELBW.loc[sorted(both_simv_and_sippv)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exclude recordings which have SIMV and SIPPV\n",
    "both_simv_and_simvpsv = no_predominant_mode - both_simv_and_sippv\n",
    "vent_modes_ELBW.loc[sorted(both_simv_and_simvpsv)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Assign those ones with >95% VG or noVG to those modes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predominantly_vg = set(vent_modes_ELBW[vent_modes_ELBW['VG_on'] > vent_modes_ELBW['total'] * 0.95].index)\n",
    "vent_modes_ELBW.loc[sorted(predominantly_vg)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(predominantly_vg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predominantly_novg = set(vent_modes_ELBW[vent_modes_ELBW['VG_on'] < vent_modes_ELBW['total'] * 0.05].index) & \\\n",
    "                     set(vent_modes_ELBW[vent_modes_ELBW['NCPAP'] < vent_modes_ELBW['total'] * 0.05].index)\n",
    "vent_modes_ELBW.loc[sorted(predominantly_novg)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(predominantly_novg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "both_vg_novg = set(vent_modes_ELBW[ (vent_modes_ELBW['VG_on'] >= vent_modes_ELBW['total'] * 0.05) & \n",
    "                  (vent_modes_ELBW['VG_on'] <= vent_modes_ELBW['total'] * 0.95)].index)\n",
    "vent_modes_ELBW.loc[sorted(both_vg_novg)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Save tables with ventilator modes into Excel file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter(os.path.join(DIR_WRITE, 'vent_modes_ELBW_categories.xlsx'))\n",
    "vent_modes_ELBW.to_excel(writer, 'all')\n",
    "vent_modes_ELBW.loc[sorted(predominantly_simv)].to_excel(writer, 'predominantly_simv')\n",
    "vent_modes_ELBW.loc[sorted(predominantly_simv_ps)].to_excel(writer, 'predominantly_simv_ps')\n",
    "vent_modes_ELBW.loc[sorted(predominantly_sippv)].to_excel(writer, 'predominantly_sippv')\n",
    "vent_modes_ELBW.loc[sorted(predominantly_ncpap)].to_excel(writer, 'predominantly_ncpap')\n",
    "vent_modes_ELBW.loc[sorted(no_predominant_mode)].to_excel(writer, 'no_predominant_mode')\n",
    "vent_modes_ELBW.loc[sorted(both_simv_and_sippv)].to_excel(writer, 'both_simv_and_sippv')\n",
    "vent_modes_ELBW.loc[sorted(both_simv_and_simvpsv)].to_excel(writer, 'both_simv_and_simvpsv')\n",
    "\n",
    "vent_modes_ELBW.loc[sorted(predominantly_vg)].to_excel(writer, 'predominantly_vg')\n",
    "vent_modes_ELBW.loc[sorted(predominantly_novg)].to_excel(writer, 'predominantly_novg')\n",
    "vent_modes_ELBW.loc[sorted(both_vg_novg)].to_excel(writer, 'both_vg_novg')\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11. Limit data to the given ventilator modes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recordings with one mode only\n",
    "\n",
    "# SIPPV only\n",
    "data_sippv = {}\n",
    "for case in predominantly_sippv:\n",
    "    data_sippv[case] = data_all[case][data_all[case]['Ventilator_mode'] == 'SIPPV']\n",
    "\n",
    "# SIMV only\n",
    "data_simv = {}\n",
    "for case in predominantly_simv:\n",
    "    data_simv[case] = data_all[case][data_all[case]['Ventilator_mode'] == 'SIMV']\n",
    "\n",
    "# NCPAP\n",
    "# None of these recordings has any other mode as the baby was not intubated\n",
    "data_ncpap = {}\n",
    "for case in predominantly_ncpap:\n",
    "    data_ncpap[case] = data_all[case][data_all[case]['Ventilator_mode']== 'NCPAP']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recordings with more than one mode\n",
    "\n",
    "# SIPPV and SIMV in the same recording\n",
    "data_both_sippv_simv = {}\n",
    "for case in both_simv_and_sippv:\n",
    "    data_both_sippv_simv[case] = data_all[case][data_all[case]['Ventilator_mode'].isin(['SIMV', 'SIPPV'])]\n",
    "    \n",
    "# SIMV and SIMVPSV in the same recording\n",
    "data_both_simv_simvpsv = {}\n",
    "for case in both_simv_and_simvpsv:\n",
    "    data_both_simv_simvpsv[case] = data_all[case][data_all[case]['Ventilator_mode'].isin(['SIMV', 'SIMVPSV'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# VG only\n",
    "data_vg = {}\n",
    "for case in predominantly_vg:\n",
    "    data_vg[case] = data_all[case][data_all[case]['VG_state'] == 'on']\n",
    "\n",
    "# no VG only\n",
    "# None of these recordings had any VG periods therefore (see above table), therefore there is nothing to remove\n",
    "data_novg = {}\n",
    "for case in predominantly_novg:\n",
    "    data_novg[case] = data_all[case]\n",
    "\n",
    "# both VG and no VG\n",
    "data_both_vg_novg = {}\n",
    "for case in both_vg_novg:\n",
    "    data_both_vg_novg[case] = data_all[case]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 12. Further process ventilator data specific to different ventilator modes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. Recordings with SIPPV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_sippv_all = pd.concat(data_sippv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For SIPPV MVresp should be n/a as it is always 100%\n",
    "data_sippv_all = data_sippv_all.drop('MVresp', axis=1)\n",
    "\n",
    "# During SIPPV VTemand_resp_kg and VTespon_resp_kg are not valid\n",
    "data_sippv_all = data_sippv_all.drop(['VTemand_resp_kg', 'VTespon_pat_kg'], axis=1)\n",
    "\n",
    "# Drop Flow_exp as no data\n",
    "data_sippv_all = data_sippv_all.drop('Flow_exp', axis=1)\n",
    "\n",
    "# RR_diff is the difference between the set and actual ventilator rate\n",
    "data_sippv_all['RR_diff'] = data_sippv_all['RRmand'] - data_sippv_all['RR_set'] \n",
    "\n",
    "data_sippv_all.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. Recordings with SIMV "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_simv_all = pd.concat(data_simv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop Flow_exp as no data\n",
    "data_simv_all = data_simv_all.drop('Flow_exp', axis=1)\n",
    "\n",
    "data_simv_all.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Recordings with both SIPPV and SIMV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_both_sippv_simv_all = pd.concat(data_both_sippv_simv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop Flow_exp as no data\n",
    "data_both_sippv_simv_all = data_both_sippv_simv_all.drop('Flow_exp', axis=1)\n",
    "\n",
    "data_both_sippv_simv_all.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### D. Recordings with both SIMV and SIMVPSV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_both_simv_simvpsv_all = pd.concat(data_both_simv_simvpsv)\n",
    "\n",
    "data_both_simv_simvpsv_all.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### E. Recordings with nasal CPAP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_ncpap_all = pd.concat(data_ncpap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_keep = ['FiO2', 'Flow_exp', 'MAP', 'FiO2_set', 'PEEP_set', 'PIP_set', 'Ventilator_mode',]\n",
    "data_ncpap_all = data_ncpap_all[columns_to_keep]\n",
    "\n",
    "data_ncpap_all.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### F. Recordings with VG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For SIMV and SIMVPSV recordings `VTemand_resp_kg` is to be used as here there can be spontaneous breaths between ventilator inflation whose tidal volume is not targeted, for SIPPV use `VTemand_kg`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for patient, dta in data_vg.items():\n",
    "    \n",
    "    # The difference between the actual and the target expired VT\n",
    "    dta['VT_diff_kg'] = dta['VTe_kg'] - dta['VG_set_kg']\n",
    "    # Calculate absolute value of the difference\n",
    "    dta['VT_diff_kg_abs'] = np.abs(dta['VT_diff_kg'])\n",
    "\n",
    "    dta['P_diff'] = dta['PIP_set'] - dta['PIP']\n",
    "    # Calculate absolute value of the difference\n",
    "    dta['P_diff_abs'] = np.abs(dta['P_diff'])\n",
    "    \n",
    "    data_vg[patient] = dta\n",
    "    data_vg_all = pd.concat(data_vg)\n",
    "\n",
    "# Drop Flow_exp as no data\n",
    "data_vg_all = data_vg_all.drop('Flow_exp', axis=1)\n",
    "\n",
    "data_vg_all.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### G. Recordings without VG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for patient, dta in data_novg.items():\n",
    "    \n",
    "    data_novg[patient] = dta\n",
    "\n",
    "data_novg_all = pd.concat(data_novg)\n",
    "\n",
    "# Drop Flow_exp as no data and 'VGset_kg' as there is no target VT here\n",
    "data_novg_all = data_novg_all.drop(['Flow_exp', 'VG_set_kg'], axis=1)\n",
    "\n",
    "data_novg_all.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### H. Recordings with both VG and non-VG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for patient, dta in data_both_vg_novg.items():\n",
    "        \n",
    "    # The difference between the actual and the target expired VT\n",
    "    dta['VT_diff_kg'] = dta['VTe_kg'] - dta['VG_set_kg']\n",
    "    # Calculate absolute value of the difference\n",
    "    dta['VT_diff_kg_abs'] = np.abs(dta['VT_diff_kg'])\n",
    "\n",
    "    dta['P_diff'] = dta['PIP_set'] - dta['PIP']\n",
    "    # Calculate absolute value of the difference\n",
    "    dta['P_diff_abs'] = np.abs(dta['P_diff'])\n",
    "    \n",
    "    data_both_vg_novg[patient] = dta\n",
    "    \n",
    "    data_both_vg_novg[patient] = dta\n",
    "    data_both_vg_novg_all = pd.concat(data_both_vg_novg)\n",
    "\n",
    "# Drop Flow_exp as no data and 'VGset_kg' as there is no target VT here\n",
    "# Also, there is no SIMV, so remove SIMV specific parameters\n",
    "columns_to_drop = ['Flow_exp', 'MVresp', 'VTemand_resp_kg','VTespon_pat_kg']\n",
    "data_both_vg_novg_all = data_both_vg_novg_all.drop(columns_to_drop, axis=1)\n",
    "\n",
    "data_both_vg_novg_all.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### I. Save descriptive statistics about the subsets with different ventilator modes in Excel file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save statistics into Excel file\n",
    "writer = pd.ExcelWriter(os.path.join(DIR_WRITE, 'descriptive_stats_vent_modes_ELBW.xlsx'))\n",
    "data_sippv_all.describe().to_excel(writer, 'sippv')\n",
    "data_simv_all.describe().to_excel(writer, 'simv')\n",
    "data_both_sippv_simv_all.describe().to_excel(writer, 'both_sippv_simv')\n",
    "data_both_simv_simvpsv_all.describe().to_excel(writer, 'both_simv_simvpsv')\n",
    "\n",
    "data_ncpap_all.describe().to_excel(writer, 'ncpap')\n",
    "\n",
    "data_vg_all.describe().to_excel(writer, 'vg')\n",
    "data_novg_all.describe().to_excel(writer, 'novg')\n",
    "data_both_vg_novg_all.describe().to_excel(writer, 'both_vg_novg')\n",
    "\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13. Compare SIPPV and SIMV in different recordings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clin_data_comput_ELBW.loc[sorted(predominantly_sippv)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clin_data_comput_ELBW.loc[sorted(predominantly_simv)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for par in ['Gestational Age (weeks)_comput', 'Birth Weight_comput']:\n",
    "    group_1 = clin_data_comput_ELBW.loc[sorted(predominantly_sippv)][par]\n",
    "    group_2 = clin_data_comput_ELBW.loc[sorted(predominantly_simv)][par]\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1, group_2,), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ventilator parameters with normal distribution\n",
    "pars_for_mean = ['PIP', 'MAP', 'PEEP', 'MV_kg', 'VTe_kg', 'VTimand_kg', 'VTemand_kg', 'RRmand', 'RR_set', ]\n",
    "# Ventilator parameters which are not normally distributed\n",
    "pars_for_median = ['FiO2_set', 'Leak', 'RR_diff',]\n",
    "\n",
    "sippv_stats_mean = round(data_sippv_all.groupby(level=0)[pars_for_mean].describe().swaplevel(axis=1)['mean'].describe().T, 1)\n",
    "sippv_stats_median = round(data_sippv_all.groupby(level=0)[pars_for_median].describe().swaplevel(axis=1)['50%'].describe().T, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sippv_stats_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sippv_stats_median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ventilator parameters with normal distribution\n",
    "pars_for_mean = ['PIP', 'MAP', 'PEEP', 'MV_kg', 'MVresp', 'VTe_kg', 'VTimand_kg', 'VTemand_kg', \n",
    "                 'VTemand_resp_kg', 'VTespon_pat_kg', 'RRmand', 'RR_set']\n",
    "# Ventilator parameters which are not normally distributed\n",
    "pars_for_median = ['FiO2_set', 'Leak',]\n",
    "\n",
    "simv_stats_mean = round(data_simv_all.groupby(level=0)[pars_for_mean].describe().swaplevel(axis=1)['mean'].describe().T, 2)\n",
    "simv_stats_median = round(data_simv_all.groupby(level=0)[pars_for_median].describe().swaplevel(axis=1)['50%'].describe().T, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simv_stats_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simv_stats_median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars = ['PIP', 'MAP', 'PEEP', 'MV_kg', 'VTe_kg', 'RRmand', 'RR_set', 'FiO2_set', 'Leak', ]\n",
    "# The aggregated values (means, medians) themselves do not show normal distribution in the groups, \n",
    "# use nonparametric measures and tests\n",
    "stts = ['count', 'min', '25%', '50%', '75%', 'max']\n",
    "\n",
    "sippv_stats = pd.concat([sippv_stats_mean[stts], sippv_stats_median[stts]])\n",
    "simv_stats = pd.concat([simv_stats_mean[stts], simv_stats_median[stts]])\n",
    "sippv_simv_stats = pd.merge(sippv_stats, simv_stats, left_index=True, right_index=True, \n",
    "                            suffixes= ['_sippv', '_simv']).sort_index(axis=1, ascending=False)\n",
    "sippv_simv_stats = sippv_simv_stats.loc[pars]\n",
    "sippv_simv_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Inferential statistics for selected parameters, SIPPV vs SIMV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mann-Whitney test\n",
    "\n",
    "for par in ['PIP', 'MAP', 'PEEP', 'MV_kg', 'VTe_kg', 'RRmand', 'RR_set']:\n",
    "    # Calculate mean of these parameters in each recording and use them in group comparison\n",
    "    group_1 = data_sippv_all.groupby(level=0)[par].mean()\n",
    "    group_2 = data_simv_all.groupby(level=0)[par].mean()\n",
    "    # The aggregated values (means, medians) themselves do not show normal distribution in the groups, \n",
    "    # use nonparametric measures and tests\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1.dropna(), group_2.dropna(),), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mann-Whitney test\n",
    "\n",
    "for par in ['FiO2_set', 'Leak',]:\n",
    "    # Calculate median of these parameters in each recording and use them in group comparison\n",
    "    group_1 = data_sippv_all.groupby(level=0)[par].median()\n",
    "    group_2 = data_simv_all.groupby(level=0)[par].median()\n",
    "    # The aggregated values (means, medians) themselves do not show normal distribution in the groups, \n",
    "    # use nonparametric measures and tests\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1.dropna(), group_2.dropna(),), '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 14. Compare SIMV and SIPPV in those recordings when both occurred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_both_sippv_simv_all = pd.concat(data_both_sippv_simv)\n",
    "data_both_sippv_simv_all.index.set_names(['recording', 'datetime'], inplace=True)\n",
    "data_both_sippv_simv_all = data_both_sippv_simv_all.reset_index().set_index(['recording', 'Ventilator_mode'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_both_sippv_simv_all.groupby(level=[0,1], observed=False).size()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- AL000089: SIMV-VG -> SIPPV-VG\n",
    "- AL000104: SIPPV-VG -> SIMV-VG\n",
    "- AL000246: SIMV-VG -> SIPPV-VG\n",
    "- AL000465: SIPPV-VG -> SIMV-VG\n",
    "- AL000687: SIMV-VG -> SIPPV-VG\n",
    "- AT000009: SIPPV-VG -> SIMV-VG\n",
    "- AT000578: SIPPV-VG -> SIMV-VG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats_sippv_simv_mean = data_both_sippv_simv_all.groupby(level=[0,1],\n",
    "            observed=False)[['PIP', 'PEEP', 'MAP', 'MV_kg', 'MVresp', 'RRmand', 'VTemand_kg', 'VTemand_resp_kg', 'VTespon_pat_kg']].mean().unstack()\n",
    "stats_sippv_simv_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats_sippv_simv_mean.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.ttest_rel(stats_sippv_simv_mean['PIP']['SIPPV'], stats_sippv_simv_mean['PIP']['SIMV'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.ttest_rel(stats_sippv_simv_mean['PEEP']['SIPPV'], stats_sippv_simv_mean['PEEP']['SIMV'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.ttest_rel(stats_sippv_simv_mean['MAP']['SIPPV'], stats_sippv_simv_mean['MAP']['SIMV'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.ttest_rel(stats_sippv_simv_mean['MV_kg']['SIPPV'], stats_sippv_simv['MV_kg']['SIMV'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.ttest_rel(stats_sippv_simv_mean['RRmand']['SIPPV'], stats_sippv_simv['RRmand']['SIMV'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.ttest_rel(stats_sippv_simv_mean['VTemand_kg']['SIPPV'], stats_sippv_simv['VTemand_kg']['SIMV'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stats_sippv_simv_median = data_both_sippv_simv_all.groupby(level=[0,1],\n",
    "            observed=False)[['FiO2_set', 'Leak']].median().unstack()\n",
    "stats_sippv_simv_median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats_sippv_simv_median.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.ttest_rel(stats_sippv_simv_median['FiO2_set']['SIPPV'], stats_sippv_simv_median['FiO2_set']['SIMV'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.ttest_rel(stats_sippv_simv_median['Leak']['SIPPV'], stats_sippv_simv_median['Leak']['SIMV'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 15. Compare VG and no VG recordings\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clin_data_comput_ELBW.loc[sorted(predominantly_vg)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clin_data_comput_ELBW.loc[sorted(predominantly_novg)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for par in ['Gestational Age (weeks)_comput', 'Birth Weight_comput']:\n",
    "    group_1 = clin_data_comput_ELBW.loc[sorted(predominantly_vg)][par]\n",
    "    group_2 = clin_data_comput_ELBW.loc[sorted(predominantly_novg)][par]\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1, group_2,), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ventilator parameters with normal distribution\n",
    "pars_for_mean = ['PIP', 'MAP', 'PEEP', 'MV_kg', 'VTe_kg', 'VG_set_kg',  'VTimand_kg', 'VTemand_kg', 'RRmand', 'RR_set']\n",
    "# Ventilator parameters which are not normally distributed\n",
    "pars_for_median = ['FiO2_set', 'Leak', 'VT_diff_kg', 'VT_diff_kg_abs', 'P_diff', 'P_diff_abs']\n",
    "\n",
    "vg_stats_mean = round(data_vg_all.groupby(level=0)[pars_for_mean].describe().swaplevel(axis=1)['mean'].describe().T, 2)\n",
    "vg_stats_median = round(data_vg_all.groupby(level=0)[pars_for_median].describe().swaplevel(axis=1)['50%'].describe().T, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vg_stats_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vg_stats_median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ventilator parameters with normal distribution\n",
    "pars_for_mean = ['PIP', 'MAP', 'PEEP', 'MV_kg', 'VTe_kg', 'VTimand_kg', 'VTemand_kg', 'RRmand', 'RR_set']\n",
    "# Ventilator parameters which are not normally distributed\n",
    "pars_for_median = ['FiO2_set', 'Leak', ]\n",
    "\n",
    "novg_stats_mean = round(data_novg_all.groupby(level=0)[pars_for_mean].describe().swaplevel(axis=1)['mean'].describe().T, 2)\n",
    "novg_stats_median = round(data_novg_all.groupby(level=0)[pars_for_median].describe().swaplevel(axis=1)['50%'].describe().T, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "novg_stats_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "novg_stats_median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars = ['PIP', 'MAP', 'PEEP', 'MV_kg', 'VTe_kg', 'RRmand', 'RR_set', 'FiO2_set', 'Leak', ]\n",
    "# The aggregated values (means, medians) themselves do not show normal distribution in the groups, \n",
    "# use nonparametric measures and tests\n",
    "stts = ['count', 'min', '25%', '50%', '75%', 'max']\n",
    "\n",
    "vg_stats = pd.concat([vg_stats_mean[stts], vg_stats_median[stts]])\n",
    "novg_stats = pd.concat([novg_stats_mean[stts], novg_stats_median[stts]])\n",
    "vg_novg_stats = pd.merge(vg_stats,  novg_stats, left_index=True, right_index=True, \n",
    "                            suffixes= ['_vg', '_novg']).sort_index(axis=1, ascending=False)\n",
    "vg_novg_stats = vg_novg_stats.loc[pars]\n",
    "vg_novg_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Inferential statistics for selected parameters, VG vs noVG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mann-Whitney test\n",
    "\n",
    "for par in ['PIP', 'MAP', 'PEEP', 'MV_kg', 'VTe_kg', 'RRmand', 'RR_set']:\n",
    "    # Calculate mean of these parameters in each recording and use them in group comparison\n",
    "    group_1 = data_vg_all.groupby(level=0)[par].mean()\n",
    "    group_2 = data_novg_all.groupby(level=0)[par].mean()\n",
    "    # The aggregated values (means, medians) themselves do not show normal distribution in the groups, \n",
    "    # use nonparametric measures and tests\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1.dropna(), group_2.dropna(),), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mann-Whitney test\n",
    "\n",
    "for par in ['FiO2_set', 'Leak',]:\n",
    "    # Calculate median of these parameters in each recording and use them in group comparison\n",
    "    group_1 = data_vg_all.groupby(level=0)[par].median()\n",
    "    group_2 = data_novg_all.groupby(level=0)[par].median()\n",
    "    # The aggregated values (means, medians) themselves do not show normal distribution in the groups, \n",
    "    # use nonparametric measures and tests\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1.dropna(), group_2.dropna(),), '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 16. Compare VG and noVG in those recordings when both occurred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- VTemand, PIP and MV and FiO2 in VG and noVG recordings\n",
    "- VTdiff and Pdiff in VG recordings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_both_vg_novg_all.index.set_names(['recording', 'datetime'], inplace=True)\n",
    "data_both_vg_novg_all = data_both_vg_novg_all.reset_index().set_index(['recording', 'VG_state'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_both_vg_novg_all.groupby(level=[0,1], observed=False).size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_both_vg_novg_all.groupby(level=[0,1], observed=False)[['MV_kg','RRmand', 'VTe_kg']].mean().unstack()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_both_vg_novg_all.groupby(level=[0,1], observed=False)[['MV_kg','RRmand', 'VTe_kg']].mean().unstack().mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- AL000110: \n",
    "- AL000115: \n",
    "- AL001019: \n",
    "- AT000437: \n",
    "- AT000613: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats_vg_novg_mean = data_both_vg_novg_all.groupby(level=[0,1],\n",
    "            observed=False)[['MAP', 'MV_kg', 'PEEP', 'PIP', 'RR',\n",
    "       'VTemand_kg', 'VTimand_kg', 'VG_set_kg', 'VTe_kg', 'RRmand',]].mean().unstack()\n",
    "stats_vg_novg_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats_vg_novg_mean.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.wilcoxon(stats_vg_novg_mean['MAP']['off'], stats_vg_novg_mean['MAP']['on'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.wilcoxon(stats_vg_novg_mean['PIP']['off'], stats_vg_novg_mean['PIP']['on'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.wilcoxon(stats_vg_novg_mean['PEEP']['off'], stats_vg_novg_mean['PEEP']['on'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.wilcoxon(stats_vg_novg_mean['MV_kg']['off'], stats_vg_novg_mean['MV_kg']['on'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.wilcoxon(stats_vg_novg_mean['RR']['off'], stats_vg_novg_mean['RR']['on'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.wilcoxon(stats_vg_novg_mean['VTe_kg']['off'], stats_vg_novg_mean['VTe_kg']['on'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.wilcoxon(stats_vg_novg_mean['VTemand_kg']['off'], stats_vg_novg_mean['VTemand_kg']['on'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.wilcoxon(stats_vg_novg_mean['VTimand_kg']['off'], stats_vg_novg_mean['VTimand_kg']['on'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats_vg_novg_median  = data_both_vg_novg_all.groupby(level=[0,1],\n",
    "            observed=False)[['Leak', 'FiO2_set', ]].median().unstack()\n",
    "stats_vg_novg_median"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "stats.wilcoxon(stats_vg_novg_median['FiO2_set']['off'], stats_vg_novg_median['FiO2_set']['on'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.wilcoxon(stats_vg_novg_median['Leak']['off'], stats_vg_novg_median['Leak']['on'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 17. Analyse blood gases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-  Blood gases at the arrival to NICU for all cases and separately for SIPPV, SIMV and VG, non-VG\n",
    "-  pCO2s at the end of transfer vs VTemand and MV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. Blood gases at the end of transfer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Blood gases were available in 47/55 and 50/50 cases, during arrival and at handover, respectively.\n",
    "None of them was arterial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All cases, arrival of Cerny\n",
    "blood_gases_arrival = blood_gases_manual['At arrival of PCAM'][['pH', 'pCO2', 'BE']]\n",
    "blood_gases_arrival.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All cases, handover on NICU\n",
    "blood_gases_handover = blood_gases_manual['At arrival to NICU'][['pH', 'pCO2', 'BE']]\n",
    "blood_gases_handover.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIPPV, arrival of Cerny\n",
    "recs_available_sippv = set(blood_gases_arrival.index) & predominantly_sippv\n",
    "blood_gases_arrival.loc[sorted(recs_available_sippv)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIPPV, handover on NICU\n",
    "recs_available_sippv = set(blood_gases_handover.index) & predominantly_sippv\n",
    "blood_gases_handover.loc[sorted(recs_available_sippv)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIMV, arrival of Cerny\n",
    "recs_available_simv = set(blood_gases_arrival.index) & predominantly_simv\n",
    "blood_gases_arrival.loc[sorted(recs_available_simv)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIMV, handover on NICU\n",
    "recs_available_simv = set(blood_gases_handover.index) & predominantly_simv\n",
    "blood_gases_handover.loc[sorted(recs_available_simv)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mann-Whitney test, blood gases at arrival of Cerny\n",
    "for par in ['pH', 'pCO2', 'BE']:\n",
    "    group_1 = blood_gases_arrival.loc[sorted(recs_available_sippv)][par]\n",
    "    group_2 = blood_gases_arrival.loc[sorted(recs_available_simv)][par]\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1.dropna(), group_2.dropna(),), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mann-Whitney test, blood gases on handover to NICU\n",
    "for par in ['pH', 'pCO2', 'BE']:\n",
    "    group_1 = blood_gases_handover.loc[sorted(recs_available_sippv)][par]\n",
    "    group_2 = blood_gases_handover.loc[sorted(recs_available_simv)][par]\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1.dropna(), group_2.dropna(),), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# VG, arrival of Cerny\n",
    "recs_available_vg = set(blood_gases_arrival.index) & predominantly_vg\n",
    "blood_gases_arrival.loc[sorted(recs_available_vg)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# VG, handover on NICU\n",
    "recs_available_vg = set(blood_gases_handover.index) & predominantly_vg\n",
    "blood_gases_handover.loc[sorted(recs_available_vg)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# No VG, arrival of Cerny\n",
    "recs_available_novg = set(blood_gases_arrival.index) & predominantly_novg\n",
    "blood_gases_arrival.loc[sorted(recs_available_novg)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# No VG, handover on NICU\n",
    "recs_available_novg = set(blood_gases_handover.index) & predominantly_novg\n",
    "blood_gases_handover.loc[sorted(recs_available_novg)].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mann-Whitney test, blood gases at arrival of Cerny\n",
    "for par in ['pH', 'pCO2', 'BE']:\n",
    "    group_1 = blood_gases_arrival.loc[sorted(recs_available_vg)][par]\n",
    "    group_2 = blood_gases_arrival.loc[sorted(recs_available_novg)][par]\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1.dropna(), group_2.dropna(),), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mann-Whitney test, blood gases at handover on NICU\n",
    "for par in ['pH', 'pCO2', 'BE']:\n",
    "    group_1 = blood_gases_handover.loc[sorted(recs_available_vg)][par]\n",
    "    group_2 = blood_gases_handover.loc[sorted(recs_available_novg)][par]\n",
    "    print(par, '\\n', stats.mannwhitneyu(group_1.dropna(), group_2.dropna(),), '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. Relationship between pCO2 at the end of transfer and VTe or MV\n",
    "\n",
    "For this, use computationally collected gases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "blood_gases_frme = pd.concat(blood_gases_comput)\n",
    "blood_gases_frme.index.names = ['case', 'datetime']\n",
    "blood_gases_frme.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "blood_gases_last_frme = blood_gases_frme.drop_duplicates(subset='case', keep='last')\n",
    "blood_gases_last_frme.set_index('case', inplace=True)\n",
    "blood_gases_last_frme.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIPPV_MV_CO2 = pd.merge(left = data_sippv_all.groupby(level=0)[['MV_kg']].mean(),\n",
    "    right = blood_gases_last_frme['pCO2'].astype('float'), how='inner', left_index=True, right_index=True).dropna(how='any')\n",
    "SIPPV_MV_CO2.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIMV_MV_CO2 = pd.merge(left = data_simv_all.groupby(level=0)[['MV_kg']].mean(),\n",
    "    right = blood_gases_last_frme['pCO2'].astype('float'), how='inner', left_index=True, right_index=True).dropna(how='any')\n",
    "\n",
    "SIMV_MV_CO2.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIPPV_SIMV_MV_CO2 = pd.concat([SIPPV_MV_CO2, SIMV_MV_CO2], keys = ['SIPPV', 'SIMV'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi=300; filetype = 'jpg'\n",
    "\n",
    "fig,ax = plt.subplots(1,1, figsize=(6,6))\n",
    "ax.scatter(SIPPV_MV_CO2['MV_kg'], SIPPV_MV_CO2['pCO2'], marker='o', facecolors='none', color ='black')\n",
    "ax.scatter(SIMV_MV_CO2['MV_kg'], SIMV_MV_CO2['pCO2'], marker = 'x', color = 'black')\n",
    "ax.set_ylim(0,100)\n",
    "ax.set_xlabel('MV [L/min/kg]')\n",
    "ax.set_ylabel('pCO$_2$ [mmHg]')\n",
    "fig.savefig(os.path.join(DIR_WRITE, f'MV_CO2.{filetype}'), dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VG_VT_CO2 = pd.merge(left = data_vg_all.groupby(level=0)[['VTemand_kg']].mean(),\n",
    "    right = blood_gases_last_frme['pCO2'].astype('float'), how='inner', left_index=True, right_index=True).dropna(how='any')\n",
    "VG_VT_CO2.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noVG_VT_CO2 = pd.merge(left = data_novg_all.groupby(level=0)[['VTemand_kg']].mean(),\n",
    "    right = blood_gases_last_frme['pCO2'].astype('float'), how='inner', left_index=True, right_index=True).dropna(how='any')\n",
    "noVG_VT_CO2.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VG_noVG_VTCO2 = pd.concat([VG_VT_CO2, noVG_VT_CO2], keys = ['VG', 'noVG'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpi=300; filetype = 'jpg'\n",
    "\n",
    "fig,ax = plt.subplots(1,1, figsize=(6,6))\n",
    "ax.scatter(VG_VT_CO2['VTemand_kg'], VG_VT_CO2['pCO2'], marker='o', facecolors='none', color ='black')\n",
    "ax.scatter(noVG_VT_CO2['VTemand_kg'], noVG_VT_CO2['pCO2'], marker = 'x', color = 'black')\n",
    "ax.hlines(37.5, 0, 10, color='black', linestyle='dashed')\n",
    "ax.hlines(60, 0, 10, color='black', linestyle='dashed')\n",
    "ax.vlines(4, 0, 100, color='black', linestyle='dashed')\n",
    "ax.vlines(6, 0, 100, color='black', linestyle='dashed')\n",
    "ax.set_xlim(0, 10)\n",
    "ax.set_ylim(0,100)\n",
    "ax.set_xlabel('VTe [mL/kg]')\n",
    "ax.set_ylabel('pCO$_2$ [mmHg]')\n",
    "fig.savefig(os.path.join(DIR_WRITE, f'VT_CO2.{filetype}'), dpi = dpi, format = filetype, bbox_inches='tight', pad_inches=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}