---
title: "Population Comparison Analysis"
subtitle: "Multi-Group Biomechanical Assessment"
author: "Research Team"
date: "`r Sys.Date()`"
output:
  html_document:
    theme: flatly
    toc: true
    toc_float: true
    code_folding: hide
    fig_width: 12
    fig_height: 8
params:
  datasets: ["group1.parquet", "group2.parquet", "group3.parquet"]
  group_names: ["Control", "Treatment A", "Treatment B"]
  primary_task: "normal_walk"
  comparison_type: "between_groups"  # or "within_group"
  demographic_data: "demographics.csv"  # optional
  alpha_level: 0.05
  multiple_comparison_correction: "bonferroni"
  include_effect_sizes: TRUE
  generate_plots: TRUE
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE,
  fig.width = 12,
  fig.height = 8,
  dpi = 300
)

# Load required libraries
library(LocomotionData)
library(ggplot2)
library(plotly)
library(dplyr)
library(tidyr)
library(DT)
library(knitr)
library(kableExtra)
library(viridis)
library(patchwork)
library(car)        # for ANOVA
library(emmeans)    # for post-hoc tests
library(effsize)    # for effect sizes
library(broom)      # for tidy statistical output
library(report)     # for automated reporting

# Set plotting theme
theme_set(theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    plot.subtitle = element_text(size = 12, color = "gray50"),
    strip.text = element_text(face = "bold"),
    legend.position = "bottom"
  ))

# Color palette for groups
group_colors <- viridis::viridis_discrete(n = length(params$group_names))
names(group_colors) <- params$group_names
```

# Study Overview

```{r load-data, include=FALSE}
# Validate inputs
if (length(params$datasets) != length(params$group_names)) {
  stop("Number of datasets must match number of group names")
}

# Load all datasets
locomotion_data <- list()
all_subjects <- list()
all_tasks <- list()
all_features <- list()

for (i in seq_along(params$datasets)) {
  dataset_path <- params$datasets[i]
  group_name <- params$group_names[i]
  
  if (!file.exists(dataset_path)) {
    stop("Dataset file not found: ", dataset_path)
  }
  
  # Load data
  loco_data <- loadLocomotionData(dataset_path)
  locomotion_data[[group_name]] <- loco_data
  
  # Extract metadata
  all_subjects[[group_name]] <- getSubjects(loco_data)
  all_tasks[[group_name]] <- getTasks(loco_data)
  all_features[[group_name]] <- getFeatures(loco_data)
}

# Find common elements across groups
common_tasks <- Reduce(intersect, all_tasks)
common_features <- Reduce(intersect, all_features)

# Set primary analysis parameters
analysis_task <- if (params$primary_task %in% common_tasks) {
  params$primary_task
} else {
  common_tasks[1]
}

# Select key biomechanical features
key_features <- intersect(c(
  "knee_flexion_angle_ipsi_rad", "knee_flexion_angle_contra_rad",
  "hip_flexion_angle_ipsi_rad", "hip_flexion_angle_contra_rad",
  "ankle_flexion_angle_ipsi_rad", "ankle_flexion_angle_contra_rad"
), common_features)

# Load demographic data if available
demographics <- NULL
if (!is.null(params$demographic_data) && 
    params$demographic_data != "" && 
    file.exists(params$demographic_data)) {
  demographics <- read.csv(params$demographic_data)
}
```

**Study Type:** `r params$comparison_type`  
**Analysis Date:** `r Sys.Date()`  
**Primary Task:** `r analysis_task`  
**Groups Compared:** `r paste(params$group_names, collapse = ", ")`

## Group Characteristics

```{r group-summary}
# Create group summary table
group_summary <- data.frame(
  Group = params$group_names,
  N_Participants = sapply(all_subjects, length),
  N_Tasks = sapply(all_tasks, length),
  N_Features = sapply(all_features, length),
  Dataset = basename(params$datasets)
)

