---
title: "Statistical Analysis for Biomechanical Data"
subtitle: "Comprehensive Statistical Methods for Gait Research"
author: "Biomechanics Research Team"
date: "`r Sys.Date()`"
output:
  html_document:
    theme: flatly
    toc: true
    toc_float: true
    code_folding: show
    fig_width: 10
    fig_height: 6
    number_sections: true
params:
  dataset: "research_data.parquet"
  control_dataset: "control_data.parquet"
  alpha_level: 0.05
  effect_size_threshold: 0.5
  correction_method: "bonferroni"
  show_advanced: TRUE
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width = 10,
  fig.height = 6,
  comment = "#>",
  collapse = TRUE
)
```

# Introduction to Statistical Analysis in Biomechanics

Welcome to advanced statistical analysis for biomechanical research! This tutorial covers the essential statistical methods used in gait analysis, from basic comparisons to sophisticated mixed-effects modeling.

## Learning Objectives

By the end of this tutorial, you'll master:

- ✅ **Descriptive Statistics:** Means, variability, and distribution assessment
- ✅ **Group Comparisons:** t-tests, ANOVA, and non-parametric alternatives
- ✅ **Effect Size Analysis:** Clinical significance beyond p-values
- ✅ **Multiple Comparisons:** Controlling Type I error in complex studies
- ✅ **Mixed-Effects Models:** Accounting for subject-level variability
- ✅ **Result Interpretation:** Translating statistics to clinical meaning

## Prerequisites

- **R Basics Tutorial:** Complete the R Basics for Biomechanics tutorial
- **Statistics Knowledge:** Basic understanding of hypothesis testing
- **Research Context:** Familiarity with biomechanical research questions

---

# Setup and Data Preparation

```{r packages}
# Load required packages
library(LocomotionData)
library(ggplot2)
library(plotly)
library(dplyr)
library(tidyr)
library(broom)
library(car)         # For ANOVA and statistical tests
library(lme4)        # For mixed-effects models
library(lmerTest)    # For p-values in mixed-effects models
library(emmeans)     # For post-hoc tests and marginal means
library(effectsize)  # For effect size calculations
library(report)      # For automated statistical reporting
library(performance) # For model diagnostics
library(see)         # For visualization of statistical models
library(parameters)  # For model parameter extraction

# Set theme for consistent plotting
theme_set(theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    strip.text = element_text(face = "bold")
  ))
```

```{r create-sample-data}
# Create comprehensive sample dataset for this tutorial
set.seed(42)

# Simulation parameters
n_subjects_per_group <- 15
n_cycles <- 8
phases <- 150
tasks <- c("normal_walk", "fast_walk", "slow_walk")
groups <- c("Control", "Patient")

# Create realistic biomechanical data
sample_data <- data.frame()

for (group in groups) {
  for (subj_num in 1:n_subjects_per_group) {
    subj_id <- paste0(group, "_S", sprintf("%02d", subj_num))
    
    # Add group-specific baseline differences
    group_effect <- ifelse(group == "Patient", 5, 0)  # Patients have 5° more knee flexion
    
    for (task in tasks) {
      # Task-specific effects
      task_effect <- switch(task,
        "slow_walk" = -3,
        "normal_walk" = 0,
        "fast_walk" = 4
      )
      
      for (cycle in 1:n_cycles) {
        phase_vec <- seq(0, 100, length.out = phases)
        
        # Generate realistic knee flexion pattern
        base_pattern <- 15 + 45 * sin(2 * pi * phase_vec/100 + pi/4)
        
        # Add individual subject variation
        subject_variation <- rnorm(1, 0, 3)
        
        # Add noise
        noise <- rnorm(phases, 0, 2)
        
        # Combine all effects
        knee_angle <- base_pattern + group_effect + task_effect + subject_variation + noise
        
        # Create cycle data
        cycle_data <- data.frame(
          subject = subj_id,
          group = group,
          task = task,
          cycle = cycle,
          phase = phase_vec,
          knee_flexion_angle_contra_rad = knee_angle * pi/180,
          age = ifelse(group == "Control", rnorm(1, 30, 5), rnorm(1, 45, 8)),
          sex = sample(c("M", "F"), 1, prob = c(0.6, 0.4))
        )
        
        sample_data <- rbind(sample_data, cycle_data)
      }
    }
  }
}

# Save and load through LocomotionData system
temp_file <- tempfile(fileext = ".csv")
write.csv(sample_data, temp_file, row.names = FALSE)

