<< 'knitr_options-simu', echo = FALSE>>= knitr::opts_knit$set(fig.pos = "HT", fig.width = 6, fig.height = 6, fig.align = "center", warnings = FALSE, echo = FALSE, format = "latex") @ % Simu: Plusieurs design, tailles etc % On sait la vérité, on peut connaitre les vrais positifs etc % Qu'est ce qu'on prend en entrées qu'est ce qu'on veut en sortie % Bien insister sur l'arbre d'entrée et l'objectif de la simu : quelle approche pour mieux détecter les gènes différentiellement exprimés. % Simulations : % \begin{itemize} % \item soit selon l'arbre des données % \item soit partir sur regarder l'impact de la taille de l'arbre etc. % \end{itemize} \subsection{Erreur de type I et puissance} Dans cette partie nous souhaitons comparer les résultats de l'ANOVA et de l'ANOVA phylogénétique classique, avec approximation de Satterthwaite et avec le \emph{Likelihood ratio test}. Pour cela nous allons simuler des données selon plusieurs modalités et évaluer l'\emph{erreur de première espèce} et la \emph{puissance} obtenue. \begin{itemize} \item Des données réparties en deux groupes au hasard par rapport à la phylogénie. \item Des données réparties en deux groupes cohérents avec la phylogénie. \end{itemize} En sélectionnant des espèces de manière aléatoire, nous cassons la structure induite par la phylogénie. Nous nous attendons donc à ce que l'ANOVA réalise de meilleurs résultats en ne prenant pas en compte l'information phylogénétique. Pour les simulations avec des groupes respectant la structure de l'arbre phylogénétique, nous nous attendons à ce que l'ANOVA phylogénétique parvienne à mieux prendre en compte l'information apportée par la phylogénie et à démêler son effet.\\ <<'modules-simulations', include = FALSE, eval=TRUE>>= necessary_simu <- c("ape", "remotes", "phylolm", "phylolimma", "phytools", "latex2exp", "here") if (any(!(necessary_simu %in% installed.packages()))){ install.packages(necessary_simu) remotes::install_github("lamho86/phylolm", quiet = TRUE) remotes::install_github("pbastide/phylolimma", quiet = TRUE) } library("ape") library("phylolm") library("phytools") library("here") library("tidyverse") library("ggplot2") library("patchwork") library("latex2exp") source(here("R","utils.R")) @ <<'Import-arbre'>>= K <- 2 nb_species <- 43 plot_group_on_tree <- function(tree, groups, ...) { plot(tree, ...) tiplabels(bg = groups, pch = 21, cex = 1.5) } tree <- read.tree(here("R","chen2019.tree")) # Normalising tree edge length taille_tree <- diag(vcv(tree))[1] tree$edge.length <- tree$edge.length / taille_tree @ <<'simus-groupes'>>= seed <- 1234 set.seed(seed) # Mus et Rat vs le reste group_mus_rat_vs_other <- sapply(1:nb_species, function(tip) { if (tip %in% getDescendants(tree = tree, 55)) { return(1) } return(2) }) rng_species <- c("chimp", "bonobo", "human", "orangutan", "marmoset", "musMusculus", "rat", "dog", "ferret", "cow", "opossum") random_groups <- rowSums(sapply(rng_species, function(spec) grepl(spec, tree$tip.label))) random_groups[random_groups == 0] <- 2 #sample(1:K, nb_species, replace = TRUE, prob = c(16, 27)) @ Pour faire nos simulations dans un contexte proche du cas réel nous allons utiliser l'arbre présenté sur la figure~\ref{fig:arbre-chen2019}. Nous choisissons de diviser les espèces en deux groupes. Pour le groupe respectant la phylogénie, on a d'un côté les espèces du genre \emph{Mus} avec les rats et les autres espèces dans un autre groupe (voir la figure~\ref{fig:simu-groupes-mus}). % DONE Choisir les espèces à la main et mettre les réplicats dans le même groupe % on attribue aléatoirement les individus à l'un des deux groupes Et pour le groupe ne respectant pas la phylogénie, nous avons sélectionnés les espèces en respectant les proportions des groupes définis avant afin de rendre les résultats comparables (voir la figure~\ref{fig:simu-groupes-prop}). Enfin pour que notre analyse soit reproductible nous fixons la graine à \Sexpr{seed}. \begin{figure}[H] \begin{subfigure}[H]{0.49\textwidth} \centering <<'plot-groupes-mus'>>= plot_group_on_tree(tree, group = group_mus_rat_vs_other) @ \caption{Groupes \emph{Mus} et rats contre les autres} \label{fig:simu-groupes-mus} \end{subfigure} \begin{subfigure}[H]{0.49\textwidth} \centering <<'plot-groupes-random'>>= plot_group_on_tree(tree, group = random_groups) @ \caption{Groupes sélectionnés sans respect de la phylogénie.} \label{fig:simu-groupes-prop} \end{subfigure} \caption{Arbre et groupes pour les simulations} \label{fig:arbres-groupes} \end{figure} <<'param-simulation'>>= # Generate data for rat&mus vs the rest N <- 500 risk_threshold <- 0.05 ## Standardized parameters total_variance <- 1.0 # sigma2_phylo + sigma2_error, fixed [as tree_height = 1] heri <- c(0.3, 0.5, 0.7, 0.9) # heritability her = sigma2_phylo / total_variance. 0 means only noise. 1 means only phylo. snr <- 1 # signal to noise ratio snr = size_effect / total_variance @ Afin d'avoir un paramètre unique à faire varier, nous re-paramétrisons le modèle, la variance totale $v_{tot}$ suit la relation $v_{tot} = \sigma^2_{phylo} + \sigma^2_{measure} = \Sexpr{total_variance} $. Nous allons faire prendre à $h$, défini comme l'héritabilité, les valeurs $h \in (\Sexpr{heri})$. L'héritabilité est liée à $\sigma^2_{phylo}$ et $\sigma^2_{phylo} = h\times v_{tot}$. Et alors $\sigma^2_{measure} = (1-h) \times v_{tot}$. Ainsi, $h = 0$ signifie qu'il y a seulement du bruit, et $h = 1$ seulement de l'information phylogénétique. Pour les valeurs quantitatives des 2 groupes, nous avons 2 valeurs différentes : \begin{align} \mu_1 = 0, & & \mu_2 = snr \times v_{tot} = \frac{taille\text{ }d'effet}{\cancel{v_{tot}}} \times \cancel{v_{tot}} = \Sexpr{snr} \label{eq:simus-base-values} \end{align} \emph{Note :} \emph{snr} signifie \emph{signal noise ratio} et comme indiqué est donc le rapport entre la taille d'effet et la variance totale. Et dans l'équation~\ref{eq:simus-base-values}, $\mu_1$ et $\mu_2$ correspondent aux $\beta_1$ et $\beta_2$ définis dans la sous-section~\ref{subsec:anova-phylogenetique}. Pour chaque valeur d'héritabilité, nous allons générer $\Sexpr{N}$ jeux de données différents sur lesquels les méthodes sont utilisées avec les valeurs définies dans l'équation~\ref{eq:simus-base-values}. <<'simus-results', warning = FALSE, fig.pos = "H", fig.cap = "Erreur de type I et puissance pour les simulations en faisant varier l'héritabilité", fig.subcap = paste0("$h = ", heri, "$"), fig.ncol = 2, out.width='.49\\linewidth'>>= plot_list <- lapply(seq_len(length(heri)), function(idx) { her <- heri[idx] filename <- here("data", "simus", paste0("real_her_", her, "_seed_", seed, ".Rds")) if (!file.exists(filename)) { sim <- N_simulation_typeI_power(N, groups_list = list(RatMus = group_mus_rat_vs_other, Sélectionnés = random_groups), tree = tree, base_values = c(0, snr * total_variance), sigma2_phylo = her * total_variance, sigma2_measure = (1 - her) * total_variance#, # REML = TRUE ) saveRDS(sim, filename) } sim <- readRDS(filename) df_sim_plot <- compute_power_typeI(df = sim) res_sim_plot <- plot_method_comparison(df_sim_plot) res_sim_plot }) for(plot in plot_list) { print(plot) } @ Sur toutes les sous-figures de la figure~\ref{fig:simus-results}, les étiquettes A présentent les erreurs de type I commises par les méthodes et les étiquettes B présentent les puissances des mêmes méthodes. L'erreur de type I est particulièrement importante à contrôler, en effet elle indique le nombre de faux positifs et l'on veut pouvoir en déterminer le seuil $\alpha$ avec comme seuil classique $0.05$. TODO Insister sur pourquoi trop de faux-positifs pour l'ANOVA classique, du fait de la structure Brownienne, deux clades peuvent être éloignés au niveau temporel beaucoup de génération. En oubliant la structure, on peut vouloir mettre un saut alors que l'écart est simplement dû à de la dérive. L'ANOVA suppose des données iid ce qui n'est pas le cas ici. TODO Important de préciser qu'il faut contrôler l'erreur de type I car les manips coûtent très cher. TODO Ajouter les commentaires sur les simulations \paragraph*{REML vs Maximum Likelihood (ML)} D'après nos simulations, les méthodes utilisant le REML contrôle toujours mieux l'erreur de première espèce que les méthodes utilisant le maximum de vraisemblance.