2_analysis.R

#*######################################################
#*######################################################
#*######################################################
#*#########    BIOCLOCK: Statistical Analysis   ########
#*######################################################
#*######################################################
#*######################################################

###################################
####           Set-up           ###
###################################

library(readxl)
library(writexl)
library(dplyr)
library(tidyr)
library(stringr)
library(arrow)
library(lme4)
library(lmerTest)
library(MuMIn)

setwd(paste0("/", file.path("data", "user_homes", "mennovd", "BIOKLOK")))

dir.create("Results/Analysis", recursive = TRUE, showWarnings = FALSE)



###################################
####        Data Loading        ###
###################################

load("Data/Objects/pheno.Rdata")

clocks <- c(
    "epigenetic_age_GrimAge",
    "epigenetic_age_Horvath"
)

bloodcells <- c(
    "B_cells_naive",
    "B_cells_memory",
    "CD4T_cells_naive",
    "CD4T_cells_memory",
    "CD8T_cells_naive",
    "CD8T_cells_memory",
    "T_regulatory_cells",
    "natural_killer_cells",
    "eosinophils",
    "basophils",
    "monocytes",
    "neutrophils"
)

predictors <- c(clocks, bloodcells)

biovars <- c(
    "vo2max_per_kg",
    "peak_power_output_per_kg",
    "hand_grip_strength",
    "body_mass_index",
    "fat_percent",
    "lean_mass",
    "bone_density",
    "nightly_heart_rate",
    "nightly_heart_rate_variability",
    "systolic_blood_pressure",  
    "diastolic_blood_pressure",
    "pulse_wave_velocity",
    "sleep_score",
    "sleep_duration",
    "diet_quality_score"
)



###################################
####        Data Summary        ###
###################################

# Calculate mean epigenetic ages & blood cell fractions for replicates
bio <- pheno[, c(
    "subject_id", 
    "exercise_timepoint", 
    predictors
    )] %>%
    group_by(subject_id, exercise_timepoint) %>%
    summarize_at(predictors, mean) %>%
    arrange(subject_id, exercise_timepoint) %>% 
    ungroup()

# Select relevant columns
pheno_subset <- pheno[, c(
        "subject_id",
        "exercise_timepoint", 
        "sex",
        "age",
        "vo2max_per_kg", 
        "peak_power_output_per_kg", 
        "hand_grip_strength",
        "weight",
        "height",
        "body_mass_index",
        "fat_percent",
        "lean_mass",
        "bone_density",
        "nightly_heart_rate", 
        "nightly_heart_rate_variability",
        "systolic_blood_pressure",
        "diastolic_blood_pressure",
        "pulse_wave_velocity",
        "sleep_duration",
        "sleep_score",
        "diet_quality_score",
        "alcohol_use",
        "smoking_status",
        "planned_training",
        "executed_training",
        "training_adherence",
        "dropout",
        "excluded"
        )
    ] %>% unique()

# Merge pheno with bio
bio <- merge(bio, pheno_subset, by = c("subject_id", "exercise_timepoint"))

# Make blood cell fractions sum to 1 for replicates
bio[,bloodcells] <- (bio[,bloodcells] / rowSums(bio[,bloodcells]))*100

# Calculate epigenetic age acceleration
GrimAge_model <- lm(epigenetic_age_GrimAge ~ age, bio)
bio$epigenetic_age_acceleration_GrimAge <-  residuals(GrimAge_model)

Horvath_model <- lm(epigenetic_age_Horvath ~ age, bio)
bio$epigenetic_age_acceleration_Horvath <- residuals(Horvath_model)

save(bio, file = "Data/Objects/bio.Rdata")

# Longitudinal changes in variables
d_bio <- bio %>%
    arrange(subject_id, exercise_timepoint) %>%
    group_by(subject_id) %>%
    summarize(across(where(is.numeric), ~ -diff(.))) %>%
    ungroup()

save(d_bio, file = "Data/Objects/d_bio.Rdata")



###################################
#####  Changes in variables   #####
###################################

load("Data/Objects/pheno.Rdata")
load("Data/Objects/bio.Rdata")
load("Data/Objects/d_bio.Rdata")

