TP53_founder_modelling.Rmd

---
title: "Demographic Modeling of the TP53 p.R337H Founder Variant"
author: "Emilia Modolo Pinto, João Carlos Degraf Muzzi, José Andres Yunes, 
Bonald C. Figueiredo, Raul C. Ribeiro, Gerard P. Zambetti"
date: "2025-04-01"
output: pdf_document
---

# Introduction

This document contains the full code and analysis pipeline used to generate the results for the article:  
**"Reconstructing the Origin and Demographic Expansion of the TP53 p.R337H Founder Variant in Brazil"**  
(Pinto EM, Muzzi J, Yunes JA, Figueiredo BC, Ribeiro RC, Zambetti GP)

For full details on data sources, dependencies, and how to reproduce the figures, please refer to the `README.md` file in this repository.

```{r libraries}
# Load required libraries for data manipulation, plotting, and geospatial analysis
library(tidyverse)
library(plotly)
library(scales)
library(cowplot)
library(geobr)
library(sf) 

to_numeric <- function(x) {
  as.numeric(gsub(" ", "", ifelse(x == "-", NA, x)))
}
```

```{r Brazilian_population_screening}
# Read the spatial data for Brazilian states
brazil_states <- read_state(year = 2020)

# 2022 IBGE Census population data
pop_data_2022 <- data.frame(
  abbrev_state = c('RO', 'AC', 'AM', 'RR', 'PA', 'AP', 'TO', 'MA', 'PI', 'CE', 'RN', 'PB', 'PE', 'AL', 'SE', 'BA', 'MG', 'ES', 'RJ', 'SP', 'PR', 'SC', 'RS', 'MS', 'MT', 'GO', 'DF'),
  population = c(1581196, 830018, 3941613, 636707, 8121025, 733759, 1511460, 6776699, 3271199, 8794957, 3302729, 3974687, 9058931, 3127683, 2210004, 14141626, 20539989, 3833712, 16055174, 44411238, 11444380, 7610361, 10882965, 2757013, 3658649, 7056495, 2817381)
)

# Create a new column with the discrete population categories
pop_data_2022 <- pop_data_2022 %>%
  mutate(
    pop_category = case_when(
      population <= 1e6 ~ "Up to 1 million",
      population > 1e6 & population <= 5e6 ~ "1-5 million",
      population > 5e6 & population <= 10e6 ~ "5-10 million",
      population > 10e6 & population <= 20e6 ~ "10-20 million",
      population > 20e6 ~ "> 20 million"
    ),
    # The 'highlight' column is useful for identifying states with established screening programs
    highlight = ifelse(abbrev_state %in% c("SP", "PR", "SC"), 
                       "States with Screening", 
                       "Other States")
  )

# Convert pop_category to a factor to ensure the legend is ordered correctly
pop_data_2022$pop_category <- factor(pop_data_2022$pop_category,
                                     levels = c("Up to 1 million",
                                                "1-5 million",
                                                "5-10 million",
                                                "10-20 million", 
                                                "> 20 million"))

# Join the population data with the spatial data
map_data <- left_join(brazil_states, pop_data_2022, by = "abbrev_state")


# --- 3. CREATE THE MAP ---

y_offset <- 0.8 
g_map = ggplot(data = map_data) +
  geom_sf(aes(fill = pop_category), color = "white", linewidth = 0.1) +
  geom_text(
    data = . %>% filter(highlight == "States with Screening"), 
    aes(x = st_coordinates(st_centroid(geom))[, 1],
        y = st_coordinates(st_centroid(geom))[, 2] - y_offset), 
    label = "*",         
    size = 8,           
    color = "grey",     
    fontface = "bold"    
  ) +
  
  scale_fill_brewer(palette = "YlGnBu", name = "Population (2022)") +
  
  labs(
    title = "Population of Brazilian States by Category",
    subtitle = "Population according to the 2022 IBGE Census",
    caption = "Source: IBGE. * States with established\nneonatal screening programs (SP, PR, SC)."
  ) +
  theme_void() +
  theme(
    legend.position = "right",
    plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
    plot.subtitle = element_text(hjust = 0.5, size = 8),
    plot.caption = element_text(hjust = 0.5, face = "italic", size = 8), legend.title = element_text(size = 8), legend.text = element_text(size = 8)
  )

g_map

ggsave(g_map, filename = "./Figure 1 map.pdf",
       device = "pdf", height = 3.5, width = 5, unit = "in")

```


