width=box_width,
axis.params=list(
axis="x", text="Adult-to-larvae\ntrophallaxis",
text.angle=0, hjust=0.5,
text.size=3, fontface="bold"
)
) +
scale_fill_manual(
values=c( "0" = "gray90", "1" = "#666600"),
labels=c("0" = "absent", "1" = "present"),
guide=guide_legend(title="Adult-to-larvae\ntrophallaxis", keywidth=0.5, keyheight=0.5, order=2)
)+
# Adult to adult trophallaxis
new_scale_fill() +
geom_fruit(
geom=geom_tile,
mapping=aes(fill=as.factor(AA_troph_freq)),
offset=offSet,
col = "white",
width=box_width,
axis.params=list(
axis="x", text="Adult-to-adult\ntrophallaxis",
text.angle=0, hjust=0.5,
text.size=3, fontface="bold"
)
) +
scale_fill_manual(
values=c( "0" = "gray90", "1" = "#FFCC66", "2" = "#CC0000"),
labels=c("0" = "never", "1" = "sometime", "2" = "frequent"),
guide=guide_legend(title="Adult-to-adult\ntrophallaxis", keywidth=0.5, keyheight=0.5, order=3)
)+
# NB of OG
geom_star(
mapping=aes(x = 225, size=n_OG),
fill = "grey20",
starshape=8,position="identity",starstroke=0
) + labs(size = "# Orthogroups\nfound in the crop") +
# species label
geom_tiplab(offset = 75, fontface = "italic", size = 3.5) +
theme(legend.position = c(0.1, 0.7))+
coord_cartesian(xlim = c(NA, 300))
p1
p1 + barplot + theme(axis.title.y = element_blank(), axis.text.y = element_blank())
barplot <- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Orthogoup Relative Abundance (iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid lines
legend.position = "none",
# Move x axis text to the top
axis.text.y = element_text(hjust = 0),  # Align text
# Move the title to the top
axis.title.y = element_text(angle = 0, vjust = 1),
# Adjust the position of the axis title
axis.title.y.left = element_text(margin = margin(r = 20))
) +
scale_fill_viridis_d(option = "plasma", guide = "none")
barplot
barplot <- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Orthogoup Relative Abundance (iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid lines
legend.position = "none",
axis.text.y = element_text(hjust = 0),  # Align text to the left
axis.title.y = element_text(angle = 0),  # Make title horizontal
axis.title.y.right = element_text(),     # Move title to the top
axis.text.y.right = element_text()       # Move text to the top
) +
scale_fill_viridis_d(option = "plasma", guide = "none")
barplot
barplot <- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Orthogoup Relative Abundance (iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid lines
legend.position = "none",
axis.text.x = element_text(hjust = 0),  # Align text to the left
axis.title.y = element_text(angle = 0),  # Make title horizontal
axis.title.y.right = element_text(),     # Move title to the top
axis.text.y.right = element_text()       # Move text to the top
) +
scale_fill_viridis_d(option = "plasma", guide = "none")
p1 + barplot + theme(axis.title.y = element_blank(), axis.text.y = element_blank())
barplot <- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Orthogoup Relative Abundance (iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid lines
legend.position = "none",
axis.text.x.bottom = element_blank(),  # Remove bottom x-axis text
axis.text.x.top = element_text(),      # Show top x-axis text
axis.title.x.bottom = element_blank(), # Remove bottom x-axis title
axis.title.x.top = element_text()      # Show top x-axis title
) +
scale_fill_viridis_d(option = "plasma", guide = "none")
p1 + barplot + theme(axis.title.y = element_blank(), axis.text.y = element_blank())
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Orthogoup Relative Abundance (iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_x_discrete(position = "top")
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Orthogoup Relative Abundance (iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_discrete(position = "top")
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Orthogoup Relative Abundance (iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none")
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Orthogoup Relative Abundance (iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_x_continuous(
position = "bottom",
guide = guide_axis(position = "top")
)
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Orthogoup Relative Abundance (iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(
position = "bottom",
guide = guide_axis(position = "top")
)
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title = element_text(size = 10)
panel.grid = element_blank(),  # Removes all grid line
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title = element_text(size = 10),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title = element_text(size = 15),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title = element_text(size = 15),
axis.text = element_text(size = 10)
panel.grid = element_blank(),  # Removes all grid line
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title = element_text(size = 15),
axis.text = element_text(size = 10),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title = element_text(size = 15),
axis.text.y = element_text(size = 10),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title = element_text(size = 15),
axis.text.y = element_text(size = 20),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title = element_text(size = 15),
axis.text.x = element_text(size = 20),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title.x = element_text(size = 15),
axis.text.x = element_text(size = 20),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title.x = element_text(size = 15),
axis.text.x = element_text(size = 13),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title.x = element_text(size = 15),
axis.text.x = element_text(size = 11),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
p1 + barplot + theme(axis.title.y = element_blank(), axis.text.y = element_blank())
# Plot
barplot<- ggplot(data_long, aes(x = species, y = Abundance, fill = OG)) +
geom_bar(stat = "identity") +
coord_flip() +  # Horizontal bars
labs(x = "Species", y = "Proportion of protein abundance (%iBAQ)", fill = "Orthogroup") +
theme_minimal() +
theme(
axis.title.x = element_text(size = 15),
axis.text.x = element_text(size = 10),
panel.grid = element_blank(),  # Removes all grid line
legend.position = "none") +
scale_fill_viridis_d(option = "plasma", guide = "none") +
scale_y_continuous(position = "bottom", guide = guide_axis(position = "top"))
barplot
p1 + barplot + theme(axis.title.y = element_blank(), axis.text.y = element_blank())
# Test model fit
test_evolutionary_model <- function(phy, trait) {
"
This function compare the model fit for nine evolutionary models
the input in a phylogentic tree with branch lengh and a continuous trait
the output is a table sumarizing the different model fits
"
require(stringr, geiger) # libraries
# prepare data
phy_pruned <- ape::drop.tip(phy = phy, setdiff(phy$tip.label, names(trait))) # Remove tips absent in the data
trait <- subset(trait, names(trait) %in% phy_pruned$tip.label)
# Calculate phylogenetic signals
PS.kappa <- phylosig(phy_pruned, trait, method="K", test=FALSE, nsim=100, se=NULL, start=NULL,control=list()) # the Kappa one
PS.lambda <- phylosig(phy_pruned, trait, method="lambda", test=FALSE, nsim=100, se=NULL, start=NULL,control=list()) # the Lambda one
# Run the nine models of evolution
model.BM <- geiger::fitContinuous(phy_pruned, trait, model="BM")
model.OU <- geiger::fitContinuous(phy_pruned, trait, model="OU")
model.EB <- geiger::fitContinuous(phy_pruned, trait, model="EB")
model.white <- geiger::fitContinuous(phy_pruned, trait, model="white")
model.rate_trend <- geiger::fitContinuous(phy_pruned, trait, model="rate_trend")
model.mean_trend <- geiger::fitContinuous(phy_pruned, trait, model="mean_trend")
model.lambda <- geiger::fitContinuous(phy_pruned, trait, model="lambda")
model.kappa <- geiger::fitContinuous(phy_pruned, trait, model="kappa")
model.delta <- geiger::fitContinuous(phy_pruned, trait, model="delta")
#Compile AIC and logL information
AIC <- setNames(c(model.BM$opt$aic,model.OU$opt$aic,model.EB$opt$aic,model.white$opt$aic,model.rate_trend$opt$aic,model.mean_trend$opt$aic,model.lambda$opt$aic,model.kappa$opt$aic,model.delta$opt$aic),
c('BM', 'OU', 'EB', 'WN', 'rate_trend', 'mean_trend', 'lambda', 'kappa', 'delta'))
LnL <- setNames(c(model.BM$opt$lnL,model.OU$opt$lnL,model.EB$opt$lnL,model.white$opt$lnL,model.rate_trend$opt$lnL,model.mean_trend$opt$lnL,model.lambda$opt$lnL,model.kappa$opt$lnL,model.delta$opt$lnL),
c('BM', 'OU', 'EB', 'WN', 'rate_trend', 'mean_trend', 'lambda', 'kappa', 'delta'))
# Calculate difference from the best model
AIC_comparison <- round(AIC - min(AIC),0)
LnL_comparison <- round(LnL - max(LnL),0)
model_comparison <- sprintf("%s(%s)", LnL_comparison, AIC_comparison)
# Get the sigsq in brownian motion model
sigsq <- model.BM$opt$sigsq
# Compile all data in a table
table <- setNames(c(model_comparison, round(PS.lambda$lambda,3), round(PS.kappa[1],3), sigsq, length(trait)),
c('BM', 'OU', 'EB', 'WN', 'rate_trend', 'mean_trend', 'lambda', 'kappa', 'delta', 'PS.lambda', 'PS.kappa', 'sigsq', 'n'))
result <- list(phy_pruned = phy_pruned, trait = trait, table = table)
return(result)
}
# Load data
source("C:/Users/matte/Desktop/cuticle_evolution/R_scripts/resources/load_data.R")
# Import the tree
backbone_tree <- read.tree("C:/Users/matte/Desktop/ant_utilities/phylogenetical_data/Economo et al 2018/Dryad_archive/Dryad_archive/backbone_trees/backbone_NCuniform_mcc.tre")
# Prepare the tree
backbone_tree$tip.label <- str_sub(backbone_tree$tip.label, 2, nchar(backbone_tree$tip.label) - 1) # Adjust tip labels format
backbone_tree$tip.label <- gsub("\\d$", "", backbone_tree$tip.label) # Remove numeric suffix from tip labels
backbone_tree <- update_phylo_names(backbone_tree, name) # Update names of tips according to Antcats taxonomic revision
foo(backbone_tree)
phy <- backbone_tree
trait <- cuticle_ratio_w
trait <- data_w_mean$ratio
# Filter rows based on presence in backbone_tree_pruned$tip.label
cuticle_ratio_w <- setNames(data_w_mean$ratio , data_w_mean$genus.species)
cuticle_ratio_q <- setNames(data_q_mean$ratio , data_q_mean$genus.species)
cuticle_ratio_m <- setNames(data_m_mean$ratio , data_m_mean$genus.species)
# Filter rows based on presence in backbone_tree_pruned$tip.label
cuticle_ratio_w <- setNames(data_w_mean$ratio , data_w_mean$genus.species)
cuticle_ratio_q <- setNames(data_q_mean$ratio , data_q_mean$genus.species)
cuticle_ratio_m <- setNames(data_m_mean$ratio , data_m_mean$genus.species)
cuticle_evolution_w <- test_evolutionary_model(backbone_tree, cuticle_ratio_w)
cuticle_evolution_q <- test_evolutionary_model(backbone_tree, cuticle_ratio_q)
cuticle_evolution_m <- test_evolutionary_model(backbone_tree, cuticle_ratio_m)
#compile all tables in one final
final <- rbind(cuticle_evolution_w$table,
cuticle_evolution_q$table,
cuticle_evolution_m$table)
print(final)
## Combine results from model 1 and 2
#####################################
pg1 <- summarise_pgls(pgls.r1)
pg2 <- summarise_pgls(pgls.r2)
# Usefull links
#https://www.mpcm-evolution.com/OPM/Chapter5_OPM/OPM_chap5.pdf
#https://lukejharmon.github.io/ilhabela/instruction/2015/07/03/PGLS/
#https://dwbapst.github.io/PaleoSoc_phylo_short_course_2019/articles/module_09_worked_PCM_example/module_08.html
################################
#### SET UP THE ENVIRONMENT ####
################################
# load packages
library(MASS)
library(nlme)
library(ape)
library(geiger)
library(phytools)
library(readxl)
library(tidyr)
library(dplyr)
library(tidyverse)
# Load personal function library and data
source("C:/Users/matte/Desktop/cuticle_evolution/R_scripts/resources/cuticle_00_personal_library.R")
source("C:/Users/matte/Desktop/cuticle_evolution/R_scripts/resources/load_data.R")
# Load the trees
temp <- data_w_mean; go_two_hundred_phylo = TRUE
source("C:/Users/matte/Desktop/cuticle_evolution/R_scripts/resources/load_trees.R")
dataset <- temp
# Compute the mean covariance matrix accross all trees
cov_mat <- compute_mean_vcv(two_hundred_phy_pruned)/200
# Load the trait data.R
source("C:/Users/matte/Desktop/cuticle_evolution/R_scripts/resources/load_trait_data.R")
head(dataset)
#############################
#### RUN PGLS ON CUTICLE ####
#############################
## PGLS 1
#########
dataset$diet <- factor(dataset$diet, levels = c('omnivore', 'predator', 'herbivore', "fungivore"))
dataset$foraging <- factor(dataset$foraging, levels = c('epigaeic', 'arboreal', 'hypogaeic'))
# Set the formula
formula <- log10(cuticle_volume) ~
mean_ann_prec + mean_ann_temp +
log10(colony_size) + spine_number +
foraging + diet +
log10(body_volume)
# prepare dataset and covariance matrix
res <- process_data_for_PGLS(formula, dataset, cov_mat)
plot_subfam_balance(dataset_balancing(res$dataset_pruned))
dataset_content(res$dataset_pruned)
# Run model 1
pgls.r1 <- gls(formula, correlation = corSymm(res$cov_mat_pruned[lower.tri(res$cov_mat_pruned)], fixed = TRUE), method = "ML", data = res$dataset_pruned)
# Usefull links
#https://www.mpcm-evolution.com/OPM/Chapter5_OPM/OPM_chap5.pdf
#https://lukejharmon.github.io/ilhabela/instruction/2015/07/03/PGLS/
#https://dwbapst.github.io/PaleoSoc_phylo_short_course_2019/articles/module_09_worked_PCM_example/module_08.html
################################
#### SET UP THE ENVIRONMENT ####
################################
# load packages
library(MASS)
library(nlme)
library(ape)
library(geiger)
library(phytools)
library(readxl)
library(tidyr)
library(dplyr)
library(tidyverse)
# Load personal function library and data
source("C:/Users/matte/Desktop/cuticle_evolution/R_scripts/resources/cuticle_00_personal_library.R")
source("C:/Users/matte/Desktop/cuticle_evolution/R_scripts/resources/load_data.R")
# Load the trees
temp <- data_w_mean; go_two_hundred_phylo = TRUE
source("C:/Users/matte/Desktop/cuticle_evolution/R_scripts/resources/load_trees.R")
dataset <- temp
# Compute the mean covariance matrix accross all trees
cov_mat <- compute_mean_vcv(two_hundred_phy_pruned)/200
# Load the trait data.R
source("C:/Users/matte/Desktop/cuticle_evolution/R_scripts/resources/load_trait_data.R")
head(dataset)
dataset$diet <- factor(dataset$diet, levels = c('omnivore', 'predator', 'herbivore', "fungivore"))
dataset$foraging <- factor(dataset$foraging, levels = c('epigaeic', 'arboreal', 'hypogaeic'))
# Set the formula
formula <- log10(cuticle_volume) ~
mean_ann_prec + mean_ann_temp +
log10(colony_size) + spine_number +
foraging + diet +
log10(body_volume)
# prepare dataset and covariance matrix
res <- process_data_for_PGLS(formula, dataset, cov_mat)
plot_subfam_balance(dataset_balancing(res$dataset_pruned))
dataset_content(res$dataset_pruned)
# Run model 1
pgls.r1 <- gls(formula, correlation = corSymm(res$cov_mat_pruned[lower.tri(res$cov_mat_pruned)], fixed = TRUE), method = "ML", data = res$dataset_pruned)
formula
res$dataset_pruned
# Run model 1
pgls.r1 <- gls(formula,
correlation = corSymm(res$cov_mat_pruned[lower.tri(res$cov_mat_pruned)], fixed = TRUE),
method = "ML", data = res$dataset_pruned)
