scripts.qmd

---
title: "snippets"
---

# Write the synthetic data to CSV files

# Function to create the CSV files from the synthetic data
create_synthetic_datasets <- function() {
  # Create diet_data dataframe
  diet_data <- read.csv("diet_fatty_acids.csv")
  
  # Create milk_data dataframe
  milk_data <- read.csv("milk_fatty_acids.csv")
  
  return(list(diet_data = diet_data, milk_data = milk_data))
}

# Basic visualization functions

# Plot individual fatty acid trajectories
plot_fa_trajectory <- function(data_long, target_fa) {
  require(ggplot2)
  
  fa_data <- data_long %>% 
    filter(fatty_acid == target_fa)
  
  p <- ggplot(fa_data, aes(x = days_postpartum, y = abundance, 
                       group = woman_id, color = factor(woman_id))) +
    geom_line() +
    geom_point() +
    labs(title = paste0(target_fa, " Abundance Over Time"),
         x = "Days Postpartum",
         y = "Abundance (% of total FAs)",
         color = "Woman ID") +
    theme_minimal()
  
  return(p)
}

# Multivariate visualization with PCA
pca_analysis <- function(data, id_cols = c("woman_id", "timepoint")) {
  require(ggplot2)
  
  # Extract just the FA columns for PCA
  fa_cols <- data %>% select(-all_of(id_cols))
  
  # Run PCA
  pca_result <- prcomp(fa_cols, scale. = TRUE)
  
  # Extract scores
  pca_scores <- as.data.frame(pca_result$x)
  
  # Add back identifiers
  for(col in id_cols) {
    pca_scores[[col]] <- data[[col]]
  }
  
  # Create plot
  p <- ggplot(pca_scores, aes(x = PC1, y = PC2, 
                         color = factor(woman_id), 
                         shape = factor(timepoint))) +
    geom_point(size = 3) +
    labs(title = "PCA of Fatty Acid Profiles",
         color = "Woman ID",
         shape = "Timepoint") +
    theme_minimal()
  
  # Return both plot and PCA results
  return(list(plot = p, pca = pca_result, scores = pca_scores))
}

# Correlation analysis between diet and milk
correlation_analysis <- function(diet_long, milk_long) {
  require(ggplot2)
  
  # Join diet and milk data for correlation analysis
  joined_data <- diet_long %>%
    inner_join(milk_long, 
               by = c("woman_id", "timepoint", "days_postpartum", "fatty_acid"),
               suffix = c("_diet", "_milk"))
  
  # Calculate correlations for each fatty acid
  correlations <- joined_data %>%
    group_by(fatty_acid) %>%
    summarize(correlation = cor(abundance_diet, abundance_milk, method = "spearman"),
              p_value = cor.test(abundance_diet, abundance_milk, method = "spearman")$p.value)
  
  # Plot the correlations
  p <- ggplot(correlations, aes(x = reorder(fatty_acid, correlation), y = correlation)) +
    geom_bar(stat = "identity", aes(fill = p_value < 0.05)) +
    coord_flip() +
    labs(title = "Correlation Between Diet and Milk Fatty Acid Abundance",
         x = "Fatty Acid",
         y = "Spearman Correlation",
         fill = "P < 0.05") +
    theme_minimal()
  
  return(list(plot = p, correlations = correlations))
}

# Heatmap visualization by timepoint
heatmap_by_timepoint <- function(data, timepoint_val, scale = "column") {
  require(pheatmap)
  
  # Prepare data for heatmap
  heatmap_data <- data %>%
    filter(timepoint == timepoint_val) %>%
    select(-timepoint, -days_postpartum) %>%
    column_to_rownames("woman_id")
  
  # Create heatmap
  pheatmap(heatmap_data,
           scale = scale,
           clustering_distance_rows = "euclidean",
           clustering_distance_cols = "euclidean",
           main = paste0("Fatty Acid Profiles at Timepoint ", timepoint_val))
}