```{r states_data, include = F}

# Historical Population data from São Paulo, Paraná and Santa Catarina
# Source: https://seriesestatisticas.ibge.gov.br/
df <- data.frame(
                      Year = c(1782L,1808L,
                              1819L,1823L,1830L,1854L,1867L,1869L,1872L,
                              1890L,1900L,1920L,1940L,1950L,1960L,1970L,1980L,
                              1991L,2000L,2010L,2022L),
        Population = c(140874L,239095L, 342296L,330000L,650000L,677400L,
                              1220000L, 1090000L,1123878L,1918013L,2929704L,5946642L,
                              9594932L,12810472L,19224379L,27632227L,36812603L,
                              44528020L,51877510L,57955161L,63465219L)
      )

# Actual population from São Paulo (SP), Paraná (PR) and Santa Catarina (SC)
# Source: IBGE Census 2022

SP = 44411238
PR = 11444380
SC = 7609601

# Number of actual carriers, proportions retrieved from:
#  0.21% for SP, 0.306 % for PR, and 0.24% for SC
# Custódio et al., 2013; Seidinger et al., 2020; Costa et al., 2019
carriers = SP*0.0021 + PR*0.00306 + 0.0024*SC

```

## Exponential model

We applied the exponential growth model described by Labuda et al. (1996) to estimate population growth rates in Paraná (PR), São Paulo (SP), and Santa Catarina (SC) from 1782 to 2022. 

$$
P(t) = P_0 \cdot e^{rt/25}
$$

Where:

- \( P(g) \): population size after \( g \) generations  
- \( P_0 \): initial population size  
- \( r \): growth rate per generation  
- \( t \): time in years  
- \( \frac{t}{25} \): approximate number of generations, assuming 25 years per generation
- \( e \): Euler’s number (\(\approx 2.718\))

```{r exponential_model}

# Retrieve initial and final population and time period
Pg = df[nrow(df), 2]
Po = df[1, 2]
t =  df[nrow(df), 1]
to =  df[1, 1]

t_delta= t - to

# Converting time in year to generation
g = t_delta/25
# calculating growth rate
r = log(Pg/Po)/g # 0.6364981
df$t_delta= df$Year - to
df$g = df$t_delta/25

# Generating estimates with the exponential model using the calculated growth rate
df$exp = Po * exp(r*df$g)
mse_original <- mean((df$Population - df$exp)^2)
df_log <- subset(df, Population > 0 & exp > 0)

# Plot in normal scale
g1=  ggplot(df, aes(x = Year)) +
    geom_line(aes(y = Population/1000000, color = "IBGE Census"))+
    geom_line(aes(y= exp/1000000, color = "Exp. Model"))+
    annotate(
      "text", 
      x = (min(df$Year)+10), 
      y = max(df$Population/1000000),
      label = paste("r=", round(r,3)),
      hjust = 0, vjust = 1, size = 2.5
    )+
    labs(
      x = "Year",
      y = "Population in São Paulo, Santa\nCatarina, and Paraná (millions)",
      colour = "",
      title = ""
    )+
  theme_minimal() +
    scale_color_manual(
            values = c(
              "Exp. Model"     =  "#e66101",  
              "IBGE Census"    =  "#1f78b4"
            )
          )+
  theme(
    legend.position = "bottom",
    axis.title.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
  legend.title = element_text(size = 8),
  legend.text = element_text(size = 8)
    
  )

g1

# Plot in log scale
 g2= ggplot(df, aes(x = Year)) +
    geom_line(aes(y = log(Population), color = "IBGE Census"))+
    geom_line(aes(y= log(exp), color = "Exp. Model"))+
    annotate(
      "text", 
      x = (min(df$Year)+10), 
      y = max(log(df$Population)),
      label = paste("r=", round(r,3)),
      hjust = 0, vjust = 1, size = 2.5
    )+
    labs(
      x = "Year",
      y = "Population in São Paulo, Santa\nCatarina, and Paraná (Log-scale)",
      colour = "",
      #title = "Population of São Paulo (SP) and Paraná (PR) Over Time"
      title = ""
    )+
   theme_minimal() +
    scale_color_manual(
            values = c(
              "Exp. Model"     =  "#e66101",  
              "IBGE Census"    =  "#1f78b4"
            )
          )+
  theme(
    legend.position = "bottom",
    axis.title.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
  legend.title = element_text(size = 8),
  legend.text = element_text(size = 8)
  )
 


g2
```

