---
title: "Intervention Study Analysis"
subtitle: "Pre/Post Treatment 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:
  pre_data: "baseline_data.parquet"
  post_data: "followup_data.parquet"
  intervention_name: "Gait Training"
  study_id: "STUDY_001"
  analysis_group: "treatment"
  control_group: "control"
  alpha_level: 0.05
  effect_size_threshold: 0.5
  include_statistics: 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(effsize)
library(lme4)
library(emmeans)
library(broom)
library(report)

# Set 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"
  ))
```

# Study Overview

```{r load-data, include=FALSE}
# Load pre and post intervention data
if (!file.exists(params$pre_data)) {
  stop("Pre-intervention data file not found: ", params$pre_data)
}
if (!file.exists(params$post_data)) {
  stop("Post-intervention data file not found: ", params$post_data)
}

pre_loco <- loadLocomotionData(params$pre_data)
post_loco <- loadLocomotionData(params$post_data)

# Get study information
pre_subjects <- getSubjects(pre_loco)
post_subjects <- getSubjects(post_loco)
common_subjects <- intersect(pre_subjects, post_subjects)

pre_tasks <- getTasks(pre_loco)
post_tasks <- getTasks(post_loco)
common_tasks <- intersect(pre_tasks, post_tasks)

pre_features <- getFeatures(pre_loco)
post_features <- getFeatures(post_loco)
common_features <- intersect(pre_features, post_features)

# Primary analysis parameters
primary_task <- if ("normal_walk" %in% common_tasks) "normal_walk" else common_tasks[1]
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)
```

**Study ID:** `r params$study_id`  
**Intervention:** `r params$intervention_name`  
**Analysis Date:** `r Sys.Date()`  
**Primary Task:** `r primary_task`

## Participant Overview

- **Total Participants:** `r length(common_subjects)`
- **Pre-intervention Assessments:** `r length(pre_subjects)`
- **Post-intervention Assessments:** `r length(post_subjects)`
- **Complete Data (Pre + Post):** `r length(common_subjects)`
- **Tasks Analyzed:** `r paste(common_tasks, collapse = ", ")`
- **Features Analyzed:** `r length(common_features)`

```{r participant-table}
# Create participant summary table
participant_summary <- data.frame(
  Participant_ID = common_subjects,
  Pre_Assessment = common_subjects %in% pre_subjects,
  Post_Assessment = common_subjects %in% post_subjects,
  Status = "Complete"
)

# Add missing participants
missing_pre <- setdiff(post_subjects, pre_subjects)
missing_post <- setdiff(pre_subjects, post_subjects)

if (length(missing_pre) > 0) {
  missing_pre_df <- data.frame(
    Participant_ID = missing_pre,
    Pre_Assessment = FALSE,
    Post_Assessment = TRUE,
    Status = "Missing Pre"
  )
  participant_summary <- rbind(participant_summary, missing_pre_df)
}

if (length(missing_post) > 0) {
  missing_post_df <- data.frame(
    Participant_ID = missing_post,
    Pre_Assessment = TRUE,
    Post_Assessment = FALSE,
    Status = "Missing Post"
  )
  participant_summary <- rbind(participant_summary, missing_post_df)
}

DT::datatable(participant_summary, 
              options = list(pageLength = 10, dom = 't'),
              rownames = FALSE) %>%
  formatStyle("Status",
    backgroundColor = styleEqual(c("Complete", "Missing Pre", "Missing Post"),
                                c("lightgreen", "lightyellow", "lightyellow")))
```

---

# Primary Outcome Analysis

```{r primary-analysis}
# Calculate mean patterns for all participants
pre_results <- list()
post_results <- list()

for (subject in common_subjects) {
  # Pre-intervention patterns
  pre_patterns <- getMeanPatterns(pre_loco, subject, primary_task, key_features)
  pre_results[[subject]] <- pre_patterns
  
  # Post-intervention patterns
  post_patterns <- getMeanPatterns(post_loco, subject, primary_task, key_features)
  post_results[[subject]] <- post_patterns
}

