longitudinal-analysis.rmd

---
title: "Longitudinal Study Analysis"
author: "Your Name"
date: "2025-03-13"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(tidyverse)
library(broom)
library(ggpubr)
library(kableExtra)
library(patchwork)
library(lme4)
library(lmerTest)
```

## Introduction

This document demonstrates the analysis of longitudinal study data with repeated measurements. The data consists of:

- Subject metadata: Contains information about 15 subjects, including treatment group and sex
- Sample data: Contains measurements from 3 visits per subject, with features and age information

## Data Import

### Reading the Subject and Sample Data

First, let's read in the subject and sample data:

```{r read-data}
# Read the subject data
subjects <- read_csv("/home/kevin/Downloads/dummy_subjects.csv")

# Read the sample data
samples <- read_csv("/home/kevin/Downloads/dummy_samples.csv")

# Display the subject data
subjects %>%
  kable(caption = "Subject Metadata") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

# Display first few rows of the sample data
samples %>%
  head(10) %>%
  kable(caption = "First 10 Rows of Sample Data") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```

### Joining Subject and Sample Data

Now, let's join the subject metadata with the sample data to create a complete dataset:

```{r join-data}
# Join the datasets
complete_data <- samples %>%
  left_join(subjects, by = c("subject" = "subject_id"))

# Display the first few rows of the combined data
complete_data %>%
  head(10) %>%
  kable(caption = "First 10 Rows of Combined Data") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```

### Study Overview

Let's get an overview of our study design:

```{r study-overview}
# Count subjects by treatment and sex
subjects %>%
  count(treatment, sex) %>%
  pivot_wider(names_from = sex, values_from = n) %>%
  mutate(Total = F + M) %>%
  kable(caption = "Subject Distribution by Treatment and Sex") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

# Summarize number of visits per subject
samples %>%
  count(subject) %>%
  summarize(
    min_visits = min(n),
    max_visits = max(n),
    mean_visits = mean(n),
    total_samples = sum(n)
  ) %>%
  kable(caption = "Visit Summary") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```

## Exploratory Data Analysis

Let's explore the experimental data across visits and treatment groups.

### Summary Statistics

First, let's calculate summary statistics for our features by treatment and visit:

```{r summary-stats}
# Calculate summary statistics by treatment and visit for each feature
feature_summary <- complete_data %>%
  group_by(treatment, visit) %>%
  summarize(across(
    c(feature1, feature2, feature3, age),
    list(
      mean = ~mean(.x, na.rm = TRUE),
      sd = ~sd(.x, na.rm = TRUE),
      min = ~min(.x, na.rm = TRUE),
      max = ~max(.x, na.rm = TRUE),
      n = ~sum(!is.na(.x))
    )
  )) %>%
  ungroup()

# Display feature1 summary in a more readable format
feature_summary %>%
  select(treatment, visit, starts_with("feature1")) %>%
  kable(caption = "Summary Statistics for Feature 1") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

# Display feature2 summary
feature_summary %>%
  select(treatment, visit, starts_with("feature2")) %>%
  kable(caption = "Summary Statistics for Feature 2") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