kable(group_summary,
      col.names = c("Group", "Participants", "Tasks", "Features", "Dataset"),
      align = c("l", rep("c", 4))) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```

```{r demographic-summary, eval=!is.null(demographics)}
if (!is.null(demographics)) {
  cat("### Demographic Characteristics\n\n")
  
  # Merge demographic data with groups
  demo_summary <- demographics %>%
    group_by(Group) %>%
    summarise(
      N = n(),
      Age_Mean = round(mean(Age, na.rm = TRUE), 1),
      Age_SD = round(sd(Age, na.rm = TRUE), 1),
      Female_N = sum(Sex == "F", na.rm = TRUE),
      Female_Percent = round(sum(Sex == "F", na.rm = TRUE) / n() * 100, 1),
      .groups = 'drop'
    )
  
  kable(demo_summary,
        col.names = c("Group", "N", "Age Mean", "Age SD", "Female N", "Female %"),
        align = c("l", rep("c", 5))) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
}
```

**Common Elements Across Groups:**
- **Tasks:** `r paste(common_tasks, collapse = ", ")`
- **Features:** `r length(common_features)` biomechanical variables
- **Key Features for Analysis:** `r length(key_features)` selected features

---

# Statistical Analysis Framework

```{r analysis-setup}
# Prepare data for statistical analysis
analysis_data <- data.frame()

for (group_name in params$group_names) {
  loco_data <- locomotion_data[[group_name]]
  subjects <- all_subjects[[group_name]]
  
  for (subject in subjects) {
    # Calculate summary metrics for each subject
    subject_metrics <- list()
    
    for (feature in key_features) {
      # Get mean pattern and calculate ROM
      mean_pattern <- getMeanPatterns(loco_data, subject, analysis_task, feature)[[1]]
      rom_value <- max(mean_pattern) - min(mean_pattern)
      peak_value <- max(abs(mean_pattern))
      
      subject_metrics[[paste0(feature, "_ROM")]] <- rom_value * 180/pi
      subject_metrics[[paste0(feature, "_Peak")]] <- peak_value * 180/pi
    }
    
    # Create row for this subject
    subject_row <- data.frame(
      Subject = subject,
      Group = group_name,
      subject_metrics
    )
    
    analysis_data <- rbind(analysis_data, subject_row)
  }
}

# Reshape data for statistical analysis
rom_data <- analysis_data %>%
  select(Subject, Group, ends_with("_ROM")) %>%
  pivot_longer(cols = ends_with("_ROM"), 
               names_to = "Feature", 
               values_to = "ROM") %>%
  mutate(Feature = gsub("_ROM", "", Feature),
         Feature = gsub("_", " ", Feature))

peak_data <- analysis_data %>%
  select(Subject, Group, ends_with("_Peak")) %>%
  pivot_longer(cols = ends_with("_Peak"), 
               names_to = "Feature", 
               values_to = "Peak") %>%
  mutate(Feature = gsub("_Peak", "", Feature),
         Feature = gsub("_", " ", Feature))
```

## Analysis Plan

**Statistical Tests:**
- **Between-group comparisons:** One-way ANOVA with post-hoc tests
- **Effect sizes:** Cohen's d for pairwise comparisons
- **Multiple comparisons:** `r params$multiple_comparison_correction` correction
- **Alpha level:** `r params$alpha_level`

**Outcome Measures:**
- **Range of Motion (ROM):** Peak-to-peak amplitude for each joint
- **Peak Values:** Maximum absolute angle during gait cycle

---

# Range of Motion Analysis

```{r rom-anova}
# Perform ANOVA for each feature (ROM)
rom_anova_results <- list()
rom_posthoc_results <- list()

features_list <- unique(rom_data$Feature)

for (feature in features_list) {
  feature_data <- rom_data[rom_data$Feature == feature, ]
  
  # One-way ANOVA
  anova_model <- aov(ROM ~ Group, data = feature_data)
  anova_summary <- summary(anova_model)
  
  # Post-hoc tests if significant
  posthoc_results <- NULL
  if (anova_summary[[1]][["Pr(>F)"]][1] < params$alpha_level) {
    posthoc_results <- emmeans(anova_model, "Group") %>%
      pairs(adjust = params$multiple_comparison_correction) %>%
      summary()
  }
  
  rom_anova_results[[feature]] <- anova_summary
  rom_posthoc_results[[feature]] <- posthoc_results
}