# Convert categorical variables to factors
pheno$sex <- factor(pheno$sex)
pheno$exercise_timepoint <- factor(pheno$exercise_timepoint, levels = c("pre", "post"))
pheno$smoking_status <- factor(pheno$smoking_status, levels = c("never", "quit", "sometimes", "daily"), ordered = TRUE)
pheno$alcohol_use <- factor(pheno$alcohol_use, levels = c("almost never", "<1 per week", "1 per week", "2-3 per week", "4-6 per week", "1-2 per day", "3-4 per day", "5+ per day"), ordered = TRUE)
bio$sex <- factor(bio$sex)
bio$exercise_timepoint <- factor(bio$exercise_timepoint, levels = c("pre", "post"))
bio$smoking_status <- factor(bio$smoking_status, levels = c("never", "quit", "sometimes", "daily"), ordered = TRUE)
bio$alcohol_use <- factor(bio$alcohol_use, levels = c("almost never", "<1 per week", "1 per week", "2-3 per week", "4-6 per week", "1-2 per day", "3-4 per day", "5+ per day"), ordered = TRUE)

### Training adherence

# In all participants
print(bio %>%
    filter(exercise_timepoint == "post") %>%
    summarize(mean_adherence = mean(training_adherence, na.rm = TRUE),
            sd_adherence = sd(training_adherence, na.rm = TRUE)) %>%
    ungroup()
# )
#   mean_adherence sd_adherence
# 1       101.8784     34.65668

# In the 33 adhering participants
print(bio %>%
    filter(!dropout, !excluded, exercise_timepoint == "post") %>%
    summarize(mean_adherence = mean(training_adherence, na.rm = TRUE),
            sd_adherence = sd(training_adherence, na.rm = TRUE)) %>%
    ungroup()
)
#   mean_adherence sd_adherence
# 1       108.7395     29.93598

# In the 33 adhering participants, split by sex
print(bio %>%
    filter(!dropout, !excluded, exercise_timepoint == "post") %>%
    group_by(sex) %>%
    summarize(mean_adherence = mean(training_adherence, na.rm = TRUE),
            sd_adherence = sd(training_adherence, na.rm = TRUE)) %>%
    ungroup()
)
#   sex    mean_adherence sd_adherence
#   <fct>           <dbl>        <dbl>
# 1 female           104.         31.3
# 2 male             114.         28.4

# Test difference between men and women
t <- t.test(
    bio %>% filter(!dropout, !excluded, exercise_timepoint == "post", sex == "male") %>% pull(training_adherence),
    bio %>% filter(!dropout, !excluded, exercise_timepoint == "post", sex == "female") %>% pull(training_adherence)
)
print(t)
# t = 1.022, df = 30.968, p-value = 0.3147
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#  -10.57290  31.81047
# sample estimates:
# mean of x mean of y 
#  114.2098  103.5910


### Continuous variables

# Function to test for changes in continuous variables, sex-stratified
analyzeContinuous <- function(data, label = "") {

    results <- data.frame(variable = character(), mean_pre = numeric(), sd_pre = numeric(),
                        N_pre = numeric(), mean_post = numeric(), sd_post = numeric(),
                        N_post = numeric(), perc_change = numeric(), p_value = numeric(),
                        stringsAsFactors = FALSE)

    for (variable in colnames(data)) {
        if (is.numeric(data[[variable]])) {
            pre <- data %>% filter(exercise_timepoint == "pre") %>% pull(variable)
            post <- data %>% filter(exercise_timepoint == "post") %>% pull(variable)
            m_pre <- mean(pre, na.rm = TRUE)
            sd_pre <- sd(pre, na.rm = TRUE)
            N_pre <- sum(!is.na(pre))
            m_post <- mean(post, na.rm = TRUE)
            sd_post <- sd(post, na.rm = TRUE)
            N_post <- sum(!is.na(post))
            perc_change <- (m_post - m_pre) / m_pre * 100

            formula <- as.formula(paste(variable, "~ (1|subject_id) + factor(exercise_timepoint)"))
            model <- lmer(formula, data = data)
            p <- summary(model)$coefficients["factor(exercise_timepoint)post", "Pr(>|t|)"]

            results <- rbind(results, data.frame(
                variable = paste0(variable, label),
                mean_pre = m_pre, sd_pre = sd_pre, N_pre = N_pre,
                mean_post = m_post, sd_post = sd_post, N_post = N_post,
                perc_change = perc_change, p_value = p,
                stringsAsFactors = FALSE
            ))
        }
    }
    return(results)
}

## For 33 adhering participants

# Run analyses
table_data <- bio %>% filter(!dropout, !excluded) %>% select(all_of(c("subject_id", "exercise_timepoint", biovars, "sex")))
results_all <- analyzeContinuous(data = table_data, label = "")
results_M <- analyzeContinuous(data = filter(table_data, sex == "male"), label = "_M")
results_F <- analyzeContinuous(data = filter(table_data, sex == "female"), label = "_F")

# Combine and order
combined_results <- rbind(results_all, results_M, results_F)
ix <- match(paste0(rep(biovars, each = 3), c("", "_M", "_F")), combined_results$variable)
combined_results <- combined_results[ix, ]

# Write to file
write_xlsx(combined_results, "Results/Analysis/Table_Variable_Changes_SexStratified.xlsx")

## For all 38 non-dropout participants

table_data <- bio %>% filter(!dropout) %>% select(all_of(c("subject_id", "exercise_timepoint", biovars, "sex")))
results_all <- analyzeContinuous(data = table_data, label = "")
results_M <- analyzeContinuous(data = filter(table_data, sex == "male"), label = "_M")
results_F <- analyzeContinuous(data = filter(table_data, sex == "female"), label = "_F")
combined_results <- rbind(results_all, results_M, results_F)
ix <- match(paste0(rep(biovars, each = 3), c("", "_M", "_F")), combined_results$variable)
combined_results <- combined_results[ix, ]
write_xlsx(combined_results, "Results/Analysis/Table_Variable_Changes_SexStratified_all_38.xlsx")


### Categorical variables

# Function to test for changes in categorical variables, sex-stratified
analyzeCategorical <- function(data, variable) {
    var_sym <- rlang::sym(variable)
    
    # Keep only non-dropouts with valid values
    df <- data %>%
        group_by(subject_id) %>%
        filter(!any(is.na(!!var_sym))) %>%
        ungroup() %>%
        select(subject_id, sex, exercise_timepoint, !!var_sym) %>%
        arrange(subject_id, exercise_timepoint)
    
    # Function for one subset
    run_subset <- function(df_sub, label) {
        pre <- df_sub %>% filter(exercise_timepoint == "pre") %>% pull(!!var_sym) %>% as.numeric()
        post <- df_sub %>% filter(exercise_timepoint == "post") %>% pull(!!var_sym) %>% as.numeric()
        test <- wilcox.test(pre, post, paired = TRUE)
        
        data.frame(
        variable = paste0(variable, label),
        N = length(pre),
        V = test$statistic,
        p_value = test$p.value,
        direction_neg = sum(post - pre < 0),
        direction_same = sum(post - pre == 0),
        direction_pos = sum(post - pre > 0),
        stringsAsFactors = FALSE
        )
    }
    
    # Stratify by sex
    results_all <- run_subset(df, "")
    results_M   <- run_subset(df %>% filter(sex == "male"), "_M")
    results_F   <- run_subset(df %>% filter(sex == "female"), "_F")

    t <- table(
        table_data %>% pull(!!var_sym), 
        table_data %>% pull(exercise_timepoint)
    )
    t_M <- table(
        table_data %>% filter(sex == "male") %>% pull(!!var_sym), 
        table_data %>% filter(sex == "male") %>% pull(exercise_timepoint)
    )
    rownames(t_M) <- paste(rownames(t_M), "_M", sep = "")
    t_F <- table(
        table_data %>% filter(sex == "female") %>% pull(!!var_sym), 
        table_data %>% filter(sex == "female") %>% pull(exercise_timepoint)
    )
    rownames(t_F) <- paste(rownames(t_F), "_F", sep = "")
    table <- rbind(t, t_M, t_F)

    return(list(stats = rbind(results_all, results_M, results_F), table = table))
}

table_data <- bio %>% filter(!dropout, !excluded) %>% select(all_of(c("subject_id", "exercise_timepoint", "sex", "smoking_status", "alcohol_use")))

# Smoking status

print(analyzeCategorical(table_data, "smoking_status"))
# $stats
#            variable  N   V p_value direction_neg direction_same direction_pos
# V    smoking_status 33 1.5       1             1             31             1
# V1 smoking_status_M 16 0.0       1             0             15             1
# V2 smoking_status_F 17 1.0       1             1             16             0

# $table
#             pre post
# never        21   21
# quit         11   11
# sometimes     0    0
# daily         1    1
# never_M       8    7
# quit_M        7    8
# sometimes_M   0    0
# daily_M       1    1
# never_F      13   14
# quit_F        4    3
# sometimes_F   0    0
# daily_F       0    0

# Alcohol use

print(analyzeCategorical(table_data, "alcohol_use"))
# $stats
#         variable  N    V    p_value direction_neg direction_same direction_pos
# V    alcohol_use 33 32.0 0.04183059             7             25             1
# V1 alcohol_use_M 16 10.0 0.08897301             4             12             0
# V2 alcohol_use_F 17  7.5 0.42371080             3             13             1

# $table
#                pre post
# almost never    10   11
# <1 per week      3    6
# 1 per week      11    9
# 2-3 per week     6    4
# 4-6 per week     2    2
# 1-2 per day      1    1
# 3-4 per day      0    0
# 5+ per day       0    0
# almost never_M   3    4
# <1 per week_M    1    2
# 1 per week_M     6    6
# 2-3 per week_M   5    3
# 4-6 per week_M   1    1
# 1-2 per day_M    0    0
# 3-4 per day_M    0    0
# 5+ per day_M     0    0
# almost never_F   7    7
# <1 per week_F    2    4
# 1 per week_F     5    3
# 2-3 per week_F   1    1
# 4-6 per week_F   1    1
# 1-2 per day_F    1    1
# 3-4 per day_F    0    0
# 5+ per day_F     0    0      

## For all 38 non-dropout participants

table_data <- bio %>% filter(!dropout) %>% select(all_of(c("subject_id", "exercise_timepoint", "sex", "smoking_status", "alcohol_use")))

# Smoking status

print(analyzeCategorical(table_data, "smoking_status"))
# $stats
#            variable  N   V p_value direction_neg direction_same direction_pos
# V    smoking_status 38 1.5       1             1             36             1
# V1 smoking_status_M 18 0.0       1             0             17             1
# V2 smoking_status_F 20 1.0       1             1             19             0

# $table
#             pre post
# never        24   24
# quit         13   13
# sometimes     0    0
# daily         1    1
# never_M       9    8
# quit_M        8    9
# sometimes_M   0    0
# daily_M       1    1
# never_F      15   16
# quit_F        5    4
# sometimes_F   0    0
# daily_F       0    0
# Alcohol use

print(analyzeCategorical(table_data, "alcohol_use"))
# $stats
#         variable  N    V    p_value direction_neg direction_same direction_pos
# V    alcohol_use 38 50.0 0.01465960             9             28             1
# V1 alcohol_use_M 18 10.0 0.08897301             4             14             0
# V2 alcohol_use_F 20 17.5 0.12943067             5             14             1

# $table
#                pre post
# almost never    12   14
# <1 per week      4    7
# 1 per week      12    9
# 2-3 per week     7    5
# 4-6 per week     2    2
# 1-2 per day      1    1
# 3-4 per day      0    0
# 5+ per day       0    0
# almost never_M   4    5
# <1 per week_M    1    2
# 1 per week_M     6    6
# 2-3 per week_M   6    4
# 4-6 per week_M   1    1
# 1-2 per day_M    0    0
# 3-4 per day_M    0    0
# 5+ per day_M     0    0
# almost never_F   8    9
# <1 per week_F    3    5
# 1 per week_F     6    3
# 2-3 per week_F   1    1
# 4-6 per week_F   1    1
# 1-2 per day_F    1    1
# 3-4 per day_F    0    0
# 5+ per day_F     0    0



###################################
#######   Clock statistics  #######
###################################


### Calculate partial marginal variance of EA explained by CA (including replicates)

##  GrimAge clock

# Build the full and reduced mixed models
model <- lmer(epigenetic_age_GrimAge ~ (1|subject_id) + exercise_timepoint + sex + age, data = pheno)
reduced_model <- lmer(epigenetic_age_GrimAge ~ (1|subject_id) + exercise_timepoint + sex, data = pheno)

# Perform a likelihood ratio test to compare the models (p-val for variance explained by age)
anova_result <- anova(model, reduced_model)
p_value <- anova_result$`Pr(>Chisq)`[2]
print(p_value)
# 1.217268e-21

# Calculate partial marginal R² for age (MuMin package)
r2_full <- r.squaredGLMM(model)
print(r2_full)
r2_reduced <- r.squaredGLMM(reduced_model)
print(r2_reduced)
partial_r2_Age <- r2_full[1] - r2_reduced[1]
print(partial_r2_Age)
# 0.8563092

## Horvath clock

# Build the full and reduced mixed models
model <- lmer(epigenetic_age_Horvath ~ (1|subject_id) + exercise_timepoint + sex + age, data = pheno)
reduced_model <- lmer(epigenetic_age_Horvath ~ (1|subject_id) + exercise_timepoint + sex, data = pheno)

# Perform a likelihood ratio test to compare the models (p-val for variance explained by age)
anova_result <- anova(model, reduced_model)
p_value <- anova_result$`Pr(>Chisq)`[2]
print(p_value)
# 1.764362e-13

# Calculate partial marginal R² for age (MuMin package)
r2_full <- r.squaredGLMM(model)
print(r2_full)
r2_reduced <- r.squaredGLMM(reduced_model)
print(r2_reduced)
partial_r2_Age <- r2_full[1] - r2_reduced[1]
print(partial_r2_Age)
# 0.7070314


### Difference in EAA in males vs females, pre-EET

## GrimAge clock
t <- t.test(
    bio %>% filter(sex == "male", exercise_timepoint == "pre") %>% pull(epigenetic_age_acceleration_GrimAge), 
    bio %>% filter(sex == "female", exercise_timepoint == "pre") %>% pull(epigenetic_age_acceleration_GrimAge),
    var.equal = FALSE
)
print(t$p.value)
# 0.0008406405
print(t$estimate[1] - t$estimate[2])
# 2.730406

## Horvath clock
t <- t.test(
    bio %>% filter(sex == "male", exercise_timepoint == "pre") %>% pull(epigenetic_age_acceleration_Horvath), 
    bio %>% filter(sex == "female", exercise_timepoint == "pre") %>% pull(epigenetic_age_acceleration_Horvath),
    var.equal = FALSE
)
print(t$p.value)
# 0.1087038
print(t$estimate[1] - t$estimate[2])
# 1.972114


### Correlation of d_EAA between both clocks
model <- lm(epigenetic_age_acceleration_Horvath ~ epigenetic_age_acceleration_GrimAge, d_bio)
model <- lm(epigenetic_age_acceleration_GrimAge ~ epigenetic_age_acceleration_Horvath, d_bio)
print(summary(model))
# Residual standard error: 1.149 on 36 degrees of freedom
# Multiple R-squared:  0.2388,    Adjusted R-squared:  0.2177 
# F-statistic:  11.3 on 1 and 36 DF,  p-value: 0.00185


### Difference in EAA pre- and post-EET (in months)

# Function for testing the difference in EAA, adjusted and unadjusted for leukocyte composition, stratified by sex
EAA_change <- function(data, clocks) {

    sexes <- list(All = data, M = filter(data, sex == "male"), F = filter(data, sex == "female"))

    run_model <- function(data, var, covar = NULL) {
        f <- as.formula(paste0(var, "*12 ~ (1|subject_id) + exercise_timepoint", if (!is.null(covar)) paste0(" + ", covar) else ""))
        m <- lmer(f, data = data)
        res <- list(
            estimate = summary(m)$coefficients["exercise_timepointpost", "Estimate"],
            p = summary(m)$coefficients["exercise_timepointpost", "Pr(>|t|)"]
        )
        return(res)
    }
    
    results <- do.call(rbind, lapply(names(clocks), function(clock) {
            covar <- clocks[[clock]]
            var <- paste0("epigenetic_age_acceleration_", clock)
            do.call(rbind, lapply(names(sexes), function(s) {
                d <- sexes[[s]]
                rbind(
                    cbind(clock, adj = "unadj", sex = s, t(run_model(d, var, NULL))),
                    cbind(clock, adj = "adj",   sex = s, t(run_model(d, var, covar)))
                )
            }))
        }))

    results <- as.data.frame(results) %>% arrange(clock, adj, sex)

    return(as.data.frame(results))
}

results <- EAA_change(bio %>% filter(!dropout, !excluded), clocks = list(GrimAge = "neutrophils", Horvath = "basophils"))
print(results)
#      clock   adj sex  estimate            p
# 1  GrimAge unadj All -7.446353   0.01210277
# 2  GrimAge unadj   M -5.140452    0.2383759
# 3  GrimAge unadj   F -9.616611   0.02190524
# 4  GrimAge   adj All -5.943747 5.890311e-05
# 5  GrimAge   adj   M -7.451078  0.002297967
# 6  GrimAge   adj   F -4.949702   0.01090589
# 7  Horvath unadj All  -5.52984    0.0230796
# 8  Horvath unadj   M -2.229527    0.4570118
# 9  Horvath unadj   F -8.636017   0.02437446
# 10 Horvath   adj All -5.057068  0.002673217
# 11 Horvath   adj   M -5.480436    0.0122493
# 12 Horvath   adj   F -4.472946    0.1339921



################################################
####  Baseline Associations EAA ~ variable  ####
################################################


# Function for testing the associations between epigenetic age acceleration and other variables at baseline (pre-EET), while adjusting for sex and with or without adjusting for other variables
crosssectionalAssociationsTest <- function(data, outcome, predictors, adjust = NULL) {
    # Scale numeric columns
    numeric_columns <- sapply(data, is.numeric)
    data_scaled <- scale(data[, numeric_columns])
    data[,numeric_columns] <- data_scaled
    data_scaled <- data

    # Helper function to test a single predictor without adjust variable
    modelTest <- function(data_scaled, outcome, predictor) {
        # Unadjusted model
        unadjusted_formula <- as.formula(paste0(outcome, " ~ ", predictor, " + sex")) # always adjust for sex when looking cross-sectionally
        unadjusted_model <- lm(unadjusted_formula, data = as.data.frame(data_scaled), na.action = na.omit)
        unadjusted_summary <- summary(unadjusted_model)

        # Return results
        list(
            Predictor = predictor,
            Coefficient = unadjusted_summary$coefficients[2, 1],
            P_value = unadjusted_summary$coefficients[2, 4]
        )
    }

    # Helper function to test a single predictor with adjust variable
    adjustedModelTest <- function(data_scaled, outcome, predictor, adjust = NULL) {
        # Unadjusted model
        unadjusted_formula <- as.formula(paste0(outcome, " ~ ", predictor, " + sex")) # always adjust for sex when looking cross-sectionally
        unadjusted_model <- lm(unadjusted_formula, data = as.data.frame(data_scaled), na.action = na.omit)
        unadjusted_summary <- summary(unadjusted_model)

        # Adjusted model
        adjust <- paste(adjust, collapse = " + ")
        adjusted_formula <- as.formula(paste0(outcome, " ~ ", predictor, " + sex + ", adjust))
        adjusted_model <- lm(adjusted_formula, data = as.data.frame(data_scaled), na.action = na.omit)
        adjusted_summary <- summary(adjusted_model)

        # Return results
        if (is.numeric(data_scaled[, predictor])) {
            # return results for continuous predictors
            return(list(
                Predictor = predictor,
                Coefficient_Unadjusted = unadjusted_summary$coefficients[2, 1],
                P_value_Unadjusted = unadjusted_summary$coefficients[2, 4],
                Coefficient_Adjusted = adjusted_summary$coefficients[2, 1],
                P_value_Adjusted = adjusted_summary$coefficients[2, 4]
            ))
        } else {
            # return results for categorical predictors
            predictor_levels <- rownames(unadjusted_summary$coefficients)[grepl(predictor, rownames(unadjusted_summary$coefficients))]
            return(list(
                Predictor = predictor_levels,
                Coefficient_Unadjusted = unadjusted_summary$coefficients[predictor_levels, 1],
                P_value_Unadjusted = unadjusted_summary$coefficients[predictor_levels, 4],
                Coefficient_Adjusted = adjusted_summary$coefficients[predictor_levels, 1],
                P_value_Adjusted = adjusted_summary$coefficients[predictor_levels, 4]
            ))
        }
    }

    # Iterate over predictors and apply either modelTest or adjustedModelTest
    if (is.null(adjust)) {
        results <- lapply(predictors, function(predictor) {
            modelTest(data_scaled, outcome, predictor)
        })
    } else {
        results <- lapply(predictors, function(predictor) {
            adjustedModelTest(data_scaled, outcome, predictor, adjust)
        })
    }

    # Combine results into a data frame
    results_df <- do.call(rbind, lapply(results, as.data.frame))
    rownames(results_df) <- NULL
    return(results_df)
}


### Associations between EAA and leukocytes at baseline

# GrimAge clock
GrimAge_pre_BC_associations <- crosssectionalAssociationsTest(
    data = bio %>% filter(exercise_timepoint == "pre"),
    outcome = "epigenetic_age_acceleration_GrimAge", 
    predictors = bloodcells
)
write_xlsx(GrimAge_pre_BC_associations, "Results/Analysis/Baseline_BloodCell_Associations_GrimAge.xlsx")

# Horvath clock
Horvath_pre_BC_associations <- crosssectionalAssociationsTest(
    data = bio %>% filter(exercise_timepoint == "pre"),
    outcome = "epigenetic_age_acceleration_Horvath", 
    predictors = bloodcells
)
write_xlsx(Horvath_pre_BC_associations, "Results/Analysis/Baseline_BloodCell_Associations_Horvath.xlsx")


### Associations between EAA and other variables, unadjusted and adjusted for leukocytes

# Convert ordinal alcohol use to numeric for analysis
bio$alcohol_use <- factor(bio$alcohol_use, levels = c("almost never", "<1 per week", "1 per week", "2-3 per week", "4-6 per week", "1-2 per day", "3-4 per day", "5+ per day"), ordered = TRUE) %>% as.numeric()

# Convert ordinal smoking status to unordered factor for analysis
bio$smoking_status <- factor(bio$smoking_status, levels = c("never", "quit", "sometimes", "daily"), ordered = FALSE)

# GrimAge clock
GrimAge_pre_associations <- crosssectionalAssociationsTest(
    data = bio %>% filter(exercise_timepoint == "pre"),
    outcome = "epigenetic_age_acceleration_GrimAge", 
    predictors = c(biovars, "alcohol_use", "smoking_status"),
    adjust = "neutrophils"
)
write_xlsx(GrimAge_pre_associations, "Results/Analysis/Baseline_Associations_GrimAge.xlsx")

# Horvath clock
Horvath_pre_associations <- crosssectionalAssociationsTest(
    data = bio %>% filter(exercise_timepoint == "pre"),
    outcome = "epigenetic_age_acceleration_Horvath", 
    predictors = c(biovars, "alcohol_use", "smoking_status"),
    adjust = "basophils"
)
write_xlsx(Horvath_pre_associations, "Results/Analysis/Baseline_Associations_Horvath.xlsx")


########################################################
####  Longitudinal Associations d_EAA ~ d_variable  ####
########################################################

# Function to test longitudinal associations between changes in epigenetic age acceleration and changes in other variables, with or without adjusting for other variables)
longitudinalAssociationsTest <- function(data, outcome, predictors, adjust = NULL) {
    
    # Scale numeric columns
    numeric_columns <- sapply(data, is.numeric)
    data_scaled <- scale(data[, numeric_columns])
    data[,numeric_columns] <- data_scaled
    data_scaled <- data

    # Helper function to test a single predictor without adjust variable
    modelTest <- function(data_scaled, outcome, predictor) {
        # Unadjusted model
        unadjusted_formula <- as.formula(paste0(outcome, " ~ ", predictor))
        unadjusted_model <- lm(unadjusted_formula, data = as.data.frame(data_scaled), na.action = na.omit)
        unadjusted_summary <- summary(unadjusted_model)

        # Return results
        list(
            Predictor = predictor,
            Coefficient = unadjusted_summary$coefficients[2, 1],
            P_value = unadjusted_summary$coefficients[2, 4],
            R_squared = unadjusted_summary$r.squared
        )
    }

    # Helper function to test a single predictor with adjust variable
    adjustedModelTest <- function(data_scaled, outcome, predictor, adjust) {
        # Unadjusted model
        unadjusted_formula <- as.formula(paste0(outcome, "  ~ ", predictor))
        unadjusted_model <- lm(unadjusted_formula, data = as.data.frame(data_scaled), na.action = na.omit)
        unadjusted_summary <- summary(unadjusted_model)

        # Adjusted model
        adjust <- paste(adjust, collapse = " + ")
        adjusted_formula <- as.formula(paste0(outcome, " ~ ", predictor, " + ", adjust))
        adjusted_model <- lm(adjusted_formula, data = as.data.frame(data_scaled), na.action = na.omit)
        adjusted_summary <- summary(adjusted_model)

        # Return results
        list(
            Predictor = predictor,
            Coefficient_Unadjusted = unadjusted_summary$coefficients[2, 1],
            P_value_Unadjusted = unadjusted_summary$coefficients[2, 4],
            R_squared_Unadjusted = unadjusted_summary$r.squared,
            Coefficient_Adjusted = adjusted_summary$coefficients[2, 1],
            P_value_Adjusted = adjusted_summary$coefficients[2, 4],
            R_squared_Adjusted = adjusted_summary$r.squared
        )
    }

    # Iterate over predictors and apply either modelTest or adjustedModelTest
    if (is.null(adjust)) {
        results <- lapply(predictors, function(predictor) {
            modelTest(data_scaled, outcome, predictor)
        })
    } else {
        results <- lapply(predictors, function(predictor) {
            adjustedModelTest(data_scaled, outcome, predictor, adjust)
        })
    }

    # Combine results into a data frame
    results_df <- do.call(rbind, lapply(results, as.data.frame))
    rownames(results_df) <- NULL
    return(results_df)
}


### Associations between EAA and leukocytes

# GrimAge clock
GrimAge_long_BC_associations <- longitudinalAssociationsTest(
    data = d_bio,
    outcome = "epigenetic_age_acceleration_GrimAge", 
    predictors = bloodcells
)
write_xlsx(GrimAge_long_BC_associations, "Results/Analysis/Longitudinal_BloodCell_Associations_GrimAge.xlsx")

# Horvath clock
Horvath_long_BC_associations <- longitudinalAssociationsTest(
    data = d_bio,
    outcome = "epigenetic_age_acceleration_Horvath", 
    predictors = bloodcells
)
write_xlsx(Horvath_long_BC_associations, "Results/Analysis/Longitudinal_BloodCell_Associations_Horvath.xlsx")


### Associations between EAA and other variables, unadjusted and adjusted for a single leukocyte fraction

# GrimAge clock
GrimAge_long_associations <- longitudinalAssociationsTest(
    data = d_bio,
    outcome = "epigenetic_age_acceleration_GrimAge", 
    predictors = c(biovars, "executed_training", "training_adherence"),
    adjust = "neutrophils"
)
write_xlsx(GrimAge_long_associations, "Results/Analysis/Longitudinal_Associations_GrimAge.xlsx")

# Horvath clock
Horvath_long_associations <- longitudinalAssociationsTest(
    data = d_bio,
    outcome = "epigenetic_age_acceleration_Horvath", 
    predictors = c(biovars, "executed_training", "training_adherence"),
    adjust = "basophils"
)
write_xlsx(Horvath_long_associations, "Results/Analysis/Longitudinal_Associations_Horvath.xlsx")


### Associations between EAA and other variables, unadjusted and adjusted for all leukocyte fractions

# GrimAge clock
GrimAge_long_associations_allBC <- longitudinalAssociationsTest(
    data = d_bio,
    outcome = "epigenetic_age_acceleration_GrimAge", 
    predictors = c(biovars, "executed_training", "training_adherence"),
    adjust = bloodcells[bloodcells != "eosinophils"] # adjust for all except eosinophils
)
write_xlsx(GrimAge_long_associations_allBC, "Results/Analysis/Longitudinal_Associations_GrimAge_AllBC.xlsx")

# Horvath clock
Horvath_long_associations_allBC <- longitudinalAssociationsTest(
    data = d_bio,
    outcome = "epigenetic_age_acceleration_Horvath", 
    predictors = c(biovars, "executed_training", "training_adherence"),
    adjust = bloodcells[bloodcells != "eosinophils"] # adjust for all except eosinophils
)
write_xlsx(Horvath_long_associations_allBC, "Results/Analysis/Longitudinal_Associations_Horvath_AllBC.xlsx")