# Display feature3 summary
feature_summary %>%
  select(treatment, visit, starts_with("feature3")) %>%
  kable(caption = "Summary Statistics for Feature 3") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```

### Age Distribution by Treatment

Let's examine the age distribution across treatment groups:

```{r age-distribution}
# Plot age distribution by treatment
ggplot(complete_data, aes(x = treatment, y = age, fill = treatment)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  labs(
    title = "Age Distribution by Treatment Group",
    x = "Treatment",
    y = "Age (years)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# Statistical test for age differences between treatment groups
age_aov <- aov(age ~ treatment, data = complete_data)
tidy(age_aov) %>%
  kable(caption = "ANOVA for Age Differences Between Treatment Groups") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```

### Feature Trends Across Visits

Let's visualize how features change across visits by treatment group:

```{r feature-trends}
# Function to create trend plots for each feature
p1 <- complete_data %>%
    group_by(treatment, visit) %>%
    summarize(
      mean_value = mean(get("feature1"), na.rm = TRUE),
      se = sd(get("feature1"), na.rm = TRUE) / sqrt(n()),
      .groups = "drop"
    ) %>%
    ggplot(aes(x = factor(visit), y = mean_value, color = treatment, group = treatment)) +
    geom_line(size = 1) +
    geom_point(size = 3) +
    geom_errorbar(aes(ymin = mean_value - se, ymax = mean_value + se), width = 0.2) +
    labs(
      title = paste("Mean", "feature1", "by Visit and Treatment"),
      x = "Visit",
      y = paste("Mean", "feature1"),
      color = "Treatment"
    ) +
    theme_minimal()


p2 <- complete_data %>%
    group_by(treatment, visit) %>%
    summarize(
      mean_value = mean(get("feature2"), na.rm = TRUE),
      se = sd(get("feature1"), na.rm = TRUE) / sqrt(n()),
      .groups = "drop"
    ) %>%
    ggplot(aes(x = factor(visit), y = mean_value, color = treatment, group = treatment)) +
    geom_line(size = 1) +
    geom_point(size = 3) +
    geom_errorbar(aes(ymin = mean_value - se, ymax = mean_value + se), width = 0.2) +
    labs(
      title = paste("Mean", "feature2", "by Visit and Treatment"),
      x = "Visit",
      y = paste("Mean", "feature2"),
      color = "Treatment"
    ) +
    theme_minimal()

p3 <- complete_data %>%
    group_by(treatment, visit) %>%
    summarize(
      mean_value = mean(get("feature3"), na.rm = TRUE),
      se = sd(get("feature1"), na.rm = TRUE) / sqrt(n()),
      .groups = "drop"
    ) %>%
    ggplot(aes(x = factor(visit), y = mean_value, color = treatment, group = treatment)) +
    geom_line(size = 1) +
    geom_point(size = 3) +
    geom_errorbar(aes(ymin = mean_value - se, ymax = mean_value + se), width = 0.2) +
    labs(
      title = paste("Mean", "feature3", "by Visit and Treatment"),
      x = "Visit",
      y = paste("Mean", "feature3"),
      color = "Treatment"
    ) +
    theme_minimal()

# Display the plots
p1 / p2 / p3
```

### Individual Subject Trajectories

Let's visualize individual subject trajectories for each feature:

```{r individual-trajectories}

ggplot(complete_data, aes(x = factor(visit), y = get("feature1"), 
                        group = subject, color = treatment)) +
  geom_line(alpha = 0.6) +
  geom_point(size = 2, alpha = 0.7) +
  facet_wrap(~treatment) +
  labs(
    title = paste("Individual Trajectories -", "feature1"),
    x = "Visit",
    y = "feature1"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

```

### Feature Distributions by Treatment and Visit

Let's examine the distribution of features across treatments and visits:

```{r feature-distributions}
# Create box plots for feature1 by treatment and visit
ggplot(complete_data, aes(x = factor(visit), y = feature1, fill = treatment)) +
  geom_boxplot(alpha = 0.7) +
  labs(
    title = "Feature 1 Distribution by Treatment and Visit",
    x = "Visit",
    y = "Feature 1",
    fill = "Treatment"
  ) +
  theme_minimal()

# Create box plots for feature2 by treatment and visit
ggplot(complete_data, aes(x = factor(visit), y = feature2, fill = treatment)) +
  geom_boxplot(alpha = 0.7) +
  labs(
    title = "Feature 2 Distribution by Treatment and Visit",
    x = "Visit",
    y = "Feature 2",
    fill = "Treatment"
  ) +
  theme_minimal()

# Create box plots for feature3 by treatment and visit
ggplot(complete_data, aes(x = factor(visit), y = feature3, fill = treatment)) +
  geom_boxplot(alpha = 0.7) +
  labs(
    title = "Feature 3 Distribution by Treatment and Visit",
    x = "Visit",
    y = "Feature 3",
    fill = "Treatment"
  ) +
  theme_minimal()
```

### Correlation Analysis

Let's examine the correlations between numerical features:

```{r correlation}
# Select only numeric features for correlation analysis
numeric_data <- complete_data %>%
  select(feature1, feature2, feature3, age)

# Calculate correlation matrix
cor_matrix <- cor(numeric_data, use = "pairwise.complete.obs")

# Create a correlation heatmap
corrplot::corrplot(
  cor_matrix,
  method = "color",
  type = "upper",
  order = "hclust",
  tl.col = "black",
  tl.srt = 45,
  addCoef.col = "black",
  number.cex = 0.7
)

# Create scatter plots for feature pairs
ggplot(complete_data, aes(x = feature1, y = feature2, color = treatment)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~visit) +
  labs(
    title = "Relationship between Feature 1 and Feature 2",
    x = "Feature 1",
    y = "Feature 2",
    color = "Treatment"
  ) +
  theme_minimal()
```

### Sex Differences in Features

Let's explore whether there are differences in features between male and female subjects:

```{r sex-differences}
# Create box plots for feature1 by sex and treatment
ggplot(complete_data, aes(x = treatment, y = feature1, fill = sex)) +
  geom_boxplot(alpha = 0.7) +
  labs(
    title = "Feature 1 by Sex and Treatment",
    x = "Treatment",
    y = "Feature 1",
    fill = "Sex"
  ) +
  theme_minimal()

# Create box plots for feature2 by sex and treatment
ggplot(complete_data, aes(x = treatment, y = feature2, fill = sex)) +
  geom_boxplot(alpha = 0.7) +
  labs(
    title = "Feature 2 by Sex and Treatment",
    x = "Treatment",
    y = "Feature 2",
    fill = "Sex"
  ) +
  theme_minimal()

# Create box plots for feature3 by sex and treatment
ggplot(complete_data, aes(x = treatment, y = feature3, fill = sex)) +
  geom_boxplot(alpha = 0.7) +
  labs(
    title = "Feature 3 by Sex and Treatment",
    x = "Treatment",
    y = "Feature 3",
    fill = "Sex"
  ) +
  theme_minimal()
```

## Statistical Analysis

### Linear Mixed-Effects Models

For longitudinal data with repeated measures, linear mixed-effects models are appropriate to account for within-subject correlations:

```{r mixed-models}
# Convert visit to a factor for proper modeling
complete_data$visit_f <- factor(complete_data$visit)

# Linear mixed model for feature1
model_feature1 <- lmer(feature1 ~ treatment * visit_f + sex + age + (1|subject), 
                       data = complete_data)

# Display model summary
summary(model_feature1)

# Tidy model output for feature1
tidy_model1 <- broom.mixed::tidy(model_feature1, conf.int = TRUE)
tidy_model1 %>%
  filter(effect == "fixed") %>%
  select(-effect) %>%
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>%
  kable(caption = "Linear Mixed Model Results for Feature 1") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

# ANOVA-style tests for fixed effects
anova(model_feature1) %>%
  tidy() %>%
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>%
  kable(caption = "ANOVA Table for Feature 1 Mixed Model") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

# Linear mixed model for feature2
model_feature2 <- lmer(feature2 ~ treatment * visit_f + sex + age + (1|subject), 
                       data = complete_data)

# ANOVA-style tests for fixed effects (feature2)
anova(model_feature2) %>%
  tidy() %>%
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>%
  kable(caption = "ANOVA Table for Feature 2 Mixed Model") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

# Linear mixed model for feature3
model_feature3 <- lmer(feature3 ~ treatment * visit_f + sex + age + (1|subject), 
                       data = complete_data)

# ANOVA-style tests for fixed effects (feature3)
anova(model_feature3) %>%
  tidy() %>%
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>%
  kable(caption = "ANOVA Table for Feature 3 Mixed Model") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```

### Pairwise Comparisons

Let's perform pairwise comparisons between treatment groups for feature1:

```{r pairwise-comparisons}
# Pairwise comparisons for feature1 between treatment groups
emmeans::emmeans(model_feature1, pairwise ~ treatment) %>%
  .$contrasts %>%
  broom::tidy() %>%
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>%
  kable(caption = "Pairwise Comparisons Between Treatment Groups for Feature 1") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

# Pairwise comparisons for feature1 between visits
emmeans::emmeans(model_feature1, pairwise ~ visit_f) %>%
  .$contrasts %>%
  broom::tidy() %>%
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>%
  kable(caption = "Pairwise Comparisons Between Visits for Feature 1") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

# Treatment x Visit interaction comparisons
emmeans::emmeans(model_feature1, pairwise ~ treatment | visit_f) %>%
  .$contrasts %>%
  broom::tidy() %>%
  mutate(across(where(is.numeric), ~round(.x, digits = 3))) %>%
  kable(caption = "Pairwise Comparisons Between Treatment Groups at Each Visit for Feature 1") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```

### Traditional Repeated Measures ANOVA

For comparison, let's also run a traditional repeated measures ANOVA:

```{r rmanova}
# Reshape data to wide format for repeated measures ANOVA
wide_data <- complete_data %>%
  select(subject, treatment, sex, visit, feature1) %>%
  pivot_wider(
    id_cols = c(subject, treatment, sex),
    names_from = visit,
    values_from = feature1,
    names_prefix = "visit"
  )

# Repeated measures ANOVA using ez package
if (requireNamespace("ez", quietly = TRUE)) {
  wide_data_ez <- complete_data %>%
    mutate(visit_f = factor(visit)) %>%
    rename(id = subject)
  
  feature1_aov <- ez::ezANOVA(
    data = wide_data_ez,
    dv = feature1,
    wid = id,
    within = visit_f,
    between = c(treatment, sex)
  )
  
  print(feature1_aov)
}
```

## Advanced Visualizations

### Spaghetti Plots with Fitted Trajectories

Let's create spaghetti plots showing individual trajectories with fitted group means:

```{r spaghetti-plots}
# Function to create spaghetti plots
create_spaghetti_plot <- function(feature_name) {
  # Calculate group means
  group_means <- complete_data %>%
    group_by(treatment, visit) %>%
    summarize(mean_value = mean(get(feature_name), na.rm = TRUE),
              .groups = "drop")
  
  # Create plot
  ggplot() +
    # Individual trajectories
    geom_line(data = complete_data, 
              aes(x = visit, y = get(feature_name), group = subject, color = treatment),
              alpha = 0.3) +
    # Group mean trajectories
    geom_line(data = group_means,
              aes(x = visit, y = mean_value, group = treatment, color = treatment),
              size = 1.5) +
    geom_point(data = group_means,
               aes(x = visit, y = mean_value, color = treatment),
               size = 3) +
    facet_wrap(~treatment) +
    labs(
      title = paste("Individual and Mean Trajectories -", feature_name),
      x = "Visit",
      y = feature_name
    ) +
    theme_minimal() +
    theme(legend.position = "none")
}

# Create spaghetti plots for each feature
create_spaghetti_plot("feature1")
create_spaghetti_plot("feature2")
create_spaghetti_plot("feature3")
```


### Change from Baseline

Let's analyze the change from baseline for each feature:

```{r baseline-change}
# Calculate change from baseline (visit 1)
baseline_data <- complete_data %>%
  filter(visit == 1) %>%
  select(subject, feature1, feature2, feature3) %>%
  rename(
    baseline_feature1 = feature1,
    baseline_feature2 = feature2,
    baseline_feature3 = feature3
  )

# Join with original data and calculate changes
change_data <- complete_data %>%
  left_join(baseline_data, by = "subject") %>%
  mutate(
    change_feature1 = feature1 - baseline_feature1,
    change_feature2 = feature2 - baseline_feature2,
    change_feature3 = feature3 - baseline_feature3
  )

# Plot change from baseline for feature1
ggplot(change_data %>% filter(visit > 1), 
       aes(x = factor(visit), y = change_feature1, fill = treatment)) +
  geom_boxplot(alpha = 0.7) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    title = "Change from Baseline - Feature 1",
    x = "Visit",
    y = "Change from Baseline",
    fill = "Treatment"
  ) +
  theme_minimal()

# Statistical test for change from baseline differences
visit2_changes <- change_data %>% 
  filter(visit == 2) %>%
  aov(change_feature1 ~ treatment, data = .)

tidy(visit2_changes) %>%
  kable(caption = "ANOVA for Feature 1 Change at Visit 2") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

visit3_changes <- change_data %>% 
  filter(visit == 3) %>%
  aov(change_feature1 ~ treatment, data = .)

tidy(visit3_changes) %>%
  kable(caption = "ANOVA for Feature 1 Change at Visit 3") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```