BZea Flowering Time Analysis

Load Libraries

Purpose: Load all required packages for data manipulation, visualization, and analysis.

Approach: Load tidyverse ecosystem packages plus specialized plotting libraries.

Expected outcome: All dependencies available for downstream analysis.

library(dplyr)
library(lubridate)
library(ggplot2)
library(ggpubr)
library(tidyr)
library(ggbeeswarm)
library(forcats)

Data Import

Purpose: Import flowering time data and population metadata.

Approach: Read CSV files containing days to silking/anthesis measurements and line metadata.

Expected outcome: Two dataframes ready for cleaning and merging.

data <- read.csv("BZea DTS_DTA  - Sheet1.csv")
meta <- read.csv("B73xTEO-Pops-Metadata.csv")

# Diagnostics
{
  cat("Data dimensions:", nrow(data), "rows x", ncol(data), "cols\n")
  cat("Metadata dimensions:", nrow(meta), "rows x", ncol(meta), "cols\n")
  cat("Data columns:", paste(colnames(data), collapse = ", "), "\n")
}

## Data dimensions: 4432 rows x 6 cols
## Metadata dimensions: 175 rows x 8 cols
## Data columns: CLY23.D4, Description, Female.genotype, DTS, DTA, Planting.Date

Data Cleaning and Transformation

Purpose: Rename columns, calculate flowering time metrics, and extract donor population information.

Approach: Rename measurement columns to standard names, convert dates to days after planting, parse genotype strings to extract donor population and line identifiers using regex.

Expected outcome: Clean dataframe with DTA (days to anthesis), DTS (days to silking), ASI (anthesis-silking interval), and donor identifiers.

# Rename and calculate date differences
data <- data %>%
  rename(
    silking = "DTS",
    anthesis = "DTA",
    Replicate = "Description"
  ) %>%
  mutate(
    DTS = mdy(silking) - mdy(Planting.Date),
    DTA = mdy(anthesis) - mdy(Planting.Date),
    donor_pop = gsub("-.*", "", Female.genotype, perl = TRUE),
    donor_line = gsub("_P.*?$", "", Female.genotype, perl = TRUE)
  )

# Clean donor_line strings
data$donor_line <- gsub("^.*?-", "", data$donor_line, perl = TRUE)
data$donor_line <- gsub("-+", "_", data$donor_line, perl = TRUE)
data$donor_line <- gsub("\\(.*", "", data$donor_line, perl = TRUE)
data$donor_line <- gsub("-", "_", data$donor_line, perl = TRUE)
data$donor_line <- gsub(" +$", "", data$donor_line, perl = TRUE)
data$donor_line <- factor(data$donor_line)

# Convert to numeric and calculate ASI, clean replicate names, filter outliers
data <- data %>%
  mutate(
    DTA = as.numeric(DTA),
    DTS = as.numeric(DTS),
    ASI = as.numeric(DTS - DTA),
    Replicate = gsub("BZea-Rep", "", Replicate)
  ) %>%
  filter(DTA < 140)  # Remove typos/outliers

# Diagnostics
{
  cat("After cleaning:", nrow(data), "observations\n")
  cat("DTA range:", range(data$DTA, na.rm = TRUE), "\n")
  cat("DTS range:", range(data$DTS, na.rm = TRUE), "\n")
  cat("Replicates:", paste(unique(data$Replicate), collapse = ", "), "\n")
  }

## After cleaning: 4335 observations
## DTA range: 78 100 
## DTS range: 78 100 
## Replicates: 1, 2, 3

Species Classification

Purpose: Assign species labels based on donor population and filter out controls.

Approach: Map donor population codes to Zea species names, remove Purple Check controls.

Expected outcome: Data with species classification for downstream analysis and visualization.

data <- data %>%
  filter(donor_pop != "Purple Check") %>%
  droplevels() %>%
  mutate(
    spp = case_when(
      donor_pop %in% c("Dura", "Nabo", "Mesa", "Chal") ~ "mexicana",
      donor_pop == "B73" ~ "mays",
      donor_pop == "Bals" ~ "parviglumis",
      donor_pop == "Hueh" ~ "huehuetenagensis",
      donor_pop == "Zdip" ~ "diploperennis",
      donor_pop == "Zlux" ~ "luxurians"
    )
  )