loco_data <- loadLocomotionData(temp_file, 
                               subject_col = "subject",
                               task_col = "task",
                               phase_col = "phase")

cat("✅ Sample dataset created with realistic group differences\n")
cat("📊 Dataset: ", nrow(sample_data), "observations\n")
cat("👥 Groups: ", paste(groups, collapse = ", "), "\n")
cat("🚶 Tasks: ", paste(tasks, collapse = ", "), "\n")
```

---

# Descriptive Statistics and Data Exploration

Before any hypothesis testing, we must thoroughly understand our data distribution and characteristics.

## Summary Statistics by Group

```{r descriptive-stats}
cat("📊 DESCRIPTIVE STATISTICS ANALYSIS\n")
cat("=====================================\n\n")

# Calculate ROM for each subject-task combination
rom_data <- data.frame()
subjects <- getSubjects(loco_data)

for (subj in subjects) {
  # Extract group from subject ID
  group <- ifelse(grepl("Control", subj), "Control", "Patient")
  
  for (task in tasks) {
    rom_result <- calculateROM(loco_data, subj, task, "knee_flexion_angle_contra_rad", by_cycle = FALSE)
    
    if (length(rom_result) > 0) {
      rom_row <- data.frame(
        Subject = subj,
        Group = group,
        Task = task,
        ROM_degrees = rom_result[[1]] * 180/pi
      )
      rom_data <- rbind(rom_data, rom_row)
    }
  }
}

# Calculate descriptive statistics
descriptive_stats <- rom_data %>%
  group_by(Group, Task) %>%
  summarise(
    N = n(),
    Mean = round(mean(ROM_degrees, na.rm = TRUE), 2),
    SD = round(sd(ROM_degrees, na.rm = TRUE), 2),
    SE = round(sd(ROM_degrees, na.rm = TRUE) / sqrt(n()), 2),
    Min = round(min(ROM_degrees, na.rm = TRUE), 2),
    Max = round(max(ROM_degrees, na.rm = TRUE), 2),
    Q25 = round(quantile(ROM_degrees, 0.25, na.rm = TRUE), 2),
    Median = round(median(ROM_degrees, na.rm = TRUE), 2),
    Q75 = round(quantile(ROM_degrees, 0.75, na.rm = TRUE), 2),
    .groups = 'drop'
  )

# Display descriptive statistics
knitr::kable(descriptive_stats, 
             caption = "Descriptive Statistics: Knee ROM by Group and Task",
             align = c("l", "l", rep("c", 8)))
```

## Distribution Assessment

```{r distribution-plots, fig.height=8}
# Create comprehensive distribution plots
p1 <- ggplot(rom_data, aes(x = ROM_degrees, fill = Group)) +
  geom_histogram(alpha = 0.7, bins = 15, position = "identity") +
  facet_wrap(~ Task, ncol = 3) +
  labs(
    title = "Distribution of Knee ROM by Group and Task",
    x = "ROM (degrees)",
    y = "Frequency",
    fill = "Group"
  ) +
  scale_fill_manual(values = c("Control" = "#2E86AB", "Patient" = "#A23B72"))

p2 <- ggplot(rom_data, aes(x = Group, y = ROM_degrees, fill = Group)) +
  geom_boxplot(alpha = 0.7, outlier.alpha = 0.6) +
  geom_jitter(width = 0.2, alpha = 0.4) +
  facet_wrap(~ Task, ncol = 3) +
  labs(
    title = "ROM Distribution: Boxplots by Group and Task",
    x = "Group",
    y = "ROM (degrees)",
    fill = "Group"
  ) +
  scale_fill_manual(values = c("Control" = "#2E86AB", "Patient" = "#A23B72")) +
  theme(legend.position = "none")

# Combine plots
print(p1)
print(p2)
```

## Normality Testing

```{r normality-tests}
cat("🔍 NORMALITY ASSESSMENT\n")
cat("======================\n\n")

# Test normality for each group-task combination
normality_results <- rom_data %>%
  group_by(Group, Task) %>%
  summarise(
    N = n(),
    Shapiro_W = round(shapiro.test(ROM_degrees)$statistic, 4),
    Shapiro_p = round(shapiro.test(ROM_degrees)$p.value, 4),
    Normal = ifelse(shapiro.test(ROM_degrees)$p.value > 0.05, "Yes", "No"),
    .groups = 'drop'
  )

knitr::kable(normality_results,
             caption = "Normality Tests (Shapiro-Wilk) by Group and Task",
             align = c("l", "l", "c", "c", "c", "c"))

