library(readxl)
library(tidyverse)
library(readr)
rel <- read_excel("FPPA_data.xlsx", sheet = "Reliability")Example code
Reliability
Reliability data named “rel”
head(rel)# A tibble: 6 × 3
Original Remeasure Difference
<dbl> <dbl> <dbl>
1 21.4 22.4 1.06
2 12.5 12.1 0.41
3 7.99 7.36 0.63
4 16.2 15.8 0.5
5 10.8 12.0 1.19
6 7.78 8.21 0.43ggplot(rel) +
aes(x = Original, y = Remeasure) +
geom_point(color = "grey2") +
geom_smooth(method = "lm", se = TRUE) +
labs(x = "FPPA (test)", y = "FPPA (re-test)") + theme_minimal()`geom_smooth()` using formula = 'y ~ x'
rel2<- rel %>%subset(select= -c(Difference))
cor<-cor.test(x = rel$Original, y = rel$Remeasure)
library(psych)
Attaching package: 'psych'
The following objects are masked from 'package:ggplot2':
%+%, alphaICC(rel2)Call: ICC(x = rel2)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.99 133 29 30 1.1e-24 0.97 0.99
Single_random_raters ICC2 0.99 138 29 29 3.2e-24 0.97 0.99
Single_fixed_raters ICC3 0.99 138 29 29 3.2e-24 0.97 0.99
Average_raters_absolute ICC1k 0.99 133 29 30 1.1e-24 0.98 1.00
Average_random_raters ICC2k 0.99 138 29 29 3.2e-24 0.98 1.00
Average_fixed_raters ICC3k 0.99 138 29 29 3.2e-24 0.98 1.00
Number of subjects = 30 Number of Judges = 2
See the help file for a discussion of the other 4 McGraw and Wong estimates,rel2<-rel2 %>% mutate(dif = (Remeasure - Original))
TE=sd(rel2$dif)/sqrt(2)
TE[1] 0.5865561# Calculate the degrees of freedom
n <- nrow(rel2)
df <- n - 1
# Calculate the lower and upper bounds of the 95% confidence interval for TE. through chi-squared distribution
lower_CI <- sqrt((df * TE^2) / qchisq(0.975, df ))
upper_CI <- sqrt((df * TE^2) / qchisq(0.025, df ))
# Print the results
cat("95% Confidence Interval for Typical Error (TE): [", lower_CI, ", ", upper_CI, "]\n")95% Confidence Interval for Typical Error (TE): [ 0.4671373 , 0.7885163 ]Frontal plane projection angle data
“Pre” and “Post” test data for each files imported and cleaned
# A tibble: 6 × 8
ID Limb `Attempt 1` `Attempt 2` `Attempt 3` `Attempt 4` `Attempt 5`
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 RTC02 Right 12.9 11.9 10.0 5.34 6.9
2 RTC02 Left 10.3 10.4 7.9 8.16 5.84
3 RTC62 Right 20.3 25.7 17.6 22.2 23.7
4 RTC62 Left 12.2 11.8 11.2 6.81 11.0
5 RTC29 Right 5.53 7.95 10.9 7.57 2.68
6 RTC29 Left 5.59 7.01 5 5.28 5.6
# ℹ 1 more variable: Average <dbl>Individual trial data cleaned.
pre_long<-pivot_longer(Pre, cols = c(3:7), names_to = "Trial", values_to = "FPPA")
post_long<-pivot_longer(Post, cols = c(3:7), names_to = "Trial", values_to = "FPPA")
head(pre_long)# A tibble: 6 × 5
ID Limb Average Trial FPPA
<chr> <chr> <dbl> <chr> <dbl>
1 RTC02 Right 9.42 Attempt 1 12.9
2 RTC02 Right 9.42 Attempt 2 11.9
3 RTC02 Right 9.42 Attempt 3 10.0
4 RTC02 Right 9.42 Attempt 4 5.34
5 RTC02 Right 9.42 Attempt 5 6.9
6 RTC02 Left 8.52 Attempt 1 10.3 meanFFPA <-mean(post_long$FPPA)
sd_FFPA <-sd(post_long$FPPA)
pre_meanFFPA <-mean(pre_long$FPPA)
pre_sd_FFPA <-sd(pre_long$FPPA)
# Threshold for outliers as 2 standard deviations from the mean
threshold <- 2 * sd_FFPA
# Remove rows containing outliers
post_long_clean <- post_long[abs(post_long$FPPA - meanFFPA) < threshold, ]
threshold <- 2 * pre_sd_FFPA
pre_long_clean <- pre_long[abs(pre_long$FPPA - pre_meanFFPA) < threshold, ]Calculate mean of remaining trials and sort by “limb dominance”. Information on Limb dominance in data frame (“ID”)
# A tibble: 6 × 3
ID Leg_dominance Group
<chr> <chr> <chr>
1 RTC29 left Land
2 RTC84 left Land
3 ETC03 left Land
4 ETC15 left Land
5 RTC85 right Land
6 RTC89 right Land pre_left <-pre_long_clean %>%
subset(select = -c(Average)) %>%
filter(Limb == "Left") %>%
group_by(ID) %>% summarise( mean(FPPA))
post_left <-post_long_clean %>%
subset(select = -c(Average)) %>%
filter(Limb == "Left") %>%
group_by(ID) %>% summarise( mean(FPPA))
left<- left_join(pre_left, post_left,
by = c("ID")) %>% rename(preFPPA= "mean(FPPA).x", postFPPA = "mean(FPPA).y")%>% filter(!is.na(postFPPA))
left<- left_join(left, ID,
by = c("ID")) %>% mutate(dif = postFPPA - preFPPA) %>% mutate(limb_dom = ifelse( Leg_dominance == "left", "dom", "non_dom"), limb = "left")
head(left)# A tibble: 6 × 8
ID preFPPA postFPPA Leg_dominance Group dif limb_dom limb
<chr> <dbl> <dbl> <chr> <chr> <dbl> <chr> <chr>
1 ETC02 4.59 2.82 right Sand -1.77 non_dom left
2 ETC03 12.7 7.59 left Land -5.14 dom left
3 ETC05 8.21 5.8 right Sand -2.41 non_dom left
4 ETC06 12.0 7.68 right Land -4.34 non_dom left
5 ETC07 8.41 5.56 right Land -2.85 non_dom left
6 ETC09 9.08 6.16 right Land -2.92 non_dom left Code repeated for right limb and bound.
dom<-rbind(right, left) %>%
filter(limb_dom == "dom")
non_dom<-rbind(right, left) %>%
filter(limb_dom == "non_dom")
head(dom)# A tibble: 6 × 8
ID preFPPA postFPPA Leg_dominance Group dif limb_dom limb
<chr> <dbl> <dbl> <chr> <chr> <dbl> <chr> <chr>
1 ETC02 12.2 6.22 right Sand -5.96 dom right
2 ETC05 7.81 4.47 right Sand -3.34 dom right
3 ETC06 12.4 10.6 right Land -1.82 dom right
4 ETC07 9.94 9.94 right Land -0.00600 dom right
5 ETC09 7.74 5.71 right Land -2.04 dom right
6 ETC14 13.2 6.42 right Sand -6.74 dom rightSummary statistics
tapply(dom$preFPPA, list(dom$Group), mean) Land Sand
13.38808 13.51518 tapply(dom$preFPPA, list(dom$Group), sd) Land Sand
3.640343 4.009116 tapply(non_dom$preFPPA, list(dom$Group), mean) Land Sand
10.06077 10.88787 tapply(non_dom$preFPPA, list(dom$Group), sd) Land Sand
3.883249 4.471404 Inferential statistics
library(lme4)
library(emmeans)
library(performance)
library(ggpubr)
# Fit the linear model
m1 <- lm(dif ~ Group + preFPPA, data = non_dom)
# Obtain estimated marginal means
emm <- emmeans(m1, pairwise ~ Group)
# Function to calculate confidence intervals at various levels
calc_conf_intervals_emm <- function(emm, conf_levels = seq(0.01,0.50, 0.99, by = 0.01)) {
emm_summary <- summary(emm)
conf_intervals <- lapply(conf_levels, function(conf_level) {
conf <- confint(emm, level = conf_level)
data.frame(
conf_level = conf_level,
conf_low = conf$contrasts$lower.CL,
conf_high = conf$contrasts$upper.CL,
estimate = conf$contrasts$estimate
)
})
do.call(rbind, conf_intervals)
}
# Calculate confidence intervals at various levels
conf_intervals <- calc_conf_intervals_emm(emm, conf_levels = seq(0.01, 0.99, by = 0.001))
# Calculate p-values from the confidence intervals
conf_intervals$p_value <- 1 - conf_intervals$conf_level
conf_intervals$conf_level <- conf_intervals$conf_level*100
conf_intervals$suprisal <- -log2(conf_intervals$p_value)
# Filter the data based on your condition (e.g., confidence interval level)
filtered_intervals <- conf_intervals %>%
filter(conf_level %in% c(75, 95)) # Adjust this condition as needed
# Plot the confidence intervals and p-values with horizontal line based on filtered data
a<-ggplot(conf_intervals, aes(x = estimate, y = p_value)) +
geom_ribbon(aes(xmin = conf_low, xmax = conf_high), fill = "lightblue", alpha = 0.4) +
geom_line(aes(x = conf_low), linetype = "solid", color = "darkblue") +
geom_line(aes(x = conf_high), linetype = "solid", color = "darkblue") +
geom_segment(data = filtered_intervals, aes(x = conf_low, y = p_value, xend = conf_high, yend = p_value),
color = "black", linetype = "solid") +
geom_point(data = filtered_intervals, aes(x = conf_low, y = p_value), color = "black") +
geom_point(data = filtered_intervals, aes(x = conf_high, y = p_value), color = "black") +
scale_y_continuous(
breaks = seq(0, 1, by = 0.1),
sec.axis = sec_axis(~ (1 - .) * 100, breaks = seq(0, 100, by = 10), name = "Confidence Interval Level (%)")
) +
geom_vline(xintercept = 0, color = "darkgrey", linetype = "dashed")+
geom_vline(xintercept = 0.79, color = "red", linetype = "dashed")+
geom_vline(xintercept = -0.79, color = "red", linetype = "dashed")+
labs(title = "Compatibility Curves for Between-Group Differences",
x = "Difference in FPPA (°)",
y = "P-value") +
theme_minimal()
b<-ggplot(conf_intervals, aes(x = estimate, y = suprisal)) +
geom_ribbon(aes(xmin = conf_low, xmax = conf_high), fill = "lightblue", alpha = 0.4) +
geom_line(aes(x = conf_low), linetype = "solid", color = "darkblue") +
geom_line(aes(x = conf_high), linetype = "solid", color = "darkblue") +
geom_segment(data = filtered_intervals, aes(x = conf_low, y = suprisal, xend = conf_high, yend = suprisal),
color = "black", linetype = "solid") +
geom_point(data = filtered_intervals, aes(x = conf_low, y = suprisal), color = "black") +
geom_point(data = filtered_intervals, aes(x = conf_high, y = suprisal), color = "black") +
geom_vline(xintercept = 0, color = "darkgrey", linetype = "dashed")+
geom_vline(xintercept = 0.79, color = "red", linetype = "dashed")+
geom_vline(xintercept = -0.79, color = "red", linetype = "dashed")+
labs(
x = "Difference in FPPA (°)",
y = "s-value") +
theme_minimal()+
scale_y_continuous(breaks = seq(0, 10, by = 1))
fig3<-ggarrange(a,b)
fig3
Visualizing the behavior of the residuals
residuals <- residuals(m1)
qqnorm(residuals)
qqline(residuals)
plot(fitted(m1), residuals, ylab = "Residuals", xlab = "Fitted Values")
abline(h = 0, col = "red") # Add a horizontal line at y = 0
hist(residuals, freq = FALSE, main = "Residuals Histogram and Normal Distribution",
xlab = "Residuals", ylab = "Density", xlim = c(min(residuals), max(residuals)))
# Add density plot for normal distribution
curve(dnorm(x, mean = mean(residuals), sd = sd(residuals)),
col = "blue", lwd = 2, add = TRUE, yaxt = "n")
# Normality checks
check_heteroscedasticity(m1)OK: Error variance appears to be homoscedastic (p = 0.705).check_normality(m1)OK: residuals appear as normally distributed (p = 0.656).residuals <- residuals(m1)