# Order donor_pop by mean DTA
data$donor_pop <- fct_reorder(data$donor_pop, data$DTA, .fun = mean, .na_rm = TRUE)

# Diagnostics
{
  cat("Species counts:\n")
  print(table(data$spp))
  cat("\nDonor populations:", nlevels(data$donor_pop), "\n")
}

## Species counts:
## 
##    diploperennis huehuetenagensis        luxurians             mays 
##              712              160              322              392 
##         mexicana      parviglumis 
##             1681             1067 
## 
## Donor populations: 9

Merge Metadata

Purpose: Join population metadata (coordinates, elevation) to flowering data.

Approach: Left join metadata by donor line ID, reorder columns to place metadata first.

Expected outcome: Complete dataset with geographic and flowering information.

data <- data %>%
  left_join(meta, by = c(donor_line = "Line.ID")) %>%
  dplyr::select(all_of(colnames(meta)[-1]), everything())

# Diagnostics
{
  cat("Final dataset:", nrow(data), "rows x", ncol(data), "cols\n")
  cat("Lines with metadata:", sum(!is.na(data$Latitude)), "\n")
  cat("Replicate distribution:\n")
  print(table(data$Replicate))
  }

## Final dataset: 4334 rows x 19 cols
## Lines with metadata: 3900 
## Replicate distribution:
## 
##    1    2    3 
## 1444 1438 1452

Replicate Effects Visualization

Purpose: Visualize flowering time distributions across replicates to assess experimental consistency.

Approach: Create beeswarm plots showing individual observations with mean ± CI overlays for each replicate.

Expected outcome: Visual assessment of replicate-to-replicate variation in DTA, DTS, and ASI.

data %>%
  filter(DTA < 140) %>%
  ggplot(aes(y = DTA, x = Replicate, group = Replicate, col = Replicate)) +
  ggtitle("Days to Anthesis by Replicate") +
  geom_quasirandom() +
  stat_summary(fun.data = mean_cl_normal, geom = "pointrange", col = "black") +
  theme_classic2(base_size = 15) +
  theme(legend.position = "none")

Days to anthesis by replicate

# Helper function for y-axis breaks
equal_breaks <- function(n = 10, s = 1, ...) {
  function(x) {
    seq <- seq(min(x), max(x), s)
    round(seq)
  }
}

# Pivot data for faceted plotting
long_data <- data %>%
  filter(DTA < 140, donor_pop != "Purple Check") %>%
  droplevels() %>%
  pivot_longer(
    cols = c("DTA", "DTS", "ASI"),
    names_to = "var",
    values_to = "Days"
  ) %>%
  group_by(var) %>%
  mutate(mean = mean(Days, na.rm = TRUE)) %>%
  ungroup()

# Set logical variable order
long_data$var <- factor(long_data$var, levels = c("DTA", "DTS", "ASI"))

long_data %>%
  ggplot(aes(y = Days, x = Replicate, group = Replicate, col = Replicate)) +
  scale_y_continuous(breaks = equal_breaks(20)) +
  geom_quasirandom() +
  facet_wrap(~var, scales = "free_y") +
  geom_hline(aes(yintercept = mean, group = var), colour = "grey", linetype = 2) +
  stat_summary(fun.data = mean_cl_normal, geom = "pointrange", col = "black") +
  theme_classic2(base_size = 15) +
  theme(
    legend.position = "none",
    strip.background = element_rect(colour = "white", fill = "white")
  )

Flowering time variables by replicate

Linear Models: Replicate Effects

Purpose: Quantify replicate effects on flowering time using linear regression.

Approach: Fit linear models with DTA and DTS as response variables, replicate and donor population as predictors.

Expected outcome: Estimates of replicate effects in days, significance tests.

# DTA model
dta_mod <- lm(DTA ~ Replicate + donor_pop, data = data)

# DTS model
dts_mod <- lm(DTS ~ Replicate + donor_pop, data = data)

# Diagnostics
{
  cat("=== DTA Model (Days to Anthesis) ===\n")
  print(summary(dta_mod))
  cat("\n=== DTS Model (Days to Silking) ===\n")
  print(summary(dts_mod))
}