# Overall assessment
normal_count <- sum(normality_results$Normal == "Yes")
total_tests <- nrow(normality_results)

cat("\n📈 Normality Summary:\n")
cat("   Normal distributions:", normal_count, "out of", total_tests, "conditions\n")
cat("   Percentage normal:", round(normal_count/total_tests * 100, 1), "%\n\n")

if (normal_count/total_tests >= 0.8) {
  cat("✅ Recommendation: Parametric tests appropriate\n")
} else {
  cat("⚠️  Recommendation: Consider non-parametric tests\n")
}
```

---

# Two-Group Comparisons

Let's start with basic two-group comparisons between Control and Patient groups.

## Independent t-tests

```{r t-tests}
cat("📊 INDEPENDENT T-TESTS\n")
cat("=====================\n\n")

# Perform t-tests for each task
t_test_results <- data.frame()

for (task in tasks) {
  task_data <- rom_data[rom_data$Task == task, ]
  control_data <- task_data$ROM_degrees[task_data$Group == "Control"]
  patient_data <- task_data$ROM_degrees[task_data$Group == "Patient"]
  
  # Perform t-test
  t_result <- t.test(patient_data, control_data, var.equal = FALSE)
  
  # Calculate effect size (Cohen's d)
  cohens_d <- effectsize::cohens_d(patient_data, control_data)
  
  # Store results
  result_row <- data.frame(
    Task = task,
    Control_Mean = round(mean(control_data), 2),
    Patient_Mean = round(mean(patient_data), 2),
    Mean_Diff = round(mean(patient_data) - mean(control_data), 2),
    t_statistic = round(t_result$statistic, 3),
    df = round(t_result$parameter, 1),
    p_value = round(t_result$p.value, 4),
    CI_Lower = round(t_result$conf.int[1], 2),
    CI_Upper = round(t_result$conf.int[2], 2),
    Cohens_d = round(cohens_d$Cohens_d, 3),
    Effect_Size = cohens_d$magnitude
  )
  
  t_test_results <- rbind(t_test_results, result_row)
}

# Display results
knitr::kable(t_test_results,
             caption = "Independent t-test Results: Patient vs Control",
             align = c("l", rep("c", 10)))
```

## Multiple Comparisons Correction

```{r multiple-comparisons}
cat("🔄 MULTIPLE COMPARISONS CORRECTION\n")
cat("=================================\n\n")

# Apply multiple comparisons correction
corrected_results <- t_test_results %>%
  mutate(
    p_raw = p_value,
    p_bonferroni = p.adjust(p_value, method = "bonferroni"),
    p_fdr = p.adjust(p_value, method = "fdr"),
    Sig_Raw = ifelse(p_raw < params$alpha_level, "Yes", "No"),
    Sig_Bonferroni = ifelse(p_bonferroni < params$alpha_level, "Yes", "No"),
    Sig_FDR = ifelse(p_fdr < params$alpha_level, "Yes", "No")
  ) %>%
  select(Task, Mean_Diff, Cohens_d, p_raw, p_bonferroni, p_fdr, 
         Sig_Raw, Sig_Bonferroni, Sig_FDR)

knitr::kable(corrected_results,
             caption = "Multiple Comparisons Corrections",
             align = c("l", rep("c", 8))) %>%
  kableExtra::kable_styling() %>%
  kableExtra::column_spec(7:9, background = ifelse(
    corrected_results[, 7:9] == "Yes", "lightgreen", "white"))

cat("\n📝 Correction Methods:\n")
cat("   - Raw p-values: No correction applied\n")
cat("   - Bonferroni: Conservative family-wise error rate control\n")
cat("   - FDR: False discovery rate control (less conservative)\n")
```

## Effect Size Interpretation

```{r effect-size-viz, fig.height=6}
# Visualize effect sizes
effect_plot_data <- t_test_results %>%
  mutate(
    Significance = ifelse(p_value < params$alpha_level, "Significant", "Not Significant"),
    Effect_Category = case_when(
      abs(Cohens_d) < 0.2 ~ "Negligible",
      abs(Cohens_d) < 0.5 ~ "Small", 
      abs(Cohens_d) < 0.8 ~ "Medium",
      TRUE ~ "Large"
    )
  )