# Create ANOVA summary table
anova_summary_df <- data.frame()

for (feature in features_list) {
  anova_result <- rom_anova_results[[feature]]
  f_stat <- anova_result[[1]][["F value"]][1]
  p_value <- anova_result[[1]][["Pr(>F)"]][1]
  
  anova_row <- data.frame(
    Feature = feature,
    F_statistic = round(f_stat, 3),
    p_value = ifelse(p_value < 0.001, "<0.001", round(p_value, 3)),
    Significant = ifelse(p_value < params$alpha_level, "Yes", "No")
  )
  
  anova_summary_df <- rbind(anova_summary_df, anova_row)
}

cat("## ANOVA Results (Range of Motion)\n\n")

kable(anova_summary_df,
      col.names = c("Feature", "F", "p-value", "Significant"),
      align = c("l", "c", "c", "c")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
  column_spec(4, color = ifelse(anova_summary_df$Significant == "Yes", "green", "red"))
```

```{r rom-posthoc}
# Display post-hoc results for significant features
significant_features <- anova_summary_df$Feature[anova_summary_df$Significant == "Yes"]

if (length(significant_features) > 0) {
  cat("\n### Post-hoc Comparisons (Significant Features)\n\n")
  
  for (feature in significant_features) {
    posthoc_result <- rom_posthoc_results[[feature]]
    
    if (!is.null(posthoc_result)) {
      cat("**", feature, "**\n\n")
      
      posthoc_df <- data.frame(
        Comparison = paste(posthoc_result$contrast),
        Estimate = round(posthoc_result$estimate, 2),
        SE = round(posthoc_result$SE, 2),
        p_value = ifelse(posthoc_result$p.value < 0.001, "<0.001", 
                        round(posthoc_result$p.value, 3)),
        Significant = ifelse(posthoc_result$p.value < params$alpha_level, "Yes", "No")
      )
      
      print(kable(posthoc_df,
            col.names = c("Comparison", "Difference (°)", "SE", "p-value", "Significant"),
            align = c("l", rep("c", 4))) %>%
        kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
        column_spec(5, color = ifelse(posthoc_df$Significant == "Yes", "green", "red")))
      
      cat("\n")
    }
  }
}
```

---

# Effect Size Analysis

```{r effect-sizes, eval=params$include_effect_sizes}
if (params$include_effect_sizes && length(params$group_names) == 2) {
  
  cat("## Effect Sizes (Cohen's d)\n\n")
  
  # Calculate effect sizes for two-group comparison
  group1 <- params$group_names[1]
  group2 <- params$group_names[2]
  
  effect_sizes_df <- data.frame()
  
  for (feature in features_list) {
    feature_data <- rom_data[rom_data$Feature == feature, ]
    
    group1_data <- feature_data$ROM[feature_data$Group == group1]
    group2_data <- feature_data$ROM[feature_data$Group == group2]
    
    # Calculate Cohen's d
    cohens_d <- cohen.d(group1_data, group2_data)
    
    effect_row <- data.frame(
      Feature = feature,
      Group1_Mean = round(mean(group1_data), 2),
      Group2_Mean = round(mean(group2_data), 2),
      Cohens_d = round(cohens_d$estimate, 3),
      Effect_Size = cohens_d$magnitude,
      CI_Lower = round(cohens_d$conf.int[1], 3),
      CI_Upper = round(cohens_d$conf.int[2], 3)
    )
    
    effect_sizes_df <- rbind(effect_sizes_df, effect_row)
  }
  
  kable(effect_sizes_df,
        col.names = c("Feature", paste(group1, "Mean (°)"), paste(group2, "Mean (°)"),
                      "Cohen's d", "Effect Size", "95% CI Lower", "95% CI Upper"),
        align = c("l", rep("c", 6))) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
    column_spec(5, background = ifelse(effect_sizes_df$Effect_Size %in% c("large", "medium"), 
                                     "lightblue", "white"))

} else if (params$include_effect_sizes && length(params$group_names) > 2) {
  cat("## Effect Sizes\n\n")
  cat("*Effect size calculation for multiple groups requires pairwise comparisons.*\n")
  cat("*Consider using eta-squared from ANOVA results for overall effect sizes.*\n\n")
}
```

---

# Visualization

```{r visualization-prep, eval=params$generate_plots}
if (params$generate_plots) {
  
  # Prepare group mean patterns for visualization
  group_patterns <- list()
  
  for (group_name in params$group_names) {
    loco_data <- locomotion_data[[group_name]]
    subjects <- all_subjects[[group_name]]
    
    # Calculate group means for each feature
    group_means <- list()
    
    for (feature in key_features) {
      # Get all subject patterns
      subject_patterns <- list()
      
      for (subject in subjects) {
        pattern <- getMeanPatterns(loco_data, subject, analysis_task, feature)[[1]]
        subject_patterns[[subject]] <- pattern
      }
      
      # Calculate group mean
      pattern_matrix <- do.call(rbind, subject_patterns)
      group_mean <- colMeans(pattern_matrix, na.rm = TRUE)
      group_means[[feature]] <- group_mean
    }
    
    group_patterns[[group_name]] <- group_means
  }
  
  phase_vector <- seq(0, 100, length.out = 150)
}
```

## Group Comparison Plots

```{r group-comparison-plots, eval=params$generate_plots, fig.height=12}
if (params$generate_plots && length(key_features) > 0) {
  
  # Create comparison plots for each feature
  comparison_plots <- list()
  
  for (i in seq_along(key_features)) {
    feature <- key_features[i]
    
    # Prepare plot data
    plot_data <- data.frame()
    
    for (group_name in params$group_names) {
      group_data <- data.frame(
        Phase = phase_vector,
        Angle = group_patterns[[group_name]][[feature]] * 180/pi,
        Group = group_name
      )
      plot_data <- rbind(plot_data, group_data)
    }
    
    # Create plot
    p <- ggplot(plot_data, aes(x = Phase, y = Angle, color = Group)) +
      geom_line(size = 1.5, alpha = 0.8) +
      labs(
        title = gsub("_", " ", feature),
        x = "Gait Cycle (%)",
        y = "Angle (degrees)",
        color = "Group"
      ) +
      scale_color_manual(values = group_colors) +
      theme(legend.position = "bottom")
    
    comparison_plots[[i]] <- p
  }
  
  # Arrange plots
  if (length(comparison_plots) >= 6) {
    combined_plot <- (comparison_plots[[1]] | comparison_plots[[2]] | comparison_plots[[3]]) /
                    (comparison_plots[[4]] | comparison_plots[[5]] | comparison_plots[[6]])
  } else if (length(comparison_plots) >= 3) {
    combined_plot <- comparison_plots[[1]] | comparison_plots[[2]] | comparison_plots[[3]]
  } else if (length(comparison_plots) >= 2) {
    combined_plot <- comparison_plots[[1]] | comparison_plots[[2]]
  } else {
    combined_plot <- comparison_plots[[1]]
  }
  
  print(combined_plot)
}
```

## Interactive Group Analysis

```{r interactive-group-plot, eval=params$generate_plots}
if (params$generate_plots && length(key_features) > 0) {
  
  # Create interactive plot for primary feature
  primary_feature <- key_features[1]
  
  # Prepare interactive data
  interactive_data <- data.frame()
  
  for (group_name in params$group_names) {
    group_data <- data.frame(
      Phase = phase_vector,
      Angle = group_patterns[[group_name]][[primary_feature]] * 180/pi,
      Group = group_name
    )
    interactive_data <- rbind(interactive_data, group_data)
  }
  
  # Create interactive plot
  p_interactive <- ggplot(interactive_data, aes(x = Phase, y = Angle, color = Group)) +
    geom_line(size = 1.5, 
              aes(text = paste("Group:", Group, "<br>Phase:", Phase, "%<br>Angle:", round(Angle, 1), "°"))) +
    labs(
      title = paste("Interactive Group Comparison:", gsub("_", " ", primary_feature)),
      subtitle = "Hover for detailed values | Compare group patterns",
      x = "Gait Cycle (%)",
      y = "Angle (degrees)",
      color = "Group"
    ) +
    scale_color_manual(values = group_colors)
  
  ggplotly(p_interactive, tooltip = "text") %>%
    layout(hovermode = "x unified")
}
```

## Range of Motion Box Plots

```{r rom-boxplots, eval=params$generate_plots, fig.height=8}
if (params$generate_plots) {
  
  # Create box plots for ROM comparison
  p_rom <- ggplot(rom_data, aes(x = Group, y = ROM, fill = Group)) +
    geom_boxplot(alpha = 0.7, outlier.alpha = 0.6) +
    geom_jitter(width = 0.2, alpha = 0.4, size = 1) +
    facet_wrap(~ Feature, scales = "free_y", ncol = 3) +
    labs(
      title = "Range of Motion Comparison Across Groups",
      subtitle = paste("Task:", analysis_task),
      x = "Group",
      y = "ROM (degrees)",
      fill = "Group"
    ) +
    scale_fill_manual(values = group_colors) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      legend.position = "none"
    )
  
  print(p_rom)
}
```

---

# Clinical Interpretation

```{r clinical-interpretation, results='asis'}
# Generate clinical interpretation
cat("## Statistical Summary\n\n")