# Calculate group means
group_pre_means <- list()
group_post_means <- list()

for (feature in key_features) {
  # Extract feature data across subjects
  pre_feature_data <- sapply(pre_results, function(x) x[[feature]])
  post_feature_data <- sapply(post_results, function(x) x[[feature]])
  
  # Calculate group means
  group_pre_means[[feature]] <- rowMeans(pre_feature_data, na.rm = TRUE)
  group_post_means[[feature]] <- rowMeans(post_feature_data, na.rm = TRUE)
}
```

## Feature-Specific Changes

```{r change-analysis}
# Calculate change metrics for each feature
change_summary <- data.frame()

for (feature in key_features) {
  pre_data <- group_pre_means[[feature]]
  post_data <- group_post_means[[feature]]
  
  # Calculate ROM for pre and post
  pre_rom <- max(pre_data) - min(pre_data)
  post_rom <- max(post_data) - min(post_data)
  rom_change <- post_rom - pre_rom
  rom_change_percent <- (rom_change / pre_rom) * 100
  
  # Calculate peak values
  pre_peak <- max(abs(pre_data))
  post_peak <- max(abs(post_data))
  peak_change <- post_peak - pre_peak
  peak_change_percent <- (peak_change / pre_peak) * 100
  
  # Store results
  feature_summary <- data.frame(
    Feature = gsub("_", " ", feature),
    Pre_ROM = round(pre_rom * 180/pi, 1),
    Post_ROM = round(post_rom * 180/pi, 1),
    ROM_Change = round(rom_change * 180/pi, 1),
    ROM_Change_Percent = round(rom_change_percent, 1),
    Pre_Peak = round(pre_peak * 180/pi, 1),
    Post_Peak = round(post_peak * 180/pi, 1),
    Peak_Change = round(peak_change * 180/pi, 1),
    Peak_Change_Percent = round(peak_change_percent, 1)
  )
  
  change_summary <- rbind(change_summary, feature_summary)
}

kable(change_summary,
      col.names = c("Feature", "Pre ROM (°)", "Post ROM (°)", "Δ ROM (°)", "Δ ROM (%)",
                    "Pre Peak (°)", "Post Peak (°)", "Δ Peak (°)", "Δ Peak (%)"),
      align = c("l", rep("c", 8))) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
  column_spec(c(4, 5, 8, 9), 
              color = ifelse(c(change_summary$ROM_Change, change_summary$Peak_Change) > 0, 
                           "blue", "red"))
