1 Packages Initialization
In this project, the open-source R programming language is used to examine the relationship between RPE and strength decline. R is maintained by an international team of developers who make the language available at The Comprehensive R Archive Network. Readers interested in reusing our code and reproducing our results should have R installed locally on their machines. R can be installed on a number of different operating systems (see Windows, Mac, and Linux for the installation instructions for these systems). We also recommend using the RStudio interface for R. The reader can download RStudio for free by following the instructions at the link. For non-R users, we recommend the Hands-on Programming with R for a brief overview of the software’s functionality. Hereafter, we assume that the reader has an introductory understanding of the R programming language.
In the code chunk below, we load the packages used to support our analysis. Note that the code of this and any of the code chunks can be hidden by clicking on the ‘Hide’ button to facilitate the navigation. The reader can hide all code and/or download the Rmd file associated with this document by clicking on the Code button on the top right corner of this document. Our input data file can be accessed/ downloaded from blindedResearch/corr_analysis.
rm(list = ls()) # clear global environment
graphics.off() # close all graphics
#p_load is equivalent to combining both install.packages() and library()
pacman::p_load(stringr, readxl, ggplot2, dplyr, rje, signal, stats, tibble, kableExtra, magrittr, polyreg)2 Load Files and Functions
In this section we first load the “Mixed_Data_Full”, the data frame including all RPE scores and Strength decline data for the whole experiment, which involves three work periods and two rest periods. Since the focus of this report is on the first two work periods, the other periods are removed for the later analysis.
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
load('Mixed_Data_Full.RData')
Mixed_Data_Full$Unique_code = paste(Mixed_Data_Full$Unique_code, Mixed_Data_Full$Period)
data_to_remove = which(Mixed_Data_Full$Period == '1st_break' | Mixed_Data_Full$Period == '2nd_break'| Mixed_Data_Full$Period == '3rd_45min')
Mixed_Data_Full = Mixed_Data_Full[-data_to_remove,]
Mixed_Data_List = split(Mixed_Data_Full, f = Mixed_Data_Full$Unique_code)
Mixed_Data = as.data.frame(do.call(rbind, Mixed_Data_List))3 Results
This section includes the process of data preparation, plotting, and correlation analysis of RPE scores and strength decline.
3.1 Period 1
First, the process starts with interpolating RPE scores for the strength test’s time points using the cubic fitting for the first work period. Since the RPE scores are reported every 5 minutes and strength tests are done every 9 minutes as well as the begin and end of each period, this interpolation is needed to have RPE scores at timepoints that they are not originally reported. This interpolation is done using the function “model = lm(Y ~ stats::poly(X, 3))” to fit to current data and then using “Y_fit = predict(model, data.frame(X))” to find the fitted RPE scores for the missing data points. It can be replaced with “model = approx(X, Y, xout = X)” and “Y_fit = do.call(cbind, model)” to have the linear interpolation of the RPE scores. In this report we have proceeded with the cubic fit, however, the results of the linear fit is also provided in the paper and its results does not differ significantly from the results shown in the next steps.
Subjects = unique(Mixed_Data$ID)
Conditions = unique(Mixed_Data$Condition)
Full_Data1 = NULL
for (i in 1:length(Subjects))
{
for (j in 1:length(Conditions))
{
Subject = Subjects[i]
Condition = Conditions[j]
Data = Mixed_Data[Mixed_Data$ID == Subject &
unlist(Mixed_Data$Condition) == Condition & Mixed_Data$Period == '1st_45min',]
to_remove = which((is.na(Data$Change_in_Max_Load))& (is.na(Data$RPE)))
if (length(to_remove)>0)
{
Data = Data[-to_remove,]
}
if (dim(Data)[1]>0)
{
tryCatch({
X = Data$Time - min(Data$Time)
Y = Data$RPE
#Fitting the cubic polynomial to the RPE scores
model = lm(Y ~ stats::poly(X, 3))
Y_fit = predict(model, data.frame(X))
Data$Fitted_RPE = Y_fit
Full_Data1 = rbind(Full_Data1, Data)
}, error = function(e) {})
}
}
}3.2 Period 2
The same process is followed to obtain the interpolated values of RPE scores for the second work period of the experiment. Similarly, a cubic fit is used in this section.
Full_Data2 = NULL
for (i in 1:length(Subjects))
{
for (j in 1:length(Conditions))
{
Subject = Subjects[i]
Condition = Conditions[j]
Data = Mixed_Data[Mixed_Data$ID == Subject &
unlist(Mixed_Data$Condition) == Condition & Mixed_Data$Period == '2nd_45min',]
to_remove = which((is.na(Data$Change_in_Max_Load))& (is.na(Data$RPE)))
if (length(to_remove)>0)
{
Data = Data[-to_remove,]
}
if (dim(Data)[1]>0)
{
tryCatch({
X = Data$Time - min(Data$Time)
Y = Data$RPE
#Fitting the cubic polynomial to the RPE scores
model = lm(Y ~ stats::poly(X, 3))
Y_fit = predict(model, data.frame(X))
Data$Fitted_RPE = Y_fit
Full_Data2 = rbind(Full_Data2, Data)
}, error = function(e) {})
}
}
}4 Plotting-1st Work Period
In this section, the fitted RPE scores are plotted together with the original reported RPE scores along with the original strength decline ratios. As demonstrated in these plots, the interpolation has provided us with the fitted RPE scores in the time-points where it is not originally reported (in the strength test time points). Then the correlation analysis of RPE scores and the strength decline is calculated which provides the correlation coefficients.
Corr_Matrix = matrix(data = NA, nrow = length(Conditions),
ncol = length(Subjects))
rownames(Corr_Matrix) = Conditions
colnames(Corr_Matrix) = Subjects
Full_Data = NULL
for (i in 1:length(Subjects))
{
cat('##',Subjects[i], "{.tabset .tabset-fade}", '\n')
for (j in 1:length(Conditions))
{
cat('###', as.character(Conditions[j]), "{.tabset .tabset-fade}",'\n')
Subject = Subjects[i]
Condition = as.character(Conditions[j])
Data = Mixed_Data[which(Mixed_Data$ID == Subject &
Mixed_Data$Condition == Condition & Mixed_Data$Period == "1st_45min"),]
to_remove = which((is.na(Data$Change_in_Max_Load))&
(is.na(Data$RPE)))
if(length(to_remove)==0)({
Data = Data
}) else {
Data = Data[-to_remove,]}
if (dim(Data)[1]>0)
{
tryCatch({
X = Data$Time
Y = Data$RPE
model = lm(Y ~ stats::poly(X, 3))
Y_fit = predict(model, data.frame(X))
Y_fit[Y_fit<0] = c(0)
Y_fit[Y_fit>10] = c(10)
Data$Fitted_RPE = Y_fit
Full_Data = rbind(Full_Data, Data)
#Calculate the correlation
Corr_Matrix[j, i] = cor(Data$Change_in_Max_Load,
Data$Fitted_RPE, use = "complete.obs")
gg_Data1 = Data[,c(5,6,7,9)]
colnames(gg_Data1) = c('Time', 'Y', 'Y', 'Y')
gg_Data1 = rbind(gg_Data1[,c(1,2)], gg_Data1[,c(1,3)], gg_Data1[,c(1,4)])
gg_Data1$Variable = c(rep('Strength', dim(Data)[1]),
rep('RPE', 2*dim(Data)[1]))
gg_Data1$Type = c(rep('Reported', 2*dim(Data)[1]),
rep('Fitted', dim(Data)[1]))
gg_Data1 = na.omit(gg_Data1)
gg_Data1$Y = as.numeric(gg_Data1$Y)
p = ggplot(gg_Data1, aes(Time, Y, colour = Type)) + geom_point()
# Use vars() to supply faceting variables:
p = p + facet_wrap(vars(Variable), nrow = 2, scales = "free_y")+
scale_x_continuous(breaks = scales::pretty_breaks(n=5), limits = c(0, 50))
print(p)
}, error = function(e) {})
}
cat('\n \n')
}
cat('\n \n')
}4.1 Sub01
4.1.1 Condition 2.5-15
4.1.2 Condition 2.5-10
4.1.3 Condition 2.5-5
4.1.4 Condition 1.5-15
4.2 Sub02
4.2.1 Condition 2.5-15
4.2.2 Condition 2.5-10
4.2.3 Condition 2.5-5
4.2.4 Condition 1.5-15
4.3 Sub03
4.3.1 Condition 2.5-15
4.3.2 Condition 2.5-10
4.3.3 Condition 2.5-5
4.3.4 Condition 1.5-15
4.4 Sub04
4.4.1 Condition 2.5-15
4.4.2 Condition 2.5-10
4.4.3 Condition 2.5-5
4.4.4 Condition 1.5-15
4.5 Sub05
4.5.1 Condition 2.5-15
4.5.2 Condition 2.5-10
4.5.3 Condition 2.5-5
4.5.4 Condition 1.5-15
4.6 Sub06
4.6.1 Condition 2.5-15
4.6.2 Condition 2.5-10
4.6.3 Condition 2.5-5
4.6.4 Condition 1.5-15
4.7 Sub07
4.7.1 Condition 2.5-15
4.7.2 Condition 2.5-10
4.7.3 Condition 2.5-5
4.7.4 Condition 1.5-15
4.8 Sub08
4.8.1 Condition 2.5-15
4.8.2 Condition 2.5-10
4.8.3 Condition 2.5-5
4.8.4 Condition 1.5-15
4.9 Sub09
4.9.1 Condition 2.5-15
4.9.2 Condition 2.5-10
4.9.3 Condition 2.5-5
4.9.4 Condition 1.5-15
4.10 Sub10
4.10.1 Condition 2.5-15
4.10.2 Condition 2.5-10
4.10.3 Condition 2.5-5
4.10.4 Condition 1.5-15
4.11 Sub13
4.11.1 Condition 2.5-15
4.11.2 Condition 2.5-10
4.11.3 Condition 2.5-5
4.11.4 Condition 1.5-15
4.12 Sub15
4.12.1 Condition 2.5-15
4.12.2 Condition 2.5-10
4.12.3 Condition 2.5-5
4.12.4 Condition 1.5-15
4.13 Sub16
4.13.1 Condition 2.5-15
4.13.2 Condition 2.5-10
4.13.3 Condition 2.5-5
4.13.4 Condition 1.5-15
4.14 Sub17
4.14.1 Condition 2.5-15
4.14.2 Condition 2.5-10
4.14.3 Condition 2.5-5
4.14.4 Condition 1.5-15
5 Plotting-2nd Work Period
The same approach as the previous section is followed to plot the fitted RPE scores for the second work period of the experiment and the correlation coefficients are also calculated in a similar manner.
Corr_Matrix2 = matrix(data = NA, nrow = length(Conditions),
ncol = length(Subjects))
rownames(Corr_Matrix2) = Conditions
colnames(Corr_Matrix2) = Subjects
Full_Data2 = NULL
for (i in 1:length(Subjects))
{
cat('##',Subjects[i], "{.tabset .tabset-fade}", '\n')
for (j in 1:length(Conditions))
{
cat('###', as.character(Conditions[j]), "{.tabset .tabset-fade}",'\n')
Subject = Subjects[i]
Condition = as.character(Conditions[j])
Data2 = Mixed_Data[which(Mixed_Data$ID == Subject &
Mixed_Data$Condition == Condition & Mixed_Data$Period == "2nd_45min"),]
to_remove = which((is.na(Data2$Change_in_Max_Load))&
(is.na(Data2$RPE)))
if(length(to_remove)==0)({
Data2 = Data2
}) else {
Data2 = Data2[-to_remove,]}
if (dim(Data2)[1]>0)
{
tryCatch({
X = Data2$Time
Y = Data2$RPE
model = lm(Y ~ stats::poly(X, 3))
Y_fit = predict(model, data.frame(X))
Y_fit[Y_fit<0] = c(0)
Y_fit[Y_fit>10] = c(10)
Data2$Fitted_RPE = Y_fit
Full_Data2 = rbind(Full_Data2, Data2)
#Calculate the correlation
Corr_Matrix2[j, i] = cor(Data2$Change_in_Max_Load,
Data2$Fitted_RPE, use = "complete.obs")
gg_Data1 = Data2[,c(5,6,7,9)]
colnames(gg_Data1) = c('Time', 'Y', 'Y', 'Y')
gg_Data1 = rbind(gg_Data1[,c(1,2)], gg_Data1[,c(1,3)], gg_Data1[,c(1,4)])
gg_Data1$Variable = c(rep('Strength', dim(Data2)[1]),
rep('RPE', 2*dim(Data2)[1]))
gg_Data1$Type = c(rep('Reported', 2*dim(Data2)[1]),
rep('Fitted', dim(Data2)[1]))
gg_Data1 = na.omit(gg_Data1)
gg_Data1$Y = as.numeric(gg_Data1$Y)
p = ggplot(gg_Data1, aes(Time, Y, colour = Type)) + geom_point()
# Use vars() to supply faceting variables:
p = p + facet_wrap(vars(Variable), nrow = 2, scales = "free_y")+
scale_x_continuous(breaks = scales::pretty_breaks(n=5), limits = c(min(gg_Data1$Time), max(gg_Data1$Time)))
print(p)
}, error = function(e) {})
}
cat('\n \n')
}
cat('\n \n')
}5.1 Sub01
5.1.1 Condition 2.5-15
5.1.2 Condition 2.5-10
5.1.3 Condition 2.5-5
5.1.4 Condition 1.5-15
5.2 Sub02
5.2.1 Condition 2.5-15
5.2.2 Condition 2.5-10
5.2.3 Condition 2.5-5
5.2.4 Condition 1.5-15
5.3 Sub03
5.3.1 Condition 2.5-15
5.3.2 Condition 2.5-10
5.3.3 Condition 2.5-5
5.3.4 Condition 1.5-15
5.4 Sub04
5.4.1 Condition 2.5-15
5.4.2 Condition 2.5-10
5.4.3 Condition 2.5-5
5.4.4 Condition 1.5-15
5.5 Sub05
5.5.1 Condition 2.5-15
5.5.2 Condition 2.5-10
5.5.3 Condition 2.5-5
5.5.4 Condition 1.5-15
5.6 Sub06
5.6.1 Condition 2.5-15
5.6.2 Condition 2.5-10
5.6.3 Condition 2.5-5
5.6.4 Condition 1.5-15
5.7 Sub07
5.7.1 Condition 2.5-15
5.7.2 Condition 2.5-10
5.7.3 Condition 2.5-5
5.7.4 Condition 1.5-15
5.8 Sub08
5.8.1 Condition 2.5-15
5.8.2 Condition 2.5-10
5.8.3 Condition 2.5-5
5.8.4 Condition 1.5-15
5.9 Sub09
5.9.1 Condition 2.5-15
5.9.2 Condition 2.5-10
5.9.3 Condition 2.5-5
5.9.4 Condition 1.5-15
5.10 Sub10
5.10.1 Condition 2.5-15
5.10.2 Condition 2.5-10
5.10.3 Condition 2.5-5
5.10.4 Condition 1.5-15
5.11 Sub13
5.11.1 Condition 2.5-15
5.11.2 Condition 2.5-10
5.11.3 Condition 2.5-5
5.11.4 Condition 1.5-15
5.12 Sub15
5.12.1 Condition 2.5-15
5.12.2 Condition 2.5-10
5.12.3 Condition 2.5-5
5.12.4 Condition 1.5-15
5.13 Sub16
5.13.1 Condition 2.5-15
5.13.2 Condition 2.5-10
5.13.3 Condition 2.5-5
5.13.4 Condition 1.5-15
5.14 Sub17
5.14.1 Condition 2.5-15
5.14.2 Condition 2.5-10
5.14.3 Condition 2.5-5
5.14.4 Condition 1.5-15
6 Correlation Results_1st Work Period
In this section, the values of the Pearson’s correlation coefficients are reported in a table for the first work period. In the first 45 minute work period, the average correlation between RPE and strength change, across all conditions, is -0.62. When considering each individual’s correlation, in many cases the correlations exceed -0.8 for the first period. This indicates a somewhat strong negative linear relationship between RPE and strength decline. The negative sign is expected since as a participant gets fatigued, one expects the RPE value to increase and strength to decrease/decline.
# Corr_Matrix = data.frame(Corr_Matrix)
Corr_Matrix %>%
kable(format = "html", col.names = colnames(Corr_Matrix)) %>%
kable_styling() %>%
kableExtra::scroll_box(width = "100%", height = "300px")| Sub01 | Sub02 | Sub03 | Sub04 | Sub05 | Sub06 | Sub07 | Sub08 | Sub09 | Sub10 | Sub13 | Sub15 | Sub16 | Sub17 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Condition 2.5-15 | -0.6803455 | -0.0470419 | -0.9917536 | -0.7335948 | 0.5569103 | -0.8016492 | -0.5306839 | -0.7817394 | -0.9212997 | -0.8048402 | -0.9542397 | -0.8017957 | -0.9451933 | 0.6200871 |
| Condition 2.5-10 | -0.9207853 | -0.4409184 | -0.5007067 | -0.8156772 | -0.3920030 | -0.5754893 | -0.0195457 | -0.8657196 | -0.6834956 | -0.8720612 | -0.8386104 | -0.9132669 | -0.5565792 | -0.6428274 |
| Condition 2.5-5 | -0.3633217 | -0.8366559 | -0.9129173 | -0.8592127 | -0.4760247 | -0.2293649 | -0.7431481 | -0.4529402 | -0.3027686 | -0.0339077 | -0.9640506 | -0.7636875 | -0.7820622 | -0.8152096 |
| Condition 1.5-15 | -0.3350294 | -0.9204215 | -0.9061243 | -0.7140332 | -0.6381683 | -0.9446565 | -0.9229072 | -0.7302866 | -0.8297332 | 0.2679715 | -0.9989992 | -0.5138187 | -0.8148924 | -0.7354459 |
Average_Corr_1st = mean(na.omit(Corr_Matrix))
Average_Corr_1st[1] -0.6272622
7 Correlation Results_2nd Work Period
Similar to the previous section, the values of the Pearson’s correlation coefficients are reported in a table for the second work period. In the second period, the average correlation between RPE and strength change, across all conditions, is -0.43, which shows the negative relationship between RPE and strength decline, however, it is less stronger than the first work period.
# Corr_Matrix = data.frame(Corr_Matrix)
Corr_Matrix2 %>%
kable(format = "html", col.names = colnames(Corr_Matrix2)) %>%
kable_styling() %>%
kableExtra::scroll_box(width = "100%", height = "300px")| Sub01 | Sub02 | Sub03 | Sub04 | Sub05 | Sub06 | Sub07 | Sub08 | Sub09 | Sub10 | Sub13 | Sub15 | Sub16 | Sub17 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Condition 2.5-15 | 0.3582046 | -0.5990801 | -0.9523749 | -0.6112975 | -0.2032299 | -0.9950105 | -0.2239527 | -0.7561728 | -0.9977609 | -0.5933673 | -0.8860376 | -0.6849213 | -0.8483693 | -0.6888065 |
| Condition 2.5-10 | -0.3149191 | -0.2125776 | -0.7896680 | -0.8680324 | -0.8056649 | -0.9781610 | 0.8236219 | -0.5626364 | -0.2578060 | -0.8891650 | -0.5598962 | -0.6438110 | -0.7273062 | 0.4103357 |
| Condition 2.5-5 | 0.6249696 | -0.3193794 | -0.7233387 | -0.0087296 | -0.1922238 | -0.7439048 | -0.9143338 | -0.1691954 | 0.9252212 | -0.5140405 | -0.7553018 | -0.9679565 | 0.4597668 | 0.2086452 |
| Condition 1.5-15 | -0.0186620 | 0.0337413 | -0.6188169 | 0.1016194 | 0.7543394 | -0.9317877 | 0.5833613 | NA | -0.1777579 | 0.2265232 | -0.7941124 | -0.9140903 | -0.4027939 | -0.7283759 |
Average_Corr_2nd = mean(na.omit(Corr_Matrix2))
Average_Corr_2nd[1] -0.4320872
8 Box-plots_1st Work Period
In this section the box plots of correlation coefficients of RPE scores and strength decline are plotted for each condition type (combination of work loads and task paces) for the first work period.The conditions with higher paces (Condition 2.5-15, and 1.5-15) indicate stronger correlation between the RPE scores and Strength decline.
gg_Data2 = data.frame(as.vector(Corr_Matrix), rep(unlist(Conditions), length(Subjects)))
colnames(gg_Data2) = c('Correlation', 'Condition')
p = ggplot(gg_Data2, aes(x= Condition, y= Correlation)) +
geom_boxplot()
print(p)
9 Box-plots_2nd Work Period
Similar to the previous section, the box plots of correlation coefficients of RPE scores and strength decline are plotted for each condition type (combination of work loads and task paces) for the second work period. The condition with highest pace and load (condition 2.5-15), has the strongest correlation coefficients among other conditions.
gg_Data3 = data.frame(as.vector(Corr_Matrix2), rep(unlist(Conditions), length(Subjects)))
colnames(gg_Data3) = c('Correlation', 'Condition')
p = ggplot(gg_Data3, aes(x= Condition, y= Correlation)) +
geom_boxplot()
print(p)