significant_count <- sum(anova_summary_df$Significant == "Yes")
total_features <- nrow(anova_summary_df)

cat("**Overall Results:**\n")
cat("- **Significant differences found:** ", significant_count, " out of ", total_features, " features\n")
cat("- **Proportion of significant findings:** ", round(significant_count/total_features * 100, 1), "%\n\n")

if (significant_count > 0) {
  cat("**Features with significant group differences:**\n")
  significant_features <- anova_summary_df$Feature[anova_summary_df$Significant == "Yes"]
  cat(paste("-", significant_features, collapse = "\n"))
  cat("\n\n")
}

# Effect size interpretation (for two-group studies)
if (params$include_effect_sizes && length(params$group_names) == 2) {
  large_effects <- effect_sizes_df$Feature[effect_sizes_df$Effect_Size == "large"]
  medium_effects <- effect_sizes_df$Feature[effect_sizes_df$Effect_Size == "medium"]
  
  if (length(large_effects) > 0) {
    cat("**Large effect sizes (d > 0.8):**\n")
    cat(paste("-", large_effects, collapse = "\n"))
    cat("\n\n")
  }
  
  if (length(medium_effects) > 0) {
    cat("**Medium effect sizes (0.5 < d < 0.8):**\n")
    cat(paste("-", medium_effects, collapse = "\n"))
    cat("\n\n")
  }
}