```

---

# Statistical Analysis

```{r statistical-tests, eval=params$include_statistics}
if (params$include_statistics) {
  
  # Prepare data for statistical analysis
  stats_results <- list()
  
  for (feature in key_features) {
    # Extract individual subject ROM values
    pre_roms <- sapply(common_subjects, function(subj) {
      data <- pre_results[[subj]][[feature]]
      max(data) - min(data)
    })
    
    post_roms <- sapply(common_subjects, function(subj) {
      data <- post_results[[subj]][[feature]]
      max(data) - min(data)
    })
    
    # Paired t-test
    t_test <- t.test(post_roms, pre_roms, paired = TRUE)
    
    # Effect size (Cohen's d)
    effect_size <- cohen.d(post_roms, pre_roms, paired = TRUE)
    
    # Store results
    stats_results[[feature]] <- list(
      t_statistic = t_test$statistic,
      p_value = t_test$p.value,
      ci_lower = t_test$conf.int[1] * 180/pi,
      ci_upper = t_test$conf.int[2] * 180/pi,
      effect_size = effect_size$estimate,
      effect_magnitude = effect_size$magnitude
    )
  }
  
  # Create statistical summary table
  stats_df <- data.frame(
    Feature = gsub("_", " ", names(stats_results)),
    t_statistic = round(sapply(stats_results, function(x) x$t_statistic), 3),
    p_value = sapply(stats_results, function(x) 
      ifelse(x$p_value < 0.001, "<0.001", round(x$p_value, 3))),
    CI_Lower = round(sapply(stats_results, function(x) x$ci_lower), 2),
    CI_Upper = round(sapply(stats_results, function(x) x$ci_upper), 2),
    Cohens_d = round(sapply(stats_results, function(x) x$effect_size), 3),
    Effect_Magnitude = sapply(stats_results, function(x) x$effect_magnitude),
    Significant = sapply(stats_results, function(x) 
      ifelse(x$p_value < params$alpha_level, "Yes", "No"))
  )
  
  cat("## Statistical Test Results\n\n")
  cat("**Test:** Paired t-test (Pre vs Post)\n")
  cat("**Alpha Level:** ", params$alpha_level, "\n")
  cat("**Effect Size Threshold:** ", params$effect_size_threshold, "\n\n")
  
  print(kable(stats_df,
        col.names = c("Feature", "t", "p-value", "95% CI Lower (°)", "95% CI Upper (°)", 
                      "Cohen's d", "Effect Size", "Significant"),
        align = c("l", rep("c", 7))) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
    column_spec(8, color = ifelse(stats_df$Significant == "Yes", "green", "red")) %>%
    column_spec(7, background = ifelse(stats_df$Effect_Magnitude %in% c("large", "medium"), 
                                     "lightblue", "white")))
}
```

```{r effect-size-interpretation, eval=params$include_statistics, results='asis'}
if (params$include_statistics) {
  cat("\n### Effect Size Interpretation\n\n")
  
  significant_features <- stats_df$Feature[stats_df$Significant == "Yes"]
  large_effects <- stats_df$Feature[stats_df$Effect_Magnitude == "large"]
  medium_effects <- stats_df$Feature[stats_df$Effect_Magnitude == "medium"]
  
  if (length(significant_features) > 0) {
    cat("**Statistically Significant Changes:**\n")
    cat(paste("-", significant_features, collapse = "\n"))
    cat("\n\n")
  }
  
  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")
  }
}
```

---

# Visualization of Changes

```{r visualization-setup, eval=params$generate_plots}
if (params$generate_plots) {
  phase_vector <- seq(0, 100, length.out = 150)
}
```

## Pre vs Post Comparison

```{r pre-post-comparison, eval=params$generate_plots, fig.height=10}
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 data
    plot_data <- data.frame(
      Phase = rep(phase_vector, 2),
      Angle = c(group_pre_means[[feature]] * 180/pi, 
                group_post_means[[feature]] * 180/pi),
      Timepoint = rep(c("Pre", "Post"), each = 150)
    )
    
    # Create plot
    p <- ggplot(plot_data, aes(x = Phase, y = Angle, color = Timepoint)) +
      geom_line(size = 1.2) +
      labs(
        title = gsub("_", " ", feature),
        x = "Gait Cycle (%)",
        y = "Angle (degrees)",
        color = "Timepoint"
      ) +
      scale_color_manual(values = c("Pre" = "#E31A1C", "Post" = "#1F78B4")) +
      theme(legend.position = "bottom")
    
    comparison_plots[[i]] <- p
  }
  
  # Arrange plots
  if (length(comparison_plots) >= 4) {
    combined_plot <- (comparison_plots[[1]] | comparison_plots[[2]]) /
                    (comparison_plots[[3]] | comparison_plots[[4]])
  } else if (length(comparison_plots) >= 2) {
    combined_plot <- comparison_plots[[1]] | comparison_plots[[2]]
  } else {
    combined_plot <- comparison_plots[[1]]
  }
  
  print(combined_plot)
}
```

## Interactive Change Visualization

```{r interactive-change-plot, eval=params$generate_plots}
if (params$generate_plots && length(key_features) > 0) {
  
  # Create interactive plot for primary feature
  primary_feature <- key_features[1]
  
  change_data <- data.frame(
    Phase = phase_vector,
    Pre = group_pre_means[[primary_feature]] * 180/pi,
    Post = group_post_means[[primary_feature]] * 180/pi
  )
  
  change_data$Difference <- change_data$Post - change_data$Pre
  
  # Create interactive plot
  p_interactive <- ggplot(change_data, aes(x = Phase)) +
    geom_ribbon(aes(ymin = pmin(Pre, Post), ymax = pmax(Pre, Post)),
                fill = "lightblue", alpha = 0.3) +
    geom_line(aes(y = Pre, text = paste("Phase:", Phase, "%<br>Pre:", round(Pre, 1), "°")), 
              color = "#E31A1C", size = 1.2) +
    geom_line(aes(y = Post, text = paste("Phase:", Phase, "%<br>Post:", round(Post, 1), "°<br>Change:", round(Difference, 1), "°")), 
              color = "#1F78B4", size = 1.2) +
    labs(
      title = paste("Interactive Change Analysis:", gsub("_", " ", primary_feature)),
      subtitle = "Hover for detailed values | Blue ribbon shows change magnitude",
      x = "Gait Cycle (%)",
      y = "Angle (degrees)"
    )
  
  ggplotly(p_interactive, tooltip = "text") %>%
    layout(hovermode = "x unified")
}
```

---

# Individual Participant Analysis

```{r individual-analysis}
# Calculate individual changes
individual_changes <- data.frame()

