polarolouis/anova-phylogenetique-projet-msv
1
0
Fork 0
mirror of https://github.com/Polarolouis/anova-phylogenetique-projet-msv.git synced 2026-06-17 18:25:25 +02:00
Code Issues Projects Releases Packages Wiki Activity

Added a plot to compare the different cases

Browse source
This commit is contained in:
Louis Lacoste 2024-01-04 12:12:34 +01:00
parent 28377d6339
commit 074d3d9dc1
1 changed files with 57 additions and 4 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

59
simulations/test-anova.phylo.R
View file

@ -2,6 +2,7 @@ library(phylolm)
library(phylotools)
library(phytools)
library(ape)
library(ggplot2)
set.seed(1234)
N <- 100 # Number of different simulations
@ -28,9 +29,61 @@ source("./simulations/functions-anova.R")
phylo_groups <- as.factor(sapply(1:n, get_group))
non_phylo_groups <- as.factor(sample(c(1, 2, 3), n, replace = TRUE))
# N répétitions pour les 2 groupes générés
phylo_groups_results <- do.call("rbind", lapply(1:N, function(id) simulate_ANOVAs(sim_id = id, groups = phylo_groups, tree = tree)))
non_phylo_groups_results <- do.call("rbind", lapply(1:N, function(id) simulate_ANOVAs(sim_id = id, groups = non_phylo_groups, tree = tree)))
calcul_puissance <- function(data, test_method) {
mean(data[which(data$tested_method == test_method), ]$has_selected_correctly)
}
# Mu tous différents
mu_vect <- c(1, 5, 10)
# N répétitions pour les 2 groupes générés
mu_diff_phylo_groups_results <- do.call("rbind", lapply(1:N, function(id) simulate_ANOVAs(sim_id = id, groups = phylo_groups, tree = tree, mu_vect = mu_vect)))
mu_diff_non_phylo_groups_results <- do.call("rbind", lapply(1:N, function(id) simulate_ANOVAs(sim_id = id, groups = non_phylo_groups, tree = tree, mu_vect = mu_vect)))
puissance_mu_diff_phylo_for_phylo_groups <- calcul_puissance(mu_diff_phylo_groups_results, "ANOVA-Phylo")
puissance_mu_diff_classic_for_phylo_groups <- calcul_puissance(mu_diff_phylo_groups_results, "ANOVA")
puissance_mu_diff_phylo_for_non_phylo_groups <- calcul_puissance(mu_diff_non_phylo_groups_results, "ANOVA-Phylo")
puissance_mu_diff_classic_for_non_phylo_groups <- calcul_puissance(mu_diff_non_phylo_groups_results, "ANOVA")
# Mu égaux
mu_vect <- rep(1, K)
# N répétitions pour les 2 groupes générés
mu_equals_phylo_groups_results <- do.call("rbind", lapply(1:N, function(id) simulate_ANOVAs(sim_id = id, groups = phylo_groups, tree = tree, mu_vect = mu_vect)))
mu_equals_non_phylo_groups_results <- do.call("rbind", lapply(1:N, function(id) simulate_ANOVAs(sim_id = id, groups = non_phylo_groups, tree = tree, mu_vect = mu_vect)))
# Calcul de puissance
puissance_mu_equals_phylo_for_phylo_groups <- calcul_puissance(mu_equals_phylo_groups_results, "ANOVA-Phylo")
puissance_mu_equals_classic_for_phylo_groups <- calcul_puissance(mu_equals_phylo_groups_results, "ANOVA")
puissance_mu_equals_phylo_for_non_phylo_groups <- calcul_puissance(mu_equals_non_phylo_groups_results, "ANOVA-Phylo")
puissance_mu_equals_classic_for_non_phylo_groups <- calcul_puissance(mu_equals_non_phylo_groups_results, "ANOVA")
# Graphiques
puissances <- c(
puissance_mu_diff_phylo_for_phylo_groups,
puissance_mu_diff_classic_for_phylo_groups,
puissance_mu_equals_phylo_for_phylo_groups,
puissance_mu_equals_classic_for_phylo_groups,
puissance_mu_diff_phylo_for_non_phylo_groups,
puissance_mu_diff_classic_for_non_phylo_groups,
puissance_mu_equals_phylo_for_non_phylo_groups,
puissance_mu_equals_classic_for_non_phylo_groups
)
plot_data <- data.frame(
puissance = puissances,
tested_method = rep(c("ANOVA-Phylo", "ANOVA"), 4),
group_type = rep(c("phylo", "non_phylo"), each = 4),
mu_type = rep(rep(c("different", "equals"), each = 2), 2)
)
ggplot(plot_data, aes(x = tested_method, y = puissance, fill = interaction(group_type, mu_type))) +
geom_bar(stat = "identity", position = "dodge") +
labs(
title = "Puissance vs Tested Method",
x = "Tested Method",
y = "Puissance"
) +
theme_minimal()
# TODO : Regarder la notice de lmertest pour l'implémentation de Satterthwaite
# TODO : En utilisant l'arbre étoile, on obtient un modele mixte classique donc on peut appliquer lmerTest