p_effect <- ggplot(effect_plot_data, aes(x = Task, y = Cohens_d, fill = Effect_Category)) +
  geom_col(alpha = 0.8) +
  geom_hline(yintercept = c(-0.8, -0.5, -0.2, 0.2, 0.5, 0.8), 
             linetype = "dashed", alpha = 0.5) +
  geom_text(aes(label = paste0("d = ", Cohens_d)), 
            vjust = ifelse(effect_plot_data$Cohens_d > 0, -0.5, 1.5),
            fontface = "bold") +
  labs(
    title = "Effect Sizes: Patient vs Control Differences",
    subtitle = "Cohen's d with interpretation thresholds",
    x = "Walking Task",
    y = "Cohen's d (Effect Size)",
    fill = "Effect Magnitude"
  ) +
  scale_fill_manual(values = c("Negligible" = "gray80", "Small" = "lightblue", 
                              "Medium" = "orange", "Large" = "red")) +
  annotate("text", x = 0.7, y = 0.8, label = "Large (0.8)", size = 3, angle = 90) +
  annotate("text", x = 0.7, y = 0.5, label = "Medium (0.5)", size = 3, angle = 90) +
  annotate("text", x = 0.7, y = 0.2, label = "Small (0.2)", size = 3, angle = 90)

print(p_effect)
```

<div class="alert alert-info">
**💡 Effect Size Interpretation:**
- **Small (0.2):** Noticeable to trained observers
- **Medium (0.5):** Apparent to casual observers  
- **Large (0.8):** Grossly perceptible differences
- **Clinical relevance:** Consider both statistical significance and effect magnitude
</div>

---

# ANOVA for Multiple Group Comparisons

When comparing more than two conditions, ANOVA is the appropriate statistical approach.

## Two-Way ANOVA: Group × Task

```{r two-way-anova}
cat("📊 TWO-WAY ANOVA: GROUP × TASK\n")
cat("=============================\n\n")

# Perform two-way ANOVA
anova_model <- aov(ROM_degrees ~ Group * Task, data = rom_data)
anova_summary <- summary(anova_model)

# Display ANOVA table
print(anova_summary)

# Calculate effect sizes (eta-squared)
eta_squared <- effectsize::eta_squared(anova_model)
print(eta_squared)

# Check assumptions
cat("\n🔍 ASSUMPTION CHECKING:\n")
cat("======================\n")

# Homogeneity of variance (Levene's test)
levene_test <- car::leveneTest(ROM_degrees ~ Group * Task, data = rom_data)
cat("Levene's Test (Homogeneity of Variance):\n")
cat("   F =", round(levene_test$`F value`[1], 3), 
    ", p =", round(levene_test$`Pr(>F)`[1], 4), "\n")

# Normality of residuals
residuals_shapiro <- shapiro.test(residuals(anova_model))
cat("Shapiro-Wilk Test (Residual Normality):\n")
cat("   W =", round(residuals_shapiro$statistic, 4),
    ", p =", round(residuals_shapiro$p.value, 4), "\n\n")

# Assumption summary
cat("📝 Assumption Assessment:\n")
homogeneity_ok <- levene_test$`Pr(>F)`[1] > 0.05
normality_ok <- residuals_shapiro$p.value > 0.05

cat("   Homogeneity of variance:", ifelse(homogeneity_ok, "✅ Met", "❌ Violated"), "\n")
cat("   Normality of residuals:", ifelse(normality_ok, "✅ Met", "❌ Violated"), "\n")
```

## Post-hoc Tests

```{r post-hoc-tests}
cat("🔄 POST-HOC ANALYSIS\n")
cat("===================\n\n")

# Estimated marginal means
emm_group <- emmeans(anova_model, ~ Group)
emm_task <- emmeans(anova_model, ~ Task)
emm_interaction <- emmeans(anova_model, ~ Group * Task)

cat("📊 Estimated Marginal Means:\n\n")

# Group comparisons
cat("GROUP DIFFERENCES:\n")
group_contrasts <- contrast(emm_group, "pairwise", adjust = params$correction_method)
print(group_contrasts)

cat("\nTASK DIFFERENCES:\n")
task_contrasts <- contrast(emm_task, "pairwise", adjust = params$correction_method)
print(task_contrasts)

