Polar Verity Sense: Guide to data extraction and analysis

^ A picture of the polar verity sense heart rate monitor ^

The Polar Verity Sense is a $89 armband heart rate sensor that has often been reviewed as a viable alternative to ECG chest-strap style monitors. Portable and comfortable, it can be worn for exercises indoors, outdoors, and even in the pool.

Beyond its accuracy as a heart rate monitor, the Verity Sense’s utility lies in its interface with Polar’s official app, which helps keep track of data collected from the sensor as well as providing plenty of useful visualizations and analytics.

In the following sections, we will present what we have learned from using the Verity Sense and how to extract its data for visualizations and statistical analyses. With the wearipedia library, this process will only require your username and password!

We will be able to extract the following parameters (in bold are parameters we will be using in this notebook). Per Activity means one measurement through a single session of exercise. Per second means one measurement every second throughout an exercise.

Parameter Name	Sampling Frequency
Calories	Per Activity
Duration	Per Activity
Sessions	Per day
Continuous Heart Rate	Per second

In this guide, we sequentially cover the following nine topics:

1. Set up

2. Authentication/Authorization

   • Requires only username and password, no OAuth.

3. Data extraction

• We get data via wearipedia in a couple lines of code.

4. Data Exporting

• We export all of this data to file formats compatible by R, Excel, and MatLab.

5. Adherence

   • We simulate non-adherence by dynamically removing datapoints from our simulated data.

6. Visualization

   • We create a simple plot to visualize our data.

7. Data visualization & analysis

• 7.1 Heart Rate Zones during an Exercise Session

• 7.2 Visualizing Heart Rate Averages and Sessions over a Month

• 7.3 Visualizing Continuous Heart Rate data

8. Outlier Detection and Data Cleaning

   • We detect outliers in our data and filter them out.

9. Data Analysis

• 9.1 Analyzing correlation between Average Heart Rate and Calories Burned per Minute

• 9.2 Analyzing correlation between Duration of Workouts and Calories Burned

• 9.3 Analyzing how activity intensity compares in first half and second half of a session

Note that we are not making any scientific claims here as our sample size is small and the data collection process was not rigorously vetted (it is our own data), only demonstrating that this code could potentially be used to perform rigorous analyses in the future.

Disclaimer: this notebook is purely for educational purposes. All of the data currently stored in this notebook is purely synthetic, meaning randomly generated according to rules we created. Despite this, the end-to-end data extraction pipeline has been tested on our own data, meaning that if you enter your own email and password on your own Colab instance, you can visualize your own real data. That being said, we were unable to thoroughly test the timezone functionality, though, since we only have one account, so beware.

1. Set up

Participant Setup

Dear Participant,

Once you unbox your Polar Verity Sense, please set up the device by following the video:

• Video Guide: https://www.youtube.com/watch?v=wiA_ucJJV7Y&t

Make sure that your phone is paired to it using the Polar Flow login credentials (email and password) given to you by the study coordinator.

Best,

Wearipedia

Data Receiver Setup

Please follow the below steps:

1. Create an email address for the participant, for example foo@email.com

2. Create a Polar Flow account with the email foo@email.com and some random password.

3. Keep foo@email.com and password stored somewhere safe.

4. Distribute the device to the participant and instruct them to follow the participant setup letter above.

5. Install the wearipedia Python package to easily extract data from this device.

!pip install git+https://TrafficCop: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@github.com/TrafficCop/wearipedia.git

Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting git+https://TrafficCop:****@github.com/TrafficCop/wearipedia.git
  Cloning https://TrafficCop:****@github.com/TrafficCop/wearipedia.git to /tmp/pip-req-build-3248x7v1
  Running command git clone --filter=blob:none --quiet 'https://TrafficCop:****@github.com/TrafficCop/wearipedia.git' /tmp/pip-req-build-3248x7v1
  Resolved https://TrafficCop:****@github.com/TrafficCop/wearipedia.git to commit b570a805866b7c6f8ea7b57b85b37e1b29106f97
  Installing build dependencies ... done
  Getting requirements to build wheel ... done
  Preparing metadata (pyproject.toml) ... done
Collecting garminconnect<0.2.0,>=0.1.48
  Downloading garminconnect-0.1.50.tar.gz (17 kB)
  Preparing metadata (setup.py) ... done
Collecting typer[all]<0.7.0,>=0.6.1
  Downloading typer-0.6.1-py3-none-any.whl (38 kB)