## === DTA Model (Days to Anthesis) ===
## 
## Call:
## lm(formula = DTA ~ Replicate + donor_pop, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3921 -1.9548 -0.1374  1.5707 12.9716 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    85.5539     0.1454 588.430  < 2e-16 ***
## Replicate2      0.5627     0.1005   5.601 2.26e-08 ***
## Replicate3      0.4009     0.1002   4.000 6.43e-05 ***
## donor_popNabo   0.2518     0.1896   1.328  0.18427    
## donor_popMesa   0.4745     0.1775   2.672  0.00756 ** 
## donor_popChal   0.8381     0.1974   4.246 2.23e-05 ***
## donor_popB73    1.0128     0.1906   5.314 1.12e-07 ***
## donor_popBals   1.4246     0.1568   9.085  < 2e-16 ***
## donor_popZdip   1.5835     0.1673   9.465  < 2e-16 ***
## donor_popHueh   2.4352     0.2514   9.685  < 2e-16 ***
## donor_popZlux   2.4743     0.2009  12.317  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.696 on 4323 degrees of freedom
## Multiple R-squared:  0.0717, Adjusted R-squared:  0.06955 
## F-statistic: 33.39 on 10 and 4323 DF,  p-value: < 2.2e-16
## 
## 
## === DTS Model (Days to Silking) ===
## 
## Call:
## lm(formula = DTS ~ Replicate + donor_pop, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.4062 -1.9630 -0.0509  1.6647 12.3974 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    86.4730     0.1559 554.793  < 2e-16 ***
## Replicate2      0.4901     0.1078   4.547 5.59e-06 ***
## Replicate3      0.1964     0.1076   1.826  0.06788 .  
## donor_popNabo   0.5401     0.2034   2.656  0.00794 ** 
## donor_popMesa   0.4838     0.1904   2.542  0.01107 *  
## donor_popChal   0.9332     0.2116   4.411 1.06e-05 ***
## donor_popB73    2.0813     0.2043  10.189  < 2e-16 ***
## donor_popBals   1.5779     0.1682   9.381  < 2e-16 ***
## donor_popZdip   1.6659     0.1794   9.286  < 2e-16 ***
## donor_popHueh   2.5233     0.2701   9.342  < 2e-16 ***
## donor_popZlux   2.5989     0.2153  12.071  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.89 on 4312 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.06946,    Adjusted R-squared:  0.06731 
## F-statistic: 32.19 on 10 and 4312 DF,  p-value: < 2.2e-16

Combined Flowering Time Distribution

Purpose: Visualize overall distributions of flowering time variables.

Approach: Beeswarm plots with all replicates combined, showing global means.

Expected outcome: Publication-ready figure summarizing flowering time phenotypes.

p_combined <- long_data %>%
  ggplot(aes(y = Days, x = 1, group = var, col = var)) +
  scale_y_continuous(breaks = equal_breaks(20)) +
  geom_quasirandom() +
  facet_wrap(~var, scales = "free_y") +
  xlab("") +
  geom_hline(aes(yintercept = mean, group = var), colour = "grey", linetype = 2) +
  stat_summary(fun.data = mean_cl_normal, geom = "pointrange", col = "black") +
  theme_classic2(base_size = 15) +
  theme(
    legend.position = "none",
    strip.background = element_rect(colour = "white", fill = "white"),
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  )

print(p_combined)

Combined flowering time distributions

# Save to PDF (run interactively)
pdf(file = "flowering_plot.pdf")
print(p_combined)
dev.off()

Flowering Time by Donor Population

Purpose: Compare flowering time across donor populations and species.

Approach: Faceted beeswarm plots colored by species with population-level summaries.

Expected outcome: Visual comparison of flowering time variation among teosinte populations.

