Phenotypic_data_clonal_bank
Thomas Francisco
2024-09-16
1 Introduction
The goal of this script is to estimate BLUP values for the phenotypic data, which will be used to evaluate genomic offset predictions. BLUP estimates allow us to generate population-level phenotypic values for each trait while accounting for covariates such as genotype for clones, tree ID when multiple measurements per tree exist, and tree size to account for age in growth-related traits.
Depending on the trait, we will select a specific period for climate
data extraction. For all traits, the initial year is 1992 (the year of
establishment of the clonal bank), with the following end years:
- 2021 for leaf thickness
- 2012 for shoot volume, stem length, spring elongation, and total open
male strobili
2 Phenotypic data
We loaded the phenotypic traits associated with:
Growth -> Shoot volume Growth_pheno -> spring elongation Repro_pheno -> total open strobili leaf_traits -> leaf thickness
The first step is to load the functions used in the analyses.
clean.MCMC <- function(x) {
sols <- summary(x)$solutions ## pull out relevant info from model summary
Gcovs <- summary(x)$Gcovariances
Rcovs <- summary(x)$Rcovariances
fixed <- data.frame(row.names(sols), sols, row.names = NULL) ## convert to dataframes with the row.names as the first col
random <- data.frame(row.names(Gcovs), Gcovs, row.names = NULL)
residual <- data.frame(row.names(Rcovs), Rcovs, row.names = NULL)
names(fixed)[names(fixed) == "row.names.sols."] <- "variable" ## change the columns names to variable, so they all match
names(random)[names(random) == "row.names.Gcovs."] <- "variable"
names(residual)[names(residual) == "row.names.Rcovs."] <- "variable"
fixed$effect <- "fixed" ## add ID column for type of effect (fixed, random, residual)
random$effect <- "random"
residual$effect <- "residual"
modelTerms <- as.data.frame(bind_rows(fixed, random, residual)) # merge it all together
}
plot.estimates <- function(x) {
if (class(x) != "summary.mcmc")
x <- summary(x)
n <- dim(x$statistics)[1]
par(mar=c(2, 9, 4, 1))
plot(x$statistics[,1], n:1,
yaxt="n", ylab="",
xlim=range(x$quantiles)*1.2,
pch=19,
main="Posterior means and 95% credible intervals")
grid()
axis(2, at=n:1, rownames(x$statistics), las=2)
arrows(x$quantiles[,1], n:1, x$quantiles[,5], n:1, code=0)
abline(v=0, lty=2)
}
myFun_h2Qst1 <- function(mcmc){
out <- c(
posterior.mode((mcmc$VCV[,2])/(rowSums(mcmc$VCV))),
HPDinterval((mcmc$VCV[,2])/(rowSums(mcmc$VCV)))[1],
HPDinterval((mcmc$VCV[,2])/(rowSums(mcmc$VCV)))[2],
posterior.mode((mcmc$VCV[,2])/((mcmc$VCV[,2])+(mcmc$VCV[,3]))),
HPDinterval((mcmc$VCV[,2])/((mcmc$VCV[,2])+(mcmc$VCV[,3])))[1],
HPDinterval((mcmc$VCV[,2])/((mcmc$VCV[,2])+(mcmc$VCV[,3])))[2],
posterior.mode((mcmc$VCV[,1])/(rowSums(mcmc$VCV))),
HPDinterval((mcmc$VCV[,1])/(rowSums(mcmc$VCV)))[1],
HPDinterval((mcmc$VCV[,1])/(rowSums(mcmc$VCV)))[2],
posterior.mode((mcmc$VCV[,1])/(2*(mcmc$VCV[,2])+(mcmc$VCV[,1]))),
HPDinterval((mcmc$VCV[,1])/(2*(mcmc$VCV[,2])+(mcmc$VCV[,1])))[1],
HPDinterval((mcmc$VCV[,1])/(2*(mcmc$VCV[,2])+(mcmc$VCV[,1])))[2]
)
names(out) <- c("H2c_all","lo.ci.H2c_all","up.ci.H2c_all" ,"H2c","lo.ci.H2c","up.ci.H2c" ,"H2p","lo.ci.H2p","up.ci.H2p" ,"Qst","lo.ci.Qst","up.ci.Qst")
return(out)}
# for the logit link, which is used when family=categorical
myFun_h2Qst2 <- function(mcmc){
out <- c(
posterior.mode((mcmc$VCV[,2])/(rowSums(mcmc$VCV)+(pi^2)/3)),
HPDinterval((mcmc$VCV[,2])/(rowSums(mcmc$VCV)+(pi^2)/3))[1],
HPDinterval((mcmc$VCV[,2])/(rowSums(mcmc$VCV)+(pi^2)/3))[2],
posterior.mode((mcmc$VCV[,2])/((mcmc$VCV[,2])+(mcmc$VCV[,3])+(pi^2)/3)),
HPDinterval((mcmc$VCV[,2])/((mcmc$VCV[,2])+(mcmc$VCV[,3])+(pi^2)/3))[1],
HPDinterval((mcmc$VCV[,2])/((mcmc$VCV[,2])+(mcmc$VCV[,3])+(pi^2)/3))[2],
posterior.mode(((mcmc$VCV[,1]))/(rowSums(mcmc$VCV)+(pi^2)/3)),
HPDinterval(((mcmc$VCV[,1]))/(rowSums(mcmc$VCV)+(pi^2)/3))[1],
HPDinterval(((mcmc$VCV[,1]))/(rowSums(mcmc$VCV)+(pi^2)/3))[2],
posterior.mode((mcmc$VCV[,1])/(2*(mcmc$VCV[,2])+(mcmc$VCV[,1]))),
HPDinterval((mcmc$VCV[,1])/(2*(mcmc$VCV[,2])+(mcmc$VCV[,1])))[1],
HPDinterval((mcmc$VCV[,1])/(2*(mcmc$VCV[,2])+(mcmc$VCV[,1])))[2]
)
names(out) <- c("H2c_all","lo.ci.H2c_all","up.ci.H2c_all" ,"H2c","lo.ci.H2c","up.ci.H2c" ,"H2pc","lo.ci.H2pc","up.ci.H2pc" ,"Qst","lo.ci.Qst","up.ci.Qst")
return(out)}2.1 Data
list_traits <- c("Shoot_Growth","Basal_area","Shoot_Elongation","Male_Strobili","leaf_traits")
list_names <- c("growth","growth_bis","growth_pheno","repro_pheno","leaf_traits")
for(x in 1: length(list_traits)){
trait <- list_traits[x]
name <- list_names[x]
path <- paste0("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/phenotypic_data/common_garden/Individual_level_data/Taxus_",trait,".xlsx")
phenotypic_data <- read_excel(path, na="NA")
assign(paste0("phenotypic_data_",name),phenotypic_data)
}2.2 Selection of the period
For each trait, we calculated average values across years, as trait values in a given year are influenced by conditions and events from previous years. We retained the maximum number of sampling years while also maximizing the number of trees.
#Growth
##Shootvol
#Identify the number of measurement per plant
phenotypic_data_growth_table <- data.frame(table(phenotypic_data_growth$Plant))
#Filter the number of measurement to keep the maximum number of individuals
phenotypic_data_growth_filtered <- phenotypic_data_growth %>%
filter(Plant %in% phenotypic_data_growth_table$Var1[phenotypic_data_growth_table$Freq == 4])
#Mean values of the traits per plant
mean_traits_shootvol <- phenotypic_data_growth_filtered %>%
group_by(Plant) %>%
summarize(mean_AverageShootVolume = mean(AverageShootVolume.mm3, na.rm = TRUE),
mean_StemLength = mean(StemLength.cm, na.rm = TRUE))
mean_traits_shootvol_merge <- merge(mean_traits_shootvol,phenotypic_data_growth,"Plant")
#Keep only one value per plant (ex: 2012, it's not only the values of 2012)
mean_traits_growth_f <- mean_traits_shootvol_merge %>%
filter(Year == 2012)
##Basal area
mean_traits_growth_bis_f <- data.frame(phenotypic_data_growth_bis) %>%
filter(across(c(SumAB.2021.cm2, StemLength.cm), ~ !is.na(.)))# Growth phenology
##Spring elong
#identify the number of measurement per plant
phenotypic_data_phenology_table <- data.frame(table(phenotypic_data_growth_pheno$Plant))
#filter the number of measurement to keep the maximum number of individuals
phenotypic_data_pheno_filtered <- phenotypic_data_growth_pheno %>%
filter(Plant %in% phenotypic_data_phenology_table$Var1[phenotypic_data_phenology_table$Freq == 3])
#mean values of the traits per plant
mean_traits_pheno <- phenotypic_data_pheno_filtered %>%
group_by(Plant) %>%
summarize(mean_springelong = mean(ProportionSpringElongation, na.rm = TRUE))
mean_traits_pheno_merge <- merge(mean_traits_pheno,phenotypic_data_growth_pheno,"Plant")
#keep only one value per plant (ex: 2012, it's not only the values of 2012)
mean_traits_growth_pheno_f <- mean_traits_pheno_merge %>%
filter(Year == 2012)# Reproductive phenology
##Reproduction
#identify the number of measurement per plant
phenotypic_data_repro_table <- data.frame(table(phenotypic_data_repro_pheno$Plant))
#filter the number of measurement to keep the maximum number of individuals
phenotypic_data_repro_filtered <- phenotypic_data_repro_pheno %>%
filter(Plant %in% phenotypic_data_repro_table$Var1[phenotypic_data_repro_table$Freq == 3])
#mean values of the traits per plant
mean_traits_repro <- phenotypic_data_repro_filtered %>%
group_by(Plant) %>%
summarize(mean_proportion_open = mean(Proportion_total_open, na.rm = TRUE))
mean_traits_repro_merge <- merge(mean_traits_repro,phenotypic_data_repro_pheno,"Plant")
#keep only one value per plant (ex: 2012, it's not only the values of 2012)
mean_traits_repro_pheno_f <- mean_traits_repro_merge %>%
filter(Year == 2012)#Leaf_traits
phenotypic_data_leaftrait_table <- data.frame(table(phenotypic_data_leaf_traits$Plant))
#table(phenotypic_data_leaftrait_table$Freq)#we retained 4
#filter the number of measurement to keep the maximum number of individuals
phenotypic_data_leaftrait_filtered <- phenotypic_data_leaf_traits %>%
filter(Plant %in% phenotypic_data_leaftrait_table$Var1[phenotypic_data_leaftrait_table$Freq == 4])
#mean values of the traits per plant
mean_traits_leaftraits <- phenotypic_data_leaftrait_filtered %>%
group_by(Plant) %>%
summarize(mean_StomDens = mean(StomDens.mm2, na.rm = TRUE),
mean_leafthickness = mean(LeafThickness.mm, na.rm = TRUE),
mean_LMA = mean(LMA.gm2,na.rm=T),
mean_D13C = mean(D13C,na.rm=T))
mean_traits_leaftraits_merge <- merge(mean_traits_leaftraits,phenotypic_data_leaf_traits[,c(1:5)],"Plant")
#keep only one value per plant (ex: 2012, it's not only the values of 2012)
mean_traits_leaftraits_f <- mean_traits_leaftraits_merge %>%
filter(Year == "Y2021")2.3 Selection of the populations with enought trees
We also checked the number of measured individuals per population, retaining only populations with more than two individuals.
table(mean_traits_growth_f$Population)##
## ACEBEDO AGRES ALCOY ARAGUES ASPONTES
## 2 16 5 16 5
## BARRACO BEGONTE BOCAHUERGANO CANENCIA CERVANTES
## 10 6 4 8 16
## FANLO FOLGOSO JERTE LASENIA MONTAGUT
## 9 14 6 2 8
## MONTSERRAT NUNOMORAL PUEBLALILLO QUESADA RASCAFRIA
## 5 26 5 3 16
## RASQUERA SALES_LLIERCA SILES TORLA VILLABLINO
## 5 8 9 33 5
## VIMBODI
## 5#We removed the ACEBDO, LASENIA
table(mean_traits_growth_pheno_f$Population)##
## ACEBEDO AGRES ALCOY ARAGUES ASPONTES
## 2 8 2 16 4
## BARRACO BEGONTE BOCAHUERGANO CANENCIA CERVANTES
## 10 6 4 8 16
## FANLO FOLGOSO JERTE LASENIA MONTAGUT
## 9 14 6 2 8
## MONTSERRAT NUNOMORAL PUEBLALILLO QUESADA RASCAFRIA
## 5 14 4 2 16
## RASQUERA SALES_LLIERCA SILES TORLA VILLABLINO
## 3 8 3 14 4
## VIMBODI
## 5#We removed the ACEBDO, ALCOY, LASENIA Quesada pop
table(mean_traits_repro_pheno_f$Population)##
## AGRES ARAGUES BARRACO BOCAHUERGANO CANENCIA
## 3 9 5 1 5
## FANLO FOLGOSO JERTE LASENIA MONTAGUT
## 4 9 2 2 3
## MONTSERRAT NUNOMORAL PUEBLALILLO RASCAFRIA RASQUERA
## 3 5 2 8 2
## SALES_LLIERCA TORLA VIMBODI
## 4 6 3#We removed the BOCAHUERGANO, JERTE, LASENIA,PUEBLALILLO and RASQUERA
table(mean_traits_leaftraits_f$Population)##
## ACEBEDO AGRES ALCOY ARAGUES ASPONTES
## 2 8 2 16 4
## BARRACO BEGONTE BOCAHUERGANO CANENCIA CERVANTES
## 10 6 4 8 14
## FANLO FOLGOSO JERTE LASENIA MONTAGUT
## 8 14 6 2 8
## MONTSERRAT NUNOMORAL PUEBLA_LILLO QUESADA RASCAFRIA
## 5 14 4 2 15
## RASQUERA SALES_LLIERCA SILES TORLA VILLLABLINO
## 3 8 3 15 4
## VIMBODI
## 5#we can remove: ACEBEDO, ALCOY, LASENIA, and QUESADA2.4 Graphical visualization
To select the type of distribution used in the GLMM models, we investigated the shape of the data distribution for each trait.
list_traits <- c("growth","growth_bis","growth_pheno","repro_pheno","leaftraits")
for(x in 1: length(list_traits)){
name <- list_traits[x]
df <- get(paste0("mean_traits_",name,"_f"))
df$Population<- as.factor(df$Population)
df$Plant <- as.factor(df$Plant)
df$Genet <- as.factor(df$Genet)
assign(paste0("phenotypic_data_",name),df)
}hist(phenotypic_data_growth$mean_AverageShootVolume)
#hist(phenotypic_data_growth_bis$SumAB.2021.cm2)
hist(phenotypic_data_growth_pheno$mean_springelong)
hist(phenotypic_data_repro_pheno$mean_proportion_open)
#hist(phenotypic_data_leaftraits$mean_StomDens)
hist(phenotypic_data_leaftraits$mean_leafthickness)
#hist(phenotypic_data_leaftraits$mean_LMA)
#hist(phenotypic_data_leaftraits$mean_D13C)3 MCMC GLMM models
We can now run the MCMC GLMM models for each trait. We need to specify the priors, as we are working with Bayesian models. We chose a non-informative prior due to the lack of prior information. The number of iterations was set to 1,500,000, with a burn-in period of 50,000 and a thinning interval of 500, sampling every 500 iterations. Each model was run four times per trait (four MCMC chains) to ensure stability. All traits follow a Gaussian distribution with continuous data.
3.4 Water statues
3.4.1 Leaf thickness
load("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/mean_years/leafthickness_model_wo_plantID.rda")
leafthickness_model_wo_plantID <-leafthickness_model_wo_plantID4 Evaluation of the models
To evaluate the performance of the models, we investigated several aspects:
4.1 Autocorrelation within MCMC
First, we checked the autocorrelation between values within the same MCMC chain to ensure that our thinning interval was sufficient.
list_name<- c("Population","Genet")
list_model <- c("growth","growth_pheno","repro_pheno","leafthickness")#"growth_bis","LMA","StomDens","D13C"
for(i in 1:length(list_model)){
model_pheno <- get(paste0(list_model[i],"_model_wo_plantID"))
name_pheno <- list_model[i]
for(x in 1:2){
name <- list_name[x]
acf_list <- lapply(model_pheno, function(model) {
#Extract the MCMC samples for the random effects of provenance
mcmc_samples <- mcmc(model$VCV[, name])
#Calculate ACF
acf_values <- acf(as.numeric(mcmc_samples), plot = FALSE) #Convert to numeric if needed
return(acf_values)
})
# Plot the first ACF
plot(acf_list[[1]], main=paste0("Autocorrelation for ",name," Random Effect ",name_pheno), col="black", lwd=2)
colors <- c("red", "blue", "green") #Colors for the additional chains
for (i in 2:length(acf_list)) {
lines(acf_list[[i]]$lag, acf_list[[i]]$acf, col=colors[i-1], lwd=2)
}
#Add a legend to distinguish between the chains
legend("topright", legend=paste("Chain", 1:length(model_pheno)), col=c("black", colors), lwd=2)
plot
#Save
png(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/figures_wo_plantID/checkmodel/autocorrelation/Autocorrelation_period_",name,"_",name_pheno,".png"));plot(acf_list[[1]], main=paste0("Autocorrelation for ",name," Random Effect ",name_pheno), col="black", lwd=2)
colors <- c("red", "blue", "green")
for (i in 2:length(acf_list)) {
lines(acf_list[[i]]$lag, acf_list[[i]]$acf, col=colors[i-1], lwd=2)
}
legend("topright", legend=paste("Chain", 1:length(model_pheno)), col=c("black", colors), lwd=2);dev.off()
}
}4.2 Chains convergence
Next, we ensured that the MCMC chains were converging.
for(i in 1:length(list_model)){
model_pheno <- get(paste0(list_model[i],"_model_wo_plantID"))
name_pheno <- list_model[i]
for(x in 1: length(list_name)){
name <- list_name[x]
#Create MCMC lists for overall random effects (aggregated)
mcmc_list <- mcmc.list(lapply(model_pheno, function(model) {
#Aggregate random effects for provenance
mcmc(model$VCV[, name])
}))
par(mfrow=c(2,1), mar=c(4,4,2,1)) #Two plots: one for provenance, one for plant ID
#Plot trace plot for overall provenance random effect
plot(mcmc_list, auto.layout=FALSE, main=paste0("Trace Plot for Overall ",name," Random Effect ",name_pheno), col=c("red", "blue", "green", "black"))
#Save
png(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/figures_wo_plantID/checkmodel/convergence/MCMC_period_",name,"_",name_pheno,".png"));par(mfrow=c(2,1), mar=c(4,4,2,1))
plot(mcmc_list, auto.layout=FALSE, main=paste0("Trace Plot for Overall ",name," Random Effect ",name_pheno), col=c("red", "blue", "green", "black"));dev.off()
}
}Graphically, the chains appear to converge. We calculated the Gelman-Rubin index to further assess convergence.
list_model <- c("growth","growth_pheno","repro_pheno","leafthickness")#"growth_bis","LMA","StomDens","D13C"
gelman_rubin_results <- list()
for(i in 1:length(list_model)){
#Get the model object
model_pheno <- get(paste0(list_model[i], "_model_wo_plantID"))
name_pheno <- list_model[i]
#Combine the MCMC chains into a single mcmc.list
chaine_mc <- lapply(model_pheno, function(m) m$Sol)
chaine_mc <- do.call(mcmc.list, chaine_mc)
#Gelman-Rubin (GR) criterion to assess convergence across MCMC chains
GRcriterion_index <- gelman.diag(chaine_mc)
#Store the summary of the GR criterion in the list
gelman_rubin_results[[name_pheno]] <- list(
mpsrf = GRcriterion_index$mpsrf
)
}
gelman_rubin_results## $growth
## $growth$mpsrf
## [1] 1.047963
##
##
## $growth_pheno
## $growth_pheno$mpsrf
## [1] 1.01642
##
##
## $repro_pheno
## $repro_pheno$mpsrf
## [1] 1.008652
##
##
## $leafthickness
## $leafthickness$mpsrf
## [1] 1.018937We can see that the Gelman-Rubin values are below 1.1, which is the defined threshold for convergence.
4.3 Similarity of the posterior distribution across chains
We also examined the posterior distributions to assess the convergence of the MCMC results.
for(i in 1:length(list_model)){
model_pheno <- get(paste0(list_model[i],"_model_wo_plantID"))
name_pheno <- list_model[i]
for(x in 1:length(list_name)){
name <-list_name[x]
density_list <- lapply(model_pheno, function(model) {
#Extract the MCMC samples for the random effects of provenance
mcmc_samples <- as.numeric(mcmc(model$VCV[, name]))
#Calculate density
density_values <- density(mcmc_samples)
return(density_values)
})
#Limit of the plot
max_density_value <- max(sapply(density_list, function(d) max(d$y)))
#Plot the first density
plot(density_list[[1]], main=paste0("Posterior distribution of the ",name," Random Effect value ",name_pheno), col="black", lwd=2, xlab="Value", ylab="Density",ylim=c(0, max_density_value * 1.02))
#Overlay the remaining densities
colors <- c("red", "blue", "green") #Colors for the additional chains
for (i in 2:length(density_list)) {
lines(density_list[[i]], col=colors[i-1], lwd=2)
}
#Add a legend to distinguish between the chains
legend("topright", legend=paste("Chain", 1:length(model_pheno)), col=c("black", colors), lwd=2)
#Save
png(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/figures_wo_plantID/postdistrib/Posterior_distrib_period_",name,"_",name_pheno,".png"));plot(density_list[[1]], main=paste0("Posterior distribution of the ",name," Random Effect value ",name_pheno), col="black", lwd=2, xlab="Value", ylab="Density",ylim=c(0, max_density_value * 1.02))
colors <- c("red", "blue", "green") # Colors for the additional chains
for (i in 2:length(density_list)) {
lines(density_list[[i]], col=colors[i-1], lwd=2)
}
legend("topright", legend=paste("Chain", 1:length(model_pheno)), col=c("black", colors), lwd=2);dev.off()
}
}5 Blups values
Finally, we extracted the BLUP values from the models.
To calculate the BLUPs, we used the mean BLUP value across the four MCMC runs.
#list_model <- c("growth","growth_bis","growth_pheno","repro_pheno","leafthickness","LMA","StomDens","D13C")
list_model <- c("growth","growth_pheno","repro_pheno","leafthickness")#"growth_bis","LMA","StomDens"
for(i in 1:length(list_model)){
model_pheno <- get(paste0(list_model[i],"_model_wo_plantID"))
name_pheno <- list_model[i]
#BLUPs values
for(x in 1:4){
model <- model_pheno[[x]] #Extract the chains from the model
blups<-data.frame(posterior.mode(model$Sol))
assign(paste0("blups_",name_pheno,"_",x),blups)
}
}
for(i in 1:length(list_model)){
name_pheno <- list_model[i]
MC1 <- get(paste0("blups_",name_pheno,"_1"))
MC2 <- get(paste0("blups_",name_pheno,"_2"))
MC3 <- get(paste0("blups_",name_pheno,"_3"))
MC4 <- get(paste0("blups_",name_pheno,"_4"))
Blups_df <- data.frame(MC1,MC2,MC3,MC4)
Blups_df$mean_blup <- rowMeans(Blups_df)
write_xlsx(data.frame(row.names(Blups_df),Blups_df[,c(5)]),paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/BLUPs_mean/Blups_",name_pheno,"final_period.xlsx"))
}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
## ape * 5.7-1 2023-03-13 [1] CRAN (R 4.3.2)
## boot 1.3-30 2024-02-26 [2] CRAN (R 4.3.2)
## 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)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.3.2)
## 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)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.1)
## corpcor 1.6.10 2021-09-16 [1] CRAN (R 4.3.1)
## cubature 2.1.1 2024-07-14 [1] CRAN (R 4.3.3)
## 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)
## 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)
## evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.2)
## 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)
## 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)
## 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)
## 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)
## 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)
## lme4 * 1.1-37 2025-03-26 [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)
## MCMCglmm * 2.36 2024-05-06 [1] CRAN (R 4.3.3)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.1)
## 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)
## minqa 1.2.8 2024-08-17 [1] CRAN (R 4.3.3)
## munsell 0.5.1 2024-04-01 [1] CRAN (R 4.3.3)
## nlme 3.1-164 2023-11-27 [2] CRAN (R 4.3.2)
## nloptr 2.2.1 2025-03-17 [1] CRAN (R 4.3.3)
## pastecs * 1.4.2 2024-02-01 [1] CRAN (R 4.3.3)
## 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)
## 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)
## QGglmm * 0.7.4 2020-01-07 [1] CRAN (R 4.3.3)
## R6 2.6.1 2025-02-15 [1] CRAN (R 4.3.3)
## rbibutils 2.3 2024-10-04 [1] CRAN (R 4.3.3)
## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.1)
## Rdpack 2.6.4 2025-04-09 [1] CRAN (R 4.3.2)
## readxl * 1.4.3 2023-07-06 [1] CRAN (R 4.3.2)
## reformulas 0.4.1 2025-04-30 [1] CRAN (R 4.3.2)
## remotes 2.5.0 2024-03-17 [1] CRAN (R 4.3.3)
## 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)
## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.2)
## 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)
## stargazer * 5.2.3 2022-03-04 [1] CRAN (R 4.3.1)
## 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)
## tensorA 0.36.2.1 2023-12-13 [1] CRAN (R 4.3.2)
## tibble 3.2.1 2023-03-20 [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)
## 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)
## writexl * 1.5.0 2024-02-09 [1] CRAN (R 4.3.2)
## 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)
##
## [1] C:/Users/tfrancisco/AppData/Local/R/win-library/4.3
## [2] C:/Users/tfrancisco/AppData/Local/Programs/R/R-4.3.2/library
##
## ──────────────────────────────────────────────────────────────────────────────6 Draft
df_growth <- read_excel("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/BLUPs_mean/Blups_growthfinal_period.xlsx")
df_growth <- df_growth[c(2:27),]
df_growth_V2 <- read_excel("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/phenotypic_data/BLUPs_mean/Blups_growth_V2final_period.xlsx")
df_growth_V2 <- df_growth_V2[c(2:27),]
cor(df_growth$Blups_df...c.5..,df_growth_V2$Blups_df...c.5..)## [1] 0.3034145We retained the full posterior distribution information rather than using a single point estimate, as suggested by Reviewer 2 to keep error information.
6.1 Mean value of CFI after accounting for uncertainty in trait values
6.2 Phenotypic values not using BLUPs (not using population as a random factor), as suggested by Reviewer 2
As the reviewer suggested, “Is BLUP for population i better than a simple mean for a given population? If yes, please explain how BLUP improves further inferences.” We will test whether not using BLUPs provides similar results. To do this, we will run a model for each phenotypic trait accounting for covariates, then calculate the mean value for each individual grouped by population to obtain the phenotypic mean per population.
6.5 Water statues
6.5.1 Leaf thickness
Traits:
Basal area
LMA
StomDens
D13C