cat("## Clinical Significance\n\n")

# Determine overall clinical interpretation
if (significant_count >= total_features * 0.5) {
  cat("**Interpretation:** Strong evidence of biomechanical differences between groups. ")
  cat("Multiple joint parameters show significant variations, suggesting meaningful ")
  cat("physiological or pathological distinctions.\n\n")
} else if (significant_count >= total_features * 0.25) {
  cat("**Interpretation:** Moderate evidence of group differences. Some biomechanical ")
  cat("parameters show significant variations, indicating potential differences that ")
  cat("warrant further investigation.\n\n")
} else {
  cat("**Interpretation:** Limited evidence of biomechanical differences between groups. ")
  cat("Most parameters show similar patterns, suggesting comparable gait characteristics ")
  cat("across the studied populations.\n\n")
}

cat("## Recommendations\n\n")
cat("### Clinical Practice\n")
cat("- **Assessment protocols:** Consider group-specific reference values\n")
cat("- **Intervention planning:** Target features showing significant differences\n")
cat("- **Monitoring:** Focus on parameters with large effect sizes\n\n")

cat("### Future Research\n")
cat("- **Mechanistic studies:** Investigate underlying causes of observed differences\n")
cat("- **Longitudinal follow-up:** Track changes over time within groups\n")
cat("- **Intervention studies:** Test targeted treatments for identified differences\n")
cat("- **Sample size:** Consider larger samples for features showing trends\n\n")
```

---

# Study Quality and Limitations

```{r study-quality, results='asis'}
cat("## Data Quality Assessment\n\n")

