diff --git a/simulations/functions-anova.R b/simulations/functions-anova.R index 0332598..4942b56 100644 --- a/simulations/functions-anova.R +++ b/simulations/functions-anova.R @@ -60,7 +60,12 @@ simulate_ANOVAs <- function( ## ANOVAs fit_ANOVA <- lm(y ~ groups) - fitphy_ANOVA <- phylolm(y ~ groups, phy = tree, model = stoch_process) + if (stoch_process == "OU"){ + model = "OUfixedRoot" + } else { + model = "BM" + } + fitphy_ANOVA <- phylolm(y ~ groups, phy = tree, model = model) ## DONE refaire avec ces modalités et évaluer les erreurs de type 1 et erreurs de type 2 ## faire scénario H_0: mu egaux -> ANOVA se plante car dep entre les indivs diff --git a/simulations/test-anova.phylo.R b/simulations/test-anova.phylo_MB.R similarity index 97% rename from simulations/test-anova.phylo.R rename to simulations/test-anova.phylo_MB.R index 5ee59bb..53f81f5 100644 --- a/simulations/test-anova.phylo.R +++ b/simulations/test-anova.phylo_MB.R @@ -79,9 +79,9 @@ plot_data <- data.frame( 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", + title = "Proportions correctes vs Tested Method (BM)", x = "Tested Method", - y = "Puissance" + y = "Proportions correctes" ) + theme_minimal() # TODO : Regarder la notice de lmertest pour l'implémentation de Satterthwaite diff --git a/simulations/test-anova.phylo_OU.R b/simulations/test-anova.phylo_OU.R new file mode 100644 index 0000000..0332aac --- /dev/null +++ b/simulations/test-anova.phylo_OU.R @@ -0,0 +1,89 @@ +library(phylolm) +library(phylotools) +library(phytools) +library(ape) +library(ggplot2) +set.seed(1234) + +N <- 100 # Number of different simulations +n <- 100 + +# Arbre +tree <- rphylo(n, 0.1, 0) + +## Groupes +K <- 3 +get_group <- function(tip) { + if (tip %in% getDescendants(tree, 105)) { + return(2) + } + if (tip %in% getDescendants(tree, 110)) { + return(3) + } + return(1) +} + +source("./simulations/functions-anova.R") + +# Computing groups +phylo_groups <- as.factor(sapply(1:n, get_group)) +non_phylo_groups <- as.factor(sample(c(1, 2, 3), n, replace = TRUE)) + +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, stoch_process = "OU"))) +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, stoch_process = "OU"))) + +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 = "Proportions correctes vs Tested Method (OU)", + x = "Tested Method", + y = "Proportions correctes" + ) + + 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 + diff --git a/sources/test_satterthwaite_star_tree.R b/sources/test_satterthwaite_star_tree.R index 49af5df..f6c43f2 100644 --- a/sources/test_satterthwaite_star_tree.R +++ b/sources/test_satterthwaite_star_tree.R @@ -1,7 +1,7 @@ library(lmerTest) library(phylolm) -source("test_satterthwaite_utils.R") +source("./sources/test_satterthwaite_utils.R") set.seed(12891289)