# Simple effects (if interaction is significant)
if (anova_summary[[1]][["Pr(>F)"]][3] < params$alpha_level) {
  cat("\n🔄 SIMPLE EFFECTS (Interaction detected):\n")
  
  # Group differences within each task
  simple_group <- contrast(emm_interaction, "pairwise", by = "Task", adjust = params$correction_method)
  cat("\nGroup differences within each task:\n")
  print(simple_group)
  
  # Task differences within each group
  simple_task <- contrast(emm_interaction, "pairwise", by = "Group", adjust = params$correction_method)
  cat("\nTask differences within each group:\n")
  print(simple_task)
}
```

## ANOVA Visualization

```{r anova-visualization, fig.height=8}
# Create comprehensive ANOVA visualization
p_means <- ggplot(rom_data, aes(x = Task, y = ROM_degrees, color = Group)) +
  stat_summary(fun = mean, geom = "point", size = 4, 
               position = position_dodge(width = 0.2)) +
  stat_summary(fun = mean, geom = "line", aes(group = Group), size = 1.5,
               position = position_dodge(width = 0.2)) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1,
               position = position_dodge(width = 0.2)) +
  geom_point(alpha = 0.3, position = position_jitterdodge(dodge.width = 0.2, jitter.width = 0.1)) +
  labs(
    title = "Two-Way ANOVA: Group × Task Interaction",
    subtitle = "Mean ± SE with individual data points",
    x = "Walking Task",
    y = "Knee ROM (degrees)",
    color = "Group"
  ) +
  scale_color_manual(values = c("Control" = "#2E86AB", "Patient" = "#A23B72")) +
  theme(legend.position = "bottom")

# Model diagnostics plots
p_residuals <- ggplot(data.frame(fitted = fitted(anova_model), 
                                residuals = residuals(anova_model)),
                     aes(x = fitted, y = residuals)) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(se = FALSE, color = "blue") +
  labs(
    title = "Residuals vs Fitted Values",
    subtitle = "Check for homoscedasticity",
    x = "Fitted Values",
    y = "Residuals"
  )

p_qq <- ggplot(data.frame(residuals = residuals(anova_model)),
               aes(sample = residuals)) +
  stat_qq(alpha = 0.6) +
  stat_qq_line(color = "red") +
  labs(
    title = "Q-Q Plot of Residuals",
    subtitle = "Check for normality",
    x = "Theoretical Quantiles",
    y = "Sample Quantiles"
  )

# Display plots
print(p_means)
print(p_residuals)
print(p_qq)
```

---

# Mixed-Effects Modeling

Mixed-effects models account for the hierarchical structure of biomechanical data (multiple cycles within subjects).

```{r mixed-effects-prep}
cat("🔄 MIXED-EFFECTS MODELING\n")
cat("========================\n\n")

# Prepare data for mixed-effects analysis
# Include cycle-level data for proper modeling
mixed_data <- data.frame()

for (subj in getSubjects(loco_data)) {
  group <- ifelse(grepl("Control", subj), "Control", "Patient")
  
  for (task in tasks) {
    # Get individual cycle ROMs
    cycles_result <- getCycles(loco_data, subj, task, "knee_flexion_angle_contra_rad")
    
    if (!is.null(cycles_result)) {
      data_3d <- cycles_result$data_3d
      
      # Calculate ROM for each cycle
      for (cycle in 1:dim(data_3d)[1]) {
        cycle_angles <- data_3d[cycle, , 1]
        cycle_rom <- (max(cycle_angles) - min(cycle_angles)) * 180/pi
        
        cycle_row <- data.frame(
          Subject = subj,
          Group = group,
          Task = task,
          Cycle = cycle,
          ROM_degrees = cycle_rom
        )
        
        mixed_data <- rbind(mixed_data, cycle_row)
      }
    }
  }
}

cat("📊 Mixed-effects dataset prepared:\n")
cat("   Observations:", nrow(mixed_data), "\n")
cat("   Subjects:", length(unique(mixed_data$Subject)), "\n")
cat("   Average cycles per subject-task:", 
    round(nrow(mixed_data) / (length(unique(mixed_data$Subject)) * length(tasks)), 1), "\n\n")
```

## Fit Mixed-Effects Models

```{r mixed-effects-models}
# Model 1: Random intercept only
model1 <- lmer(ROM_degrees ~ Group * Task + (1 | Subject), data = mixed_data)

# Model 2: Random intercept and slope for Task
model2 <- lmer(ROM_degrees ~ Group * Task + (Task | Subject), data = mixed_data)

# Model comparison
model_comparison <- anova(model1, model2)
print(model_comparison)

# Select best model based on AIC
best_model <- if (AIC(model2) < AIC(model1)) model2 else model1
cat("\n📊 Selected Model:", ifelse(identical(best_model, model2), "Random intercept + slope", "Random intercept only"), "\n\n")

# Model summary
print(summary(best_model))

