### Libraries and functions ----
library(rstan)
library(here)
library(bayesplot)
library(stringr)
library(here)
library(rstan)
library(ggplot2)
library(dplyr)
library(kableExtra)
library(tidyr)
library(knitr)
library(RColorBrewer)
library(readxl)
library(visreg)
library(emmeans)
library(car)
options(mc.cores = parallel::detectCores())
completeFun <- function(data, desiredCols) {
completeVec <- complete.cases(data[, desiredCols])
return(data[completeVec, ])
}
###########################################
########### Trait Divergence ##############
###########################################
Datos2<-read_excel(here("Functional_traits_analysis","Traits_Data.xlsx"))
Datos2$PLOIDY <- as.factor(Datos2$PLOIDY)
Datos2$TIME <- as.factor(Datos2$TIME)
Datos2$TREATMENT <- as.factor(Datos2$TREATMENT)
Datos2$TIPO <- as.factor(Datos2$TIPO)
Datos2$Fluorimetro<- as.numeric(Datos2$Fluorimetro)
str(Datos2)
contrasts(Datos2$PLOIDY)<-contr.sum
contrasts(Datos2$TIPO)<-contr.sum
contrasts(Datos2$TREATMENT)<-contr.sum
# ANOVA and Tukey test for each variable
A2 <- aov(AN ~ TIPO + TREATMENT + PLOIDY + TIPO*TREATMENT + TREATMENT*PLOIDY + PLOIDY*TIPO, data = Datos2, na.action = na.omit)
summary(A2)
C2 <- aov(gsw ~ TIPO + TREATMENT + PLOIDY + TIPO*TREATMENT + TREATMENT*PLOIDY + PLOIDY*TIPO, data = Datos2, na.action = na.omit)
summary(C2)
H2 <- aov(SLA~ TIPO + TREATMENT + PLOIDY + TIPO*TREATMENT + TREATMENT*PLOIDY + PLOIDY*TIPO, na.action= na.omit,data = Datos2)
Anova(H2, type=2)
I2 <- aov(LDMC ~ TIPO + TREATMENT + PLOIDY + TIPO*TREATMENT + TREATMENT*PLOIDY + PLOIDY*TIPO, data = Datos2, na.action = na.omit)
Anova(I2, type=2)
J2 <- aov(Fluorimetro ~ TIPO + TREATMENT + PLOIDY + TIPO*TREATMENT + TREATMENT*PLOIDY + PLOIDY*TIPO, data = Datos2, na.action = na.omit)
Anova(J2, type=2)
K2 <- aov(WUEi ~ TIPO + TREATMENT + PLOIDY + TIPO*TREATMENT + TREATMENT*PLOIDY + PLOIDY*TIPO, data = Datos2, na.action = na.omit)
Anova(K2, type=2)
# Extract ANOVA summaries without the residual row.
summary_A <- summary(A2)[[1]][-7,]
summary_gsw <- summary(C2)[[1]][-7,]
summary_SLA <- summary(H2)[[1]][-7,]
summary_LDMC <- summary(I2)[[1]][-7,]
summary_FvFm <- summary(J2)[[1]][-7,]
summary_WUEi <- summary(K2)[[1]][-7,]
summary_list <- list(summary_A, summary_gsw, summary_SLA, summary_LDMC, summary_FvFm, summary_WUEi)
# Factor and variable names
factor_names <- c("TYPE", "TREATMENT", "PLOIDY", "TYPE:TREATMENT", "TREATMENT:PLOIDY", "TYPE:PLOIDY")
column_names <- c("A", "gsw", "SLA", "LDMC", "FvFm", "WUEi")
# Create data frame
combinations <- expand.grid(Factor = factor_names, Variable = column_names)
f_values_df <- data.frame(combinations, F_value = NA, Df = NA, stringsAsFactors = FALSE)
for (i in seq_along(summary_list)) {
for (j in seq_along(factor_names)) {
f_value <- summary_list[[i]]$"F value"[j]
p_value <- summary_list[[i]]$'Pr(>F)'[j]
df_value <- summary_list[[i]]$Df[j] # Extraer los grados de libertad
if (is.numeric(p_value) && length(p_value) > 0 && !is.na(p_value)) {
if (p_value < 0.05) {
f_values_df[(i - 1) * length(factor_names) + j, 3] <- paste0(format(round(f_value, 4), nsmall = 4), "*")
} else {
f_values_df[(i - 1) * length(factor_names) + j, 3] <- format(round(f_value, 4), nsmall = 4)
}
f_values_df[(i - 1) * length(factor_names) + j, 4] <- df_value # A?adir los grados de libertad
} else {
f_values_df[(i - 1) * length(factor_names) + j, 3] <- NA
f_values_df[(i - 1) * length(factor_names) + j, 4] <- NA
}
}
}
f_values_df_transformed <- spread(f_values_df, Variable, F_value)
# Create final table
table <- kable(f_values_df_transformed, format = "html", align = "c") %>%
kable_styling(bootstrap_options = "striped", latex_options = "hold_position", full_width = FALSE)
print(table)
###########################################
## MODEL STAN: "AN" WITH WATER TREATMENT ###
###########################################
datacomp<-read.table(here("Functional_traits_analysis", "input","Traits_data.txt"),header = T, stringsAsFactors = T)
#Prepare data
datacomp$PLOIDY<-as.numeric(datacomp$PLOIDY)
datacomp$PLOIDY<-datacomp$PLOIDY*-1
datacomp$PLOIDY[datacomp$PLOIDY==-1]<-4
datacomp$PLOIDY[datacomp$PLOIDY==-2]<-1
datacomp$PLOIDY[datacomp$PLOIDY==-3]<-2
datacomp$PLOIDY[datacomp$PLOIDY==-4]<-3
datacomp$TREATMENT<-as.numeric(datacomp$TREATMENT)
datacomp$TREATMENT[datacomp$TREATMENT==2]<-0
datacomp2<-scale(datacomp[,-c(1,2,3,4,5,6,7,8,9,11,12)])
datacomp3<-cbind(datacomp[,c(1,2,3,4,5,6,7,8,9,11,12)],datacomp2)
#View(datacomp3)
datacomp3$PLOIDY<-as.numeric(datacomp3$PLOIDY)
datacomp3<- completeFun(datacomp3, c("Bwf","Bwc", "A_f", "A_c"))
Biomasa<-datacomp3$Bwf
tf<-datacomp3$A_f
tc<-datacomp3$A_c
species<-datacomp3$PLOIDY
C<-datacomp3$C
N_species <- 4
Bc<-datacomp3$Bwc
water_status<-datacomp3$TREATMENT
N<- nrow(datacomp3) #Number of observations
# Create a list
data_list <- list(N = N,
Biomasa = Biomasa,
water_status = water_status,
tf = tf,
tc = tc,
Bc = Bc,
species = species,
C = C,
N_species = N_species
)
# Compile the model
stan_model_water <- stan_model(here("Functional_traits_analysis","script", "stan_models",
"Waterstatus_model.stan"))
# Adjusts the model using the generated data
fitA_water <- sampling(stan_model_water, data = data_list,
init = function() list(m0 = 0, m1 = 0,
sigma = 1, alpha_t = 0, alpha_e = 0, alpha_d = 0),
iter = 10000, chains = 5)
##saveRDS(fitA_water, file = here("Functional_traits_analysis","output", "stan_results", "fitA_water.rds"))
#Results
A1<-print(fitA_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ))
plot(fitA_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ), ci_level = 0.95, include = T)
###########################################
# MODEL STAN: "gsw" WITH WATER TREATMENT ##
###########################################
datacomp3<- completeFun(datacomp3, c("Bwf","Bwc", "gsw_f", "gsw_c"))
Biomasa<-datacomp3$Bwf
tf<-datacomp3$gsw_f
tc<-datacomp3$gsw_c
species<-datacomp3$PLOIDY
C<-datacomp3$C
N_species <- 4
Bc<-datacomp3$Bwc
water_status<-datacomp3$TREATMENT
N<- nrow(datacomp3) #Number of observations
# Create a list
data_list <- list(N = N,
Biomasa = Biomasa,
water_status = water_status,
tf = tf,
tc = tc,
Bc = Bc,
species = species,
C = C,
N_species = N_species
)
# Compile the stan model
stan_model_water <- stan_model(here("Functional_traits_analysis","script", "stan_models",
"Waterstatus_model.stan"))
# Adjust the model using the generated data
fitgsw_water <- sampling(stan_model_water, data = data_list,
init = function() list(m0 = 0, m1 = 0,
sigma = 1, alpha_t = 0, alpha_e = 0, alpha_d = 0),
iter = 10000, chains = 5)
#saveRDS(fitgsw_water, file = here("Functional_traits_analysis","output", "stan_results", "fitgsw_water.rds"))
# Results
A1<-print(fitgsw_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ))
plot(fitgsw_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ), ci_level = 0.95, include = T)
###########################################
# MODEL STAN: "WUE" WITH WATER TREATMENT ##
###########################################
datacomp3<- completeFun(datacomp3, c("Bwf","Bwc", "WUEi_f", "WUEi_c"))
Biomasa<-datacomp3$Bwf
tf<-datacomp3$WUEi_f
tc<-datacomp3$WUEi_c
species<-datacomp3$PLOIDY
C<-datacomp3$C
N_species <- 4
Bc<-datacomp3$Bwc
water_status<-datacomp3$TREATMENT
N<- nrow(datacomp3) # Observations
# Create a list
data_list <- list(N = N,
Biomasa = Biomasa,
water_status = water_status,
tf = tf,
tc = tc,
Bc = Bc,
species = species,
C = C,
N_species = N_species
)
# Compile the stan model
stan_model_water <- stan_model(here("Functional_traits_analysis","script", "stan_models",
"Waterstatus_model.stan"))
# Adjust the model using the generated data
fitWUEi_water <- sampling(stan_model_water, data = data_list,
init = function() list(m0 = 0, m1 = 0,
sigma = 1, alpha_t = 0, alpha_e = 0, alpha_d = 0),
iter = 10000, chains = 8)
#saveRDS(fitWUEi_water, file = here("Functional_traits_analysis","output", "stan_results", "fitWUEi_water.rds"))
# Results
A1<-print(fitWUEi_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ))
plot(fitWUEi_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ), ci_level = 0.95, include = T)
###########################################
# MODEL STAN: "FvFm" WITH WATER TREATMENT #
###########################################
datacomp3<- completeFun(datacomp3, c("Bwf","Bwc", "Fluorimetro_f", "Fluorimetro_c"))
Biomasa<-datacomp3$Bwf
tf<-datacomp3$Fluorimetro_f
tc<-datacomp3$Fluorimetro_c
species<-datacomp3$PLOIDY
C<-datacomp3$C
N_species <- 4
Bc<-datacomp3$Bwc
water_status<-datacomp3$TREATMENT
N<- nrow(datacomp3) # Observations
# Create a list
data_list <- list(N = N,
Biomasa = Biomasa,
water_status = water_status,
tf = tf,
tc = tc,
Bc = Bc,
species = species,
C = C,
N_species = N_species
)
# Compile the stan model
stan_model_water <- stan_model(here("Functional_traits_analysis","script", "stan_models",
"Waterstatus_model.stan"))
# Adjust the model using the generated data
fitfvfm_water <- sampling(stan_model_water, data = data_list,
init = function() list(m0 = 0, m1 = 0,
sigma = 1, alpha_t = 0, alpha_e = 0, alpha_d = 0),
iter = 10000, chains = 4)
#saveRDS(fitfvfm_water, file = here("Functional_traits_analysis","output", "stan_results", "fitfvfm_water.rds"))
# Results
A1<-print(fitfvfm_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ))
plot(fitfvfm_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ), ci_level = 0.95, include = T)
###########################################
# MODEL STAN: "LDMC" WITH WATER TREATMENT #
###########################################
datacomp3<- completeFun(datacomp3, c("Bwf","Bwc", "LDMC_f", "LDMC_c"))
Biomasa<-datacomp3$Bwf
tf<-datacomp3$LDMC_f
tc<-datacomp3$LDMC_c
species<-datacomp3$PLOIDY
C<-datacomp3$C
N_species <- 4
Bc<-datacomp3$Bwc
water_status<-datacomp3$TREATMENT
N<- nrow(datacomp3) # Observations
# Create a list
data_list <- list(N = N,
Biomasa = Biomasa,
water_status = water_status,
tf = tf,
tc = tc,
Bc = Bc,
species = species,
C = C,
N_species = N_species
)
# Compile the stan model
stan_model_water <- stan_model(here("Functional_traits_analysis","script", "stan_models",
"Waterstatus_model.stan"))
# Adjust the model using the generated data
fitLDMC_water <- sampling(stan_model_water, data = data_list,
init = function() list(m0 = 0, m1 = 0,
sigma = 1, alpha_t = 0, alpha_e = 0, alpha_d = 0),
iter = 10000, chains = 8)
#saveRDS(fitLDMC_water, file = here("Functional_traits_analysis","output", "stan_results", "fitLDMC_water.rds"))
# Results
A1<-print(fitLDMC_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ))
plot(fitLDMC_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ), ci_level = 0.95, include = T)
###########################################
## MODEL STAN: "SLA" WITH WATER TREATMENT #
###########################################
datacomp3<- completeFun(datacomp3, c("Bwf","Bwc", "SLA_f", "SLA_c"))
Biomasa<-datacomp3$Bwf
tf<-datacomp3$SLA_f
tc<-datacomp3$SLA_c
species<-datacomp3$PLOIDY
C<-datacomp3$C
N_species <- 4
Bc<-datacomp3$Bwc
water_status<-datacomp3$TREATMENT
N<- nrow(datacomp3) # Observations
# Create a list
data_list <- list(N = N,
Biomasa = Biomasa,
water_status = water_status,
tf = tf,
tc = tc,
Bc = Bc,
species = species,
C = C,
N_species = N_species
)
# Compile the stan model
stan_model_water <- stan_model(here("Functional_traits_analysis","script", "stan_models",
"Waterstatus_model.stan"))
# Adjust the model using the generated data
fitSLA_water <- sampling(stan_model_water, data = data_list,
init = function() list(m0 = 0, m1 = 0,
sigma = 1, alpha_t = 0, alpha_e = 0, alpha_d = 0),
iter = 10000, chains = 8)
#saveRDS(fitSLA_water, file = here("Functional_traits_analysis","output", "stan_results", "fitSLA_water.rds"))
# Results
A1<-print(fitSLA_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ))
plot(fitSLA_water, pars=c("m1_dry","m1_wet",
"alpha_t_dry", "alpha_t_wet",
"alpha_e_dry","alpha_e_wet",
"alpha_d_dry","alpha_d_wet" ), ci_level = 0.95, include = T)
###########################################
########## MERGE STAN RESULTS #############
###########################################
substrRight <- function(x, n){
substr(x, nchar(x)-n+1, nchar(x))
}
archivos_rds <- list.files(path = here("Functional_traits_analysis","output","stan_results"),pattern = "\\.rds$")
archivos_rds_path <- list.files(path = here("Functional_traits_analysis","output","stan_results"),pattern = "\\.rds$",full.names = T)
nombre_variable <- sub("\\.rds$", "", archivos_rds)
for(i in 1:length(archivos_rds_path)) {
nombre_variable <- sub("\\.rds$", "", archivos_rds)
datos <- readRDS(archivos_rds_path[i])
assign(nombre_variable[i], datos)
}
list_fit<-list(fitA,fitgsw, fitWUEi, fitfvfm, fitSLA,fitLDMC)
list_fit_water<-list(fitA_water, fitgsw_water , fitWUEi_water,
fitfvfm_water, fitSLA_water,fitLDMC_water)
traits<-c("A","gsw","WUEi","FvFm","SLA","LDMC")
# all models
reslistc<-list()
for (i in 1:length(list_fit)) {
reslistc[[i]] <- as.data.frame(summary(list_fit[[i]], pars=c("m1","alpha_t","alpha_e","alpha_d"),
probs=c(0.025, 0.50, 0.975))[[1]][,4:6])
reslistc[[i]]$parameter<-str_replace(str_replace(rownames(reslistc[[i]]),
pattern = "_dry", replacement = ""), pattern = "_wet", replacement = "")
reslistc[[i]]$water_status <- "all"
reslistc[[i]]$trait <- traits[i]
}
df_post <- do.call("rbind",reslistc)
# water_status models
reslist<-list()
for (i in 1:length(list_fit_water)) {
reslist[[i]] <- as.data.frame(summary(list_fit_water[[i]], pars=c("m1_dry", "m1_wet",
"alpha_t_dry", "alpha_t_wet", "alpha_e_dry",
"alpha_e_wet", "alpha_d_dry",
"alpha_d_wet"), probs=c(0.025, 0.50, 0.975))[[1]][,4:6])
reslist[[i]]$parameter<-str_replace(str_replace(rownames(reslist[[i]]),
pattern = "_dry", replacement = ""), pattern = "_wet", replacement = "")
reslist[[i]]$water_status <- substrRight(rownames(reslist[[i]]), n=3)
reslist[[i]]$trait <- traits[i]
}
df_water <- do.call("rbind",reslist)
res<-rbind(df_post,df_water)
###########################################
################## GRAPH ##################
###########################################
colores <- c("red", "blue", "black")
# we assign these colors to the water_status labels
names(res)[1:3]<-c("min", "mean","max")
res$water_status <- factor(res$water_status, levels = c("dry", "wet", "all"))
col_pal<-c("black",brewer.pal(n = 3, name = "Set1")[1:2])
# Create graph
p <- ggplot(res, aes(x = parameter, y = mean, ymin = min, ymax = max, color = water_status)) +
geom_hline(yintercept = 0, linetype = "dashed", color="grey55") +
geom_pointrange(position = position_dodge(width = 0.6), fatten = 5) +
coord_flip() +
# ylim(-1,1) +
facet_wrap(~trait, nrow=1,scales = "free_x") +
theme_bw() +
theme(axis.text.y = element_text(size=15),
axis.text.x = element_text(angle = 0, hjust = 0.5,size=9),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 24),
strip.text = element_text(size = 15,hjust=0.75,face= "bold"),
strip.background = element_blank()) +
labs(x = "", y = "Standardized coefficients") +
scale_color_manual(values = c("wet" = col_pal[3], "dry" = col_pal[2], "all" = col_pal[1]))+
guides(color = guide_legend(reverse=TRUE))
print(p)