Banister Impulse~Response Model in R – Part 1
Before you start reading this post, please read EXCELLENT paper by Clark and Skiba, especially on the topic of Banister impulse-response model.
I decided to write code in R, but also allow for multivariate analysis (where impulse can be multiple variables, as is the case in sports) which can speed the thing up. Might end up being a package after I tune it up a bit and add model diagnostics.
I had troubles with optimizer function (optim) so I need to sort that one out, especially because with some types of random data it doesn’t produce results.
Taking all these into account I won’t be showing the code of this potential R package and I will use data presented in Excel workbook that accompany above mentioned paper.
Here is the data:
## Date BikeScore X5.min.Test.Power
## 1 2004-01-29 0.00 NA
## 2 2004-01-30 46.46 NA
## 3 2004-01-31 26.70 NA
## 4 2004-02-01 177.37 NA
## 5 2004-02-02 0.00 NA
## 6 2004-02-03 36.64 NA
## 7 2004-02-04 62.52 NA
## 8 2004-02-05 52.46 NA
## 9 2004-02-06 41.00 NA
## 10 2004-02-07 0.00 NA
## 11 2004-02-08 29.98 NA
## 12 2004-02-09 0.00 NA
## 13 2004-02-10 38.17 NA
## 14 2004-02-11 75.23 NA
## 15 2004-02-12 35.73 NA
## 16 2004-02-13 0.00 NA
## 17 2004-02-14 39.03 NA
## 18 2004-02-15 0.00 NA
## 19 2004-02-16 51.51 NA
## 20 2004-02-17 61.12 NA
## 21 2004-02-18 31.32 NA
## 22 2004-02-19 0.00 NA
## 23 2004-02-20 14.37 NA
## 24 2004-02-21 38.40 NA
## 25 2004-02-22 0.00 NA
## 26 2004-02-23 0.00 NA
## 27 2004-02-24 0.00 NA
## 28 2004-02-25 0.00 NA
## 29 2004-02-26 60.31 NA
## 30 2004-02-27 38.45 NA
## 31 2004-02-28 0.00 NA
## 32 2004-02-29 98.57 NA
## 33 2004-03-01 45.64 NA
## 34 2004-03-02 0.00 NA
## 35 2004-03-03 21.03 NA
## 36 2004-03-04 23.54 NA
## 37 2004-03-05 55.53 NA
## 38 2004-03-06 43.31 NA
## 39 2004-03-07 0.00 NA
## 40 2004-03-08 88.48 NA
## 41 2004-03-09 35.77 NA
## 42 2004-03-10 0.00 NA
## 43 2004-03-11 0.00 NA
## 44 2004-03-12 0.00 NA
## 45 2004-03-13 46.83 NA
## 46 2004-03-14 88.10 271
## 47 2004-03-15 0.00 NA
## 48 2004-03-16 32.09 NA
## 49 2004-03-17 65.53 274
## 50 2004-03-18 0.00 NA
## 51 2004-03-19 32.09 NA
## 52 2004-03-20 46.55 NA
## 53 2004-03-21 51.06 NA
## 54 2004-03-22 0.00 NA
## 55 2004-03-23 96.66 NA
## 56 2004-03-24 0.00 NA
## 57 2004-03-25 53.62 282
## 58 2004-03-26 0.00 NA
## 59 2004-03-27 0.00 NA
## 60 2004-03-28 0.00 NA
## 61 2004-03-29 50.95 NA
## 62 2004-03-30 59.19 NA
## 63 2004-03-31 50.53 NA
## 64 2004-04-01 39.33 NA
## 65 2004-04-02 0.00 NA
## 66 2004-04-03 0.00 NA
## 67 2004-04-04 0.00 NA
## 68 2004-04-05 55.38 NA
## 69 2004-04-06 0.00 NA
## 70 2004-04-07 25.19 NA
## 71 2004-04-08 85.53 290
## 72 2004-04-09 0.00 NA
## 73 2004-04-10 0.00 NA
## 74 2004-04-11 0.00 NA
## 75 2004-04-12 61.28 NA
## 76 2004-04-13 0.00 NA
## 77 2004-04-14 41.69 NA
## 78 2004-04-15 90.10 292
## 79 2004-04-16 0.00 NA
## 80 2004-04-17 47.40 NA
## 81 2004-04-18 52.16 NA
## 82 2004-04-19 0.00 NA
## 83 2004-04-20 0.00 NA
## 84 2004-04-21 49.22 296
## 85 2004-04-22 0.00 NA
## 86 2004-04-23 0.00 NA
## 87 2004-04-24 48.26 NA
## 88 2004-04-25 0.00 NA
## 89 2004-04-26 57.94 NA
## 90 2004-04-27 0.00 NA
## 91 2004-04-28 75.17 NA
## 92 2004-04-29 0.00 NA
## 93 2004-04-30 26.70 NA
## 94 2004-05-01 0.00 NA
## 95 2004-05-02 51.99 NA
## 96 2004-05-03 0.00 NA
## 97 2004-05-04 86.60 NA
## 98 2004-05-05 54.09 NA
## 99 2004-05-06 0.00 NA
## 100 2004-05-07 43.19 NA
## 101 2004-05-08 0.00 NA
## 102 2004-05-09 100.79 NA
## 103 2004-05-10 40.49 NA
## 104 2004-05-11 0.00 NA
## 105 2004-05-12 0.00 NA
## 106 2004-05-13 62.46 299
## 107 2004-05-14 71.87 NA
## 108 2004-05-15 0.00 NA
## 109 2004-05-16 36.21 NA
## 110 2004-05-17 46.51 NA
## 111 2004-05-18 0.00 NA
## 112 2004-05-19 0.00 NA
## 113 2004-05-20 72.89 NA
## 114 2004-05-21 69.89 NA
## 115 2004-05-22 0.00 NA
## 116 2004-05-23 167.60 296
## 117 2004-05-24 39.47 NA
## 118 2004-05-25 0.00 NA
## 119 2004-05-26 43.38 NA
## 120 2004-05-27 17.63 NA
## 121 2004-05-28 0.00 NA
## 122 2004-05-29 60.49 NA
## 123 2004-05-30 98.95 NA
## 124 2004-05-31 50.46 NA
## 125 2004-06-01 0.00 NA
## 126 2004-06-02 25.00 NA
## 127 2004-06-03 0.00 NA
## 128 2004-06-04 85.41 NA
## 129 2004-06-05 0.00 NA
## 130 2004-06-06 136.27 NA
## 131 2004-06-07 0.00 NA
## 132 2004-06-08 0.00 NA
## 133 2004-06-09 51.17 302
## 134 2004-06-10 0.00 NA
## 135 2004-06-11 99.37 NA
## 136 2004-06-12 0.00 NA
## 137 2004-06-13 139.22 NA
## 138 2004-06-14 0.00 NA
## 139 2004-06-15 0.00 NA
## 140 2004-06-16 0.00 NA
## 141 2004-06-17 0.00 NA
## 142 2004-06-18 0.00 NA
## 143 2004-06-19 0.00 NA
## 144 2004-06-20 0.00 NA
## 145 2004-06-21 0.00 NA
## 146 2004-06-22 0.00 NA
## 147 2004-06-23 0.00 NA
## 148 2004-06-24 0.00 NA
## 149 2004-06-25 0.00 NA
## 150 2004-06-26 0.00 NA
## 151 2004-06-27 0.00 NA
## 152 2004-06-28 0.00 NA
## 153 2004-06-29 0.00 NA
## 154 2004-06-30 0.00 NA
## 155 2004-07-01 0.00 NA
## 156 2004-07-02 0.00 NA
## 157 2004-07-03 0.00 NA
## 158 2004-07-04 0.00 NA
## 159 2004-07-05 0.00 NA
## 160 2004-07-06 0.00 NA
## 161 2004-07-07 0.00 NA
## 162 2004-07-08 0.00 NA
## 163 2004-07-09 0.00 NA
## 164 2004-07-10 0.00 NA
## 165 2004-07-11 0.00 NA
What are we basically trying to do is to predict 5 min power on the bike (which is Response or training effect) measured couple of times during this period, by using Bike Score (which is Impulse or training load) measure daily. The Banister model is not a “black box”, but rather give us the parameters for positive training effect (PTE) and negative training effect (NTE). Those parameters are tau1 and k1 for PTE and tau2 and k2 for NTE. This is essentially two-factor model of preparedness. Pretty intuitive for coaches as well – because they know that every training load cause both PTE and NTE. Highly recommend checking the article for more info.
So I have created a simple functions in R that do just that. Let’s create a model
skiba.model <- train.Banister(data=Skiba,
impulse="BikeScore",
response = "X5.min.Test.Power",
dates = Skiba$Date)
str(skiba.model)## List of 4
## $ Impulse.Variables : chr "BikeScore"
## $ Response.Variables: chr "X5.min.Test.Power"
## $ Impulse :'data.frame': 165 obs. of 1 variable:
## ..$ BikeScore: num [1:165] 0 46.5 26.7 177.4 0 ...
## $ X5.min.Test.Power :List of 5
## ..$ p0 : num 261
## ..$ BikeScore :'data.frame': 1 obs. of 4 variables:
## .. ..$ tau1: num 28.5
## .. ..$ k1 : num 1.29
## .. ..$ tau2: num 26.3
## .. ..$ k2 : num 1.34
## ..$ Response :'data.frame': 165 obs. of 3 variables:
## .. ..$ date : POSIXct[1:165], format: "2004-01-29" ...
## .. ..$ X5.min.Test.Power: int [1:165] NA NA NA NA NA NA NA NA NA NA ...
## .. ..$ predicted : num [1:165] 261 258 257 250 251 ...
## ..$ SumOfSquares : num 18.3
## ..$ Variance.Explained: num 0.981When we plot the model we get pretty nice model with great variance explained (R squared)
plot.banister(skiba.model)
The grey line is the response (5min bike power) predicted by our model. As you can see it is pretty accurate. We can use model parameters to make informed decisions about the taper duration and generally about the design of the training cycle. This is what data-driven decision making actually means.
Let’s create our own data (just simple random data generation): in this case bench press 1RM testing [response] and bench press 1RM tonnage [impulse]. But why tonnage? Why not volume over 80%? That is a GREAT question – the point is that we can select one or more metrics that gives us the best prediction without overfitting (this is a term from machine learning, basically model is very sensitive to noise in the data). So maybe we do collect multiple different impulse metrics, but over time we select only few for a given individual that have the best prediction.
Ok, let’s create our data for bench press
number.of.days <- 60
Date <- dmy("1/1/2014") + ddays(1:number.of.days)
# Response
benchpress.1RM <- 100 + (1:number.of.days)/3 + runif(min=0, max=5, n=number.of.days)
benchpress.1RM[runif(min=1, max=60, n=number.of.days * 1.5)] <- NA # Add a lot of NAs
# Impulse
benchpress.tonnage <- runif(min=0, max=3000, n=number.of.days)
benchpress.tonnage[benchpress.tonnage<500] <- 0
training.data <- data.frame(Date, benchpress.tonnage, benchpress.1RM)
print(training.data)## Date benchpress.tonnage benchpress.1RM
## 1 2014-01-02 1843.9 NA
## 2 2014-01-03 1671.5 102.5
## 3 2014-01-04 986.3 NA
## 4 2014-01-05 1359.4 NA
## 5 2014-01-06 1501.3 NA
## 6 2014-01-07 542.6 106.5
## 7 2014-01-08 1588.9 NA
## 8 2014-01-09 0.0 NA
## 9 2014-01-10 833.3 NA
## 10 2014-01-11 638.1 103.6
## 11 2014-01-12 854.4 NA
## 12 2014-01-13 2685.3 NA
## 13 2014-01-14 1338.7 NA
## 14 2014-01-15 2340.0 NA
## 15 2014-01-16 2641.9 NA
## 16 2014-01-17 1239.4 NA
## 17 2014-01-18 0.0 109.3
## 18 2014-01-19 1006.5 NA
## 19 2014-01-20 2171.2 108.2
## 20 2014-01-21 1012.8 NA
## 21 2014-01-22 1891.2 NA
## 22 2014-01-23 2521.8 NA
## 23 2014-01-24 2568.4 110.9
## 24 2014-01-25 1174.1 NA
## 25 2014-01-26 1141.5 NA
## 26 2014-01-27 2686.3 NA
## 27 2014-01-28 1932.9 NA
## 28 2014-01-29 2223.2 NA
## 29 2014-01-30 1815.9 NA
## 30 2014-01-31 2709.2 NA
## 31 2014-02-01 881.2 NA
## 32 2014-02-02 573.8 NA
## 33 2014-02-03 2659.4 113.5
## 34 2014-02-04 1510.0 NA
## 35 2014-02-05 2631.2 115.8
## 36 2014-02-06 567.6 NA
## 37 2014-02-07 2274.3 116.3
## 38 2014-02-08 2173.5 NA
## 39 2014-02-09 2831.2 NA
## 40 2014-02-10 1642.9 115.4
## 41 2014-02-11 2135.2 NA
## 42 2014-02-12 1166.7 NA
## 43 2014-02-13 0.0 NA
## 44 2014-02-14 2781.9 NA
## 45 2014-02-15 849.7 NA
## 46 2014-02-16 1771.7 NA
## 47 2014-02-17 0.0 115.8
## 48 2014-02-18 2521.5 NA
## 49 2014-02-19 953.9 120.0
## 50 2014-02-20 2348.6 NA
## 51 2014-02-21 802.5 NA
## 52 2014-02-22 655.9 NA
## 53 2014-02-23 1550.4 NA
## 54 2014-02-24 806.9 NA
## 55 2014-02-25 543.5 NA
## 56 2014-02-26 1555.7 119.2
## 57 2014-02-27 1688.3 NA
## 58 2014-02-28 0.0 NA
## 59 2014-03-01 769.1 NA
## 60 2014-03-02 2153.8 122.0Now we create a model from this data and plot it
strength.model <- train.Banister(data=training.data,
impulse="benchpress.tonnage",
response = "benchpress.1RM",
dates = training.data$Date)
plot.banister(strength.model)
We got decent fit. What about if I did another training that might affect my adaptation? For example I have also trained squats, so I add the tonnage of the squats in there as well.
squat.tonnage <- runif(min=0, max=5000, n=number.of.days)
squat.tonnage[squat.tonnage<500] <- 0
training.data <- data.frame(training.data, squat.tonnage)Now we can add the info from this impulse as well and see if our predictions are better
strength.model <- train.Banister(data=training.data,
impulse=c("benchpress.tonnage", "squat.tonnage"),
response = "benchpress.1RM",
dates = training.data$Date)
plot.banister(strength.model)
As we can see from the graph taking into account squat tonnage gave us more predictive power. The question that needs to be answered is which one is more predictive of bench 1RM and does doing squats give positive or negative transfer to bench 1RM. We can do that by comparing each impulse tau and k scores. Who is better predictor might need to involve normalization of the scores first. I will cover these concepts in the future.
What if I’ve tracked squat 1RMs during that period also? I can plug that into the formula as well and generate a model for that response variable using both bench and squat tonnage. Using the Velocity Based Strength and 1RM estimates one can probably estimate 1RMs daily. Taken together with this model, I think that is very potentially useful approach.
squat.1RM <- 160 + (1:number.of.days)/3 + runif(min=0, max=10, n=number.of.days)
squat.1RM[runif(min=1, max=60, n=number.of.days * 1.5)] <- NA #Add a lot of NAs
training.data <- data.frame(training.data, squat.1RM)
strength.model <- train.Banister(data=training.data,
impulse=c("benchpress.tonnage", "squat.tonnage"),
response = c("benchpress.1RM", "squat.1RM"),
dates = training.data$Date)
plot.banister(strength.model)
This model could be made much more complex by taking into account (which is additional optimization parameter) the saturation effect or even dose-response of he given impulse on PTE. I might put this in the model some other time.
Sine we have measured a lot of variables, one needs to decide which one to include in the model. We can even plug in metrics of recovery (e.g. time spent in ice water) and wellness. Model checking might gives us information what variable gives the best prediction, when models become too complex (AIC criterion variable, etc) and when it overfits. These are all very interesting topics that I am working on at the moment.
Anyone interested in testing their own data from the field give me a head up on the forum or on email.
Responses