```{r}
# Calculate R² of log scale
log_obs <- log(df$Population)
log_pred <- log(df$exp)

# Linear regression
r_squared <- summary(lm(log_obs ~ log_pred))$r.squared
cat("R² (log-scale):", r_squared, "\n")

# Adjustment metrics on normal scale
mean(abs(df$Population - df$exp))  # MAE 3875792
sqrt(mean((df$Population - df$exp)^2))  # RMSE 6924844
mean(abs((df$Population - df$exp) / df$Population)) * 100 #MAPE 19.70762

# Adjustment metrics on log scale
mean(abs(log_obs - log_pred))  # MAE 0.2179778
sqrt(mean((log_obs - log_pred)^2))  # RMSE 0.2664648
mean(abs((log_obs - log_pred) / log_obs)) * 100 #MAPE 1.39116


```

The logarithmic scale improves visualization of the overall fit over time, highlighting that, on average, the exponential curve adequately captures the general pattern of population growth.



```{r year_founder_event}
# Estimating the year of the founder event using the fixed growth rate
g_test = log(carriers/1)/r
t_test = g_test*25
2022-t_test # 1554
```


```{r monte_carlo}
# Monte Carlo simulation to estimate TMRCA distribution under varying founder scenarios
set.seed(123)

# São Paulo carriers data  (Seidinger et al., 2022)
k_sp <- 75
n_sp <- 34344

# Paraná carriers data (Costa et al., 2019)
k_pr <- 461+159
n_pr <- 42438 + 171649

# Santa Catarina carriers data  (Costa et al., 2019)
k_sc <- 125
n_sc <- 50115

# Applying Labuda's exponential growth method (anchored on earliest and latest population estimates)
N_final_pop <- tail(df$Population, 1)
N_initial_pop <- df$Population[1]
g_total <- tail(df$g, 1)

r_fix <- log(N_final_pop / N_initial_pop) / g_total

cat("--- Mean Growth Model ---\n")
cat("Fixed Generational Growth Rate (r):", r_fix, "\n\n")

# --- BOOTSTRAP SIMULATION SETUP ---
n_sims <- 10000

# Parameters for generation time uncertainty (in years)
# Assuming a normal distribution based on literature data.
mean_gen_time <- 25
sd_gen_time <- 3 # Reasonable standard deviation, can be adjusted.

# Number of founders
for(j in c(2,4,6,1)){
  print(j)
  n_founders <- j
  print(paste0("Number of founders = ", n_founders))
  # --- 4. MONTE CARLO SIMULATION ---
  
  tmrca_results <- numeric(n_sims)
  
  for (i in 1:n_sims) {
    # Draw carrier frequencies and generation time
    freq_sp_i <- rbeta(1, k_sp + 1, n_sp - k_sp + 1)
    freq_sc_i <- rbeta(1, k_sc + 1, n_sc - k_sc + 1)
    freq_pr_i <- rbeta(1, k_pr + 1, n_pr - k_pr + 1)
    gen_time_i <- rnorm(1, mean = mean_gen_time, sd = sd_gen_time)
    
    # Calculate total number of carriers
    total_carriers_i <- (SP * freq_sp_i) + (SC * freq_sc_i) + (PR * freq_pr_i)
    
    # Calculate number of generations using FIXED 'r'
    number_generations <- log(total_carriers_i / n_founders) / r_fix
    
    # Calculate TMRCA in years
    tmrca_results[i] <- number_generations * gen_time_i
  }
  
  # --- RESULTS AND VISUALIZATION OF TMRCA UNCERTAINTY ---
  
  # Remove NAs and calculate median and CI
  tmrca_results <- tmrca_results[!is.na(tmrca_results)]
  median_tmrca <- median(tmrca_results)
  ci_tmrca <- quantile(tmrca_results, probs = c(0.025, 0.975))
  
  # Print final results
  cat("--- Final Simulation Results (Mean Growth Model) ---\n")
  cat("Estimated Median TMRCA:", round(median_tmrca, 1), "years\n")
  cat("95% Confidence Interval:", round(ci_tmrca[1], 1), "-", round(ci_tmrca[2], 1), "years\n")
  cat("Median Year of Founder Event:", 2022 - round(median_tmrca, 0), "\n")
  cat("Founder Event Year Interval (95% CI):", (2022 - round(ci_tmrca[2], 0)), "-", (2022 - round(ci_tmrca[1], 0)), "\n")
  
  # Plot to visualize TMRCA UNCERTAINTY (histogram)
  results_df <- data.frame(TMRCA = tmrca_results)
  g_mc = ggplot(results_df, aes(x = TMRCA)) +
    geom_histogram(bins = 50, fill = "#08519c", alpha = 0.8, color = "white") +
    geom_vline(xintercept = median_tmrca, color = "red", 
               linetype = "dashed", linewidth = 1) +
    geom_vline(xintercept = ci_tmrca[1], color = "black", 
               linetype = "dotted", linewidth = 1) +
    geom_vline(xintercept = ci_tmrca[2], color = "black", 
               linetype = "dotted", linewidth = 1) +
    labs(
      title = paste0("Estimated TMRCA Distribution (Mean Growth Model) - ", 
                     n_founders, 
                     " founder(s)"),
      subtitle = paste0("Median = ", round(median_tmrca, 1), 
                        " years | 95% CI: [", round(ci_tmrca[1], 1), ", ",
                        round(ci_tmrca[2], 1), "]"),
      x = "Time to the Most Recent Common Ancestor (years)",
      y = "Frequency"
    ) +
    theme_minimal()+
    theme(
      legend.position = "bottom",
      axis.title.y = element_text(size = 8),
      axis.title.x = element_text(size = 8),
      plot.title = element_text(size = 10),
      plot.subtitle = element_text(size = 8),
    legend.title = element_text(size = 8),
    legend.text = element_text(size = 8)
    )
  g_mc
  
  ggsave(g_mc, filename = paste0("./Monte Carlo_", n_founders, ".pdf"),
         device = "pdf", height = 3.5, width = 7, unit = "in")
}
```