# Assess data quality for each group
quality_summary <- data.frame()

for (group_name in params$group_names) {
  loco_data <- locomotion_data[[group_name]]
  subjects <- all_subjects[[group_name]]
  
  # Calculate validation rates
  total_cycles <- 0
  valid_cycles <- 0
  
  for (subject in subjects[1:min(5, length(subjects))]) {  # Sample subset for efficiency
    valid_mask <- validateCycles(loco_data, subject, analysis_task)
    total_cycles <- total_cycles + length(valid_mask)
    valid_cycles <- valid_cycles + sum(valid_mask)
  }
  
  validity_rate <- if (total_cycles > 0) valid_cycles / total_cycles * 100 else 0
  
  quality_row <- data.frame(
    Group = group_name,
    N_Subjects = length(subjects),
    Validity_Rate = round(validity_rate, 1),
    Data_Quality = ifelse(validity_rate >= 80, "Excellent",
                         ifelse(validity_rate >= 60, "Good", "Fair"))
  )
  
  quality_summary <- rbind(quality_summary, quality_row)
}

kable(quality_summary,
      col.names = c("Group", "N Subjects", "Validity Rate (%)", "Data Quality"),
      align = c("l", "c", "c", "c")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

cat("\n## Study Limitations\n\n")

limitations <- c(
  "**Cross-sectional design:** Cannot establish causal relationships",
  "**Sample characteristics:** Consider demographic and clinical heterogeneity",
  "**Multiple comparisons:** Risk of Type I error despite corrections",
  "**Effect size interpretation:** Clinical significance may differ from statistical significance"
)

# Add specific limitations based on analysis
if (min(quality_summary$Validity_Rate) < 70) {
  limitations <- c(limitations, "**Data quality:** Some groups show reduced data validity")
}

if (length(common_features) < 10) {
  limitations <- c(limitations, "**Feature coverage:** Limited biomechanical parameter set")
}

cat(paste("-", limitations, collapse = "\n"))
cat("\n\n")

cat("## Statistical Considerations\n\n")
cat("- **Power analysis:** Post-hoc power calculation recommended\n")
cat("- **Assumption checking:** Verify ANOVA assumptions (normality, homoscedasticity)\n")
cat("- **Missing data:** Consider impact of incomplete datasets\n")
cat("- **Confounding variables:** Account for potential demographic differences\n\n")
```

---

# Technical Appendix

```{r technical-appendix}
cat("## Analysis Configuration\n\n")
cat("**Datasets:**\n")
for (i in seq_along(params$datasets)) {
  cat("-", params$group_names[i], ":", basename(params$datasets[i]), "\n")
}

cat("\n**Statistical Parameters:**\n")
cat("- **Alpha level:** ", params$alpha_level, "\n")
cat("- **Multiple comparison correction:** ", params$multiple_comparison_correction, "\n")
cat("- **Primary task:** ", analysis_task, "\n")
cat("- **Analysis date:** ", as.character(Sys.Date()), "\n\n")

cat("**Features Analyzed:**\n")
for (i in seq_along(key_features)) {
  cat(i, ". ", gsub("_", " ", key_features[i]), "\n")
}

cat("\n**Software Information:**\n")
cat("- **Package:** LocomotionData R Package\n")
cat("- **Statistical tests:** One-way ANOVA, post-hoc comparisons\n")
cat("- **Effect sizes:** Cohen's d\n")
cat("- **Visualization:** ggplot2, plotly\n\n")

# Display sample sizes by group
cat("**Sample Sizes:**\n")
for (group_name in params$group_names) {
  cat("- ", group_name, ": n = ", length(all_subjects[[group_name]]), "\n")
}

cat("\n## R Session Information\n\n")
print(sessionInfo())
```

---

*This population comparison analysis was generated using the LocomotionData R package. Results should be interpreted within the context of study design, sample characteristics, and clinical expertise.*