# Calculate fatty acid composition by class
calculate_fa_classes <- function(fa_long_data) {
  # Create a function to classify fatty acids
  classify_fa <- function(fa_name) {
    if (grepl("C[0-9]+:0$", fa_name)) {
      return("Saturated")
    } else if (grepl("C[0-9]+:1", fa_name)) {
      if (grepl("t$", fa_name)) {
        return("Trans")
      } else {
        return("Monounsaturated")
      }
    } else if (grepl("n3$", fa_name)) {
      return("Omega-3")
    } else if (grepl("n6$", fa_name)) {
      return("Omega-6")
    } else if (grepl("t$", fa_name)) {
      return("Trans")
    } else {
      return("Other")
    }
  }
  
  # Add classifications to the data
  fa_long_data$fa_class <- sapply(fa_long_data$fatty_acid, classify_fa)
  
  # Calculate totals by class
  class_totals <- fa_long_data %>%
    group_by(woman_id, timepoint, days_postpartum, fa_class) %>%
    summarize(total_abundance = sum(abundance), .groups = "drop")
  
  return(list(classified_data = fa_long_data, class_totals = class_totals))
}

# Create a mixed effects model for a specific fatty acid
mixed_model_analysis <- function(long_data, target_fa) {
  require(lme4)
  require(lmerTest)
  
  # Filter for the target fatty acid
  fa_data <- long_data %>% 
    filter(fatty_acid == target_fa)
  
  # Linear mixed model with random intercept by woman
  model <- lmer(abundance ~ days_postpartum + (1|woman_id), data = fa_data)
  
  # Visualize individual trajectories with overall trend
  p <- ggplot(fa_data, aes(x = days_postpartum, y = abundance, group = woman_id)) +
    geom_line(alpha = 0.3) +
    geom_smooth(aes(group = 1), method = "lm", color = "red") +
    labs(title = paste0(target_fa, " Abundance Over Time with Overall Trend"),
         x = "Days Postpartum",
         y = paste0(target_fa, " Abundance (%)")) +
    theme_minimal()
  
  return(list(model = model, plot = p))
}

# Complete analysis pipeline
run_fatty_acid_analysis <- function(diet_data, milk_data) {
  require(tidyverse)
  
  # Convert to long format
  diet_long <- diet_data %>%
    pivot_longer(cols = -c(woman_id, timepoint, days_postpartum),
                 names_to = "fatty_acid", 
                 values_to = "abundance")
  
  milk_long <- milk_data %>%
    pivot_longer(cols = -c(woman_id, timepoint, days_postpartum),
                 names_to = "fatty_acid", 
                 values_to = "abundance")
  
  # List of key fatty acids to analyze
  key_fas <- c("C16:0", "C18:0", "C18:1n9", "C18:2n6", "C18:3n3", "C20:4n6", "C22:6n3")
  
  # Run analyses
  results <- list()
  
  # 1. Individual trajectories for key FAs
  results$trajectories <- lapply(key_fas, function(fa) {
    plot_fa_trajectory(milk_long, fa)
  })
  names(results$trajectories) <- key_fas
  
  # 2. PCA analysis
  results$pca <- pca_analysis(
    milk_data %>% select(-days_postpartum), 
    id_cols = c("woman_id", "timepoint")
  )
  
  # 3. Diet-milk correlations
  results$correlations <- correlation_analysis(diet_long, milk_long)
  
  # 4. Heatmaps by timepoint
  results$heatmaps <- lapply(1:6, function(tp) {
    heatmap_by_timepoint(milk_data, tp)
  })
  
  # 5. FA class analysis
  milk_classes <- calculate_fa_classes(milk_long)
  diet_classes <- calculate_fa_classes(diet_long)
  
  results$class_analysis <- list(
    milk = milk_classes,
    diet = diet_classes
  )
  
  # 6. Mixed models for key FAs
  results$mixed_models <- lapply(key_fas, function(fa) {
    mixed_model_analysis(milk_long, fa)
  })
  names(results$mixed_models) <- key_fas
  
  return(results)
}

# Example usage:
# data <- create_synthetic_datasets()
# analysis_results <- run_fatty_acid_analysis(data$diet_data, data$milk_data)
# analysis_results$trajectories$`C22:6n3`  # View DHA trajectory
# analysis_results$pca$plot  # View PCA plot