# Calculate R-squared
r_squared <- performance::r2(best_model)
cat("\n📈 Model Fit (R-squared):\n")
print(r_squared)
```

## Mixed-Effects Results Interpretation

```{r mixed-effects-interpretation}
# Extract fixed effects with confidence intervals
fixed_effects <- broom.mixed::tidy(best_model, conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  mutate(
    estimate = round(estimate, 3),
    std.error = round(std.error, 3),
    statistic = round(statistic, 3),
    p.value = round(p.value, 4),
    conf.low = round(conf.low, 3),
    conf.high = round(conf.high, 3)
  )

cat("📊 FIXED EFFECTS RESULTS:\n")
cat("========================\n\n")

knitr::kable(fixed_effects,
             caption = "Mixed-Effects Model: Fixed Effects",
             align = c("l", rep("c", 7)))

# Extract random effects
random_effects <- as.data.frame(VarCorr(best_model))
cat("\n🔄 RANDOM EFFECTS VARIANCE:\n")
cat("==========================\n")
print(random_effects)

# Calculate intraclass correlation
icc_value <- performance::icc(best_model)
cat("\n📈 Intraclass Correlation Coefficient (ICC):\n")
print(icc_value)
```

---

# Advanced Statistical Reporting

Let's create a comprehensive statistical report that could be included in a research paper.

## Automated Statistical Report

```{r automated-report}
cat("📝 AUTOMATED STATISTICAL REPORT\n")
cat("==============================\n\n")

# Generate automated report using the 'report' package
anova_report <- report::report(anova_model)
mixed_report <- report::report(best_model)

cat("🔍 ANOVA RESULTS:\n")
cat(as.character(anova_report), "\n\n")

cat("🔄 MIXED-EFFECTS RESULTS:\n")  
cat(as.character(mixed_report), "\n\n")

# Effect size summary
effect_summary <- t_test_results %>%
  mutate(
    Clinical_Significance = case_when(
      abs(Cohens_d) >= 0.8 ~ "Large clinical difference",
      abs(Cohens_d) >= 0.5 ~ "Moderate clinical difference", 
      abs(Cohens_d) >= 0.2 ~ "Small clinical difference",
      TRUE ~ "Negligible difference"
    )
  ) %>%
  select(Task, Mean_Diff, Cohens_d, p_value, Clinical_Significance)

cat("📊 CLINICAL SIGNIFICANCE SUMMARY:\n")
cat("=================================\n")
knitr::kable(effect_summary,
             caption = "Clinical Significance of Group Differences",
             col.names = c("Task", "Mean Diff (°)", "Cohen's d", "p-value", "Clinical Interpretation"))
```

## Results Visualization Summary

```{r results-summary, fig.height=10}
# Create comprehensive results summary plot
library(patchwork)

# Plot 1: Group means with confidence intervals
p_summary <- ggplot(rom_data, aes(x = Task, y = ROM_degrees, fill = Group)) +
  stat_summary(fun.data = mean_cl_normal, geom = "col", 
               position = position_dodge(width = 0.8), alpha = 0.7) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar",
               position = position_dodge(width = 0.8), width = 0.2) +
  geom_text(data = effect_plot_data, 
            aes(x = Task, y = max(rom_data$ROM_degrees), 
                label = paste0("d=", Cohens_d, "\np=", round(p_value, 3))),
            position = position_dodge(width = 0.8), vjust = 0, size = 3,
            inherit.aes = FALSE) +
  labs(
    title = "A) Group Differences with Effect Sizes",
    x = "Walking Task",
    y = "Knee ROM (degrees)",
    fill = "Group"
  ) +
  scale_fill_manual(values = c("Control" = "#2E86AB", "Patient" = "#A23B72"))

# Plot 2: Individual data with model predictions
mixed_predictions <- data.frame(
  Task = rep(tasks, 2),
  Group = rep(c("Control", "Patient"), each = 3),
  Predicted = predict(best_model, newdata = expand.grid(
    Task = tasks, 
    Group = c("Control", "Patient"),
    Subject = mixed_data$Subject[1]  # Use first subject as reference
  ), re.form = NA)
)

p_model <- ggplot(mixed_data, aes(x = Task, y = ROM_degrees, color = Group)) +
  geom_jitter(alpha = 0.3, width = 0.2) +
  geom_point(data = mixed_predictions, aes(y = Predicted), 
             size = 4, shape = 17) +
  geom_line(data = mixed_predictions, aes(y = Predicted, group = Group), 
            size = 1.5) +
  labs(
    title = "B) Mixed-Effects Model Predictions",
    subtitle = "Triangles show model predictions",
    x = "Walking Task", 
    y = "Knee ROM (degrees)",
    color = "Group"
  ) +
  scale_color_manual(values = c("Control" = "#2E86AB", "Patient" = "#A23B72"))

