THIS IS A DRAFT! It will appear on the main page once finished!

Investigating effect of spinning power output on HR

I was curious how effort exerted during the exercise impacts heart rate and whether that correlates strongly with subjective fatigue and exhaustion.

It seemed that easiest way to experiment would be spinning machine. Running only has single variable, speed. One could f change incline, but that would be harder to predict.

For elliptical machine or rowing machine it would be even more complicated!

With stationary bike, there are two variables that potentially impact power you need to exert:

• revolutions per minute, or angular velocity, that has linear effect on power (no air resistance as it's stationary!)
• resistance, which isn't exactly specified, but one would expect it to be proportional to the force you have to apply (i.e. torque)

The spinning machine I was using displays the pace (i.e. revolutions per minute), so all you have to do is maintain it. In addition it's displaying current resistance level and reports power in watts.

During exercise, I'm using a chest HR tracker, so simplest thing to do would be take whatever power spinning machine reports and try to find correlation with moving total/average of HR.

However, being paranoid and lacking any documentation for the machine, I decided no to trust its numbers blindly and check them instead. Technogym's website doesn't help in clarifying how power is computed. They have some sensible information like:

The power meter must be accurate.

The power measurement must be precise and repeatable. A 3-5 watt error is not significant, but if a system is not reliable there may be deviations of many tens of watts, i.e. equal to or greater than the amount of power that is gained from one year's training.

Let's see how accurate is their power meter!

Initial measurements (2019-11-09)¶

Throughout different exercise sessions, I've taken bunch of measurements of RPM, resistance and power:

Measurements (click to expand)
In [1]:
```# TODO could only collapse inputs?
datas = """
58 10 129
56 10 127
56 10 127
56 8  94
57 8  98
56 8  94
58 10 133
56 10 126
56 8  93
55 8  91
56 8  94
56 10 128
55 10 124
54 10 119
53 8  87
55 8  93
55 8  90
95 8  240

70 10 198
55 8  85
95 8  226
95 8  229
95 8  228
95 8  227

97 8  236
95 8  227
95 8  227
95 8  230
60 10 156
61 10 154
62 10 162
61 10 156
55 10 125
56 10 128
57 8  89
56 8  87
57 8  90
57 8  91
60 8  101
56 10 129
57 10 131
"""

import pandas as pd
def make_df(datas: str) -> pd.DataFrame:
df = pd.DataFrame(
(map(int, l.split()) for l in datas.splitlines() if len(l.strip()) > 0),
columns=['rpm', 'resistance', 'watts'],
)
df.index = df.apply(lambda row: f"{row['rpm']}_{row['resistance']}_{row['watts']}", axis='columns')
return df

df = make_df(datas)
old_df = df.copy()
# btw, if anyone knows a more elegant way of converting such a table in dataframe, I'd be happy to know!
```
In [2]:
```display(df.sample(10, random_state=0))
```
rpm resistance watts
95_8_227 95 8 227
56_8_87 56 8 87
61_10_154 61 10 154
57_8_98 57 8 98
56_8_94 56 8 94
61_10_156 61 10 156
95_8_230 95 8 230
56_10_128 56 10 128
57_8_90 57 8 90
62_10_162 62 10 162

It's reasonable to assume that power depends linearly both on RPM and resistance, so we conjecture `watts = rpm x resistance`. Let's see if it holds against what the spinning machine reports:

In [3]:
```%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
import matplotlib
matplotlib.rc('font', size=17, weight='regular')

def plot_df(df):
dff = df.copy()
dff['rpm x resistance'] = df['rpm'] * df['resistance']

fig, ax = plt.subplots(figsize=(10, 10))
plt.xlim((0, max(dff['rpm x resistance']) * 1.1))
plt.ylim((0, max(dff['watts']) * 1.1))
plt.grid(True)

g = sns.regplot(
data=dff,
x='rpm x resistance', y='watts',
ax=ax,
)

plt.xlabel('Power, theoretical, angular velocity multiplied by resistance')
plt.ylabel('Power, watts as reported by the machine')

plt.show()

plot_df(df)
```
In [4]:
```import statsmodels.api as sm
def print_stats(df):
dff = df.copy()
dff['rpm x resistance'] = df['rpm'] * df['resistance']
X = dff[['rpm x resistance']]
res = sm.OLS(dff['watts'], X).fit()
# TODO ugh. remove Date/Time from summary...
# https://www.statsmodels.org/stable/_modules/statsmodels/iolib/summary.html#Summary
print(res.summary())
print_stats(df)
```
```                            OLS Regression Results
==============================================================================
Dep. Variable:                  watts   R-squared:                       0.982
Method:                 Least Squares   F-statistic:                     2109.
Date:                Sat, 07 Dec 2019   Prob (F-statistic):           1.43e-35
Time:                        23:26:17   Log-Likelihood:                -138.70
No. Observations:                  41   AIC:                             281.4
Df Residuals:                      39   BIC:                             284.8
Df Model:                           1
Covariance Type:            nonrobust
====================================================================================
coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------
const             -110.6016      5.605    -19.732      0.000    -121.939     -99.264
rpm x resistance     0.4406      0.010     45.929      0.000       0.421       0.460
==============================================================================
Omnibus:                        1.196   Durbin-Watson:                   1.017
Prob(Omnibus):                  0.550   Jarque-Bera (JB):                1.049
Skew:                          -0.190   Prob(JB):                        0.592
Kurtosis:                       2.315   Cond. No.                     2.87e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.87e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
```
```/home/karlicos/.local/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
return ptp(axis=axis, out=out, **kwargs)
```

Wow, that doesn't look so well.

Free parameter is about `-100` watts, which is quite a lot considering my high intensity intervals are `250` watts. I don't think it can be explained by friction either: if anything, friction would shift the plot up and make the free coefficient positive.

At this point, I'm not sure what it means. I guess I'll try to make more measurements at really low resistances and speeds to make the model more complete, but I would be too surprised if either watts or resistance reported by the machine are just made up.

More data (2019-11-14)¶

I collected more data corresponding to different resistances/velocities. It's actually quite hard to consistently spin under low resistance setting, so I think I might need one more round of data collection to complete the picture!

More measurements (click to expand)
In [5]:
```new_df = make_df("""
96  4  66
99  4  69
101 6  146
103 6  149
110 6  166
111 6  170
50  13 186
36  13 107
36  13 105
31  10 41
30  10 44
28  10 39
117 8  323
116 8  320
116 8  322
48  6  40
49  6  37
60  2  24
59  2  23
86  2  40
106 2  48
62  5  44
61  5  44
81  5  70
81  5  70
93  5  90
97  5  97
35  12 87
33  12 81
25  12 51
26  12 55
27  1  50
39  8  46
39  8  44
30  8  29
30  8  31
32  8  31
32  8  31
29  8  29
""")
```
In [6]:
```df = pd.concat([old_df, new_df])
plot_df(df)
print_stats(df)
```
```                            OLS Regression Results
==============================================================================
Dep. Variable:                  watts   R-squared:                       0.924
Method:                 Least Squares   F-statistic:                     952.2
Date:                Sat, 07 Dec 2019   Prob (F-statistic):           1.82e-45
Time:                        23:26:17   Log-Likelihood:                -352.79
No. Observations:                  80   AIC:                             709.6
Df Residuals:                      78   BIC:                             714.3
Df Model:                           1
Covariance Type:            nonrobust
====================================================================================
coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------
const              -60.3329      6.129     -9.843      0.000     -72.535     -48.130
rpm x resistance     0.3614      0.012     30.857      0.000       0.338       0.385
==============================================================================
Omnibus:                       57.648   Durbin-Watson:                   1.094
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              238.330
Skew:                           2.326   Prob(JB):                     1.77e-52
Kurtosis:                      10.062   Cond. No.                     1.42e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.42e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
```

Ok, it clearly started diverging from the nice linear dependency, especially at lower values of theoretical power. It's time to try to break it down and see what is to blame: e.g., resistance or speed component, or perhaps some individual measurements.

Analysing data and looking at outliers (2019-11-24)¶

In [7]:
```import statsmodels.api as sm
def plot_influence(df):
# TODO FIXME use it in prev section
res = sm.formula.ols("watts ~ rpm * resistance - resistance - rpm", data=df).fit()
fig, ax = plt.subplots(figsize=(10, 10))
sm.graphics.influence_plot(res, ax=ax, criterion="cooks")
plt.show()
plot_influence(df)
# sm.graphics.plot_partregress_grid(res, fig=fig)
```

TODO hmm, 0.10 is not that high leverage right? Although depends on residual too, and here the residual is very high, so it would have high influence.. https://www.statsmodels.org/dev/examples/notebooks/generated/regression_plots.html

'27_1_50' seems like a typo that should have been '27 12 50', let's drop it and ignore. We also ignore few other points that seem to be outliers.

In [8]:
```fdf = df.copy() # TODO Do i really need to copy?
fdf = fdf.drop([
'27_1_50',
'86_2_40',
'59_2_23',
'60_2_24',
'116_8_322',
'116_8_320',
'117_8_323',
'106_2_48',
'95_8_240',
])
# TODO I wonder what does it do to multiple observations with exactly same input/output params? Is it a set or multiset?
plot_influence(fdf)
```

Ok, that's somewhat better at least in terms of outliers. Let's see if that helps:

In [9]:
```plot_df(fdf)
print_stats(fdf)
```
```                            OLS Regression Results
==============================================================================
Dep. Variable:                  watts   R-squared:                       0.976
Method:                 Least Squares   F-statistic:                     2773.
Date:                Sat, 07 Dec 2019   Prob (F-statistic):           1.90e-57
Time:                        23:26:18   Log-Likelihood:                -257.57
No. Observations:                  71   AIC:                             519.1
Df Residuals:                      69   BIC:                             523.7
Df Model:                           1
Covariance Type:            nonrobust
====================================================================================
coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------
const              -75.1847      3.695    -20.348      0.000     -82.556     -67.813
rpm x resistance     0.3798      0.007     52.655      0.000       0.365       0.394
==============================================================================
Omnibus:                       15.050   Durbin-Watson:                   0.523
Prob(Omnibus):                  0.001   Jarque-Bera (JB):                4.881
Skew:                           0.325   Prob(JB):                       0.0871
Kurtosis:                       1.892   Cond. No.                     1.73e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.73e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
```

So, on one hand that did make fit look more linear. On the other hand we've had to filter out all the low-resistance observations to achieve that.

I guess I'll collect more observations to be absolutely sure.