Relationship_fitness_genomic_offset
Thomas Francisco
2024-09-18
1 Introduction
To evaluate the relevance of the genomic offset in providing insight
into potential maladaptation to climate across studied T.
baccata populations, we investigated the potential relationship
between fitness proxies and genomic offset predictions using an
independent dataset of populations not used to build the gene-climate
models. Genomic offset values for each population were calculated as the
Euclidean distance between the predicted optimal genomic compositions
under the population’s reference climate and the clonal bank’s climate
from the planting to the measurement dates. We used two methods to
calculate the genomic offset:
- RDA method with 7 sets of SNPs
- Gradient forest with 7 sets of SNPs
For the phenotypic data, we aimed to retain the most fitness-related
traits. We selected four traits: - Growth: Shoot volume
- Growth phenology: Spring elongation
- Reproductive phenology: proportion of male open strobili
- Ecophysiological: leaf thickness
To search for an association between the phenotypic data and the
genomic offset predictions, we used several analyses: - Pearson
correlation
- linear model: inferential and Bayesian approaches
- non linear model: inferential and Bayesian approaches
2 Data
meta_data <- read.csv("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Populations/Populations_common_garden.csv",h=T,sep=";",dec=",")2.1 BLUPs data
First, we loaded the phenotypic BLUP data at the population level.
Only the BLUPs at the population level were kept.
For growth, growth, and reproductive phenology:
#list_traits <-c("growth","growth_pheno","repro_pheno","leafthickness","LMA","StomDens","D13C")
list_traits <-c("growth","growth_pheno","repro_pheno")#"growth_newformulationfinal"
#selection of the right coeff
list_x <- c(2,2,2)#,3
list_y <- c(27,27,19)#,28
for(i in 1:length(list_traits)){
trait <- list_traits[i]
X <- list_x[i]
Y <- list_y[i]
data <- read.csv(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/Blups_mean/Blups_",trait,"_period.csv"),h=T,sep=";",dec=",")
data_pop <- data[X:Y,] #keep only the Blup at the population level
colnames(data_pop) <- c("Population","Blups")
data_pop$Population <- str_remove(data_pop$Population, "Population\\.")
assign(paste0("df_",trait),data_pop)
}For leaf thickness:
list_traits <-c("leafthickness")#"growth_bis","StomDens","LMA"
#selection of the right coeff
list_x <- c(2,2,2,2)
list_y <- c(27,27,27,27)
for(i in 1:length(list_traits)){
trait <- list_traits[i]
X <- list_x[i]
Y <- list_y[i]
data <- read.csv(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/Blups_mean/Blups_",trait,"final_period.csv"),h=T,sep=";",dec=",")
data_pop <- data[X:Y,] #keep only the Blup at the population level
colnames(data_pop) <- c("Population","Blups")
data_pop$Population <- str_remove(data_pop$Population, "Population\\.")
assign(paste0("df_",trait),data_pop)
}We removed populations that had too few individuals measured.
#Shootvol
#We remove the ACEBDO, LASENIA pop
df_growth_final <- df_growth %>%
filter(!(Population %in% c("ACEBEDO", "LASENIA")))
#springelong
#We remove the ACEBEDO, ALCOY, LASENIA and QUESADA pop
df_growth_pheno_final <- df_growth_pheno %>%
filter(!(Population %in% c("ACEBEDO", "ALCOY", "LASENIA","QUESADA")))
#totalstrobopen
#We remove the BOCA, JERTE, LASENIA, PUBLA and RASQUERA
df_repro_pheno_final <- df_repro_pheno %>%
filter(!(Population %in% c("BOCAHUERGANO", "JERTE", "LASENIA"," PUEBLALILLO","RASQUERA")))
#leaves traits
#remove ACEBEDO, ALCOY, LASENIA and QUESADA
##leafthickness
df_leafthickness_final <- df_leafthickness %>%
filter(!(Population %in% c("ACEBEDO", "ALCOY", "LASENIA"," QUESADA")))
#Shootvol new
#We remove the ACEBDO, LASENIA pop
#df_growth_newformulation_final <- df_growth_newformulationfinal %>%
#filter(!(Population %in% c("ACEBEDO", "LASENIA")))2.2 Correlation phenotypic traits
We investigated the correlation of the phenotypic traits of populations:
colnames(df_growth_final) <- c("Population","Growth")
colnames(df_growth_pheno_final) <- c("Population","Growth_pheno")
colnames(df_repro_pheno_final) <- c("Population","Repro_pheno")
colnames(df_leafthickness_final) <- c("Population","Leafthickness")
merged_df <- df_growth_final %>%
full_join(df_growth_pheno_final, by = "Population") %>%
full_join(df_repro_pheno_final, by = "Population") %>%
full_join(df_leafthickness_final, by = "Population")correlation_phenotypic_traits <- cor(merged_df[,-1], use="pairwise.complete.obs")
corrplot(correlation_phenotypic_traits, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6)
pdf("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/figures/Correlation/correlation_phenotypic_traits_BLUPs.pdf");corrplot(correlation_phenotypic_traits, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6);dev.off()## png
## 2We can see that the phenotypic traits do not appear to be closely related, as the correlations are quite low. Could this indicate that they are responding to different selective pressures? Or is there a trade-off? Does that also mean that they provide complementary information about fitness?
2.3 Composite fitness index
We also attempted to combine these various traits into a single index, as fitness is multi-faceted, particularly for long-lived organisms such as trees. The idea is to organise populations based on their trait values and what appears to be adaptive in the common garden. For instance, high growth is advantageous, so the population with the highest growth rate out of the 26 populations will have the highest value. These values will be scaled to give equal weight to each trait, since we do not have the same number of populations for each trait (using rank values would give more weight to traits with more populations). To standardise the values, we are applying the following formula: Standardized trait= (trait - min_value_trait)/(max_value_trait - min_value_trait)
list_traits <- c("growth","growth_pheno","repro_pheno","leafthickness")#"growth_newformulation",
for(x in 1: length(list_traits)){
trait <- list_traits[x]
df <- get(paste0("df_", trait,"_final"))
df_standardized <- df %>%
mutate(
trait_scaled = (df[,2] - min(df[,2])) / (max(df[,2]) - min(df[,2])))
colnames(df_standardized) <- c("Population","Blups",paste0(trait,"_standardized"))
assign(paste0("df_standardized_", trait), df_standardized)
}final_tot_standardized_traits <- reduce(list(df_standardized_growth[,-2], df_standardized_growth_pheno[,-2], df_standardized_repro_pheno[,-2], df_standardized_leafthickness[,-2]), full_join, by = "Population")#df_standardized_growth_newformulationTrait_calculation <- final_tot_standardized_traits %>%
rowwise() %>%
dplyr::mutate(
# Calculate the sum of standardized traits, ignoring NAs
valid_traits_sum = sum(c_across(ends_with("_standardized")), na.rm = TRUE),
# Count non-NA standardized traits
valid_traits_count = sum(!is.na(c_across(ends_with("_standardized")))),
# Calculate the trait index only if there is at least one valid trait
trait_index = if_else(valid_traits_count > 0,
valid_traits_sum / valid_traits_count,
NA_real_)
) %>%
dplyr::select(-valid_traits_sum, -valid_traits_count) %>% # Optionally remove helper columns
ungroup()
# View the combined data frame with trait index
print(Trait_calculation)## # A tibble: 26 × 6
## Population growth_standardized growth_pheno_standar…¹ repro_pheno_standard…²
## <chr> <dbl> <dbl> <dbl>
## 1 AGRES 0 0.767 0.276
## 2 ALCOY 0.264 NA NA
## 3 ARAGUES 0.206 0.269 0.297
## 4 ASPONTES 0.186 0.438 NA
## 5 BARRACO 0.213 0.474 1
## 6 BEGONTE 0.293 0.554 NA
## 7 BOCAHUERGA… 0.279 0.789 NA
## 8 CANENCIA 0.588 0.755 0.441
## 9 CERVANTES 0.318 0.755 NA
## 10 FANLO 0.591 0.576 0.0433
## # ℹ 16 more rows
## # ℹ abbreviated names: ¹growth_pheno_standardized, ²repro_pheno_standardized
## # ℹ 2 more variables: leafthickness_standardized <dbl>, trait_index <dbl>df_composite_fitness_final <- data.frame(Population=Trait_calculation[,1],Fitness_index=Trait_calculation$trait_index)2.4 GO predictions
We loaded the genomic offset predictions for the 2 genomic offset methods:
#RDA standard GO
list_RDA_standard <- c("all","random","random_same_AF","LC","CG","random_AF_V2","all_outliers")
list_period <- c("CG2012","CG2021")
for(i in 1: length(list_period)){
period <- list_period[i]
for(x in 1:length(list_RDA_standard)){
name <- list_RDA_standard[x]
load(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/standard_GO/RDA/data/CG_period/GO_RDA_",name,"_",period,"_Taxus.Rdata"))
}
}
#GF GO
list_GF_standard <- c("all","random","random_same_AF","LC","CG","random_AF_V2","all_outliers")
list_period <- c("2012","2021")
for(i in 1: length(list_period)){
period <- list_period[i]
for(x in 1:length(list_GF_standard)){
name <- list_GF_standard[x]
load(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/Standard_GO/GF/data/CG_period/GO_standard_GF_",name,"_",period,".Rdata"))
}
}We can investigate the correlation between the predictions across methods and models:
df_GO <- data.frame(RDA_all=GO_RDA_all_CG2012_Taxus$GO_all_CG2012,RDA_ran=GO_RDA_random_CG2012_Taxus$GO_random_CG2012,RDA_rand_AF=GO_RDA_random_same_AF_CG2012_Taxus$GO_random_same_AF_CG2012,RDA_LC=GO_RDA_LC_CG2012_Taxus$GO_LC_CG2012,RDA_CG=GO_RDA_CG_CG2012_Taxus$GO_CG_CG2012,RDA_ran_AF_V2=GO_RDA_random_AF_V2_CG2012_Taxus$GO_random_AF_V2_CG2012,RDA_all_outliers=GO_RDA_all_outliers_CG2012_Taxus$GO_all_outliers_CG2012,GF_ran=GO_standard_GF_random_2012$genomic_offset_GF_random,GF_ran_AF=GO_standard_GF_random_same_AF_2012$genomic_offset_GF_random_same_AF,GF_LC=GO_standard_GF_LC_2012$genomic_offset_GF_LC,GF_CG=GO_standard_GF_CG_2012$genomic_offset_GF_CG,GF_ran_AF_V2=GO_standard_GF_random_AF_V2_2012$genomic_offset_GF_random_AF_V2,GF_all_outliers=GO_standard_GF_all_outliers_2012$genomic_offset_GF_all_outliers)
df <- cor(df_GO)
ggcorrplot(df, hc.order = F, type = "lower", lab = TRUE, title = "Correlation Matrix between Genomic Offsets across models")
3 Relationship fitness ~ genomic offset
3.1 Pearson correlations
# List of RDA genomic offsets for 2012 and 2021
list_GO_2012_RDA <- c("GO_RDA_all_CG2012_Taxus","GO_RDA_random_CG2012_Taxus", "GO_RDA_random_same_AF_CG2012_Taxus",
"GO_RDA_LC_CG2012_Taxus","GO_RDA_random_AF_V2_CG2012_Taxus","GO_RDA_all_outliers_CG2012_Taxus")
list_GO_2021_RDA <- c("GO_RDA_all_CG2021_Taxus","GO_RDA_random_CG2021_Taxus", "GO_RDA_random_same_AF_CG2021_Taxus",
"GO_RDA_LC_CG2021_Taxus","GO_RDA_random_AF_V2_CG2021_Taxus","GO_RDA_all_outliers_CG2021_Taxus")
# List of Gradient Forest genomic offsets for 2012 and 2021
list_GO_2012_GF <- c("GO_standard_GF_all_2012","GO_standard_GF_random_2012", "GO_standard_GF_random_same_AF_2012",
"GO_standard_GF_LC_2012","GO_standard_GF_random_AF_V2_2012","GO_standard_GF_all_outliers_2012")
list_GO_2021_GF <- c("GO_standard_GF_all_2021","GO_standard_GF_random_2021", "GO_standard_GF_random_same_AF_2021",
"GO_standard_GF_LC_2021","GO_standard_GF_random_AF_V2_2021","GO_standard_GF_all_outliers_2021")
# Combine lists into one for easier processing
list_GO_2012 <- c(list_GO_2012_RDA, list_GO_2012_GF)
list_GO_2021 <- c(list_GO_2021_RDA, list_GO_2021_GF)
# Corresponding methods and models
list_method <- c(rep("RDA", 6), rep("GF", 6)) # 4 RDA models and 4 GF models
list_model <- c("all","random", "random_same_AF", "LC", "all","random", "random_same_AF", "LC")
list_model <- c("all","random", "random_same_AF", "outlier","random_same_AF_V2","all_outliers", "all","random", "random_same_AF", "outlier","random_same_AF_V2","all_outliers")
# Traits considered
#list_traits <- c("growth","growth_pheno","repro_pheno","leafthickness","LMA","StomDens","D13C")
list_traits <- c("growth","growth_pheno","repro_pheno","leafthickness","composite_fitness")
# Create an empty dataframe to store correlations
correlation_df <- data.frame(
Method = character(),
Model = character(),
Trait = character(),
Correlation = numeric(),
P_value = numeric(),
stringsAsFactors = FALSE
)
# Loop over each Genomic Offset for both RDA and GF models
for(i in 1:12) { # Now we have 8 models (4 RDA + 4 GF)
GO_pred_2012 <- get(list_GO_2012[i])
GO_pred_2021 <- get(list_GO_2021[i])
method <- list_method[i]
model <- list_model[i]
# Loop over each trait
for(x in 1:length(list_traits)) {
trait <- list_traits[x]
# Use CG2021 data for Height, and CG2012 for other traits
if (trait == "Leafthickness" || trait=="Composite_fitness") {
df_trait <- get(paste0("df_", trait,"_final"))
# Merge data frames based on population for CG2021
merged_df <- merge(df_trait, GO_pred_2021, by = "Population", suffixes = c(".trait", ".GO"))
} else {
df_trait <- get(paste0("df_", trait,"_final"))
# Merge data frames based on population for CG2012
merged_df <- merge(df_trait, GO_pred_2012, by = "Population", suffixes = c(".trait", ".GO"))
}
# Calculate the correlation
correlation <- cor(merged_df[, 2], merged_df[, 3])
pval <- cor.test(merged_df[, 2],merged_df[, 3])
P_value <- pval$p.value
# Store the correlation in the dataframe
correlation_df <- rbind(correlation_df, data.frame(Method = method, Model = model, Trait = trait, Correlation = correlation, P_value= P_value))
}
}
# Print the correlation data frame
#print(correlation_df)Corrplot phenotypic traits, genomic offset
# Reshape the data for corrplot if necessary
cor_matrix <- reshape(correlation_df[,-5],
timevar = "Trait",
idvar = c("Method", "Model"),
direction = "wide")
# Set rownames based on a combination of 'Method' and 'Model'
rownames(cor_matrix) <- paste(cor_matrix$Method, cor_matrix$Model, sep = "_")
# Remove the 'Method' and 'Model' columns since they are now in the rownames
cor_matrix <- cor_matrix[ , !(names(cor_matrix) %in% c("Method", "Model"))]
colnames(cor_matrix) <- list_traits
# create a correlation plot
cor_matrix_fitness_index <- t(cor_matrix)#to have the traits on the abscissa
ggcorrplot(cor_matrix_fitness_index, hc.order = F, type = "full", lab = TRUE, title = "Correlation Matrix between Traits
and Genomic Offsets")
We want to plot the correlation results in a format similar to Camus et al 2024.
# Reshape the correlation matrix for ggplot
cor_long <- melt(cor_matrix_fitness_index, varnames = c("Trait", "Method_Model"), value.name = "Correlation")
#we remove the random correlation
cor_long_reduce <- cor_long[-c(6:10,26:30),]
cor_long_reduce <- cor_long
#change the name of the random same AF correlation to simply random
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "RDA_random_same_AF"] <- "RDA_random"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "GF_random_same_AF"] <- "GF_random"
# Capitalize the first letter of each trait name
cor_long_reduce$Trait <- str_to_title(cor_long_reduce$Trait)
# Separate 'Method_Model' back into 'Method' and 'Model' for better grouping
#cor_long_reduce <- cor_long_reduce %>%
# tidyr::separate(Method_Model, into = c("Method", "Model"), sep = "_")
cor_long_reduce <- cor_long_reduce %>%
tidyr::separate(Method_Model, into = c("Method", "Model"), sep = "_",extra="merge")
cor_long_reduce$Trait[cor_long_reduce$Trait == "Growth_pheno"] <- "Growth Phenology"
cor_long_reduce$Trait[cor_long_reduce$Trait == "Repro_pheno"] <- "Reproductive Phenology"
cor_long_reduce$Trait[cor_long_reduce$Trait == "Leafthickness"] <- "Leaf Thickness"
cor_long_reduce$Trait[cor_long_reduce$Trait == "Composite_fitness"] <- "Composite index"
#Change the order of the model and pheno in the plots
order_model <- c("outlier", "random", "all")
cor_long_reduce$Model <- factor(cor_long_reduce$Model, levels = order_model)
order_pheno <- c("Growth", "Growth Phenology", "Reproductive Phenology","Leaf Thickness","Composite index")
cor_long_reduce$Trait <- factor(cor_long_reduce$Trait, levels = order_pheno)
# Plot using ggplot2
plot <- ggplot(cor_long_reduce, aes(x = Model, y = Correlation, color = Method, group = Method)) +
geom_point(size=3.5) +
facet_wrap(~ Trait, scales = "fixed", nrow = 2) +
scale_y_reverse(
breaks = seq(0, -0.6, by = -0.1), # Custom y-axis breaks
limits = c(0, -0.6) # Ensure consistent limits across facets
) +
labs(
x = "SNP set",
y = "Pearson Correlation"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "bottom"
) +
theme_bw(base_size = 11)
plot
#pdf("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Final_figures/Validation/Graph_nice_correlation_GO_fitness_v2.pdf");print(plot);dev.off()# Reshape the correlation matrix for ggplot
cor_long <- melt(cor_matrix_fitness_index, varnames = c("Trait", "Method_Model"), value.name = "Correlation")
#we remove the random correlation
#cor_long_reduce <- cor_long[-c(6:10,26:30),]
cor_long_reduce <- cor_long
#change the name of the random same AF correlation to simply random
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "RDA_random_same_AF"] <- "RDA_AF"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "GF_random_same_AF"] <- "GF_AF"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "RDA_random_same_AF_V2"] <- "RDA_AF2"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "GF_random_same_AF_V2"] <- "GF_AF2"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "RDA_outlier"] <- "RDA_out"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "GF_outlier"] <- "GF_out"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "RDA_random"] <- "RDA_r"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "GF_random"] <- "GF_r"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "RDA_all_outliers"] <- "RDA_alout"
cor_long_reduce$Method_Model[cor_long_reduce$Method_Model == "GF_all_outliers"] <- "GF_alout"
cor_long_reduce <- cor_long_reduce %>%
mutate(Method_Model = recode(Method_Model,
"RDA_random_same_AF" = "RDA_AF",
"GF_random_same_AF" = "GF_AF",
"RDA_random_same_AF_V2" = "RDA_AF2",
"GF_random_same_AF_V2" = "GF_AF2",
"RDA_outlier" = "RDA_out",
"GF_outlier" = "GF_out",
"RDA_random" = "RDA_r",
"GF_random" = "GF_r",
"RDA_all_outliers" = "RDA_alout",
"GF_all_outliers" = "GF_alout"
))
# Capitalize the first letter of each trait name
cor_long_reduce$Trait <- str_to_title(cor_long_reduce$Trait)
# Separate 'Method_Model' back into 'Method' and 'Model' for better grouping
#cor_long_reduce <- cor_long_reduce %>%
# tidyr::separate(Method_Model, into = c("Method", "Model"), sep = "_")
cor_long_reduce <- cor_long_reduce %>%
tidyr::separate(Method_Model, into = c("Method", "Model"), sep = "_",extra="merge")
cor_long_reduce$Trait[cor_long_reduce$Trait == "Growth_pheno"] <- "Growth Phenology"
cor_long_reduce$Trait[cor_long_reduce$Trait == "Repro_pheno"] <- "Reproductive Phenology"
cor_long_reduce$Trait[cor_long_reduce$Trait == "Leafthickness"] <- "Leaf Thickness"
cor_long_reduce$Trait[cor_long_reduce$Trait == "Composite_fitness"] <- "Composite index"
#Change the order of the model and pheno in the plots
order_pheno <- c("Growth", "Growth Phenology", "Reproductive Phenology","Leaf Thickness","Composite index")
cor_long_reduce$Trait <- factor(cor_long_reduce$Trait, levels = order_pheno)
# Plot using ggplot2
plot <- ggplot(cor_long_reduce, aes(x = Model, y = Correlation, color = Method, group = Method)) +
geom_point(size=3.5) +
facet_wrap(~ Trait, scales = "fixed", nrow = 2) +
scale_y_reverse(
breaks = seq(0, -0.6, by = -0.1), # Custom y-axis breaks
limits = c(0, -0.6) # Ensure consistent limits across facets
) +
labs(
x = "SNP set",
y = "Pearson Correlation"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "bottom"
) +
theme_bw(base_size = 11)
plot
#pdf("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Final_figures/Validation/Graph_new_models_V2.pdf");print(plot);dev.off()We can see that the correlations are negative for almost all
methods and models for each trait. This is consistent with what would be
expected to be advantageous in the climate of the common garden. Higher
growth and late phenology appear to be beneficial, since the common
garden is situated in a mountainous climate sensitive to late frost
events. Additionally, the annual precipitation is relatively low
compared to that experienced by the populations, suggesting that
drought-tolerant traits, such as high leaf thickness, could also be
advantageous.
Depending on the trait, the strength of the correlation varies across
models, indicating that no single model performs significantly better
than the others. Regarding the traits, the composite index appears to be
more closely associated with genomic offset, suggesting that this index
could serve as a better fitness proxy if genomic offset is indeed
negatively correlated with fitness.
3.2 Linear models
To investigate this further, we built linear models to explore and test for relationships between fitness and genomic offset. We used two types of model: inferential and Bayesian, the latter of which is known to perform better with small sample sizes.
3.2.1 Inferential model
list_GO_2012_RDA <- c("GO_RDA_random_AF_V2_CG2012_Taxus","GO_RDA_all_outliers_CG2012_Taxus")#"GO_RDA_all_CG2012_Taxus","GO_RDA_random_CG2012_Taxus", "GO_RDA_random_same_AF_CG2012_Taxus","GO_RDA_LC_CG2012_Taxus",
list_GO_2021_RDA <- c("GO_RDA_random_AF_V2_CG2021_Taxus","GO_RDA_all_outliers_CG2021_Taxus")#"GO_RDA_all_CG2021_Taxus","GO_RDA_random_CG2021_Taxus", "GO_RDA_random_same_AF_CG2021_Taxus","GO_RDA_LC_CG2021_Taxus",
# List of Gradient Forest genomic offsets for 2012 and 2021
list_GO_2012_GF <- c("GO_standard_GF_random_AF_V2_2012","GO_standard_GF_all_outliers_2012")#"GO_standard_GF_all_2012","GO_standard_GF_random_2012", "GO_standard_GF_random_same_AF_2012","GO_standard_GF_LC_2012",
list_GO_2021_GF <- c("GO_standard_GF_random_AF_V2_2021","GO_standard_GF_all_outliers_2021")#"GO_standard_GF_all_2021","GO_standard_GF_random_2021", "GO_standard_GF_random_same_AF_2021","GO_standard_GF_LC_2021",
# Combine lists into one for easier processing
list_GO_2012 <- c(list_GO_2012_RDA, list_GO_2012_GF)
list_GO_2021 <- c(list_GO_2021_RDA, list_GO_2021_GF)
# Corresponding methods and models
list_method <- c(rep("RDA", 2), rep("GF", 2)) # 4 RDA models and 4 GF models
list_name_models <- c("Random_2","all_outlier","Random_2","all_outlier")#"All", "Random","Outlier","All", "Random","Outlier"
list_model <- c("random_same_AF_V2","all_outliers","random_same_AF_V2","all_outliers")#"all","random", "random_same_AF", "outlier","all","random", "random_same_AF", "outlier",
# Traits considered
list_traits <- c("growth","growth_pheno","repro_pheno","leafthickness","composite_fitness")
list_name_final_figure <- c("Growth","Growth Pheno","Reproductive Pheno","Leaf Thickness","Composite index")
# Create an empty dataframe to store correlations
lm_inf_df <- data.frame(
Method = character(),
Model = character(),
Trait = character(),
R2 = numeric(),
Intercept = numeric(),
Coeff = numeric(),
Pvalue = numeric(),
stringsAsFactors = FALSE
)
# Loop over each Genomic Offset for both RDA and GF models
for(i in 1:4) { # Now we have 8 models (4 RDA + 4 GF)
GO_pred_2012 <- get(list_GO_2012[i])
GO_pred_2021 <- get(list_GO_2021[i])
method <- list_method[i]
model <- list_model[i]
model_name <- list_name_models[i]
# Loop over each trait
for(x in 1:length(list_traits)) {
trait <- list_traits[x]
trait_name <- list_name_final_figure[x]
# Use CG2021 data for Height, and CG2012 for other traits
if (trait == "leafthickness" || trait== "composite_fitness") {
df_trait <- get(paste0("df_", trait,"_final"))
# Merge data frames based on population for CG2021
merged_df <- merge(df_trait, GO_pred_2021, by = "Population", suffixes = c(".trait", ".GO"))
} else {
df_trait <- get(paste0("df_", trait,"_final"))
# Merge data frames based on population for CG2012
merged_df <- merge(df_trait, GO_pred_2012, by = "Population", suffixes = c(".trait", ".GO"))
}
#linear model
lm_model <- lm(merged_df[, 2] ~ merged_df[, 3])
#assumptions:
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/all_models_period/Linear_models/Inferential/plot_model/",method,"/",model,"/diagnostic_plots/",method,"_",model,"_",trait,"_diagnostics.pdf"))
# Plot all four diagnostic plots at once in a 2x2 grid
par(mfrow = c(2, 2))
plot(lm_model)
dev.off()
# Reset plot settings
par(mfrow = c(1, 1))
coeff <- data.frame(coef(lm_model))[2,1]
result_model <- data.frame(anova(lm_model))
r_squared <- summary(lm_model)$r.squared
pvalues <- result_model[1,5]
Intercept <- data.frame(coef(lm_model))[1,1]
# Store the correlation in the dataframe
lm_inf_df <- rbind(lm_inf_df, data.frame(Method=method,Model = model_name, Trait = trait_name,R2=round(r_squared,3),Intercept= round(Intercept,3),Coeff = round(coeff,3),Pvalue=round(pvalues,3)))
# Create the data frame for plotting
plot_data <- data.frame(
GO_values = merged_df[, 3],
Blup = merged_df[, 2],
fitted = lm_model$fitted.values
)
# Generate the ggplot
p <- ggplot(plot_data, aes(x = GO_values, y = Blup)) +
geom_point(color = "blue") +
geom_smooth(method = "lm", se = TRUE, level = 0.95, color = "lightgrey", alpha = 0.2) + # Add confidence interval bands
geom_line(aes(y = fitted), color = "red") +
labs(title = paste("Linear Model:", trait, "~", method, model),
x = paste0("Genomic offset ",method," ",model),
y = paste0("proxy fitness ",trait)) +
theme_minimal()
# Print the plot (or use ggsave() to save each plot to a file)
#print(p)
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/all_models_period/Linear_models/Inferential/plot_model/",method,"/",model,"/plot_lm_",method,"_",model,"_",trait,"_Taxus.pdf"));print(p);dev.off()
}
}
# Print the correlation data frame
print(lm_inf_df)##
## Call:
## lm(formula = merged_df[, 2] ~ merged_df[, 3])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.22415 -0.07153 0.01335 0.07365 0.26696
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.62108 0.06604 9.405 3.64e-09 ***
## merged_df[, 3] -8.16882 3.05005 -2.678 0.0137 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1322 on 22 degrees of freedom
## Multiple R-squared: 0.2459, Adjusted R-squared: 0.2116
## F-statistic: 7.173 on 1 and 22 DF, p-value: 0.01373
We observed that for some traits and models, there is a linear negative relationship. For example, for the composite index and the genomic offset derived from the Gradient Forest (GF) model, using the random set with the same allele frequencies, we found a distinct negative relationship. However, the size of our dataset could be a limitation.
3.2.2 Bayesian model
Based on the observations from the inferential linear model, we also applied a Bayesian linear model. Indeed, the Bayesian framework is generally more effective at handling smaller datasets compared to frequentist approaches.
We found similar results with the Bayesian approach. Overall, the best associations align with the findings from the Pearson correlation: both the RDA and GF models appear to be negatively associated with all the traits, with the strength of the association varying depending on the model and the trait. Furthermore, the composite index seems to be more strongly negatively associated with genomic offset compared to other traits.
3.3 Non linear models
Due to the complexity of the relationship between the traits and the genomic offset predictions, we may not expect a linear relationship. To test this, we calculated non-linear models with quadratic parameters to capture any potential non-linear responses. We employed both inferential and Bayesian approaches for this analysis.
3.3.1 Inferential model
We performed the model using the nls function from stats package:
list_GO_2012_RDA <- c("GO_RDA_all_CG2012_Taxus","GO_RDA_random_CG2012_Taxus", "GO_RDA_random_same_AF_CG2012_Taxus",
"GO_RDA_LC_CG2012_Taxus")#, "GO_RDA_CG_CG2012_Taxus"
list_GO_2021_RDA <- c("GO_RDA_all_CG2021_Taxus","GO_RDA_random_CG2021_Taxus", "GO_RDA_random_same_AF_CG2021_Taxus",
"GO_RDA_LC_CG2021_Taxus")#, "GO_RDA_CG_CG2021_Taxus"
# List of Gradient Forest genomic offsets for 2012 and 2021
list_GO_2012_GF <- c("GO_standard_GF_all_2012","GO_standard_GF_random_2012", "GO_standard_GF_random_same_AF_2012",
"GO_standard_GF_LC_2012")#, "GO_standard_GF_CG_2012"
list_GO_2021_GF <- c("GO_standard_GF_all_2021","GO_standard_GF_random_2021", "GO_standard_GF_random_same_AF_2021",
"GO_standard_GF_LC_2021")#, "GO_standard_GF_CG_2021"
# Combine lists into one for easier processing
list_GO_2012 <- c(list_GO_2012_RDA, list_GO_2012_GF)
list_GO_2021 <- c(list_GO_2021_RDA, list_GO_2021_GF)
# Corresponding methods and models
list_method <- c(rep("RDA", 3), rep("GF", 3)) # 4 RDA models and 4 GF models
list_model <- c("all", "random_same_AF", "LC", "all", "random_same_AF", "LC")#, "CG", "CG",,"random","random",
list_name_models <- c("All", "Random","Outlier","All", "Random","Outlier")
# Traits considered
list_traits <- c("growth","growth_pheno","repro_pheno","leafthickness","composite_fitness")
list_name_final_figure <- c("Growth","Growth Pheno","Reproductive Pheno","Leaf Thickness","Composite index")
# Create an empty dataframe to store correlations
nl_infer_df_taxus <- data.frame(
Method = character(),
Model = character(),
Trait = character(),
R2 = numeric(),
coeff_b1 = numeric(),
b1_Pvalue = numeric(),
coeff_b2 = numeric(),
b2_Pvalue = numeric(),
stringsAsFactors = FALSE
)
# Loop over each Genomic Offset for both RDA and GF models
for(i in 1:6) { # Now we have 8 models (4 RDA + 4 GF)
GO_pred_2012 <- get(list_GO_2012[i])
GO_pred_2021 <- get(list_GO_2021[i])
method <- list_method[i]
model <- list_model[i]
model_name <- list_name_models[i]
# Loop over each trait
for(x in 1:length(list_traits)) {
trait <- list_traits[x]
trait_name <- list_name_final_figure[x]
# Use CG2021 data for Height, and CG2012 for other traits
if (trait == "leafthickness" || trait== "composite_fitness" ) {
df_trait <- get(paste0("df_", trait))
# Merge data frames based on population for CG2021
merged_df <- merge(df_trait, GO_pred_2021, by = "Population", suffixes = c(".trait", ".GO"))
} else {
df_trait <- get(paste0("df_", trait))
# Merge data frames based on population for CG2012
merged_df <- merge(df_trait, GO_pred_2012, by = "Population", suffixes = c(".trait", ".GO"))
}
y_data <- merged_df[,2]
x_data <- merged_df[,3]
df_nl_model <- data.frame(y=y_data,x=x_data)
infer_model <- stats::nls(y ~ b0+b1*x+b2*(x^2),data=df_nl_model,
start = list(b0=0.1,b1 = 0.1, b2 = .1),control = nls.control(minFactor = 1e-10, maxiter = 500))
#save the model:
assign(paste0("nl_infer_quadra_",method,"_",model,"_",trait),infer_model)
save(list = paste0("nl_infer_quadra_",method,"_",model,"_",trait),file=paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/all_models_period/Non_linear/Inferential/models/",method,"/",model,"/nl_infer_quadra_",method,"_",model,"_",trait,".Rdata"))
# Extract the estimates for b1 and b2
coef_b1 <- data.frame(coef(infer_model)["b1"])
coef_b2 <- data.frame(coef(infer_model)["b2"])
R2_model <- modelr::rsquare(infer_model, df_nl_model)
#pvalues
model_summary <- summary(infer_model)
p_values_b1 <- coef(model_summary)["b1", "Pr(>|t|)"]
p_values_b2 <- coef(model_summary)["b2", "Pr(>|t|)"]
nl_infer_df_taxus <- rbind(nl_infer_df_taxus, data.frame(
Method = method,
Model = model_name,
Trait = trait_name,
R2 = round(R2_model,3),
coeff_b1 = round(coef_b1[,1],3),
b1_Pvalue = round(p_values_b1,3),
coeff_b2 = round(coef_b2[,1],3),
b2_Pvalue = round(p_values_b2,3)
))
#confidence interval for the plot:
ci <- confint(infer_model)
print(ci)
# Create new data for prediction over a range of x values
new_data <- data.frame(x = seq(min(df_nl_model$x), max(df_nl_model$x), length.out = 100))
# Get the predicted values for the new data
new_data$y_pred <- predict(infer_model, newdata = new_data)
# Bootstrapping to calculate confidence intervals for predictions
# Define the prediction function for the bootstrapping
predict_nls <- function(data, indices) {
tryCatch({
fit <- nls(y ~ b1*x + b2*(x^2), data = data[indices, ], start = list(b1 = 0.1, b2 = 0.1))
predict(fit, newdata = new_data)
}, error = function(e) {
rep(NA, nrow(new_data)) # Return NA if fitting fails
})
}
boot_res <- boot::boot(
data = df_nl_model,
statistic = predict_nls,
R = 1000
)
# Filter out NA values before calculating quantiles
valid_bootstrap <- boot_res$t[complete.cases(boot_res$t), ]
# Calculate confidence intervals from bootstrapped predictions
ci_low <- apply(valid_bootstrap, 2, quantile, probs = 0.025, na.rm = TRUE)
ci_high <- apply(valid_bootstrap, 2, quantile, probs = 0.975, na.rm = TRUE)
# Add confidence intervals to the data
new_data$ci_low <- ci_low
new_data$ci_high <- ci_high
# Plotting the model with confidence intervals using ggplot2
plot <- ggplot(df_nl_model, aes(x = x, y = y)) +
geom_point(color = "blue") + # Observed data points
geom_line(data = new_data, aes(x = x, y = y_pred), color = "red") + # Fitted curve
geom_ribbon(data = new_data, aes(y = x, ymin = ci_low, ymax = ci_high), alpha = 0.2, fill = "lightgrey") + # Confidence interval band
labs(title = "Non-Linear Model with Confidence Intervals",
x = "x",
y = "y") +
theme_minimal()
plot
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/all_models_period/Non_linear/Inferential/figures/",method,"/",model,"/plot_nl_inf_quadra_",method,"_",model,"_",trait,"_Taxus.pdf"));plot(plot);dev.off()
}
}
# Print the correlation data frame
print(nl_infer_df_taxus)## 2.5% 97.5%
## b0 0.1564119 0.6605815
## b1 -3.5252380 9.2337728
## b2 -57.7261773 12.5217681
3.3.2 Bayesian model
list_GO_2012_RDA <- c("GO_RDA_all_CG2012_Taxus","GO_RDA_random_CG2012_Taxus", "GO_RDA_random_same_AF_CG2012_Taxus",
"GO_RDA_LC_CG2012_Taxus")#, "GO_RDA_CG_CG2012_Taxus"
list_GO_2021_RDA <- c("GO_RDA_all_CG2021_Taxus","GO_RDA_random_CG2021_Taxus", "GO_RDA_random_same_AF_CG2021_Taxus",
"GO_RDA_LC_CG2021_Taxus")#, "GO_RDA_CG_CG2021_Taxus"
# List of Gradient Forest genomic offsets for 2012 and 2021
list_GO_2012_GF <- c("GO_standard_GF_all_2012","GO_standard_GF_random_2012", "GO_standard_GF_random_same_AF_2012",
"GO_standard_GF_LC_2012")#, "GO_standard_GF_CG_2012"
list_GO_2021_GF <- c("GO_standard_GF_all_2021","GO_standard_GF_random_2021", "GO_standard_GF_random_same_AF_2021",
"GO_standard_GF_LC_2021")#, "GO_standard_GF_CG_2021"
# Combine lists into one for easier processing
list_GO_2012 <- c(list_GO_2012_RDA, list_GO_2012_GF)
list_GO_2021 <- c(list_GO_2021_RDA, list_GO_2021_GF)
# Corresponding methods and models
list_method <- c(rep("RDA", 3), rep("GF", 3)) # 4 RDA models and 4 GF models
list_model <- c("all","random_same_AF", "LC", "all", "random_same_AF", "LC")#, "CG", "CG","random", "random",
list_name_models <- c("All", "Random","Outlier","All", "Random","Outlier")
# Traits considered
list_traits <- c("growth","growth_pheno","repro_pheno","leafthickness","composite_fitness")
list_name_final_figure <- c("Growth","Growth Pheno","Reproductive Pheno","Leaf Thickness","Composite index")
# Create an empty dataframe to store correlations
df_nl_bayesian_taxus <- data.frame(
Method = character(),
Model = character(), # To store RDA name
Trait = character(), # To store trait name
R2_model =numeric(),
R2_ci_lower=numeric(),
R2_ci_upper=numeric(),
b1 = numeric(), # To store b1 coefficient
b2 = numeric(), # To store b2 coefficient
CI_b1 = numeric(),
CI_b2 = numeric(),
stringsAsFactors = FALSE
)
# Loop over each Genomic Offset for both RDA and GF models
for(i in 1:6) { # Now we have 8 models (4 RDA + 4 GF)
GO_pred_2012 <- get(list_GO_2012[i])
GO_pred_2021 <- get(list_GO_2021[i])
method <- list_method[i]
model <- list_model[i]
model_name <- list_name_models[i]
# Loop over each trait
for(x in 1:length(list_traits)) {
trait <- list_traits[x]
trait_name <- list_name_final_figure[x]
# Use CG2021 data for Height, and CG2012 for other traits
if (trait == "leafthickness" || trait== "composite_fitness" ) {
df_trait <- get(paste0("df_", trait,"_final"))
# Merge data frames based on population for CG2021
merged_df <- merge(df_trait, GO_pred_2021, by = "Population", suffixes = c(".trait", ".GO"))
} else {
df_trait <- get(paste0("df_", trait,"_final"))
# Merge data frames based on population for CG2012
merged_df <- merge(df_trait, GO_pred_2012, by = "Population", suffixes = c(".trait", ".GO"))
}
y_data <- merged_df[,2]
x_data <- merged_df[,3]
df_nl_model <- data.frame(y=y_data,x=x_data)
prior1 <- prior(normal(0, 10), nlpar = "b0")+prior(normal(0, 10), nlpar = "b1")+ prior(normal(0, 10), nlpar = "b2")
fit1 <- brm(bf(y ~ b0+b1*x+b2*(x^2), b0~1,b1 ~ 1,b2 ~ 1, nl = TRUE),
data = df_nl_model, prior = prior1, iter=1500000,warmup = 50000,thin=500, cores = 16)
#save the model:
#assign(paste0("nl_bayesian_quadra_",method,"_",model,"_",trait),fit1)
#save(list = paste0("nl_bayesian_quadra_",method,"_",model,"_",trait),file=paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/all_models_period/Non_linear/Bayesian/models/",method,"/nl_bayesian_quadra_",method,"_",model,"_",trait,".Rdata"))
load(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/all_models_period/Non_linear/Bayesian/models/",method,"/nl_bayesian_quadra_",method,"_",model,"_",trait,".Rdata"))
fit1 <- get(paste0("nl_bayesian_quadra_",method,"_",model,"_",trait))
#coeff
fit1_summary <- summary(fit1)
# Extract the estimates for b1 and b2
coef_b1 <- fit1_summary$fixed["b1_Intercept", "Estimate"]
coef_b2 <- fit1_summary$fixed["b2_Intercept", "Estimate"]
b1_ci_lower <- fit1_summary$fixed["b1_Intercept", "l-95% CI"]
b1_ci_upper <- fit1_summary$fixed["b1_Intercept", "u-95% CI"]
b2_ci_lower <- fit1_summary$fixed["b2_Intercept", "l-95% CI"]
b2_ci_upper <- fit1_summary$fixed["b2_Intercept", "u-95% CI"]
R2_model <- bayes_R2(fit1)["R2","Estimate"]
R2_ci_lower <- bayes_R2(fit1)["R2", "Q2.5"]
R2_ci_upper <- bayes_R2(fit1)["R2", "Q97.5"]
df_nl_bayesian_taxus <- rbind(df_nl_bayesian_taxus, data.frame(
Method = method,
Model = model_name,
Trait = trait_name,
R2= round(R2_model,3),
R2_ci_lower=round(R2_ci_lower,3),
R2_ci_upper=round(R2_ci_upper,3),
coeff_b1 = round(coef_b1,3),
CI_b1=paste0(round(b1_ci_lower,3),"_",round(b1_ci_upper,3)),
coeff_b2 = round(coef_b2,3),
CI_b2=paste0(round(b2_ci_lower,3),"_",round(b2_ci_upper,3))
))
#MCMC chains
plot(fit1)
#pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/all_models_period/Non_linear/Bayesian/figures/",method,"/",model,"/MCMC_nl_quadra_bayesian",trait,"_GO_",method,model,"_Taxus.pdf"));plot(fit1);dev.off()
#posterior distrib
pp_check(fit1)
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/all_models_period/Non_linear/Bayesian/figures/",method,"/",model,"/post_distrib_nl_quadra_bayesian",trait,"_GO_",method,model,"_Taxus.pdf"));plot(pp_check(fit1));dev.off()
plot(conditional_effects(fit1), points = TRUE)
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/Validation/GO/all_models_period/Non_linear/Bayesian/figures/",method,"/",model,"/plot_nl_quadra_",trait,"_GO_",method,model,"_Taxus.pdf"));plot(conditional_effects(fit1), points = TRUE);dev.off()
}
}
# Print the correlation data frame
print(df_nl_bayesian_taxus)Overall, we can observe that the relationships appear quite linear, as adding a quadratic parameter did not significantly change the R2 of the models. This suggests that the relationships may indeed be linear; therefore, incorporating a quadratic coefficient did not substantially impact the results as long as a linear coefficient remained in the model. Moreover, Bayesian and inferential models converged for those models with the highest negative associations. In conclusion, our dataset reveals a negative association between fitness proxies and genomic offset predictions. This suggests that genomic offset may provide valuable insights into the maladaptation of English Yew populations to climate conditions.
3.4 non-parametric tests
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.3.2 (2023-10-31 ucrt)
## os Windows 11 x64 (build 26100)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate French_France.utf8
## ctype French_France.utf8
## tz Europe/Paris
## date 2025-09-22
## pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## ! package * version date (UTC) lib source
## abind 1.4-8 2024-09-12 [1] CRAN (R 4.3.3)
## backports 1.4.1 2021-12-13 [1] CRAN (R 4.3.1)
## bayesplot 1.11.1 2024-02-15 [1] CRAN (R 4.3.3)
## boot * 1.3-30 2024-02-26 [2] CRAN (R 4.3.2)
## bridgesampling 1.1-2 2021-04-16 [1] CRAN (R 4.3.3)
## brms * 2.21.0 2024-03-20 [1] CRAN (R 4.3.3)
## Brobdingnag 1.2-9 2022-10-19 [1] CRAN (R 4.3.3)
## broom 1.0.8 2025-03-28 [1] CRAN (R 4.3.3)
## bslib 0.6.1 2023-11-28 [1] CRAN (R 4.3.2)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.1)
## checkmate 2.3.2 2024-07-29 [1] CRAN (R 4.3.3)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.1)
## coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.3.3)
## codetools 0.2-19 2023-02-01 [2] CRAN (R 4.3.2)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.1)
## corrplot * 0.92 2021-11-18 [1] CRAN (R 4.3.2)
## curl 5.2.0 2023-12-08 [1] CRAN (R 4.3.2)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.3.3)
## digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.1)
## distributional 0.5.0 2024-09-17 [1] CRAN (R 4.3.3)
## dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.2)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.1)
## emmeans 1.10.0 2024-01-23 [1] CRAN (R 4.3.2)
## estimability 1.5 2024-02-20 [1] CRAN (R 4.3.2)
## evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.2)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.3.3)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.1)
## fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.1)
## generics 0.1.4 2025-05-09 [1] CRAN (R 4.3.2)
## ggcorrplot * 0.1.4.1 2023-09-05 [1] CRAN (R 4.3.3)
## ggplot2 * 3.5.2 2025-04-09 [1] CRAN (R 4.3.2)
## glue 1.7.0 2024-01-09 [1] CRAN (R 4.3.2)
## gridExtra * 2.3 2017-09-09 [1] CRAN (R 4.3.3)
## gtable 0.3.6 2024-10-25 [1] CRAN (R 4.3.3)
## highr 0.10 2022-12-22 [1] CRAN (R 4.3.1)
## htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.3.2)
## htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.2)
## httpuv 1.6.13 2023-12-06 [1] CRAN (R 4.3.2)
## inline 0.3.19 2021-05-31 [1] CRAN (R 4.3.3)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.1)
## jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.3.2)
## knitr 1.45 2023-10-30 [1] CRAN (R 4.3.2)
## labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.1)
## later 1.3.2 2023-12-06 [1] CRAN (R 4.3.2)
## lattice 0.22-5 2023-10-24 [2] CRAN (R 4.3.2)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.2)
## loo 2.8.0 2024-07-03 [1] CRAN (R 4.3.3)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.1)
## MASS 7.3-60.0.1 2024-01-13 [2] CRAN (R 4.3.2)
## Matrix 1.6-5 2024-01-11 [1] CRAN (R 4.3.2)
## matrixStats 1.4.1 2024-09-08 [1] CRAN (R 4.3.3)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.1)
## mgcv 1.9-1 2023-12-21 [2] CRAN (R 4.3.2)
## mime 0.12 2021-09-28 [1] CRAN (R 4.3.1)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.3.2)
## modelr 0.1.11 2023-03-22 [1] CRAN (R 4.3.2)
## multcomp 1.4-25 2023-06-20 [1] CRAN (R 4.3.3)
## munsell 0.5.1 2024-04-01 [1] CRAN (R 4.3.3)
## mvtnorm 1.3-1 2024-09-03 [1] CRAN (R 4.3.3)
## nlme 3.1-164 2023-11-27 [2] CRAN (R 4.3.2)
## pillar 1.10.2 2025-04-05 [1] CRAN (R 4.3.3)
## pkgbuild 1.4.3 2023-12-10 [1] CRAN (R 4.3.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.2)
## pkgload 1.3.4 2024-01-16 [1] CRAN (R 4.3.2)
## plyr 1.8.9 2023-10-02 [1] CRAN (R 4.3.2)
## posterior 1.6.0 2024-07-03 [1] CRAN (R 4.3.3)
## profvis 0.3.8 2023-05-02 [1] CRAN (R 4.3.2)
## promises 1.2.1 2023-08-10 [1] CRAN (R 4.3.2)
## purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.3.2)
## QuickJSR 1.3.1 2024-07-14 [1] CRAN (R 4.3.3)
## R6 2.6.1 2025-02-15 [1] CRAN (R 4.3.3)
## Rcpp * 1.0.11 2023-07-06 [1] CRAN (R 4.3.1)
## D RcppParallel 5.1.7 2023-02-27 [1] CRAN (R 4.3.2)
## remotes 2.5.0 2024-03-17 [1] CRAN (R 4.3.3)
## reshape2 * 1.4.4 2020-04-09 [1] CRAN (R 4.3.2)
## rlang 1.1.3 2024-01-10 [1] CRAN (R 4.3.2)
## rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.1)
## rstan 2.32.6 2024-03-05 [1] CRAN (R 4.3.3)
## rstantools 2.4.0 2024-01-31 [1] CRAN (R 4.3.3)
## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.2)
## sandwich 3.1-0 2023-12-11 [1] CRAN (R 4.3.3)
## sass 0.4.8 2023-12-06 [1] CRAN (R 4.3.2)
## scales 1.3.0 2023-11-28 [1] CRAN (R 4.3.2)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.2)
## shiny 1.8.0 2023-11-17 [1] CRAN (R 4.3.2)
## StanHeaders 2.32.10 2024-07-15 [1] CRAN (R 4.3.3)
## stringi 1.8.3 2023-12-11 [1] CRAN (R 4.3.2)
## stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.3.2)
## survival 3.5-8 2024-02-14 [2] CRAN (R 4.3.2)
## tensorA 0.36.2.1 2023-12-13 [1] CRAN (R 4.3.2)
## TH.data 1.1-2 2023-04-17 [1] CRAN (R 4.3.3)
## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.2)
## tidyr 1.3.1 2024-01-24 [1] CRAN (R 4.3.2)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.3.3)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.3.2)
## usethis 2.2.3 2024-02-19 [1] CRAN (R 4.3.2)
## utf8 1.2.4 2023-10-22 [1] CRAN (R 4.3.2)
## V8 5.0.0 2024-08-16 [1] CRAN (R 4.3.3)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.2)
## withr 3.0.2 2024-10-28 [1] CRAN (R 4.3.3)
## xfun 0.41 2023-11-01 [1] CRAN (R 4.3.2)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.2)
## yaml 2.3.8 2023-12-11 [1] CRAN (R 4.3.2)
## zoo 1.8-12 2023-04-13 [1] CRAN (R 4.3.3)
##
## [1] C:/Users/tfrancisco/AppData/Local/R/win-library/4.3
## [2] C:/Users/tfrancisco/AppData/Local/Programs/R/R-4.3.2/library
##
## D ── DLL MD5 mismatch, broken installation.
##
## ──────────────────────────────────────────────────────────────────────────────