# Combine plots
combined_plot <- p_summary / p_model + 
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

print(combined_plot)
```

---

# Clinical Decision Making Framework

Let's translate our statistical results into clinically actionable insights.

## Clinical Significance Assessment

```{r clinical-framework}
cat("🏥 CLINICAL DECISION FRAMEWORK\n")
cat("=============================\n\n")

# Define clinical significance thresholds
clinical_thresholds <- data.frame(
  Measure = c("ROM Difference", "Effect Size", "Statistical Power"),
  Minimal_Clinically_Important = c("5 degrees", "0.5 (medium)", "0.80"),
  Substantial_Clinical_Benefit = c("10 degrees", "0.8 (large)", "0.90"),
  Source = c("Literature consensus", "Cohen 1988", "Statistical convention")
)

knitr::kable(clinical_thresholds,
             caption = "Clinical Significance Thresholds")

# Apply clinical framework to results
clinical_assessment <- t_test_results %>%
  mutate(
    Clinically_Important = abs(Mean_Diff) >= 5,
    Substantial_Benefit = abs(Mean_Diff) >= 10,
    Large_Effect = abs(Cohens_d) >= 0.8,
    Overall_Assessment = case_when(
      Substantial_Benefit & Large_Effect & p_value < 0.05 ~ "Strong clinical evidence",
      Clinically_Important & abs(Cohens_d) >= 0.5 & p_value < 0.05 ~ "Moderate clinical evidence",
      p_value < 0.05 ~ "Statistical significance only",
      abs(Cohens_d) >= 0.5 ~ "Clinical effect, not statistically significant",
      TRUE ~ "No evidence of clinical difference"
    )
  ) %>%
  select(Task, Mean_Diff, Cohens_d, p_value, Overall_Assessment)

cat("\n🎯 CLINICAL EVIDENCE ASSESSMENT:\n")
cat("===============================\n")
knitr::kable(clinical_assessment,
             caption = "Clinical Evidence Assessment by Task") %>%
  kableExtra::kable_styling() %>%
  kableExtra::column_spec(5, background = ifelse(
    grepl("Strong", clinical_assessment$Overall_Assessment), "lightgreen",
    ifelse(grepl("Moderate", clinical_assessment$Overall_Assessment), "lightyellow", "white")
  ))
```

## Recommendations for Clinical Practice

```{r clinical-recommendations, results='asis'}
cat("## Clinical Practice Recommendations\n\n")

# Generate recommendations based on results
strong_evidence <- clinical_assessment$Task[grepl("Strong", clinical_assessment$Overall_Assessment)]
moderate_evidence <- clinical_assessment$Task[grepl("Moderate", clinical_assessment$Overall_Assessment)]

if (length(strong_evidence) > 0) {
  cat("### **Strong Evidence for Clinical Differences**\n")
  for (task in strong_evidence) {
    task_result <- clinical_assessment[clinical_assessment$Task == task, ]
    cat("- **", task, ":** Mean difference =", task_result$Mean_Diff, "°, Cohen's d =", task_result$Cohens_d, "\n")
    cat("  - **Clinical action:** Consider task-specific assessment and intervention\n")
    cat("  - **Monitoring:** Use this task as a sensitive outcome measure\n\n")
  }
}

if (length(moderate_evidence) > 0) {
  cat("### **Moderate Evidence for Clinical Differences**\n")
  for (task in moderate_evidence) {
    task_result <- clinical_assessment[clinical_assessment$Task == task, ]
    cat("- **", task, ":** Mean difference =", task_result$Mean_Diff, "°, Cohen's d =", task_result$Cohens_d, "\n")
    cat("  - **Clinical action:** Consider individual patient assessment\n")
    cat("  - **Research:** Investigate in larger samples\n\n")
  }
}

cat("### **General Recommendations**\n")
cat("- **Assessment protocol:** Include tasks showing strongest group differences\n")
cat("- **Intervention targets:** Focus on biomechanical parameters with large effect sizes\n")
cat("- **Follow-up monitoring:** Use clinically significant change thresholds (≥5°)\n")
cat("- **Research priorities:** Investigate mechanisms underlying observed differences\n\n")