Requirement already satisfied: pandas<2.0,>=1.1 in /usr/local/lib/python3.8/dist-packages (from wearipedia==0.1.0) (1.3.5)
Requirement already satisfied: scipy<2.0,>=1.6 in /usr/local/lib/python3.8/dist-packages (from wearipedia==0.1.0) (1.7.3)
Requirement already satisfied: tqdm<5.0.0,>=4.64.1 in /usr/local/lib/python3.8/dist-packages (from wearipedia==0.1.0) (4.64.1)
Collecting rich<13.0.0,>=12.6.0
  Downloading rich-12.6.0-py3-none-any.whl (237 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 237.5/237.5 KB 9.5 MB/s eta 0:00:00
Collecting wget<4.0,>=3.2
  Downloading wget-3.2.zip (10 kB)
  Preparing metadata (setup.py) ... done
Requirement already satisfied: requests in /usr/local/lib/python3.8/dist-packages (from garminconnect<0.2.0,>=0.1.48->wearipedia==0.1.0) (2.25.1)
Collecting cloudscraper
  Downloading cloudscraper-1.2.68-py2.py3-none-any.whl (98 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 98.6/98.6 KB 4.2 MB/s eta 0:00:00
Requirement already satisfied: numpy>=1.17.3 in /usr/local/lib/python3.8/dist-packages (from pandas<2.0,>=1.1->wearipedia==0.1.0) (1.21.6)
Requirement already satisfied: pytz>=2017.3 in /usr/local/lib/python3.8/dist-packages (from pandas<2.0,>=1.1->wearipedia==0.1.0) (2022.7)
Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.8/dist-packages (from pandas<2.0,>=1.1->wearipedia==0.1.0) (2.8.2)
Requirement already satisfied: pygments<3.0.0,>=2.6.0 in /usr/local/lib/python3.8/dist-packages (from rich<13.0.0,>=12.6.0->wearipedia==0.1.0) (2.6.1)
Requirement already satisfied: typing-extensions<5.0,>=4.0.0 in /usr/local/lib/python3.8/dist-packages (from rich<13.0.0,>=12.6.0->wearipedia==0.1.0) (4.4.0)
Collecting commonmark<0.10.0,>=0.9.0
  Downloading commonmark-0.9.1-py2.py3-none-any.whl (51 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 51.1/51.1 KB 3.2 MB/s eta 0:00:00
Requirement already satisfied: click<9.0.0,>=7.1.1 in /usr/local/lib/python3.8/dist-packages (from typer[all]<0.7.0,>=0.6.1->wearipedia==0.1.0) (7.1.2)
Collecting shellingham<2.0.0,>=1.3.0
  Downloading shellingham-1.5.0.post1-py2.py3-none-any.whl (9.4 kB)
Collecting colorama<0.5.0,>=0.4.3
  Downloading colorama-0.4.6-py2.py3-none-any.whl (25 kB)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.8/dist-packages (from python-dateutil>=2.7.3->pandas<2.0,>=1.1->wearipedia==0.1.0) (1.15.0)
Collecting requests-toolbelt>=0.9.1
  Downloading requests_toolbelt-0.10.1-py2.py3-none-any.whl (54 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 54.5/54.5 KB 2.4 MB/s eta 0:00:00
Requirement already satisfied: pyparsing>=2.4.7 in /usr/local/lib/python3.8/dist-packages (from cloudscraper->garminconnect<0.2.0,>=0.1.48->wearipedia==0.1.0) (3.0.9)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /usr/local/lib/python3.8/dist-packages (from requests->garminconnect<0.2.0,>=0.1.48->wearipedia==0.1.0) (1.24.3)
Requirement already satisfied: idna<3,>=2.5 in /usr/local/lib/python3.8/dist-packages (from requests->garminconnect<0.2.0,>=0.1.48->wearipedia==0.1.0) (2.10)
Requirement already satisfied: chardet<5,>=3.0.2 in /usr/local/lib/python3.8/dist-packages (from requests->garminconnect<0.2.0,>=0.1.48->wearipedia==0.1.0) (4.0.0)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.8/dist-packages (from requests->garminconnect<0.2.0,>=0.1.48->wearipedia==0.1.0) (2022.12.7)
Building wheels for collected packages: wearipedia, garminconnect, wget
  Building wheel for wearipedia (pyproject.toml) ... done
  Created wheel for wearipedia: filename=wearipedia-0.1.0-py3-none-any.whl size=42684 sha256=687eb2dbca8a929f258256968d10bbe428c33415380afa1bb62acf64cefc1762
  Stored in directory: /tmp/pip-ephem-wheel-cache-umikwsiy/wheels/24/85/35/40c88fbefd67ae05ca8dd55adf3f437705bfedb5c296ddd98b
  Building wheel for garminconnect (setup.py) ... done
  Created wheel for garminconnect: filename=garminconnect-0.1.50-py3-none-any.whl size=13423 sha256=b9e9de6bacdbbedc26882dd511c22bc0070d8e6342dc83be338df9219726d3eb
  Stored in directory: /root/.cache/pip/wheels/2b/72/01/577425dabe75aecd099d4f6e4afb2b4db7e4bfe3db80c4b5a5
  Building wheel for wget (setup.py) ... done
  Created wheel for wget: filename=wget-3.2-py3-none-any.whl size=9674 sha256=96ba01d4ea9ac4e43ca5e954cfc686215449588b2a75f40ec02d37ab1c80a58b
  Stored in directory: /root/.cache/pip/wheels/bd/a8/c3/3cf2c14a1837a4e04bd98631724e81f33f462d86a1d895fae0
Successfully built wearipedia garminconnect wget
Installing collected packages: wget, commonmark, typer, shellingham, rich, colorama, requests-toolbelt, cloudscraper, garminconnect, wearipedia
  Attempting uninstall: typer
    Found existing installation: typer 0.7.0
    Uninstalling typer-0.7.0:
      Successfully uninstalled typer-0.7.0
Successfully installed cloudscraper-1.2.68 colorama-0.4.6 commonmark-0.9.1 garminconnect-0.1.50 requests-toolbelt-0.10.1 rich-12.6.0 shellingham-1.5.0.post1 typer-0.6.1 wearipedia-0.1.0 wget-3.2

For cleaning up old installations only:

!pip uninstall wearipedia

Found existing installation: wearipedia 0.1.0
Uninstalling wearipedia-0.1.0:
  Would remove:
    /usr/local/bin/wearipedia
    /usr/local/lib/python3.8/dist-packages/wearipedia-0.1.0.dist-info/*
    /usr/local/lib/python3.8/dist-packages/wearipedia/*
Proceed (Y/n)? Y
  Successfully uninstalled wearipedia-0.1.0

#2. Authentication/Authorization To obtain access to data, authorization is required. All you’ll need to do here is just put in your email and password for your Polar Verity Sense device. We’ll use this username and password to extract the data in the sections below.

#@title Enter Polar login credentials

email_address = '' #@param {type:"string"}
password = ''#@param {type:"string"}

3. Data extraction

Data can be extracted via wearipedia, our open-source Python package that unifies dozens of complex wearable device APIs into one simple, common interface.

First, we’ll set a date range and then extract all of the data within that date range. You can select whether you would like synthetic data or not with the checkbox.

#@title Enter start and end dates (in the format yyyy-mm-dd)

#set start and end dates - this will give you all the data from 2000-01-01 (January 1st, 2000) to 2100-02-03 (February 3rd, 2100), for example
start_date='2022-03-01' #@param {type:"string"}
end_date='2022-06-11' #@param {type:"string"}
synthetic = True #@param {type:"boolean"}

CHANGES TO SYNTHETIC THAT NEED TO BE DONE 1. Ability to support different range of custom dates

import wearipedia

device = wearipedia.get_device("polar/verity_sense")

if not synthetic:
  device.authenticate({"email": email_address, "password": password})
  params = {"start_date": start_date, "end_date": end_date}
  data = device.get_data("sessions", params = params)
else:
  params = {"start_date": start_date, "end_date": end_date}
  data = device.get_data("sessions", params = params)

100%|██████████| 108/108 [00:00<00:00, 8073.73it/s]
100%|██████████| 108/108 [00:06<00:00, 17.41it/s]

To get a list of the days that actually have data, we can list out the keys for the returned data. Note it might not be feasible to actually print out the entire set of returned data, as it can be quite large.

print(data.keys())

dict_keys(['2022-03-01', '2022-03-02', '2022-04-05', '2022-03-23', '2022-04-20', '2022-05-11', '2022-03-06', '2022-05-18', '2022-04-18', '2022-05-05', '2022-04-09', '2022-04-06', '2022-05-19', '2022-03-30', '2022-06-08', '2022-06-04', '2022-03-21', '2022-04-19', '2022-04-15', '2022-04-14', '2022-03-22', '2022-05-12', '2022-06-06', '2022-03-11', '2022-03-31', '2022-06-02', '2022-05-26', '2022-03-10', '2022-05-03', '2022-04-23', '2022-06-11', '2022-04-03', '2022-06-03', '2022-05-22', '2022-05-31', '2022-03-16', '2022-05-17', '2022-04-10', '2022-06-05', '2022-05-15', '2022-04-27', '2022-05-08', '2022-03-19', '2022-05-20', '2022-06-01', '2022-04-26', '2022-05-16', '2022-05-21', '2022-04-12', '2022-03-04', '2022-03-28', '2022-04-24', '2022-04-01', '2022-04-25', '2022-04-07', '2022-03-12', '2022-06-10', '2022-03-20', '2022-03-26', '2022-03-27', '2022-05-28', '2022-04-30'])

4. Data Exporting

In this section, we export all of this data to formats compatible with popular scientific computing software (R, Excel, Google Sheets, Matlab). Specifically, we will first export to JSON, which can be read by R and Matlab. Then, we will export to CSV, which can be consumed by Excel, Google Sheets, and every other popular programming language.

Exporting to JSON (R, Matlab, etc.)

Exporting to JSON is fairly simple. We export each datatype separately and also export a complete version that includes all simultaneously.

import json

hrs = [data[k]['heart_rates'] for k in data.keys()]
calories = [data[k]['calories'] for k in data.keys()]
hr_avgs = [sum(data[k]['heart_rates']) / (data[k]['minutes'] * 60) for k in data.keys()]
durations = [data[k]['minutes'] for k in data.keys()]

json.dump(list(data.keys()), open("dates.json", "w"))
json.dump(calories, open("calories.json", "w"))
json.dump(hrs, open("hrs.json", "w"))
json.dump(hr_avgs, open("hr_avgs.json", "w"))
json.dump(durations, open("durations.json", "w"))

complete = {
    "dates": list(data.keys()),
    "calories": calories,
    "hrs": hrs,
    "hr_avgs": hr_avgs,
    "durations": durations,
}


json.dump(complete, open("complete.json", "w"))

Exporting to CSV and XLSX (Excel, Google Sheets, R, Matlab, etc.)

Exporting to CSV/XLSX will require us to first process the data into a pandas dataframe format, which can then be converted into our desired file type.

Here we will export three separate files: the first will contain heart rate data (labeled by date and time), and the second will contain daily summary data such as calories, average heart rate, and minutes (labeled by date). Additionally, we will make a third dataframe containing only heart rate data for one session (of your choice) to help us with further visualizations and analysis in the following sections of this notebook.

day = "2022-03-01" #@param {type: "string"}

# set the day of interest automatically if not specified
if day == "":
  day = list(data.keys())[0]

import pandas as pd
from datetime import timedelta
from datetime import datetime
import numpy as np
import copy

# Third dataframe and list of dataframes:
hr_df = None
hr_dfs = []

# 1. First dataframe
all_df = None
for key in data.keys():
  duration = data[key]['minutes']*60
  start_time = timedelta(hours=0, minutes=0, seconds=0)
  times = [str(start_time + timedelta(seconds=i)).zfill(8) for i in range(int(duration))]
  tdf = pd.DataFrame(zip(times,data[key]['heart_rates']),columns=['time','bpm'])
  mapping = lambda x: np.datetime64(key+ " " + x)
  tdf["time"] = tdf["time"].apply(mapping)
  if key == day:
    hr_df = tdf
  if type(all_df) == type(None):
    all_df = tdf
  else:
    all_df = all_df.append(tdf)
  hr_dfs.append(copy.deepcopy(tdf))
#print("DataFrame 1: All heart rates")
#display(df1)

# 2. Second dataframe
df2 = None
for key in data.keys():
  avg_rate = sum(data[key]['heart_rates']) / (data[key]['minutes'] * 60)
  if type(df2) == type(None):
    df2 = pd.DataFrame([[key,data[key]['minutes'],data[key]['calories'], avg_rate]],columns=['day','minutes','calories', 'avg_rate'])
  else:
    tdf = pd.DataFrame([[key,data[key]['minutes'],data[key]['calories'], avg_rate]],columns=['day','minutes','calories', 'avg_rate'])
    df2 = df2.append(tdf)
#print("DataFrame 2: Summary Data")
#display(df2)

# 3. Third dataframe
print("DataFrame 3: Single session heart rate")
display(hr_df)

DataFrame 3: Single session heart rate

	time	bpm
0	2022-03-01 00:00:00	94.110535
1	2022-03-01 00:00:01	91.859208
2	2022-03-01 00:00:02	93.210171
3	2022-03-01 00:00:03	92.384613
4	2022-03-01 00:00:04	94.371856
...	...	...
3295	2022-03-01 00:54:55	137.121507
3296	2022-03-01 00:54:56	137.750669
3297	2022-03-01 00:54:57	138.724507
3298	2022-03-01 00:54:58	139.206589
3299	2022-03-01 00:54:59	138.982279

3300 rows × 2 columns

With this, we have retrieved all the data we need! Now we can start visualizing our data.

5. Adherence

The device simulator already automatically randomly deletes small chunks of the day. In this section, we will simulate non-adherence over longer periods of time from the participant (day-level and week-level).

Then, we will detect this non-adherence and give a Pandas DataFrame that concisely describes when the participant has had their device on and off throughout the entirety of the time period, allowing you to calculate how long they’ve had it on/off etc.

We will first delete a certain % of blocks either at the day level or week level, with user input.

#@title Non-adherence simulation
block_level = "day" #@param ["day","week"]
adherence_percent = 0.74 #@param {type:"slider", min:0, max:1, step:0.01}

import numpy as np
import pandas as pd
from datetime import datetime
import copy

if block_level == "day":
  block_length = 1
elif block_level == "week":
    block_length = 7

dates = np.array(list(data.keys()))

num_blocks = len(dates) // block_length

num_blocks_to_remove = int((1-adherence_percent) * num_blocks)

remove = np.random.choice(dates, replace=False, size=num_blocks_to_remove)

adhered_data = copy.deepcopy(data)
adhered_dates = copy.deepcopy(list(data.keys()))

for key in remove:
  adhered_data.pop(key)
  adhered_dates.remove(key)

And now we have significantly fewer datapoints! This will give us a more realistic situation, where participants may not use their device for days or weeks at a time.

Now let’s detect non-adherence. We will plot the days when the participant uses the watch and when they do not.

import matplotlib.pyplot as plt
from datetime import datetime
import pandas as pd

datelist = pd.date_range(start_date, end_date)
is_wearing = []
print(adhered_dates)
for date in datelist:
  if datetime.strftime(date, "%Y-%m-%d") in adhered_dates:
    is_wearing.append(1)
  else:
    is_wearing.append(0)

plt.figure(figsize=(12, 6))
plt.plot(datelist,is_wearing,drawstyle="steps-mid")

['2022-03-01', '2022-03-02', '2022-03-23', '2022-04-20', '2022-05-11', '2022-03-06', '2022-05-18', '2022-04-18', '2022-05-05', '2022-04-09', '2022-04-06', '2022-05-19', '2022-03-30', '2022-06-04', '2022-04-19', '2022-04-15', '2022-04-14', '2022-03-22', '2022-05-12', '2022-06-06', '2022-03-31', '2022-06-02', '2022-05-26', '2022-03-10', '2022-05-03', '2022-04-23', '2022-06-11', '2022-04-03', '2022-06-03', '2022-03-16', '2022-05-17', '2022-04-10', '2022-05-15', '2022-04-27', '2022-05-20', '2022-05-16', '2022-04-12', '2022-03-04', '2022-03-28', '2022-04-24', '2022-04-01', '2022-04-25', '2022-06-10', '2022-03-20', '2022-05-28', '2022-04-30']

[<matplotlib.lines.Line2D at 0x7f02b66731c0>]

6. Visualization

We’ve extracted lots of data, but what does it look like?

In this section, we will be visualizing our three kinds of data in a simple, customizable plot! This plot is intended to provide a starter example for plotting, whereas later examples emphasize deep control and aesthetics.

#@title Basic Plot
feature = "average heart rate" #@param ["heart rate", "average heart rate", "calories"]
start_date = "2022-10-21" #@param {type:"date"}
time_interval = "one week" #@param ["one day", "one week", "full time"]
smoothness = 0.05 #@param {type:"slider", min:0, max:1, step:0.01}
smooth_plot = True #@param {type:"boolean"}

from scipy.special import y0
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
import pandas as pd
from scipy.ndimage.filters import gaussian_filter1d

if time_interval == "one day":
  datelist = pd.date_range(start_date, periods=2)
elif time_interval == "one week":
  datelist = pd.date_range(start_date, periods=8)
elif time_interval == "full time":
  datelist = pd.date_range(start_date, end_date)
end = datelist[len(datelist)-1]

title_fillin = feature

if feature == "heart rate":
  tempdf = all_df[all_df["time"].map(lambda x: datetime.strptime(datetime.strftime(x,"%Y-%m-%d"),"%Y-%m-%d") < end)]
  y = tempdf["bpm"]
  x = tempdf["time"]
  sigma = 200*smoothness
elif feature == "average heart rate":
  tempdf = df2[df2["day"].map(lambda x: datetime.strptime(x,"%Y-%m-%d") < end)]
  y = tempdf["avg_rate"]
  x = tempdf["day"]
  sigma = 5*smoothness
elif feature == "calories":
  tempdf = df2[df2["day"].map(lambda x: datetime.strptime(x,"%Y-%m-%d") < end)]
  y = tempdf["calories"]
  x = tempdf["day"]
  sigma = 5*smoothness

if smooth_plot:
  y = list(gaussian_filter1d(y, sigma=sigma))
plt.figure(figsize=(16,10))
plt.plot(x, y)
plt.title(f"{title_fillin} from {start_date} to {datetime.strftime(end,'%Y-%m-%d')}",fontsize=20)
plt.xlabel("Time")
plt.ylabel(title_fillin[:-1])

Text(0, 0.5, 'average heart rat')

This plot allows you to quickly scan your data at many different time scales (day, week, and full) and for different kinds of measurements (heart rate, average heart rate, and calories), which enables easy and fast data exploration.

7. Advanced Visualization

7.1 Visualizing Participants Heart Rate Zones during Exercise

One of the graphs the you see when you open an exercise session on polar is a bar chart of heart rate zone durations, or how long your heart rate zone stayed in some range throughout your exercise. It looks something like this:

At the bottom (gray) zone is the amount of time spent with heart rate between 100-120 bpm, the next being 120-140 bpm, then 140-160 bpm, 160-180 bpm, and finally 180-200 bpm.

We will want the following data:

• The zones we want to plot

• Time in each heart rate zone

Polar uses the following zones: 100-120 bpm, 120-140 bpm, 140-160 bpm, 160-180 bpm, 180-200 bpm, which is what we will be using as well.

To get the time spent in each heart rate zone, we can use the pandas groupby function to aggregate the heart rates by zones, and then count the number of entries in each group. Since the continuous heart rate data has an entry for every second, this directly gives us the time spent in each heart rate zone.

import pandas as pd

ranges = [100,120,140,160,180,200]
z1 = hr_df.groupby(pd.cut(hr_df.bpm, ranges), as_index=False).count()
z2 = pd.to_datetime(z1['bpm'], unit='s')
z2 = pd.DataFrame(z2.dt.strftime('%H:%M:%S'))
display(z2['bpm'])

0    00:05:56
1    00:27:52
2    00:12:57
3    00:06:35
4    00:03:57
Name: bpm, dtype: object

z2 is exactly the heart rate zone data, as we wanted. Now we want to set up two more python arrays that we will use as the y-axis labels, shown here. This data is already available in our z2 series, so we can extract them as necessary.

zones = ['1','2','3','4','5']
times = [e for e in z1['bpm']][::-1]
time_labels = [e for e in z2['bpm']]

The zones array will be the left side y-axis label, and the times will be the actual data. We can place these in a pandas dataframe so it will be easier to plot using the python matplotlib library.

import pandas as pd

d = {"zones":zones, "times": times}
df1 = pd.DataFrame(data=d)
df1.set_index('zones', inplace=True)

And now we can finally plot the data. We will be using the horizontal bar plot from matplotlib, barh. Since we already have the data we need, we can use barh to directly plot the necessary bars, after which we can apply our desired styling to make the graph look like the one pictured.

from matplotlib import pyplot as plt, patches
import seaborn as sns

sns.set_theme(style="dark")
fig1, ax = plt.subplots(figsize=(16,6))
clrs = ["#c0c8c8", "#46c7ee", "#6acc2b", "#f9bf1c", "#de0f5b"][::-1] #colors array for the zones

#plot the base graph
plt. margins(y=0)
plt.xlim(0, sum(i for i in df1["times"])) #this sets the graph to have its width equal to the total session time
bar1 = ax.barh(df1.index[::-1], df1["times"][::-1], align="center", height=1., color = clrs[::-1]) #plot the bar graph
plt.xticks(color='w') #hide the x tick labels
plt.yticks(color='w')
ax.tick_params(axis='both', labelsize=25, pad=12.5)

# Adding the times for each zone
for i in range(len(time_labels)):
  plt.text(.97,0.15 + 0.15*i,time_labels[i],fontsize=18,transform=fig1.transFigure,
          horizontalalignment='center', weight=300)

#add colored backgrounds behind each zone
rect = patches.Rectangle((0, 0.805), 0.5, 0.194, color='#de0f5b')
rect2 = patches.Rectangle((0, 0.602), 0.5, 0.198, color='#f9bf1c')
rect3 = patches.Rectangle((0, 0.402), 0.5, 0.194, color='#6acc2b')
rect4 = patches.Rectangle((0, 0.202), 0.5, 0.194, color='#46c7ee')
rect5 = patches.Rectangle((0, 0.0), 0.5, 0.195, color='#c0c8c8')

lax = fig1.add_axes([0.08,0.126,1.,0.752], anchor='NE', zorder=-1)
lax.add_patch(rect)
lax.add_patch(rect2)
lax.add_patch(rect3)
lax.add_patch(rect4)
lax.add_patch(rect5)
lax.axis('off')

#add gray rectangles behind times
b1 = patches.Rectangle((0, 0.805), 0.5, 0.194, color='#eaeaf2')
b2 = patches.Rectangle((0, 0.602), 0.5, 0.194, color='#eaeaf2')
b3 = patches.Rectangle((0, 0.402), 0.5, 0.194, color='#eaeaf2')
b4 = patches.Rectangle((0, 0.202), 0.5, 0.194, color='#eaeaf2')
b5 = patches.Rectangle((0, 0.0), 0.5, 0.194, color='#eaeaf2')
rax = fig1.add_axes([.8,0.126,.42,0.752], anchor='NE', zorder=-1)
rax.add_patch(b1)
rax.add_patch(b2)
rax.add_patch(b3)
rax.add_patch(b4)
rax.add_patch(b5)
rax.axis('off')

#set grid lines separating bars
ax.set_yticks([0.5,1.5,2.5,3.5,4.5], minor=True)
ax.yaxis.grid(True, which='minor')

This plot is important because it gives us a quick way to understand the general distribution of activity levels throughout a participant’s exercise session, and can display how vigorous the exercise session was.

4.2 Visualizing Heart Rate Averages per Session over a month

Another graph that polar displays is a bar chart showing the duration of exercise sessions throughout the past thirty days along with a plot of the average heart rates per session.

To start we want to extract the necessary data. This will include:

• Dates of exercises

• Durations of each exercise

• Average heart rates of each exercise

import pandas as pd
import numpy as np

durations = []
avg_rates = []
dates = []
day_array = list(df2['day'])
rate_array = list(df2['avg_rate'])
duration_array = list(df2['minutes'])
for i in range(len(df2['day'])):

  if datetime.strptime(day_array[i],'%Y-%M-%d') < (datetime.strptime(start_date,'%Y-%M-%d') + timedelta(days=30)) and datetime.strptime(day_array[i],'%Y-%M-%d') >= datetime.strptime(start_date,'%Y-%M-%d'):
    dates.append(np.datetime64(day_array[i]))
    durations.append(duration_array[i]*60000)
    avg_rates.append(rate_array[i])

Now that we have everything we need, we can move on to plotting the graph!

To do this we will start by plotting a matplotlib bar graph, which will represent the durations of each exercise. Then we will overlay this plot with a second graph, the matplotlib line plot. This is all we need for the base graph!

All that remains after is to add the necessary styling to achieve the graph pictured above.

import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.dates as mdates

#set the style of the plot, initialize it
sns.set_style("whitegrid")
fig, ax = plt.subplots(figsize=(18,8))
ax.set_facecolor('#f2f2f2')

#before we begin, we need to add in a buffer at either end of the data so we can extend the graph
rbound = np.datetime64((pd.Timestamp(dates[len(dates)-1]) + pd.DateOffset(days=1)).strftime('%Y-%m-%d'))
lbound = np.datetime64((pd.Timestamp(dates[0]) - pd.DateOffset(days=1)).strftime('%Y-%m-%d'))
xdates = list(pd.date_range(start=lbound, end=rbound))
xdates.append(xdates[len(xdates)-1])

#process and fill our data arrays so that they match dimensions of x axis
points = {np.datetime64(datetime.strftime(pd.to_datetime(date), '%Y-%m-%d')):(avg_rate,duration) for date, avg_rate, duration in zip(dates,avg_rates,durations)}
data1 = []
data2 = {}
for i in range(len(xdates)):
  cdate = np.datetime64(datetime.strftime(xdates[i],'%Y-%m-%d'))
  if i == 0:
    added = points[np.datetime64(datetime.strftime(xdates[1],'%Y-%m-%d'))][0]
  if cdate in points:
    added = points[cdate][0]
    data1.append(added)
    data2[cdate] = points[cdate][1]/3600000
  else:
    data1.append(added)
    data2[cdate] = 0
plt.xticks(rotation=45)
ax2 = ax.twinx()

#plot the bar graph
ax.bar(data2.keys(),data2.values(), color ='#e71735', width=0.20)
ax.set_yticks(np.arange(0, max(data2.values())+0.5, 0.5))

#plot the average heart rate
plt.plot(xdates, data1, color='#e71735', linestyle='dashed')
ax.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

#clip the graph ends to get our desired look
axis = plt.axis()
plt.xlim(xdates[0], xdates[len(xdates)-2])
ax2.yaxis.grid(False)

#hide ticks
ax2.tick_params(axis='y', length=0)
ax.tick_params(axis='y', length=0)

print(points)

{numpy.datetime64('2022-03-11'): (148.60438074901958, 2880000), numpy.datetime64('2022-03-14'): (125.03413026104663, 2760000), numpy.datetime64('2022-03-27'): (90.8876123319256, 2700000), numpy.datetime64('2022-03-17'): (109.82817716867521, 2700000), numpy.datetime64('2022-05-15'): (134.32337357988095, 3000000), numpy.datetime64('2022-03-16'): (101.67126572452109, 3000000), numpy.datetime64('2022-03-28'): (89.30156282473392, 2880000), numpy.datetime64('2022-03-22'): (100.26106997019164, 3180000), numpy.datetime64('2022-03-13'): (77.30444661402656, 2760000), numpy.datetime64('2022-03-26'): (117.12185927449448, 2820000), numpy.datetime64('2022-04-30'): (100.18888458284101, 2700000), numpy.datetime64('2022-03-10'): (102.11270955167122, 3120000), numpy.datetime64('2022-04-20'): (100.23515912232084, 2940000), numpy.datetime64('2022-03-15'): (151.51392114667934, 3240000), numpy.datetime64('2022-05-12'): (122.15268124430676, 3240000), numpy.datetime64('2022-03-23'): (158.26683088729507, 2820000), numpy.datetime64('2022-04-15'): (131.07266595834122, 2760000), numpy.datetime64('2022-06-15'): (112.44969891174973, 2880000), numpy.datetime64('2022-05-14'): (133.4350159451914, 3180000), numpy.datetime64('2022-06-16'): (85.15969462663257, 2760000), numpy.datetime64('2022-05-26'): (91.9841307333696, 2700000), numpy.datetime64('2022-04-28'): (125.21753249314493, 2940000), numpy.datetime64('2022-05-27'): (107.13397021692715, 2760000), numpy.datetime64('2022-05-21'): (143.3458155760211, 3120000), numpy.datetime64('2022-04-25'): (151.4953087021972, 3360000), numpy.datetime64('2022-03-30'): (153.00516777290125, 3180000), numpy.datetime64('2022-05-29'): (103.10242645848972, 2760000), numpy.datetime64('2022-05-28'): (146.25693522751908, 3180000), numpy.datetime64('2022-04-10'): (86.58026740402435, 2940000), numpy.datetime64('2022-04-17'): (132.33100387749926, 2820000), numpy.datetime64('2022-03-29'): (119.68015580244948, 3360000), numpy.datetime64('2022-04-29'): (161.85980034725586, 2940000), numpy.datetime64('2022-03-21'): (139.88037913822555, 2820000), numpy.datetime64('2022-05-22'): (151.24360941968666, 2760000), numpy.datetime64('2022-05-16'): (126.88316898584463, 2880000), numpy.datetime64('2022-03-20'): (140.09253679690738, 2760000), numpy.datetime64('2022-04-11'): (125.91470056917359, 2760000), numpy.datetime64('2022-03-31'): (127.3760588913292, 3060000), numpy.datetime64('2022-05-23'): (147.5135913867359, 3000000), numpy.datetime64('2022-06-13'): (106.68748050611664, 2820000), numpy.datetime64('2022-04-12'): (109.3653383639486, 3360000), numpy.datetime64('2022-03-24'): (115.72448231096944, 3120000), numpy.datetime64('2022-04-19'): (83.08440428625981, 3240000), numpy.datetime64('2022-03-12'): (110.22994511488777, 3300000), numpy.datetime64('2022-05-30'): (87.58760730290345, 2880000), numpy.datetime64('2022-05-10'): (119.65725556045116, 3060000), numpy.datetime64('2022-04-16'): (126.39131742352002, 3000000), numpy.datetime64('2022-05-24'): (87.3043198208121, 2820000), numpy.datetime64('2022-03-19'): (153.789752043585, 3240000), numpy.datetime64('2022-05-25'): (114.18960726652924, 3000000), numpy.datetime64('2022-05-11'): (149.4561437472248, 3240000), numpy.datetime64('2022-05-13'): (132.18352710565964, 3240000), numpy.datetime64('2022-05-17'): (107.4544674966249, 2880000), numpy.datetime64('2022-06-11'): (110.7521333275521, 2820000), numpy.datetime64('2022-04-27'): (141.06358635016585, 2880000), numpy.datetime64('2022-06-09'): (74.70970628375426, 2940000), numpy.datetime64('2022-06-10'): (74.80748491633975, 2880000), numpy.datetime64('2022-03-25'): (120.97788470827379, 3180000), numpy.datetime64('2022-04-22'): (108.7103048662486, 2880000), numpy.datetime64('2022-05-31'): (127.00566093187098, 3120000), numpy.datetime64('2022-04-24'): (147.07278692696016, 3180000), numpy.datetime64('2022-05-20'): (157.0421130777547, 3000000), numpy.datetime64('2022-04-18'): (136.09561544189722, 3180000), numpy.datetime64('2022-06-12'): (106.69587741212065, 3300000), numpy.datetime64('2022-06-14'): (123.70767138220371, 3420000), numpy.datetime64('2022-04-23'): (111.68947949702499, 3540000), numpy.datetime64('2022-04-14'): (107.27921624551617, 3420000), numpy.datetime64('2022-04-21'): (128.67343246012445, 3240000), numpy.datetime64('2022-04-13'): (129.66215608850723, 3360000)}

This plot is important because it shows us whether the participant wore the monitor and how long they wore the monitor each day. Additionally, their average heart rate throughout the session can be a basic indicator of the intensity of their exercises.

4.3 Continuous Heart Rate Graph

One final visualization that Polar offers is a graph of the continuous heart rate data per session of exercise, shown below.

In this plot the line plot represents the heart rate (bpm) at every second of the exercise session, and the colored regions indicate the “zone” in which the heart rate (bpm) lies. We will need:

• Continuous heart rate data

import pandas as pd

heart_rts = list(hr_df["bpm"])
#we must add in '2000-01-10' as a dummy date to help with formatting concerns
times = list(hr_df["time"])

We will now graph this using the matplotlib line plot function, which will give us the base graph, after which we can focus on aesthetics.

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

sns.set_style("whitegrid")
fig, ax = plt.subplots(figsize=(18,6))
ax.set_facecolor('#f2f2f2')

ax = plt.gca()
plt.plot(times, heart_rts, color='#e71735')

#set x axis values
xts = np.arange(min(times), max(times), timedelta(minutes=10))
ax.set_xticks(xts)
ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))

#set left y axis
yts = np.arange(100, 220, 20)
ax.set_yticks(yts)

#ax.set_yscale('function', functions=(forward,inverse))
ax.axhspan(100, 120, facecolor='#e3e5e5', alpha=0.5)
ax.axhspan(120, 140, facecolor='#bee5f1', alpha=0.5)
ax.axhspan(140, 160, facecolor='#c9e7b6', alpha=0.5)
ax.axhspan(160, 180, facecolor='#f4e3b1', alpha=0.5)
ax.axhspan(180, 200, facecolor='#ecaec4', alpha=0.5)

#right side y axis
f = lambda x: (x/200)*100
i = lambda x: (x/100)*200
ax2 = ax.secondary_yaxis("right", functions=(f,i))
yticks = ax2.yaxis.get_major_ticks()
yticks[1].set_visible(False)
ax.tick_params(axis='y', colors='red')
ax2.tick_params(axis='y', colors='red', length=0)

#cut off the edges
plt.xlim(times[0], times[len(times)-1])

(738449.0, 738449.0401388889)

This plot is important because it shows us the heart rate of the participant for each second, giving us a detailed look at a particular workout session. Also, the colored zones allow us to quickly identify the level of vigor at which the user exercised.

8. Outlier Detection and Data Cleaning

NOTICE: If you are using synthetically generated data, the analyses may yield unintuitive results due to the randomly generated nature of the data

In this section, we will detect outliers in our extracted data.

Since there are currently no outliers (by construction, since it is simulated to have none), we will manually inject a couple.

import pandas as pd

outliers = {'time': ['01:11:52', '01:11:53'], 'bpm': [240, 50]}
hr = hr_df.append(pd.DataFrame(outliers),ignore_index=True)

Now we need a method to find these outliers. While there are many known methods for finding outliers, one of the quickest ways is to use something known as a z-score.

The formula for the z-score is

Z = (data point - mean of data set) / standard deviation

We will compute this formula on every data point in our data set, giving us a massive list of z-score values.

What the z score tells us is how far, in terms of standard deviations, any one data point is from the mean. In statistics, data points at increasing standard deviations away from the mean are less likely to occur. This means that the greater our z-score, the more unlikely for the data point to be real.

We can decide what distance from the mean constitutes an outlier by setting the threshold value, so if a data point’s z-score is above that threshold we might guess it to be an outlier.

One intuitive idea is that the heart rate should not vary too much from one second to the next. To capitalize on this, we can improve our analysis by first calculating the distance each heart rate sample is away from neighboring heart rate samples. We can then use the z-scores to determine which heart rates differ from their neighbors by an unusually large margin. These will likely be our outliers.

import math
import pandas as pd

nx = hr["bpm"]

distances = []

f = lambda x: x**2

#take the current count, subtract it by the points to its left and right some distance, then
#square each of these differences, add, then sqrt
for l1, l2, m, r1, r2 in zip(nx[:-4],nx[1:-3],nx[2:-2],nx[3:-1],nx[4:]):
  difference = math.sqrt(f(m-l1)+f(m-l2)+f(m-r1)+f(m-r2))
  distances.append(difference)

Now we can calculate our z-scores for the heart rate data in one training session, and identify those data points that seem like outliers. Luckily, we can save some time by using the z-score function from scipy’s statistics library, which will give us a z-score calculation for all data points.

from scipy import stats
import pandas as pd
import numpy as np

threshold = 4 #@param{type: "integer"}
z_scores = np.abs(stats.zscore(nx))
hr['z'] = z_scores
print(hr.loc[hr['z'] > threshold])

          time  bpm         z
3469  01:11:52  240  4.892948
3470  01:11:53   50  4.448598

It looks like the code was able to successfully identify the two outliers!

9. Data Analysis

9.1 Correlation between average heart rate and calories burned

We want to test the hypothesis that for a session, the average heart rate correlates with the number of calories burned.

Let us begin by extracting the necessary data. This includes:

• the average heart rate for each exercise session

• the calories burned per minute for each exercise session

To do this we will refer to the exercise_data list of exercises which we kept from section 3. In practice you can use the following code:

avg_rates = [s['heart_rate']['average'] for s in exercise_data]
calories_per_minute = [s["calories"]/(float(pd.Timedelta(s["duration"]).total_seconds()/60)) for s in exercise_data]

import pandas as pd

avg_rates = [avg_rate for avg_rate in df2['avg_rate']]
calories_per_minute = [c/(m) for c, m in zip(df2['calories'],df2['minutes'])]

And now we can actually graph the data to see the correlation. This is done by using a seaborn linear regression plot which we can use to visualize the correlation.

import seaborn as sns
import matplotlib.pyplot as plt

# plot the data
d = {'calories burned (kcal/minute)': calories_per_minute, 'average heart rate (bpm)': avg_rates}
df = pd.DataFrame(data=d)

graph = sns.lmplot(data=df, y='calories burned (kcal/minute)', x='average heart rate (bpm)')

graph.map(plt.scatter, 'average heart rate (bpm)','calories burned (kcal/minute)', edgecolor ="w").add_legend()
sns.set(context='notebook', style='whitegrid', font='sans-serif', font_scale=1, color_codes=True)

plt.show(block=True)

From the plot, we can see that the line of best fit actually fits all the data points rather tightly, and this is our first indicator that the correlation between the two factors exists.

But we can do more to confirm this correlation - we can calculate the p-value to determine statistical significance using scipy’s stats library, as shown here:

from scipy import stats

p_value = stats.linregress(calories_per_minute,avg_rates)[3]
correlation_coefficient = stats.pearsonr(calories_per_minute, avg_rates)[0]
print("P value: " + str(p_value))
print("Correlation coefficient: " + str(correlation_coefficient))

P value: 0.0033113838558418668
Correlation coefficient: 0.952639690187336

Based on the p-value 0.0023 < 0.05, we have statistical significance. This means that there is 0.2% probability that the datapoints of our dataset occurred by chance. In addition, our correlation coefficient of 0.96 implies a very strong correlation. This relationship is intuitive, because Polar uses an equation to calculate the calories burned. But looking at this data allows to see that polar’s calorie counter makes sense!

9.2 Correlation between duration of exercise and calories burned

Now we can clearly see the that average heart rate and calories burned per minute are tightly correlated. But how about the correlation between the duration of an exercise and the amount calories burned? We will first extract the necessary data through the code segment shown below:

durations = [float(pd.Timedelta(s["duration"]).total_seconds()) for s in exercise_data]
calories = [s["calories"] for s in exercise_data]

As before, our example will use manually inputted data points (which come from real sessions), but feel free to replace the code with the above segment to test out your data!

import pandas as pd

durations = [m*60 for m in df2['minutes']]
calories = [c for c in df2['calories']]

Now we will recreate the graph plotted in 5.1 this time using our new durations and calories data.

import pandas as pd
import seaborn as sns

# plot the data
d = {'calories burned (kcal)': calories, 'duration (s)': durations}
df = pd.DataFrame(data=d)

graph = sns.lmplot(data=df, y='calories burned (kcal)', x='duration (s)')

graph.map(plt.scatter, 'duration (s)','calories burned (kcal)', edgecolor ="w").add_legend()
sns.set(context='notebook', style='whitegrid', font='sans-serif', font_scale=1, color_codes=True)

plt.show(block=True)

p_value = stats.linregress(calories,durations)[3]
correlation_coefficient = stats.pearsonr(calories, durations)[0]
print("P value: " + str(p_value))
print("Correlation coefficient: " + str(correlation_coefficient))

P value: 0.02055818281629383
Correlation coefficient: 0.8805265461816261

Since our P-value is 0.017 < 0.05, meaning the we have 1.7% probability that our data occurred by chance, we again have statistical significance. Interestingly enough, our correlation coefficient for duration to calories burned is less than our correlation coefficient for average heart rate to calories burned (0.889 < 0.96). Although we can’t say from our simple analysis that average heart rate is more correlated to calories burned than duration is (further analysis has to be done to determine whether the difference in correlation coefficient values is significant, and we also should get more data before making a conclusion!), it definitely reveals a question for further investigation!

9.3. T test on activity intensity through workout

We might also be interested in determining whether the participant exercises with different overall intensity at different halves of their sessions.

Choosing on exercise session to focus on, we can use the T-test to see if the exercise intensity in the first half of the session is different from the intensity in the second half.

import pandas as pd
import numpy as np
from scipy.stats import ttest_ind

def differ(df):
  rows = len(df.index)
  f = df.iloc[range(rows//2)]["bpm"]
  s = df.iloc[range(rows//2,rows)]["bpm"]
  return ttest_ind(f,s)

print(differ(hr_df))

Ttest_indResult(statistic=-33.27175316357228, pvalue=6.677761019370051e-211)

Since the t-test returned a negative value, the sample mean in the first half of the exercise session is less than the sample mean in the second half of the session. Further, since our pvalue is small (< 0.05), we can reject the null hypothesis and conclude that the the first half of the exercise tended to have lower heart rates than the second half.

We can now repeat this process on multiple exercise sessions and visualize our findings in a plot, giving us a view of the participant’s exercise habits (i.e. whether they tend to slow down throughout a session, or whether their intensity increases).

import pandas as pd
import numpy as np
from scipy.stats import ttest_ind
import matplotlib.pyplot as plt

y = []
x = df2["day"]

for df in hr_dfs:
  v,p = differ(df)
  if p < 0.05:
    y.append(v)

plt.figure(figsize=(16,10))
plt.plot(x, y, drawstyle="steps-mid")
plt.title(f"T-test of exercise intensity between first half and second half of session from {start_date} to {datetime.strftime(end,'%Y-%m-%d')}",fontsize=20)
plt.xlabel("Time")
plt.ylabel("T-test of exercise intensity between first half and second half of session")

Text(0, 0.5, 'T-test of exercise intensity between first half and second half of session')