Skip to content

Commit

Permalink
Merge pull request #21 from MindTheGap-ERC/revision
Browse files Browse the repository at this point in the history 
Revision
  • Loading branch information
@NiklasHohmann
NiklasHohmann authored May 27, 2024
2 parents 6806fab + 417a784 commit 93bd104
Show file tree
Hide file tree
Showing 3 changed files with 29 additions and 50 deletions.
62 changes: 26 additions & 36 deletions code/make_plots.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,8 +160,8 @@ plot_AIC_strat = function(basin, simulated_mode, label){
#'
#' @title plots AIC values in boxplots for the stratigraphic domain
#'
#' @param simulated_mode: "stasis", "Brownian motion", "weak Brownian drift",
#' or "strong Brownian drift". True (= simulated) mode of evolution for which AIC is supposed to be extracted
#' @param simulated_mode: "stasis", "Brownian motion", or "Brownian drift",
#' True (= simulated) mode of evolution for which AIC is supposed to be extracted
#' @param basin: String, "A" or "B"
#' @param label: the label you want attached to the plot, String,"A" to "Z"
#'
Expand All @@ -172,18 +172,15 @@ plot_AIC_strat = function(basin, simulated_mode, label){

#Loads in the correct line thicknesses for the graph:
if(simulated_mode=="stasis"){
highlight=rep(c(0.1,0.1,0.5,0.1),5)
highlight=rep(c(0.1,0.5,0.1),5)
}else if(simulated_mode=="Brownian motion"){
highlight=rep(c(0.1,0.1,0.1,0.5),5)
}else if(simulated_mode=="strong Brownian drift"){
highlight=rep(c(0.5,0.1,0.1,0.1),5)
}else if(simulated_mode=="weak Brownian drift"){
highlight=rep(c(0.5,0.1,0.1,0.1),5)
highlight=rep(c(0.1,0.1,0.5),5)
}else if(simulated_mode=="Brownian drift"){
highlight=rep(c(0.5,0.1,0.1),5)
}

plotS1=ggplot2::ggplot(get_AIC_scenario(basin,simulated_mode), aes(x=position, y=AIC, fill=testedEvoMode)) +
geom_boxplot(outlier.shape = NA,lwd=highlight) +
geom_hline(yintercept = 0.9, linetype = "dashed", linewidth=0.1)+
labs(tag = label)+
scale_fill_brewer(palette=AIC_palette)+
theme(legend.position="none", plot.title =element_text(size = 8), plot.tag = element_text(face = "bold"), plot.tag.position = c(0.01, 0.98),
Expand All @@ -200,8 +197,8 @@ plot_AIC_time = function(no_of_sampl_loc, basin, simulated_mode, label){
#' @title plots AIC values in boxplots for the stratigraphic domain
#'
#'@param no_of_sampl_loc : string, element of noOfSamplingLoc
#' @param simulated_mode: "stasis", "Brownian motion", "weak Brownian drift",
#' or "strong Brownian drift". True (= simulated) mode of evolution for which AIC is supposed to be extracted
#' @param simulated_mode: "stasis", "Brownian motion", "Brownian drift",
#' True (= simulated) mode of evolution for which AIC is supposed to be extracted
#' @param basin: String, "A" or "B"
#' @param label: the label you want attached to the plot, String,"A" to "Z"
#'
Expand All @@ -216,19 +213,16 @@ plot_AIC_time = function(no_of_sampl_loc, basin, simulated_mode, label){

#Loads in the correct line thicknesses for the graph:
if(simulated_mode=="stasis"){
highlight=rep(c(0.1,0.1,0.5,0.1),Mult_high)
highlight=rep(c(0.1,0.5,0.1),Mult_high)
}else if(simulated_mode=="Brownian motion"){
highlight=rep(c(0.1,0.1,0.1,0.5),Mult_high)
}else if(simulated_mode=="strong Brownian drift"){
highlight=rep(c(0.5,0.1,0.1,0.1),Mult_high)
}else if(simulated_mode=="weak Brownian drift"){
highlight=rep(c(0.5,0.1,0.1,0.1),Mult_high)
highlight=rep(c(0.1,0.1,0.5),Mult_high)
}else if(simulated_mode=="Brownian drift"){
highlight=rep(c(0.5,0.1,0.1),Mult_high)
}

#Plots the plot:
plotT1= ggplot2::ggplot(plot_data, aes(x=nsp, y=AIC, fill=testedEvoMode)) +
geom_boxplot(outlier.shape = NA,linewidth=highlight) +
geom_hline(yintercept = 0.9, linetype = "dashed", linewidth=0.1)+
labs(tag = label)+
scale_fill_brewer(palette=AIC_palette)+
theme(legend.position="none", plot.title =element_text(size = 8), plot.tag = element_text(face = "bold"), plot.tag.position = c(0.01, 0.98),
Expand Down Expand Up @@ -334,7 +328,7 @@ plot_comparison_evo_modes_time_domain = function(scenario, no_of_lineages = 3, p
file_name = paste("figs/R/comparison_evo_modes_time_domain_scenario_", scenario, "_raw.pdf", sep = "")
pdf(file = file_name , width=default_width_fp_in)

combined_plot = gridExtra::grid.arrange(plot_list[[1]], plot_list[[2]], plot_list[[3]], plot_list[[4]])
combined_plot = gridExtra::grid.arrange(plot_list[[1]], plot_list[[2]], plot_list[[3]])
dev.off()
}

Expand All @@ -351,17 +345,16 @@ make_test_res_strat_plot = function(scenario){
#The stratigraphic plots:
Plot4_1=plot_AIC_strat(scenario,"stasis","A")
Plot4_2=plot_AIC_strat(scenario,"Brownian motion","B")
Plot4_3=plot_AIC_strat(scenario,"weak Brownian drift","C")
Plot4_4=plot_AIC_strat(scenario,"strong Brownian drift","D")
Plot4_3=plot_AIC_strat(scenario,"Brownian drift","C")

#The time plots:
Plot4_5=plot_AIC_time(ts_lengths,scenario,"stasis","E")
Plot4_6=plot_AIC_time(ts_lengths,scenario,"Brownian motion","F")
Plot4_7=plot_AIC_time(ts_lengths,scenario,"weak Brownian drift","G")
Plot4_8=plot_AIC_time(ts_lengths,scenario,"strong Brownian drift","H")
Plot4_5=plot_AIC_time(ts_lengths,scenario,"stasis","D")
Plot4_6=plot_AIC_time(ts_lengths,scenario,"Brownian motion","E")
Plot4_7=plot_AIC_time(ts_lengths,scenario,"Brownian drift","F")


#Combining all the plots:
combined_plot=gridExtra::grid.arrange(Plot4_1,Plot4_2,Plot4_3,Plot4_4,Plot4_5,Plot4_6,Plot4_7,Plot4_8, ncol=4)
combined_plot=gridExtra::grid.arrange(Plot4_1,Plot4_2,Plot4_3,Plot4_5,Plot4_6,Plot4_7, ncol=3)

#Printing the plots plus the legend to .pdf

Expand All @@ -385,16 +378,15 @@ make_test_res_time_plot = function(scenario){
#The plots:
Plot6_1=plot_AIC_time(ts_lengths,scenario,"stasis","A")
Plot6_2=plot_AIC_time(ts_lengths,scenario,"Brownian motion","B")
Plot6_3=plot_AIC_time(ts_lengths,scenario,"weak Brownian drift","C")
Plot6_4=plot_AIC_time(ts_lengths,scenario,"strong Brownian drift","D")
Plot6_3=plot_AIC_time(ts_lengths,scenario,"Brownian drift","C")

#Combining all the plots:
combined_plot=gridExtra::grid.arrange(Plot6_1,Plot6_2,Plot6_3,Plot6_4, ncol=2)
combined_plot=gridExtra::grid.arrange(Plot6_1,Plot6_2,Plot6_3, ncol=1)

#Printing the plots plus the legend to .pdf

file_name = paste("figs/R/test_results_time_domain_", as.character(scenario), "_raw.pdf", sep = "")
pdf(file = file_name, width = default_width_fp_in, height = 3.25)
pdf(file = file_name, width = default_width_fp_in, height = 6)
gridExtra::grid.arrange(combined_plot, Legend, nrow = 2, heights = c(10, 1)) #The multiplot
dev.off()
}
Expand Down Expand Up @@ -609,7 +601,7 @@ plot_pairwise_comparison(scenario_1 = "A",
scenario_2 = "B",
dist_1 = "2 km",
dist_2 = "6 km",
mode = "strong Brownian drift",
mode = "Brownian drift",
no_of_lineages = 3,
plot_seed = 1)

Expand Down Expand Up @@ -859,7 +851,7 @@ make_var_pres_of_modes_plot = function(pos, scenario, no_of_lineages = 3, plot_s
comb_list = list()
for (i in seq_along(simulatedEvoModes)){
comb_list[[simulatedEvoModes[i]]] = compare_evo_modes_time_strat_plot(name = simulatedEvoModes[i],
labels = LETTERS[c(i,i+4)],
labels = LETTERS[c(i,i+3)],
lwd = evo_mode_lwds[i],
no_of_lineages = 3)
}
Expand All @@ -868,12 +860,10 @@ make_var_pres_of_modes_plot = function(pos, scenario, no_of_lineages = 3, plot_s
pdf(file = file_name , width=default_width_fp_in, height = 3.25)
combined_plot = gridExtra::grid.arrange(comb_list$stasis$strat,
comb_list$`Brownian motion`$strat,
comb_list$`weak Brownian drift`$strat,
comb_list$`strong Brownian drift`$strat,
comb_list$`Brownian drift`$strat,
comb_list$stasis$time,
comb_list$`Brownian motion`$time,
comb_list$`weak Brownian drift`$time,
comb_list$`strong Brownian drift`$time,
comb_list$`Brownian drift`$time,
nrow = 2)

dev.off()
Expand Down
7 changes: 1 addition & 6 deletions code/parameters_for_tests.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,15 +25,10 @@ EvoModes[[2]] <- list(
params = c(0, 1)
) # BM with mean zero and variance 1
EvoModes[[3]] <- list(
name = "weak Brownian drift",
name = "Brownian drift",
mode = "BD",
params = c(5, 1)
) # BD with drift 5 and variance 1
EvoModes[[4]] <- list(
name = "strong Brownian drift",
mode = "BD",
params = c(10, 1)
) # BD with drift 10 and variance 1

# extract names of the simulated evolutionary modes
simulatedEvoModes <- sapply(EvoModes, function(x) x$name)
Expand Down
10 changes: 2 additions & 8 deletions code/utils.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,10 +212,9 @@ testModesOfEvolStrat <- function(distBetweenSamples,
w.grw <- paleoTS::fitSimple(fossilTimeSeries, model = "GRW") # General Random walk (Brownian drift)
w.urw <- paleoTS::fitSimple(fossilTimeSeries, model = "URW") # Unbiased Random Walk (Brownian motion)
w.stat <- paleoTS::fitSimple(fossilTimeSeries, model = "Stasis") # Stasis
w.ou <- paleoTS::fitSimple(fossilTimeSeries, model = "OU") # Ornstein-Uhlenbeck model

# This compares the models, the most likely model has the highest akaike.wt and the smallest AICc values.
compared <- paleoTS::compareModels(w.grw, w.urw, w.stat, w.ou, silent = TRUE, sort = FALSE)
compared <- paleoTS::compareModels(w.grw, w.urw, w.stat, silent = TRUE, sort = FALSE)

# store test results and input data
testRes <- list(
Expand Down Expand Up @@ -286,16 +285,11 @@ testModesOfEvolTime <- function(maxTime,
y = fossilTimeSeries,
model = "Stasis"
) # Stasis
w.ou <- paleoTS::fitSimple(
y = fossilTimeSeries,
model = "OU"
) # Ornstein-Uhlenbeck model

# This compares the models, the most likely model has the highest akaike.wt and the smallest AICc values.
compared <- paleoTS::compareModels(w.grw,
w.urw,
w.stat,
w.ou,
silent = TRUE,
sort = FALSE
)
Expand Down Expand Up @@ -409,6 +403,6 @@ get_data_from_osf = function(link){
#' @param link url to the (public) url
my_project <- osfr::osf_retrieve_node(link)
my_files <- osfr::osf_ls_files(my_project)
osfr::osf_download(my_files, recurse = TRUE, conflicts = "overwrite", progress = TRUE)
osfr::osf_download(my_files, path = "data/" , recurse = TRUE, conflicts = "overwrite", progress = TRUE)
return(invisible())
}

0 comments on commit 93bd104

Please sign in to comment.