{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":82611,"databundleVersionId":9553358,"sourceType":"competition"},{"sourceId":199742556,"sourceType":"kernelVersion"}],"dockerImageVersionId":30786,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# BrisT1D - Autocorrelations & Granger Causality\n\nHere's a check-in on some tools for analysing participant's full timeseries data.\nThe melted full-study data representation can be seen in this [**notebook**](https://www.kaggle.com/code/jagofc/brist1d-full-interactive-timeseries).","metadata":{}},{"cell_type":"code","source":"import datetime\nimport numpy as np\nimport pandas as pd\nimport plotly.express as px\n\nfrom pathlib import Path\nfrom typing import Optional\nfrom plotly.subplots import make_subplots\n\nplotly_colors = px.colors.qualitative.Plotly\n\n\ninput_path = Path('/kaggle/input')\nprocessed_data_root = input_path / 'brist1d-to-timeseries/processed' / '241006_130843'\n\nstudy_start = datetime.datetime(2000, 1, 1)\nstudy_end = datetime.datetime(2000, 4, 1)\nstudy_range = pd.date_range(study_start, study_end, freq='5min')\n\nfeatures = ['bg', 'activity', 'insulin', 'carbs', 'cals', 'steps', 'hr']\nparticipants = ['p01', 'p02', 'p03', 'p04', 'p05', 'p06', 'p10', 'p11', 'p12']\n\n\ndef compile_features(p_num: str) -> pd.DataFrame:\n    df_core = pd.DataFrame({'timestamp': study_range})\n    \n    for f in features: \n        df_feat = pd.read_parquet(processed_data_root / f'p_num={p_num}/feature={f}')\n        df_core = df_core.merge(df_feat, how='left', on='timestamp')\n    \n    return df_core\n\n# compile 5-minute granularity timeseries into a dictionary of participants\ndfs_5m = {}\nfor p in participants:\n    # get features\n    df = compile_features(p)  \n    \n    # get \"is monitored\" mask\n    mask = df.loc[:, features].notna().any(axis=1)\n    df['has_data'] = mask\n    \n    # impute core features when \"is monitored\"\n    df.loc[mask, 'is_active'] = df.loc[mask].activity.notna().astype(int)\n    df.loc[mask] = (\n        df.loc[mask]\n        .fillna({\n            'carbs': 0,\n            'cals': np.nanmin(df['cals']),\n            'steps': 0,\n            'is_active': 0\n        })\n    )\n    \n    # add dataframe to collection\n    dfs_5m[p] = df\n    \n# compile 15-minute granularity timeseries into a dictionary of participants\ndfs_15m = {}\nfor p in participants:\n    # downsample the 5-minute granularity dataframes to reduce sampling effects\n    df_15m = dfs_5m[p].resample('15min', on='timestamp', closed='right', label='right').mean(numeric_only=True)\n    \n    # add a target feature :)\n    df_15m['target'] = df_15m.bg.shift(-1, 'h')\n    df_15m = df_15m.reset_index()\n    \n    # add dataframe to collection\n    dfs_15m[p] = df_15m","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-10-09T11:53:02.218191Z","iopub.execute_input":"2024-10-09T11:53:02.218664Z","iopub.status.idle":"2024-10-09T11:53:06.449812Z","shell.execute_reply.started":"2024-10-09T11:53:02.218619Z","shell.execute_reply":"2024-10-09T11:53:06.448190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Marginal feature correlations\n\nTimeseries marginal correlation using 15-minute downsampling for each feature.\n\nTarget = \"bg+1:00\"","metadata":{}},{"cell_type":"code","source":"# Let's run with participant p01\nselected_p_num = 'p01'","metadata":{"execution":{"iopub.status.busy":"2024-10-09T11:53:06.452363Z","iopub.execute_input":"2024-10-09T11:53:06.453197Z","iopub.status.idle":"2024-10-09T11:53:06.459195Z","shell.execute_reply.started":"2024-10-09T11:53:06.453135Z","shell.execute_reply":"2024-10-09T11:53:06.457655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.imshow(\n    dfs_15m[selected_p_num].corr(numeric_only=True), \n    text_auto=True,\n    height=800,\n    color_continuous_scale='RdBu_r'\n)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-09T11:54:31.874657Z","iopub.execute_input":"2024-10-09T11:54:31.875202Z","iopub.status.idle":"2024-10-09T11:54:31.948265Z","shell.execute_reply.started":"2024-10-09T11:54:31.875153Z","shell.execute_reply":"2024-10-09T11:54:31.946676Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Autocorrelation","metadata":{}},{"cell_type":"markdown","source":"Heart-rate shows strong positive autocorrelation at a lag of ~ 288 * 5 mins = 24 hours.\n\nI.e. there's a daily pattern, as expected.","metadata":{}},{"cell_type":"code","source":"from statsmodels.tsa.stattools import acf","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2024-10-08T17:16:00.658705Z","iopub.execute_input":"2024-10-08T17:16:00.659181Z","iopub.status.idle":"2024-10-08T17:16:00.664695Z","shell.execute_reply.started":"2024-10-08T17:16:00.659135Z","shell.execute_reply":"2024-10-08T17:16:00.663244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_acf(df, feature, nlags=12*48, minutes_per_lag=5, **layout_args):\n    ac, ci = acf(\n        df.loc[:, [feature]].dropna(),\n        alpha=0.05,\n        nlags=nlags\n    )\n    \n    df_acf = pd.DataFrame({\n        'lag_minutes': minutes_per_lag * np.arange(0, len(ac)),\n        'acf': ac,\n        'ci_lower': ci[:, 0],\n        'ci_upper': ci[:, 1]\n    })\n    \n    fig = px.line(\n        df_acf,\n        x='lag_minutes',\n        y=['acf', 'ci_lower', 'ci_upper'],\n        title=f'Autocorrelation function for {feature.upper()}'\n    )\n    \n    fig.update_layout(\n        yaxis_range=[-1.1, 1.1],\n        **layout_args\n    )\n\n    \n    return fig","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-10-08T17:16:02.150270Z","iopub.execute_input":"2024-10-08T17:16:02.150720Z","iopub.status.idle":"2024-10-08T17:16:02.158727Z","shell.execute_reply.started":"2024-10-08T17:16:02.150679Z","shell.execute_reply":"2024-10-08T17:16:02.157309Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_acf(\n    dfs_5m[selected_p_num], \n    'hr', \n    nlags=12*48, \n    title=f'Autocorrelation function for HR (5-minutely) for participant {selected_p_num}',\n    width=800,\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T17:16:18.119934Z","iopub.execute_input":"2024-10-08T17:16:18.120453Z","iopub.status.idle":"2024-10-08T17:16:18.210110Z","shell.execute_reply.started":"2024-10-08T17:16:18.120404Z","shell.execute_reply":"2024-10-08T17:16:18.208866Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_acf(\n    dfs_15m[selected_p_num], \n    'hr', \n    nlags=4*48, \n    minutes_per_lag=15, \n    title=f'Autocorrelation function for HR (15-minutely) for participant {selected_p_num}',\n    width=800\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T17:17:04.165958Z","iopub.execute_input":"2024-10-08T17:17:04.166971Z","iopub.status.idle":"2024-10-08T17:17:04.248670Z","shell.execute_reply.started":"2024-10-08T17:17:04.166924Z","shell.execute_reply":"2024-10-08T17:17:04.247448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_acf(\n    dfs_15m[selected_p_num], \n    'target', \n    nlags=48*4,\n    title=f'Autocorrelation function for TARGET (15-minutely) for participant {selected_p_num}',\n    width=800\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T17:17:36.303702Z","iopub.execute_input":"2024-10-08T17:17:36.304238Z","iopub.status.idle":"2024-10-08T17:17:36.388264Z","shell.execute_reply.started":"2024-10-08T17:17:36.304189Z","shell.execute_reply":"2024-10-08T17:17:36.387154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Crosscorrelation","metadata":{}},{"cell_type":"code","source":"from statsmodels.tsa.stattools import ccf","metadata":{"execution":{"iopub.status.busy":"2024-10-08T15:40:05.741951Z","iopub.execute_input":"2024-10-08T15:40:05.742436Z","iopub.status.idle":"2024-10-08T15:40:05.747505Z","shell.execute_reply.started":"2024-10-08T15:40:05.742387Z","shell.execute_reply":"2024-10-08T15:40:05.746438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_ccf(df, feature_1, feature_2, diff=False, nlags=12*48, minutes_per_lag=5, **layout_args):\n    if not diff:\n        ccf_df = df.loc[:, [feature_1, feature_2]].dropna()\n    else:\n        ccf_df = df.loc[:, [feature_1, feature_2]].diff().dropna()\n    \n    ac, ci = ccf(\n        ccf_df.loc[:, feature_1],\n        ccf_df.loc[:, feature_2],\n        alpha=0.05,\n        nlags=nlags\n    )\n    \n    df_ccf = pd.DataFrame({\n        'lag_minutes': minutes_per_lag * np.arange(0, len(ac)),\n        'acf': ac,\n        'ci_lower': ci[:, 0],\n        'ci_upper': ci[:, 1]\n    })\n    \n    title_diff = '' if not diff else ' using 1st differences.'\n    fig = px.line(\n        df_ccf,\n        x='lag_minutes',\n        y=['acf', 'ci_lower', 'ci_upper'],\n        title=f'Crosscorrelation function for {feature_1.upper()} and {feature_2.upper()}' + title_diff\n    )\n    \n    fig.update_layout(\n        yaxis_range=[-1.1, 1.1],\n        **layout_args\n    )\n    \n    return fig","metadata":{"execution":{"iopub.status.busy":"2024-10-08T15:42:59.732586Z","iopub.execute_input":"2024-10-08T15:42:59.733910Z","iopub.status.idle":"2024-10-08T15:42:59.742772Z","shell.execute_reply.started":"2024-10-08T15:42:59.733854Z","shell.execute_reply":"2024-10-08T15:42:59.741555Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_ccf(\n    dfs_15m['p01'], \n    'hr', \n    'bg', \n    title=f'Crosscorrelation function for HR and BG (15-minutely) for participant {selected_p_num}',\n    width=800\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T17:18:12.003191Z","iopub.execute_input":"2024-10-08T17:18:12.003662Z","iopub.status.idle":"2024-10-08T17:18:12.087440Z","shell.execute_reply.started":"2024-10-08T17:18:12.003616Z","shell.execute_reply":"2024-10-08T17:18:12.086352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_ccf(\n    dfs_15m['p01'], \n    'bg',\n    'target',\n    title=f'Crosscorrelation function for BG and TARGET (15-minutely) for participant {selected_p_num}',\n    width=800,\n    minutes_per_lag=15,\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T17:19:05.370973Z","iopub.execute_input":"2024-10-08T17:19:05.372094Z","iopub.status.idle":"2024-10-08T17:19:05.454778Z","shell.execute_reply.started":"2024-10-08T17:19:05.372041Z","shell.execute_reply":"2024-10-08T17:19:05.453546Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"VERY strong cross-correlation seen above, between TARGET and BG at a lag of 60 minutes (...just a test).\n\nAlso, below, VERY strong cross-correlation between cals and steps: probably because cals are comin directly from a smartwatch ?:)","metadata":{}},{"cell_type":"code","source":"plot_ccf(\n    dfs_15m['p01'], \n    'cals', \n    'steps', \n    diff=True, \n    nlags=12*4, \n    minutes_per_lag=15, \n    title=f'Crosscorrelation function for CALS and STEPS (15-minutely) for participant {selected_p_num} using 1st differences',\n    width=1000\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T17:21:57.051065Z","iopub.execute_input":"2024-10-08T17:21:57.052078Z","iopub.status.idle":"2024-10-08T17:21:57.143137Z","shell.execute_reply.started":"2024-10-08T17:21:57.052030Z","shell.execute_reply":"2024-10-08T17:21:57.141899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Granger Causality Tests\n\nLet's test out some lag=1 Granger causality!","metadata":{}},{"cell_type":"code","source":"from statsmodels.tsa.stattools import grangercausalitytests","metadata":{"execution":{"iopub.status.busy":"2024-10-08T16:06:06.550103Z","iopub.execute_input":"2024-10-08T16:06:06.551041Z","iopub.status.idle":"2024-10-08T16:06:06.555701Z","shell.execute_reply.started":"2024-10-08T16:06:06.550976Z","shell.execute_reply":"2024-10-08T16:06:06.554616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**carbs -> hr?**\n\nLooks like there's not nothing.","metadata":{}},{"cell_type":"code","source":"_ = grangercausalitytests(\n    dfs_15m['p01'].loc[:, ['hr', 'carbs']].diff().dropna(),\n    maxlag=1\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T16:47:12.943070Z","iopub.execute_input":"2024-10-08T16:47:12.943933Z","iopub.status.idle":"2024-10-08T16:47:12.967828Z","shell.execute_reply.started":"2024-10-08T16:47:12.943886Z","shell.execute_reply":"2024-10-08T16:47:12.966393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**hr -> carbs?**\n\nLooks like there may be nothing in this direction.","metadata":{}},{"cell_type":"code","source":"_ = grangercausalitytests(\n    dfs_15m['p01'].loc[:, ['carbs', 'hr']].diff().dropna(),\n    maxlag=1\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T16:47:14.265678Z","iopub.execute_input":"2024-10-08T16:47:14.266843Z","iopub.status.idle":"2024-10-08T16:47:14.290308Z","shell.execute_reply.started":"2024-10-08T16:47:14.266791Z","shell.execute_reply":"2024-10-08T16:47:14.288802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**hr -> is_active?**\n\nNope, nothing showing here.","metadata":{}},{"cell_type":"code","source":"_ = grangercausalitytests(\n    dfs_15m['p01'].loc[:, ['is_active', 'hr']].diff().dropna(),\n    maxlag=1\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T16:47:15.011038Z","iopub.execute_input":"2024-10-08T16:47:15.011941Z","iopub.status.idle":"2024-10-08T16:47:15.035144Z","shell.execute_reply.started":"2024-10-08T16:47:15.011892Z","shell.execute_reply":"2024-10-08T16:47:15.033908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**is_active -> hr?**\n\nLooks like a pretty strong predictive causal link!","metadata":{}},{"cell_type":"code","source":"_ = grangercausalitytests(\n    dfs_15m['p01'].loc[:, ['hr', 'is_active']].diff().dropna(),\n    maxlag=1\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T16:47:15.545000Z","iopub.execute_input":"2024-10-08T16:47:15.545493Z","iopub.status.idle":"2024-10-08T16:47:15.567926Z","shell.execute_reply.started":"2024-10-08T16:47:15.545448Z","shell.execute_reply":"2024-10-08T16:47:15.566551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**bg -> insulin?**\n\nCould be something... (Should be something...)","metadata":{}},{"cell_type":"code","source":"_ = grangercausalitytests(\n    dfs_15m['p01'].loc[:, ['insulin', 'bg']].diff().dropna(),\n    maxlag=1\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T16:47:16.020171Z","iopub.execute_input":"2024-10-08T16:47:16.020589Z","iopub.status.idle":"2024-10-08T16:47:16.044466Z","shell.execute_reply.started":"2024-10-08T16:47:16.020551Z","shell.execute_reply":"2024-10-08T16:47:16.043091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**insulin -> bg?**\n\nCould be something... (Should be something...)","metadata":{}},{"cell_type":"code","source":"_ = grangercausalitytests(\n    dfs_15m['p01'].loc[:, ['bg', 'insulin']].diff().dropna(),\n    maxlag=1\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-08T16:47:19.298229Z","iopub.execute_input":"2024-10-08T16:47:19.298716Z","iopub.status.idle":"2024-10-08T16:47:19.325248Z","shell.execute_reply.started":"2024-10-08T16:47:19.298671Z","shell.execute_reply":"2024-10-08T16:47:19.323686Z"},"trusted":true},"execution_count":null,"outputs":[]}]}