for (subject in common_subjects) {
  for (feature in key_features[1:2]) {  # Limit to first 2 features for display
    pre_data <- pre_results[[subject]][[feature]]
    post_data <- post_results[[subject]][[feature]]
    
    pre_rom <- (max(pre_data) - min(pre_data)) * 180/pi
    post_rom <- (max(post_data) - min(post_data)) * 180/pi
    change <- post_rom - pre_rom
    change_percent <- (change / pre_rom) * 100
    
    individual_row <- data.frame(
      Subject = subject,
      Feature = gsub("_", " ", feature),
      Pre_ROM = round(pre_rom, 1),
      Post_ROM = round(post_rom, 1),
      Change = round(change, 1),
      Change_Percent = round(change_percent, 1),
      Responder = ifelse(abs(change_percent) > 10, "Yes", "No")
    )
    
    individual_changes <- rbind(individual_changes, individual_row)
  }
}

cat("## Individual Participant Responses\n\n")
cat("**Response Criterion:** >10% change in ROM\n\n")

DT::datatable(individual_changes,
              options = list(pageLength = 15, dom = 'tip'),
              rownames = FALSE) %>%
  formatStyle("Responder",
    backgroundColor = styleEqual(c("Yes", "No"), c("lightgreen", "lightgray"))) %>%
  formatStyle("Change",
    color = styleInterval(c(-5, 5), c("red", "black", "blue")))
```

```{r responder-analysis}
# Analyze response rates
response_summary <- individual_changes %>%
  group_by(Feature) %>%
  summarise(
    Total_Participants = n(),
    Responders = sum(Responder == "Yes"),
    Response_Rate = round((Responders / Total_Participants) * 100, 1),
    Mean_Change = round(mean(Change), 1),
    SD_Change = round(sd(Change), 1),
    .groups = 'drop'
  )

cat("\n### Response Rate Summary\n\n")

kable(response_summary,
      col.names = c("Feature", "Total N", "Responders", "Response Rate (%)", 
                    "Mean Change (°)", "SD Change (°)"),
      align = c("l", rep("c", 5))) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```

---

# Clinical Significance

```{r clinical-significance, results='asis'}
# Determine clinical significance
clinical_assessment <- list()