# Brazil and Portugal analysis

```{r Brazil}
# Brazil historical data
# Source: https://seriesestatisticas.ibge.gov.br/
df_br = data.frame(
          stringsAsFactors = FALSE,
                              Year = c(1550L,1576L,1583L,1600L,1660L,1690L,
                                       1700L,1766L,1776L,1776L,1776L,
                                       1780L,1780L,1798L,1798L,1798L,1798L,
                                       1800L,1800L,1808L,1819L,1823L,1830L,
                                       1854L,1867L,1869L,1872L,1890L,1900L,
                                       1920L,1940L,1950L,1960L,1970L,1980L,
                                       1991L,2000L,2010L,2022L
                                      ),
          Population = c(15000,17100,57000,100000,
                                       184000,242000,300000,
                                       1500000,1788480,1900000,
                                       2700000,2523000,2841000,2888078,
                                       3569000,3800000,4000000,
                                       3660000,3250000,2424463,
                                       3596132,3960866,5350000,
                                       7677800,11280000,10200000,
                                       9930478,14333915,17318556,
                                       30635605,41236315,51944397,
                                       70992343,94508583,121150573,
                                       146917459,169590693,190755799,
                                       203062512
                         )
        )

# When multiple historical records were available for the same year, their average was used.
df_br <- df_br %>%
  group_by(Year) %>%
  summarise(Population = mean(Population))
```

