From 074d3d9dc1b82b11be4189bfbc743d1e24256862 Mon Sep 17 00:00:00 2001
From: Polarolouis <louis.lacoste@hotmail.fr>
Date: Thu, 4 Jan 2024 12:12:34 +0100
Subject: [PATCH] Added a plot to compare the different cases

---
 simulations/test-anova.phylo.R | 61 +++++++++++++++++++++++++++++++---
 1 file changed, 57 insertions(+), 4 deletions(-)

diff --git a/simulations/test-anova.phylo.R b/simulations/test-anova.phylo.R
index fb5e80d..5ee59bb 100644
--- a/simulations/test-anova.phylo.R
+++ b/simulations/test-anova.phylo.R
@@ -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
\ No newline at end of file
+# TODO : En utilisant l'arbre étoile, on obtient un modele mixte classique donc on peut appliquer lmerTest
+