# Analyze effect sizes and significance
for (feature in names(stats_results)) {
  result <- stats_results[[feature]]
  
  is_significant <- result$p_value < params$alpha_level
  is_large_effect <- abs(result$effect_size) > params$effect_size_threshold
  
  if (is_significant && is_large_effect) {
    clinical_assessment[[feature]] <- "Clinically significant change"
  } else if (is_significant) {
    clinical_assessment[[feature]] <- "Statistically significant, small effect"
  } else if (is_large_effect) {
    clinical_assessment[[feature]] <- "Large effect, not statistically significant"
  } else {
    clinical_assessment[[feature]] <- "No significant change"
  }
}

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

if (params$include_statistics) {
  for (feature in names(clinical_assessment)) {
    cat("**", gsub("_", " ", feature), ":** ", clinical_assessment[[feature]], "\n")
  }
  
  cat("\n### Overall Study Conclusion\n\n")
  
  significant_count <- sum(sapply(stats_results, function(x) x$p_value < params$alpha_level))
  large_effect_count <- sum(sapply(stats_results, function(x) abs(x$effect_size) > params$effect_size_threshold))
  
  if (significant_count >= length(key_features) * 0.5) {
    cat("**Conclusion:** The intervention shows statistically significant improvements in biomechanical parameters.\n\n")
  } else if (large_effect_count >= length(key_features) * 0.5) {
    cat("**Conclusion:** The intervention shows clinically meaningful changes with large effect sizes.\n\n")
  } else {
    cat("**Conclusion:** The intervention shows limited evidence of biomechanical changes.\n\n")
  }
} else {
  cat("Statistical analysis not performed. Enable statistical analysis for detailed interpretation.\n\n")
}

# Recommendations
cat("### Clinical Recommendations\n\n")
cat("- **Follow-up:** Continue monitoring participants for sustained improvements\n")
cat("- **Dosage:** Consider intervention intensity and duration optimization\n")
cat("- **Population:** Evaluate intervention effectiveness in different populations\n")
cat("- **Mechanisms:** Investigate underlying mechanisms of observed changes\n\n")
```

---

# Study Limitations and Considerations

```{r limitations, results='asis'}
cat("## Study Limitations\n\n")

limitations <- c()

# Sample size consideration
if (length(common_subjects) < 20) {
  limitations <- c(limitations, 
    paste("**Small Sample Size:** n =", length(common_subjects), "may limit statistical power"))
}

# Missing data
missing_rate <- (length(missing_pre) + length(missing_post)) / 
                (length(pre_subjects) + length(post_subjects)) * 100
if (missing_rate > 10) {
  limitations <- c(limitations,
    paste("**Missing Data:** ", round(missing_rate, 1), "% missing data may introduce bias"))
}

# Analysis limitations
limitations <- c(limitations,
  "**Temporal Factors:** Changes may be influenced by time-related factors beyond the intervention",
  "**Individual Variability:** Response heterogeneity suggests need for personalized approaches",
  "**Control Group:** Consider matched control group for stronger causal inference")

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

cat("## Future Research Directions\n\n")
cat("- **Longitudinal Follow-up:** Assess long-term retention of improvements\n")
cat("- **Mechanistic Studies:** Investigate underlying biomechanical adaptations\n")
cat("- **Dose-Response:** Optimize intervention parameters\n")
cat("- **Predictive Factors:** Identify characteristics of intervention responders\n\n")
```

---

# Technical Appendix

```{r technical-appendix}
cat("## Analysis Parameters\n\n")
cat("**Pre-intervention Data:** ", basename(params$pre_data), "\n")
cat("**Post-intervention Data:** ", basename(params$post_data), "\n")
cat("**Analysis Date:** ", as.character(Sys.Date()), "\n")
cat("**Software:** LocomotionData R Package\n")
cat("**Statistical Methods:** Paired t-tests, Cohen's d effect size\n")
cat("**Alpha Level:** ", params$alpha_level, "\n")
cat("**Effect Size Threshold:** ", params$effect_size_threshold, "\n\n")

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

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

---

*This intervention analysis was generated using the LocomotionData R package. Results should be interpreted in the context of study design and clinical expertise.*