library(readxl) #read in excel file
library(readr) #read in .csv file
library(janitor) #tidy data (column names)
library(tidyverse) #wide range of functions including ggplot
library(ggpubr) #aranging ggplots
library(performance) #model checking
library(emmeans) # for estimated marginal means
library(WRS2) # for robust t-test (trimmed mean differences)ISO Prone versus ISO 30
Packages
The following packages were used in this analysis
Force data
Read and clean data:
force <- read_excel("Force.xlsx")
names<-names(force%>%clean_names)
colnames(force) <- names
head(force)# A tibble: 6 × 6
participant_number test l_max_force_n r_max_force_n l_avg_force_n
<dbl> <chr> <dbl> <dbl> <dbl>
1 1 ISO30 392 466. 383.
2 2 ISO30 258. 267. 219.
3 3 ISO30 391. 434 378.
4 4 ISO30 439 479. 419.
5 5 ISO30 391 376 370.
6 6 ISO30 296. 330. 269.
# ℹ 1 more variable: r_avg_force_n <dbl>Visual inspection and normality checking
a<- ggplot(force) +
aes(x =r_max_force_n, fill = test ) +
geom_density(adjust = 1L) +
scale_fill_viridis_d(option = "cividis", direction = 1) +
labs(x = "% MVC", y = "Density") +
ggtitle("Density plot") +
theme_minimal() + facet_wrap(~test) + theme(legend.position = "none")
b<- ggplot(force) +
aes(x =l_max_force_n, fill = test ) +
geom_density(adjust = 1L) +
scale_fill_viridis_d(option = "cividis", direction = 1) +
labs(x = "% MVC", y = "Density") +
theme_minimal() + facet_wrap(~test) + theme(legend.position = "none")
ggarrange(a,b)
model_l <- lm(l_max_force_n ~ test + factor(participant_number), data = force)
par(mfrow = c(2, 2)) # 2x2 layout for plots
plot(model_l)
check_normality(model_l)OK: residuals appear as normally distributed (p = 0.729).check_heteroscedasticity(model_l)OK: Error variance appears to be homoscedastic (p = 0.489).emmeans(model_l, pairwise ~ test)$emmeans
test emmean SE df lower.CL upper.CL
ISO30 351 6.85 42 337 365
PRONE 294 6.85 42 281 308
Results are averaged over the levels of: participant_number
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
ISO30 - PRONE 56.4 9.68 42 5.821 <.0001
Results are averaged over the levels of: participant_number model_r <- lm(r_max_force_n ~ test + factor(participant_number), data = force)
par(mfrow = c(2, 2)) # 2x2 layout for plots
plot(model_r)
check_normality(model_r)OK: residuals appear as normally distributed (p = 0.981).check_heteroscedasticity(model_r)OK: Error variance appears to be homoscedastic (p = 0.087).emmeans(model_r, pairwise ~ test)$emmeans
test emmean SE df lower.CL upper.CL
ISO30 367 6.73 42 353 380
PRONE 314 6.73 42 300 327
Results are averaged over the levels of: participant_number
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
ISO30 - PRONE 52.9 9.52 42 5.561 <.0001
Results are averaged over the levels of: participant_number Between-participant correlation
Filtering left and right leg data and pivoting to wide format:
left<- force %>%
subset(select = c(participant_number, test, l_max_force_n ))%>%
pivot_wider(names_from = test, values_from = l_max_force_n)
right<- force %>%
subset(select = c(participant_number, test, r_max_force_n ))%>%
pivot_wider(names_from = test, values_from = r_max_force_n)Pearson’s correlation
cor.test(x = left$ISO30, y =left$PRONE, method = "pearson")
Pearson's product-moment correlation
data: left$ISO30 and left$PRONE
t = 4.1596, df = 41, p-value = 0.0001588
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.2922375 0.7262785
sample estimates:
cor
0.5447663 model_l <- lm(ISO30 ~ PRONE, data = left)
summary(model_l)
Call:
lm(formula = ISO30 ~ PRONE, data = left)
Residuals:
Min 1Q Median 3Q Max
-149.662 -32.756 -8.021 38.542 139.884
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 140.4950 51.4195 2.732 0.009236 **
PRONE 0.7142 0.1717 4.160 0.000159 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 62.2 on 41 degrees of freedom
Multiple R-squared: 0.2968, Adjusted R-squared: 0.2796
F-statistic: 17.3 on 1 and 41 DF, p-value: 0.0001588par(mfrow = c(2, 2)) # 2x2 layout for plots
plot(model_l)
check_normality(model_l)OK: residuals appear as normally distributed (p = 0.738).check_heteroscedasticity(model_l)OK: Error variance appears to be homoscedastic (p = 0.742).cor.test(x = right$ISO30, y = right$PRONE, method = "pearson")
Pearson's product-moment correlation
data: right$ISO30 and right$PRONE
t = 5.3015, df = 41, p-value = 4.232e-06
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4173321 0.7872859
sample estimates:
cor
0.6377397 model_r <- lm(ISO30 ~ PRONE, data = right)
summary(model_r)
Call:
lm(formula = ISO30 ~ PRONE, data = right)
Residuals:
Min 1Q Median 3Q Max
-106.74 -43.06 -14.63 41.02 143.98
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 103.5845 50.5370 2.050 0.0468 *
PRONE 0.8386 0.1582 5.302 4.23e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 62.4 on 41 degrees of freedom
Multiple R-squared: 0.4067, Adjusted R-squared: 0.3922
F-statistic: 28.11 on 1 and 41 DF, p-value: 4.232e-06par(mfrow = c(2, 2)) # 2x2 layout for plots
plot(model_r)
check_normality(model_r)OK: residuals appear as normally distributed (p = 0.111).check_heteroscedasticity(model_r)OK: Error variance appears to be homoscedastic (p = 0.499).EMG data
Read in data:
raw_mean <- read_csv("sEMG_raw.csv")
head(raw_mean)# A tibble: 6 × 8
participant test gm_pk am_pk st_pk b_fl_pk b_fs_pk med_g_pk
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 P1 ISO30 2.1 2.2 3.8 1.4 2.1 1.8
2 P1 PRONE 2 2.1 3.2 1.5 1.9 1.6
3 P2 ISO30 3.2 2 1.9 3 3.4 1.7
4 P2 PRONE 2.6 1.7 1.4 2.4 2.9 1.5
5 P3 ISO30 1.9 2.6 1.8 2.4 1.8 1.6
6 P3 PRONE 2 1.3 1.5 2.1 2.6 1.6The behaviour of the residuals were visually inspected and met the assumptions of normality and homogeneity of variance”
Call:
lm(formula = gm_pk ~ test + factor(participant), data = raw_mean)
Residuals:
Min 1Q Median 3Q Max
-2.4077 -0.1923 0.0000 0.1923 2.4077
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.9077 0.7794 2.448 0.0307 *
testPRONE 0.2846 0.4166 0.683 0.5075
factor(participant)P10 -0.2000 1.0622 -0.188 0.8538
factor(participant)P11 -0.1000 1.0622 -0.094 0.9265
factor(participant)P12 -0.2500 1.0622 -0.235 0.8179
factor(participant)P13 -0.8000 1.0622 -0.753 0.4659
factor(participant)P2 0.8500 1.0622 0.800 0.4391
factor(participant)P3 -0.1000 1.0622 -0.094 0.9265
factor(participant)P4 0.7000 1.0622 0.659 0.5223
factor(participant)P5 0.2500 1.0622 0.235 0.8179
factor(participant)P6 -0.3000 1.0622 -0.282 0.7824
factor(participant)P7 0.2500 1.0622 0.235 0.8179
factor(participant)P8 1.8000 1.0622 1.695 0.1159
factor(participant)P9 -0.1000 1.0622 -0.094 0.9265
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.062 on 12 degrees of freedom
Multiple R-squared: 0.4435, Adjusted R-squared: -0.1593
F-statistic: 0.7358 on 13 and 12 DF, p-value: 0.705Analysis of Variance Table
Response: gm_pk
Df Sum Sq Mean Sq F value Pr(>F)
test 1 0.5265 0.52654 0.4667 0.5075
factor(participant) 12 10.2646 0.85538 0.7582 0.6804
Residuals 12 13.5385 1.12821 
Warning: Non-normality of residuals detected (p < .001).Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).Individual trimmed mean analysis and summary statistics:
testdata<- raw_mean %>%
subset(select = c(participant, test, gm_pk ))%>%
pivot_wider(names_from = test, values_from = gm_pk)
yuend(testdata$ISO30, testdata$PRONE, tr = 0.2)Call:
yuend(x = testdata$ISO30, y = testdata$PRONE, tr = 0.2)
Test statistic: -0.6542 (df = 8), p-value = 0.53133
Trimmed mean difference: -0.08889
95 percent confidence interval:
-0.4022 0.2244
Explanatory measure of effect size: 0.17 median(testdata$ISO30)[1] 1.9quantile(testdata$ISO30, 0.75)75%
2.2 median(testdata$PRONE)[1] 2quantile(testdata$PRONE, 0.25)25%
1.8