```{r}

# Brazilian exponential model (infer growth rate)
N_br = df_br[nrow(df_br), 2, drop = T]
N_bro = df_br[1, 2, drop = T]
t =  df_br[nrow(df_br), 1, drop = T]
to =  df_br[1, 1, drop = T]

t_delta= t - to
g = t_delta/25
# Growth rate
r_br = log(N_br/N_bro)/g # 0.5038781

# Visualizing the fit of the exponential model
df_br$t_delta= df_br$Year - to
df_br$g = df_br$t_delta/25

df_br$exp = N_bro * exp(r_br*df_br$g)
df_br_log <- subset(df_br, Population > 0 & exp > 0)

g3 =  ggplot(df_br, aes(x = Year)) +
    geom_line(aes(y = Population/1000000, color = "IBGE Census"))+
    geom_line(aes(y= exp/1000000, color = "Exp. Model"))+
    annotate(
      "text", 
      x = (min(df_br$Year)+10), 
      y = max(df_br$Population/1000000),
      label = paste("r=", round(r_br,3)),
      hjust = 0, vjust = 1, size = 2.5
    )+
    labs(
      x = "Year",
      y = "Population in Brazil (millions)",
      colour = "Model",
      title = ""
    )+
    scale_color_manual(
            values = c(
              "Exp. Model"     =  "#33a02c",  
              "IBGE Census"    = "#e31a1c" 
            )
          )+
  theme(
    legend.position = "bottom"
  )

g3
```



```{r Brazil_log_scale}

# Visualizing the model fit in log scale
log_obs <- log(df_br$Population)
log_pred <- log(df_br$exp)


 g4= ggplot(df_br, aes(x = Year)) +
    geom_line(aes(y = log(Population), color = "IBGE Census"))+
    geom_line(aes(y= log(exp), color = "Exp. Model"))+
    annotate(
      "text", 
      x = (min(df_br$Year)+10), 
      y = max(log(df_br$Population)),
      label = paste("r=", round(r_br,3)),
      hjust = 0, vjust = 1, size = 2.5
    )+
    labs(
      x = "Year",
      y = "Population in Brazil\n (Log-scale)",
      colour = "Model",
      #title = "Population of São Paulo (SP) and Paraná (PR) Over Time"
      title = ""
    )+
    scale_color_manual(
            values = c(
              "Exp. Model"     =  "#33a02c",  
              "IBGE Census"    = "#e31a1c" 
            )
          )+
  theme(
    legend.position = "bottom"
  )


g4
```


```{r}
# Portugal historical data
# source: https://www.ine.pt/
df_pt = data.frame(
  stringsAsFactors = FALSE,
       check.names = FALSE,
                       Year = c(1422L,1527L,
                               1636L,1736L,1770L,1776L,1801L,1811L,1838L,
                               1849L),
                      Pop. = c("1 043 274",
                               "1 262 376","1 100 000","2 143 368","2 850 444",
                               "3 352 310","2 931 930","2 876 602",
                               "3 200 000","3 411 454"),
                      `±%` = c("—","+21.0%",
                               "−12.9%","+94.9%","+33.0%","+17.6%","−12.5%",
                               "−1.9%","+11.2%","+6.6%")
        )


df_pt2 = data.frame(
  stringsAsFactors = FALSE,
       check.names = FALSE,
                        Year = c(1864L,1878L,
                                1890L,1900L,1911L,1920L,1930L,1940L,1950L,
                                1960L,1970L,1981L,1991L,2001L,2011L,
                                2021L),
                       Pop. = c("4 188 410",
                                "4 550 699","5 049 729","5 423 132",
                                "5 960 056","6 032 991","6 825 883","7 722 152",
                                "8 441 312","8 851 289","8 568 703","9 852 841",
                                "9 862 540","10 356 117","10 562 178",
                                "10 344 802"),
                       `±%` = c("—","+8.6%",
                                "+11.0%","+7.4%","+9.9%","+1.2%","+13.1%",
                                "+13.1%","+9.3%","+4.9%","−3.2%","+15.0%",
                                "+0.1%","+5.0%","+2.0%","−2.1%")
         )

df_pt = rbind(df_pt, df_pt2)
df_pt <- df_pt %>% filter(Year > 1500)

# Removing spaces and converting to numeric
df_pt$Pop. <- as.numeric(gsub(" ", "", df_pt$Pop.))  
df_pt <- df_pt[, c("Year", "Pop.")]
names(df_pt) <- c("Year", "Population")
df_pt$Country <- "Portugal"

# Portugal growth rate
N_pt = df_pt[nrow(df_pt), 2, drop = T]
N_pto = df_pt[1, 2, drop = T]
t =  df_pt[nrow(df_pt), 1, drop = T]
to =  df_pt[1, 1, drop = T]

t_delta= t - to
g = t_delta/25
r_pt = log(N_pt/N_pto)/g # 0.1064
```