long_data %>%
  droplevels() %>%
  ggplot(aes(y = Days, x = donor_pop, group = donor_pop, col = spp)) +
  scale_y_continuous(breaks = equal_breaks(20)) +
  geom_quasirandom() +
  geom_hline(aes(yintercept = mean, group = var), colour = "grey", linetype = 2) +
  facet_wrap(~var, scales = "free_y") +
  stat_summary(fun.data = mean_cl_normal, geom = "pointrange", col = "black") +
  theme_classic2(base_size = 15) +
  theme(
    legend.position = "top",
    strip.background = element_rect(colour = "white", fill = "white"),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

Flowering time by donor population

Elevation and Latitude Analysis

Purpose: Examine relationship between collection site geography and flowering time.

Approach: Correlation analysis and linear regression of elevation vs latitude, then flowering time vs geographic variables.

Expected outcome: Quantified relationships between geography and phenotype.

# Extract unique line-level geographic data
line_data <- data %>%
  filter(donor_pop != "B73") %>%
  dplyr::select(donor_pop, donor_line, Latitude, Altitude) %>%
  distinct() %>%
  droplevels()

r_lat_alt <- cor(line_data$Latitude, line_data$Altitude, use = "pairwise.complete.obs")

line_data %>%
  ggplot(aes(x = Latitude, y = Altitude)) +
  geom_point(aes(col = donor_pop), alpha = 1) +
  geom_smooth(method = "lm") +
  stat_regline_equation(label.y = 2500, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 2400, aes(label = ..rr.label..)) +
  theme_classic2(base_size = 20) +
  theme(legend.position = "top")

Latitude vs Altitude correlation

# Diagnostics
{
  cat("Latitude-Altitude correlation: r =", round(r_lat_alt, 3), "\n")
  cat("R-squared:", round(r_lat_alt^2, 3), "\n")
  }

## Latitude-Altitude correlation: r = 0.388 
## R-squared: 0.151

# Full latitude range
geo_mod_full <- lm(Altitude ~ Latitude, data = line_data)

# Filtered latitude range (16-21°)
lat_filtered <- line_data %>%
  filter(Latitude > 16 & Latitude < 21)

r_filtered <- cor(lat_filtered$Latitude, lat_filtered$Altitude, use = "pairwise.complete.obs")
geo_mod_filtered <- lm(Altitude ~ Latitude, data = lat_filtered)

# Diagnostics
{
  cat("=== Full Latitude Range ===\n")
  cat("Slope:", round(coef(geo_mod_full)[2], 1), "m/degree\n")
  cat("\n=== Filtered Latitude (16-21°) ===\n")
  cat("Correlation: r =", round(r_filtered, 3), "\n")
  cat("Slope:", round(coef(geo_mod_filtered)[2], 1), "m/degree\n")
}

## === Full Latitude Range ===
## Slope: 107.3 m/degree
## 
## === Filtered Latitude (16-21°) ===
## Correlation: r = 0.421 
## Slope: 283.4 m/degree

Flowering Time vs Geography

Purpose: Quantify relationship between collection site geography and flowering time.

Approach: Scatter plots with regression lines for DTA vs elevation and latitude.

Expected outcome: Effect sizes for geographic clines in flowering time.

plot_ele <- data %>%
  filter(donor_pop != "B73") %>%
  droplevels() %>%
  ggplot(aes(x = Altitude, y = DTA)) +
  geom_point(aes(col = donor_pop), alpha = 1) +
  geom_smooth(method = "lm") +
  stat_regline_equation(label.y = 102, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 101, aes(label = ..rr.label..)) +
  theme_classic2(base_size = 20)

plot_lat <- data %>%
  filter(donor_pop != "B73") %>%
  droplevels() %>%
  ggplot(aes(x = Latitude, y = DTA)) +
  geom_point(aes(col = donor_pop)) +
  geom_smooth(method = "lm") +
  stat_regline_equation(label.y = 102, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 101, aes(label = ..rr.label..)) +
  theme_classic2(base_size = 20)

ggarrange(plot_ele, plot_lat, common.legend = TRUE, legend = "top")

Flowering time vs geographic variables

# Altitude effect on DTA
mod_alt <- lm(DTA ~ Replicate + Altitude, data = data)

# Latitude effect on DTA
mod_lat <- lm(DTA ~ Replicate + Latitude, data = data)

# Donor population effect on DTS
mod_pop <- lm(DTS ~ Replicate + donor_pop, data = data)

# Diagnostics
{
  cat("=== Altitude Effect ===\n")
  cat("Effect:", round(coef(mod_alt)["Altitude"] * 1000, 2), "days/km elevation\n\n")
  cat("=== Latitude Effect ===\n")
  cat("Effect:", round(coef(mod_lat)["Latitude"], 2), "days/degree latitude\n\n")
  cat("=== Donor Population Effect (DTS) ===\n")
  print(summary(mod_pop)$coefficients[grep("donor_pop", rownames(summary(mod_pop)$coefficients)), ])
}

## === Altitude Effect ===
## Effect: -1.01 days/km elevation
## 
## === Latitude Effect ===
## Effect: -0.19 days/degree latitude
## 
## === Donor Population Effect (DTS) ===
##                Estimate Std. Error   t value     Pr(>|t|)
## donor_popNabo 0.5400888  0.2033541  2.655903 7.938752e-03
## donor_popMesa 0.4838153  0.1903646  2.541519 1.107191e-02
## donor_popChal 0.9332418  0.2115907  4.410599 1.055981e-05
## donor_popB73  2.0812898  0.2042711 10.188860 4.161729e-24
## donor_popBals 1.5778995  0.1682033  9.380907 1.028903e-20
## donor_popZdip 1.6659062  0.1793986  9.286059 2.474319e-20
## donor_popHueh 2.5232642  0.2700925  9.342223 1.473116e-20
## donor_popZlux 2.5989422  0.2153102 12.070687 5.111890e-33

Variance Components and Heritability

Purpose: Estimate broad-sense heritability (H²) for days to anthesis.

Approach: ANOVA decomposition of variance among genotypes vs residual variance.

Expected outcome: H² estimate representing proportion of phenotypic variance explained by genotype.

# Full ANOVA model
anova_full <- anova(lm(DTA ~ Latitude + Altitude + Replicate + donor_pop + donor_line + Female.genotype, data = data))

# H2 calculation from genotype model
dta_h2_mod <- lm(DTA ~ Female.genotype, data = data)
h2_anova <- anova(dta_h2_mod)
H2 <- h2_anova$`Sum Sq`[1] / sum(h2_anova$`Sum Sq`)

# Diagnostics
{
  cat("=== Full ANOVA ===\n")
  print(anova_full)
  cat("\n=== Broad-sense Heritability ===\n")
  cat("H² =", round(H2 * 100, 1), "%\n")
}

## === Full ANOVA ===
## Analysis of Variance Table
## 
## Response: DTA
##                   Df  Sum Sq Mean Sq  F value    Pr(>F)    
## Latitude           1  1397.6 1397.60 337.1574 < 2.2e-16 ***
## Altitude           1   385.9  385.90  93.0949 < 2.2e-16 ***
## Replicate          2   190.7   95.35  23.0015 1.254e-10 ***
## donor_pop          7   629.3   89.90  21.6884 < 2.2e-16 ***
## donor_line        70  1415.3   20.22   4.8776 < 2.2e-16 ***
## Female.genotype 1228 17183.0   13.99   3.3756 < 2.2e-16 ***
## Residuals       2590 10736.2    4.15                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## === Broad-sense Heritability ===
## H² = 62.3 %

Flowering Time by Donor Elevation

Purpose: Visualize flowering time ordered by collection site elevation.

Approach: Beeswarm plot with lines ordered by altitude, labeled with population and elevation.

Expected outcome: Visual pattern of elevation cline in flowering time.

# Reorder by elevation
data$donor_line <- factor(data$donor_line)
data$donor_line <- fct_reorder(data$donor_line, data$Altitude, .fun = mean, .na_rm = TRUE)

# Create labels
label_data_elev <- data %>%
  dplyr::select(donor_line, donor_pop, spp, Altitude) %>%
  unique() %>%
  mutate(pop_alt = paste(donor_pop, Altitude))

label_data_elev$donor_line <- factor(label_data_elev$donor_line, levels = levels(data$donor_line))
label_data_elev <- label_data_elev %>% arrange(donor_line)
labels_elev <- label_data_elev$pop_alt

data %>%
  ggplot(aes(y = DTA, x = donor_line, col = spp, group = donor_line)) +
  xlab("Donor Line by Elevation") +
  scale_x_discrete(labels = labels_elev) +
  scale_y_continuous(breaks = equal_breaks(20)) +
  geom_quasirandom(dodge.width = 1) +
  geom_hline(yintercept = mean(data$DTA, na.rm = TRUE), colour = "black", linetype = 2) +
  stat_summary(
    fun.data = mean_cl_normal,
    width = 0, size = 0.1,
    position = position_dodge(1),
    show.legend = FALSE,
    geom = "pointrange",
    col = "black"
  ) +
  theme_classic2(base_size = 15) +
  theme(
    legend.position = "top",
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

DTA by donor line ordered by elevation

Flowering Time by Donor Latitude

Purpose: Visualize flowering time ordered by collection site latitude.

Approach: Beeswarm plot with lines ordered by latitude, labeled with population and coordinates.

Expected outcome: Visual pattern of latitudinal cline in flowering time.

# Reorder by latitude
data$donor_line <- factor(data$donor_line)
data$donor_line <- fct_reorder(data$donor_line, data$Latitude, .fun = mean, .na_rm = TRUE)

# Create labels with formatted latitude
label_data_lat <- data %>%
  dplyr::select(donor_line, donor_pop, spp, Latitude) %>%
  unique() %>%
  mutate(
    lat_str = formatC(as.numeric(Latitude), format = "f", flag = "0", digits = 3),
    pop_lat = paste(donor_pop, lat_str)
  )

label_data_lat$donor_line <- factor(label_data_lat$donor_line, levels = levels(data$donor_line))
label_data_lat <- label_data_lat %>% arrange(donor_line)
labels_lat <- label_data_lat$pop_lat

data %>%
  ggplot(aes(y = DTA, x = donor_line, col = spp, group = donor_line)) +
  xlab("Donor Line by Latitude") +
  scale_x_discrete(labels = labels_lat) +
  scale_y_continuous(breaks = equal_breaks(20)) +
  geom_quasirandom(dodge.width = 1) +
  geom_hline(yintercept = mean(data$DTA, na.rm = TRUE), colour = "black", linetype = 2) +
  stat_summary(
    fun.data = mean_cl_normal,
    width = 0, size = 0.1,
    position = position_dodge(1),
    show.legend = FALSE,
    geom = "pointrange",
    col = "black"
  ) +
  theme_classic2(base_size = 15) +
  theme(
    legend.position = "top",
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

DTA by donor line ordered by latitude

Session Info

sessionInfo()

## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.2
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] forcats_1.0.1    ggbeeswarm_0.7.3 tidyr_1.3.2      ggpubr_0.6.2    
## [5] ggplot2_4.0.1    lubridate_1.9.4  dplyr_1.1.4     
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       beeswarm_0.4.0     xfun_0.55          bslib_0.9.0       
##  [5] htmlwidgets_1.6.4  rstatix_0.7.3      lattice_0.22-7     vctrs_0.6.5       
##  [9] tools_4.5.2        generics_0.1.4     tibble_3.3.0       cluster_2.1.8.1   
## [13] pkgconfig_2.0.3    Matrix_1.7-4       data.table_1.18.0  checkmate_2.3.3   
## [17] RColorBrewer_1.1-3 S7_0.2.1           lifecycle_1.0.4    compiler_4.5.2    
## [21] farver_2.1.2       stringr_1.6.0      carData_3.0-5      vipor_0.4.7       
## [25] htmltools_0.5.9    sass_0.4.10        yaml_2.3.12        htmlTable_2.4.3   
## [29] Formula_1.2-5      pillar_1.11.1      car_3.1-3          jquerylib_0.1.4   
## [33] cachem_1.1.0       Hmisc_5.2-4        rpart_4.1.24       abind_1.4-8       
## [37] nlme_3.1-168       tidyselect_1.2.1   digest_0.6.39      stringi_1.8.7     
## [41] purrr_1.2.0        splines_4.5.2      labeling_0.4.3     cowplot_1.2.0     
## [45] fastmap_1.2.0      grid_4.5.2         colorspace_2.1-2   cli_3.6.5         
## [49] magrittr_2.0.4     base64enc_0.1-3    dichromat_2.0-0.1  broom_1.0.11      
## [53] foreign_0.8-90     withr_3.0.2        scales_1.4.0       backports_1.5.0   
## [57] timechange_0.3.0   rmarkdown_2.30     otel_0.2.0         nnet_7.3-20       
## [61] gridExtra_2.3      ggsignif_0.6.4     evaluate_1.0.5     knitr_1.51        
## [65] mgcv_1.9-4         rlang_1.1.6        polynom_1.4-1      glue_1.8.0        
## [69] rstudioapi_0.17.1  jsonlite_2.0.0     R6_2.6.1