cat("### **Statistical Interpretation Guidelines**\n")
cat("- **p-values:** Indicate probability of observed difference due to chance\n")
cat("- **Effect sizes:** Quantify magnitude of practical significance\n")
cat("- **Confidence intervals:** Provide range of plausible true differences\n")
cat("- **Clinical thresholds:** Consider both statistical and clinical significance\n\n")
```

---

# Summary and Best Practices

## Key Statistical Concepts Mastered

```{r summary, results='asis'}
cat("## ✅ Statistical Concepts Mastered\n\n")

concepts <- data.frame(
  Concept = c(
    "Descriptive Statistics",
    "Normality Testing", 
    "Independent t-tests",
    "Effect Size Analysis",
    "Multiple Comparisons",
    "Two-way ANOVA",
    "Post-hoc Testing",
    "Mixed-Effects Models",
    "Clinical Significance"
  ),
  Application = c(
    "Mean, SD, median, quartiles for data description",
    "Shapiro-Wilk test for parametric test appropriateness",
    "Compare two independent groups",
    "Cohen's d for practical significance assessment",
    "Bonferroni/FDR correction for multiple testing",
    "Analyze group × condition interactions", 
    "Pairwise comparisons after significant ANOVA",
    "Account for subject-level clustering in data",
    "Translate statistics to clinical meaning"
  )
)

knitr::kable(concepts, align = c("l", "l"))

cat("\n## 🎯 Best Practices Checklist\n\n")

checklist <- c(
  "**Always explore data first** - plots and descriptive statistics",
  "**Check statistical assumptions** - normality, homoscedasticity, independence", 
  "**Report effect sizes** - not just p-values",
  "**Correct for multiple comparisons** - when testing multiple hypotheses",
  "**Use appropriate models** - account for data structure (repeated measures)",
  "**Include confidence intervals** - for uncertainty quantification",
  "**Consider clinical significance** - alongside statistical significance",
  "**Visualize results clearly** - plots enhance understanding",
  "**Report methods clearly** - enable reproducibility"
)

for (item in checklist) {
  cat("- ", item, "\n")
}

cat("\n## 🚀 Next Steps in Statistical Learning\n\n")

next_topics <- data.frame(
  Level = c("Intermediate", "Intermediate", "Advanced", "Advanced", "Advanced"),
  Topic = c(
    "Repeated Measures ANOVA",
    "Non-parametric Tests", 
    "Functional Data Analysis",
    "Machine Learning Methods",
    "Bayesian Statistics"
  ),
  Application = c(
    "Multiple timepoints within subjects",
    "When normality assumptions violated",
    "Analyze entire gait waveforms",
    "Pattern recognition and prediction",
    "Probabilistic interpretation of results"
  )
)

knitr::kable(next_topics, align = c("l", "l", "l"))
```

## Statistical Reporting Template

For your research papers, use this template:

```{r reporting-template, results='asis'}
cat("## 📝 Statistical Reporting Template\n\n")

cat("### **Methods Section Template:**\n")
cat('> "Statistical analyses were performed using R (version X.X.X) with the LocomotionData package. ')
cat('Descriptive statistics (mean ± SD) were calculated for all outcome measures. ')
cat('Normality was assessed using Shapiro-Wilk tests. ')
cat('Group differences were analyzed using [independent t-tests/two-way ANOVA] with ')
cat('[Bonferroni/FDR] correction for multiple comparisons. ')
cat('Effect sizes were calculated using Cohen\'s d with thresholds of 0.2, 0.5, and 0.8 ')
cat('for small, medium, and large effects, respectively. ')
cat('Mixed-effects models with random intercepts [and slopes] were used to account for ')
cat('within-subject correlation. Statistical significance was set at α = 0.05."\n\n')

cat("### **Results Section Template:**\n")
cat('> "Knee ROM differed significantly between groups (F(1,XX) = X.XX, p = 0.XXX, η² = 0.XX). ')
cat('Post-hoc analysis revealed that patients demonstrated greater ROM during [task] ')
cat('compared to controls (XX.X ± X.X° vs XX.X ± X.X°, t(XX) = X.XX, p = 0.XXX, d = X.XX). ')
cat('The difference exceeded the minimal clinically important difference of 5°, ')
cat('indicating clinical significance beyond statistical significance."\n\n')
```

---

**Congratulations!** 🎉 You've mastered advanced statistical analysis for biomechanical research. You're now equipped to conduct rigorous statistical analyses that meet the standards of peer-reviewed journals and provide clinically meaningful insights.

**Remember:** Good statistics start with good questions and good data. Always prioritize understanding your research question and data structure before selecting statistical methods.

---

*Statistical analysis completed on `r Sys.Date()`. Generated using the LocomotionData R package for biomechanical research.*