```{r}

# Merging with Brazil data
df_br <- df_br[, c("Year", "Population")]
df_br$Country <- "Brazil"

df_combined <- rbind(df_br, df_pt)

g_pt = ggplot(df_combined, aes(x = Year, y = Population / 1e6, color = Country)) +
  geom_line(size = .7) +
  #geom_point(size = .5) +
  theme_minimal() +
  theme(legend.position = "bottom")+
  scale_color_manual(values = c("Brazil" = "forestgreen", "Portugal" = "red")) +
  labs(title = "Evolution of population from Brazil and Portugal",
       x = "",
       y = "Population (millions)",
       color = "Country")+
  theme(
    legend.position = "bottom",
    axis.title.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    plot.title = element_text(size = 10),
    plot.subtitle = element_text(size = 8)
  )
g_pt

ggsave(g_pt, filename = "./Figure 2.pdf",
       device = "pdf", height = 3.5, width = 3.5, unit = "in")

```

```{r}

### log scale
g_pt2 = 
    ggplot(df_combined, aes(x = Year, y = log(Population), color = Country)) +
     geom_line(size = 1) +
     #geom_point(size = .5) +
     theme_minimal() +
     theme(legend.position = "bottom")+
     scale_color_manual(values = c("Brazil" = "forestgreen", "Portugal" = "red")) +
     labs(title = "Evolution of population from Brazil and Portugal",
          x = "",
          y = "Population (log scale)",
          color = "Country")+
     theme(
         legend.position = "bottom",
         axis.title.y = element_text(size = 8),
         axis.title.x = element_text(size = 8),
         plot.title = element_text(size = 10),
         plot.subtitle = element_text(size = 8),
          legend.title = element_text(size = 8),
          legend.text = element_text(size = 8)
     )
g_pt2
```


```{r exporting_data}
writexl::write_xlsx(list(
  Portugal = df_pt, 
  Brazil = df_br),
  "./BR and PT populations.xlsx")
```


```{r generating_final_figure}
g_final = plot_grid(g1, g2, g_mc, g_pt, labels = c("A", "B", "C", "D"))

ggsave(g_final, filename = "./Figure 2.pdf",
       device = "pdf", height = 7, width = 8, unit = "in")
g_final
```

# R session

```{r r_session}
sessionInfo()
```

R version 4.5.1 (2025-06-13)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Sao_Paulo
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] sf_1.0-21       geobr_1.9.1     cowplot_1.2.0   scales_1.4.0   
 [5] plotly_4.11.0   lubridate_1.9.4 forcats_1.0.0   stringr_1.5.1  
 [9] dplyr_1.1.4     purrr_1.1.0     readr_2.1.5     tidyr_1.3.1    
[13] tibble_3.3.0    ggplot2_3.5.2   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       xfun_0.52          htmlwidgets_1.6.4  tzdb_0.5.0        
 [5] vctrs_0.6.5        tools_4.5.1        generics_0.1.4     curl_6.4.0        
 [9] proxy_0.4-27       pkgconfig_2.0.3    KernSmooth_2.23-20 data.table_1.17.8 
[13] RColorBrewer_1.1-3 lifecycle_1.0.4    compiler_4.5.1     farver_2.1.2      
[17] textshaping_1.0.1  htmltools_0.5.8.1  class_7.3-20       yaml_2.3.10       
[21] lazyeval_0.2.2     pillar_1.11.0      crayon_1.5.3       rsconnect_1.5.0   
[25] classInt_0.4-11    wk_0.9.4           tidyselect_1.2.1   digest_0.6.37     
[29] stringi_1.8.7      labeling_0.4.3     fastmap_1.2.0      grid_4.5.1        
[33] cli_3.6.5          magrittr_2.0.3     dichromat_2.0-0.1  e1071_1.7-16      
[37] withr_3.0.2        writexl_1.5.4      timechange_0.3.0   rmarkdown_2.29    
[41] httr_1.4.7         ragg_1.4.0         hms_1.1.3          evaluate_1.0.4    
[45] knitr_1.50         viridisLite_0.4.2  s2_1.1.9           rlang_1.1.6       
[49] Rcpp_1.1.0         glue_1.8.0         DBI_1.2.3          rstudioapi_0.17.1 
[53] jsonlite_2.0.0     R6_2.6.1           fs_1.6.6           systemfonts_1.2.3 
[57] units_0.8-7