anova-phylogenetique-projet-msv/rapport.Rnw

\documentclass[a4paper, 12pt]{article}

% Packages
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[french]{babel}
\usepackage{geometry}

%Images
\usepackage{graphicx}
\graphicspath{{img/}}
\usepackage{float}
\usepackage{subcaption} %  for subfigures environments
\usepackage{wrapfig} 

% Booktabs
\usepackage{booktabs}

% Citation
\usepackage{csquotes}

\usepackage{caption}
\usepackage{subcaption}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{hyperref}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{amsthm}
\usepackage{cancel}

\usepackage[style=authoryear-comp,backend=biber]{biblatex}
%== use and define color ==%
\AtEveryCite{\color{blue}}
\addbibresource{references.bib}

% Configurations
\geometry{a4paper, margin=2.5cm}
\graphicspath{ {img/} }

% Macros utiles
\newcommand{\Normal}{\mathcal{N}}


% Titre du document
\title{Rapport de Projet : ANOVA Phylogénétique \newline Co-encadré par Mélina Gallopin et Paul Bastide}
\author{Alizée Geffroy \and Louis Lacoste}
\date{20 mars 2024}

\newtheorem*{approximation}{Approximation}

\begin{document}

\maketitle

<<'init', include=FALSE>>=
    knitr::opts_chunk$set(echo = FALSE)
    require("knitr", quietly = TRUE)
    options(knitr.table.format = "latex")
@

<<'libraries', include=FALSE>>=
# "phytools", "phylotools"
necessary_packages <- c("ape", "here")
if (any(!(necessary_packages %in% installed.packages()))) {
    install.packages(necessary_packages)
}

# require(phytools)
# require(phylotools)
require(ape)
require(here)


source(here("R","utils.R"))
@

\newpage
\tableofcontents
\newpage

\section{Introduction}
\label{chap:intro}
% Introduction au projet, contexte, objectifs.
Avec l'avènement des données massives de génomique, transcriptomique et 
protéomique, il est impératif de disposer de techniques statistiques robustes 
et adaptées à l'échelle pour mener à bien les analyses. Ces données génétiques 
fournissent généralement deux types d'informations : les mesures elles-mêmes et 
les arbres phylogénétiques. Dans certains cas, ces arbres présentent des 
ramifications intra-espèces.

Nous mesurons l'expression des gènes orthologues chez plusieurs espèces, 
souvent considérée comme une expression constitutive. Par exemple, la base de 
données BGee (\url{https://www.bgee.org}) compile les niveaux d'expression de 
gènes chez diverses espèces. En utilisant ces données inter-espèces, notre 
objectif est de détecter les gènes orthologues présentant des variations 
d'expression entre différents groupes d'espèces. En notant $Y$ le niveau 
d'expression d'un gène chez plusieurs espèces, nous modélisons ce niveau 
d'expression par une variable aléatoire ayant une moyenne $\mu_1$ dans un groupe
d'espèces (par exemple, les primates) et $\mu_2$ dans un autre groupe d'espèces 
(par exemple, toutes les espèces non primates). Nous testons alors l'hypothèse 
$H_0$ : $\mu_1 = \mu_2$ contre $H_1$ : $\mu_1 \neq \mu_2$. 

Par exemple, dans l'article 
\cite{chenQuantitativeFrameworkCharacterizing2019}, les auteurs identifient
un ensemble de gènes exprimés dans le foie qui sont sous-exprimés chez les 
primates par rapport aux autres espèces, sous-exprimés chez les rongeurs par 
rapport aux autres espèces ou sous-exprimés dans les tissus des testicules chez
les primates par rapport aux autres espèces (voir la figure~\ref{fig:chen-fig4}
extraite de l'article de
\cite{chenQuantitativeFrameworkCharacterizing2019} où les individus sont 
représentés en colonnes et les gènes en lignes, la couleur reflétant le niveau 
d'expression du gène).

\begin{figure}[H]
    \centering
    \includegraphics[width=1\textwidth]{chenFig4.png}
    \caption{Figure extraite de l'article de \cite{chenQuantitativeFrameworkCharacterizing2019}}
    \label{fig:chen-fig4}
\end{figure}

Ces données illustrent parfaitement la nécessité de techniques analytiques
robustes face à des ensembles de données complexes, combinant à la fois des
aspects écologiques et transcriptomiques.

% Format des données : arbres phylogénétiques, données génétiques 
% Arbres avec des petites branche: plusieurs individus par espèces avec chacun leurs données 
% --> problème biologique 

% Deux sujets différents écologie et transcriptomique mais une même méthode.
La figure~\ref{fig:arbre-chen2019} présente l'arbre phylogénétique :
\begin{figure}[H]
    \centering
<<'plot-arbre-chen'>>=
    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

    phytools::plotTree(tree, ftype="i")
@
    \caption{Arbre phylogénétique de \cite{chenQuantitativeFrameworkCharacterizing2019}}
    \label{fig:arbre-chen2019}
\end{figure}
La problématique qui se pose souvent est celle de l'analyse de différents gènes qu'on mesure chez plusieurs espèces. 
On cherche alors d'abord à trouver quels gènes pourraient être différentiels chez certaines espèces, en détectant un changement d'expression dans certains groupes. 
En considérant l'arbre precédent \ref{fig:arbre-chen2019}, on pourra chercher les gènes qui sont différents entre les groupes des \textit{mus} et \textit{rat} par rapport aux autres espèces. 
\newline
Le modèle le plus couramment utilisé est actuellement l'Expression Variance and Evolution modèle (EVE) présenté dans \cite{rohlfsPhylogeneticANOVAExpression2015}. 
L'EVE modèle est basé sur un Likelihood Ratio Test (LRT), une méthode statistique classique. 
Ce projet s'inscrit alors dans un questionnement plus large qui cherche à se demander si d'autres modèles classiques comme l'ANOVA, en les adaptant, pourrait produire des résultats similaires voire meilleurs que l'EVE modèle. 
En effet, avoir un bon modèle qui, en particulier, donne peu de faux positifs est important. 
On peut ensuite étudier les gènes potentiellement intéressants selon une problématique et des groupes d'espèces données. 
\newline
Au vu de la forme des données étudiées, le projet s'est tourné vers une méthode 
d'ANOVA phylogénétique. 
Celle-ci sera d'abord décrite ainsi que d'autres outils mathématiques utilisés pour affiner la fiabilité du test dans une première partie. Certains auront fait l'objet de calculs explicites en vue de leur implémentation. 
A partir de ces résultats, nous avons implémenter ces méthodes en R
d'abord appliquées à des simulations destinées à comparer et étudier la 
méthode d'ANOVA phylogénétique sur des données d'arbre simulés. 
Enfin, on a testé sur des données réelles. 

Au cours de ce projet nous avons donc eu une partie d'étude théoriques et mathématiques des modèles de l'ANOVA et de l'ANOVA phylogénétique afin de bien l'adapter à nos données. 
A partir de la formulation mathématiques des modèles
\newline
\newline
Tout le code produit est disponible sur le dépôt GitHub suivant 
\url{https://github.com/Polarolouis/anova-phylogenetique-projet-msv/}.
Ce dépôt contient le code pour implémenter la méthode, faire les 
simulations et compiler le rapport.

Nous avons au maximum indiqué le code qui n'a pas été écrit par nous, la plupart
du temps dans les commentaires du code.

\section{Méthodes}
\label{sec:methode}
Dans cette partie nous présentons les modèles statistiques d'ANOVA et sa dérivée phylogénétique. 
Après avoir posé le cadre mathématique à partir des recherches bibliographiques, nous développerons les outils mathématiques. 
En particulier pour l'approximation nous avons calculé une forme explicite afin de l'implémenter. 
Finalement, nous faisons une présentation succinte des méthodes REML et du modèle LRT. 

\subsection{L'ANOVA}

L'ANOVA est un cas classique du modèle linéaire, nous utilison ici une forme 
matricielle.

Le principe de l'ANOVA est d'expliciter le lien entre une variable quantitative
et une ou plusieurs variables qualitatives.

La forme matricielle usuelle de l'ANOVA à 1 facteur et 2 groupes de taille respectivement $n_1$ et $n_2$ est la suivante :

\begin{equation}
    Y = X\beta + u \text{,} \quad u \sim \mathcal{N}_n(0, \sigma^2I_n)
    \label{eq:ANOVA}
\end{equation} 
\[
    \text{où} \quad \mathbf{Y} = \begin{bmatrix} Y_{11} \\\vdots\\ Y_{1n_1} \\ 
        Y_{21}\\ \vdots \\ Y_{2n_2} \end{bmatrix}\text{, } 
        \mathbf{X} = \begin{bmatrix} \mathbf{1} & \mathbf{1_{n_1}} \end{bmatrix}=\begin{bmatrix} 1 & 1 \\ \vdots & \vdots\\1 & 1 \\ 1 & 0\\ \vdots & \vdots\\1 & 0 \end{bmatrix}\text{, }  \mathbf{\beta} = \begin{bmatrix} \beta_1 \\ \beta_2  \end{bmatrix} \text{, } n=n_1+n_2
\]
On noter qu'ici  $\beta_1 = \mu_1$ la moyenne du groupe 1, et $\beta_2 = \mu_2 - \mu_1$ la différence des moyennes entre les groupes dans cette paramétrisation.
\newline
\newline
Les paramètres $(\beta_1, \beta_2, \sigma^2)$ de l'ANOVA sont estimables, grâce par exemple à la méthode du maximum de vraisemblance et ont des formules bien connues.

\subsubsection*{Test statistique}
\label{subsubsec: test-ANOVA}
Dans le  cadre d'ANOVA classique nous allons rappeler les hypothèses du test et la statistique de test. 
On fait un test sur les moyennes de chaque groupe. Ce peut être la moyenne de la valeur d'un trait génétique ou bien de la valeur de la fréquance d'une séquence ou allèle. 
On testera alors les hypothèses suivantes avec $l=\begin{bmatrix}0 \\1 \end{bmatrix}$:
\[ H_0 : \beta_2 =0 \Leftrightarrow l^T\beta = \begin{bmatrix}0 \\0\end{bmatrix} \text{, les 2 groupes ont la même moyenne}\]
\[ H_1 : \beta_2\neq 0 \text{, les 2 groupes ont des moyennes différentes}\]
On a alors la statistique de test suivante : 
\begin{equation}
    \label{eq:F_ANOVA}
    F_{ANOVA}=\frac{||\hat{Y} - \bar{Y}||^2(n-2)}{||Y - \hat{Y}||^2} \underset{\mathcal{H}_0}{\sim}\mathcal{F}\text{isher} (1, n-2)
\end{equation}
\[\text{Où } \bar{Y}= \frac{1}{n} \sum_{i,j=1}^{n_1,n_2} Y_{i,j} \text{ et } \hat{Y}=X\hat{\beta}\]
\subsection{L'ANOVA phylogénétique}
\label{subsec:anova-phylogenetique}
Dans la méthode d'ANOVA classique l'information portée par l'arbre phylogénétique n'est pas prise en compte. 
Le but de cette nouvelle méthode est de ne plus mettre cette information de côté et peut être obtenir de meilleurs résultats. 
En effet on peut imaginer, en considérant des traits évolutifs ou des séquences d'ADN, que des individus d'une même espèce ou bien d'espèces proche phylogénétiquement, pourraient avoir des valeurs proches.  
Il s'agira alors de modéliser l'arbre et les informations évolutives qu'ils contient de manière à l'incorporer. 
\newline
\newline
Comme décrit dans \cite{bastideModelesEvolutionCaracteres2022} l'évolution d'un trait nécessite de décrire ses fluctuations le long de l'arbre et ses branches. 
C'est pour cela que souvent cela est le résultat d'un processus stochastique à temps continu branchant sur un arbre phylogénétique, supposé connu et fixé. 
Le processus classique est le mouvement brownien et c'est celui que nous avons utilisé. Il a cependant quelques limites qui ne font pas l'objet de ce rapport mais qui peuvent alors justifier le choix d'autres types de processus comme celui d'Ornstein-Uhlenbecks.
Le modèle de mouvement brownien va alors induire que les feuilles des arbres (nos observations) auront une distribution gausienne que l'on écrira sous la forme suivante:

\begin{equation}
    Y = X\beta + u \text{, } u \sim \mathcal{N}_n(0, \sigma^2_{phy}K)
    \label{eq:ANOVAphylo}
\end{equation} 


Les notations correspondent toujours à celles utilisées pour \eqref{eq:ANOVA}. La seule différence se trouvant dans la distribution de $u$ et la présence d'une matrice $K$.
Dans le cadre du mouvement brownien $K=(K_{i,j})_{1\leq i,j \leq n}=(t_{i,j})_{1\leq i,j \leq n}$ où $t_{i,j}$ représente le temps d’évolution commun aux espèces i et j. 
Comme on peut le voir dans l'exemple suivant, cette matrice a bien la forme attendue : deux espèces proches dans l'arbre ont un coefficient de covariance élevé (leur temps d'évolution commun est grand), alors que deux espèces éloignées sont plus faiblement corrélées.
On peut voir un exemple utilisé dans les slides de cours \cite{bastideContinuousTraitEvolution2022}:
\begin{center}
    \includegraphics[width=0.7\textwidth]{matrix_K.png}
\end{center}
<<'plot-MB', warnings = FALSE, message = FALSE, fig.cap = "Exemple d'un arbre phylogénétique dont le trait est généré selon un Mouvement Brownien", out.width = "75%", fig.height = 3.5, fig.align = "center", fig.pos = "H">>=
source(here("simulations","mouvement_brownien.R"))
set.seed(12)
df <- generate_phylo_tree(5,100, 0.1)
df$id_branche <- as.factor(df$id_branche)

ggplot(df) +
    aes(x = time, y = traj, color = as.factor(id_branche))+
    geom_line() +
    # geom_vline(data = df[df$is_spec_time,], aes(xintercept = time, color = .data$id_branche), linetype = "dashed")+
    geom_point(data = df[df$is_spec_time,], aes(x = time, y = .data$traj), color = df[df$is_spec_time,]$id_branche, size = 4)+
    labs(color = "Espèce", x = "Temps", y = "Trait") +
    coord_fixed(ratio=0.75) +
    theme_minimal()
@
A note que seules les réalisation du processus aux feuilles de l'abre (ici, à t = 10) sont observées. Le reste du processus est latent.
\subsubsection*{Test statistique}
\label{subsubsec: test-ANOVAphylo}
En considérant le même test et les mêmes hypothèses que \ref{subsubsec: test-ANOVA} mais en prenant en compte la nouvelle formule \ref{eq:ANOVAphylo}, on obtient une nouvelle statistique de test. 

\cite{bastideContinuousTraitEvolution2022} nous donne la forme de la statistique pour la méthode d'ANOVA de cette forme. 
\begin{equation}
    F_{ANOVAphylo}=\frac{||\hat{Y} - \bar{Y}||^2_{K^{-1}}(n-2)}{||Y - \hat{Y}||^2_{K^{-1}}} \underset{\mathcal{H}_0}{\sim}\mathcal{F}\text{isher} (1, n-2)
\end{equation}
\begin{align*}
&\text{Où }||Y - \hat{Y}||^2_{K^{-1}} = ||Y - X\hat{\beta}||^2_{K^{-1}}= Proj_X^{K}Y= (Y-\hat{Y})^TK^{-1}(Y-\hat{Y})\\
&\text{et }||\hat{Y} - \bar{Y}||^2_{K^{-1}}=(\hat{Y} - \bar{Y})^TK^{-1}(\hat{Y}- \bar{Y})
\end{align*}
Concernant cette statitstique, on peut dire qu'elle est toujours exacte car on connait la matrice $K$.

\subsection{ANOVA phylogénétique avec erreur de mesure}
Dans la section précedente, on a supposé que la seule source de variabilité provenait du mouvement brownien sur l'arbre.
On rajoute dans cette section une autre variabilité specifiée par $\sigma^2_{err}$ qui à partir de la formule précédente \eqref{eq:ANOVAphylo}, nous donne:
\begin{equation}
    Y = X\beta + u + \epsilon \text{, } \quad u \sim \mathcal{N}_n(0, \sigma^2_{phy}K) \text{,} \quad \epsilon \sim \mathcal{N}_n(0, \sigma^2_{err}I_n)
\label{eq:eq2err}
\end{equation} 
\begin{align}
    &\text{Alors } Y \sim \mathcal{N}_n(X\beta, \sigma^2_{phy}K + \sigma^2_{err}I_n) \notag\\
    &\text{On pose } \theta=(\sigma^2_{phy}, \sigma^2_{err}) \notag \\
    &\text{On définit pour la suite } Var_\theta(Y) = V(\theta) = \sigma^2_{phy}K + \sigma^2_{err}I_n
    \label{eq:VarTheta}
\end{align}

Comme décrit dans \cite{bastideModelesEvolutionCaracteres2022}, l'ajout de cette variance résiduelle dans notre modèle est crucial pour mieux représenter la complexité des données que nous traitons.
En effet, supposer que la seule source de variation entre les observations est le processus stochastique sur l'arbre phylogénétique (specifiée par $\sigma^2_{phy}K$) est souvent peu réaliste, surtout dans des contextes où les données sont hétérogènes ou comme on le verra plus tard, nous avons les données de plusieurs individus d'une même espèce.
C'est d'ailleurs pour ça qu'on peut parler de variation intraspécifique.
Cette hypothèse simplificatrice peut introduire des biais significatifs dans nos analyses, compromettant ainsi la validité des résultats obtenus.
En intégrant la variance résiduelle, qui capture l'effet indépendant de l'environnement sur chaque mesure, notre modèle devient plus flexible et mieux adapté pour tenir compte de la variabilité observée dans les données. 
Le modèle mixte phylogénétique résultant, combinant à la fois la variance phylogénétique et la variance résiduelle, nous permet de distinguer les effets héritables des effets non héritables, offrant ainsi une approche plus nuancée et réaliste de l'analyse comparative des données évolutives.
\newline
\newline
En posant $\lambda = \frac{\sigma^2_{phy}}{\sigma^2_{err}}$ et $E=u+\epsilon$, on peut obtenir une nouvelle forme pour $Y$
\begin{align}
    \label{eq:V_lambda}
    &Y = X\beta + E \text{, où } Var(E)=V(\theta)=\sigma^2_{phy}(K-\lambda I_n)=\sigma^2_{phy}V_\lambda \\
    &E \sim \mathcal{N}_n(0, V_\lambda) \notag
\end{align}

\subsubsection*{Test statistique}
Le test statistique et ses hypothèses sont conservés et on obtient de la même façon que dans la section précédente une nouvelle statistique de test en lien avec l'équation \ref{eq:V_lambda}.
\begin{equation}
    F_{ANOVAphylo-error}=\frac{||\hat{Y} - \bar{Y}||^2_{V_\lambda^{-1}}(n-2)}{||Y - \hat{Y}||^2_{V_\lambda^{-1}}} \underset{\mathcal{H}_0}{\sim}\mathcal{F}\text{isher} (1, n-2)
\end{equation}
\begin{align*}
&\text{Où }||Y - \hat{Y}||^2_{V_\lambda^{-1}} = ||Y - X\hat{\beta}||^2_{V_\lambda^{-1}}= Proj_X^{V_\lambda}Y= (Y-\hat{Y})^TV^{-1}_\lambda(Y-\hat{Y})\\
&\text{et }||\hat{Y} - \bar{Y}||^2_{V_\lambda^{-1}}=(\hat{Y} - \bar{Y})^TV^{-1}_\lambda(\hat{Y}- \bar{Y})
\end{align*}

Il est important de noter que lorsque le paramètre $\lambda$ est connu, l'ANOVA phylogénétique est exacte.
Cependant, dans la pratique, $\lambda$ est généralement inconnu et doit être estimé à partir des données.
 Dans ce cas, l'approximation de la distribution de F par une distribution de Fisher ne tient plus, et il est nécessaire d'utiliser des méthodes alternatives telles que la méthode de Satterthwaite pour estimer les degrés de liberté.
Cette méthode tient compte de l'incertitude associée à l'estimation de $\lambda$ et fournit une approximation plus précise de la distribution de la statistique de test.

\subsection{Approximation de Satterthwaite}
On s'est basé sur la documentation du package \texttt{lmerTest} \cite{kuznetsovaLmerTestPackageTests2017} pour calculer les formules explicites de l'approximation dans notre cadre.
En effet il existe des formules explicite dans le cadre du modèle mixte. 
Dans notre cas, on peut voir l'équation \eqref{eq:eq2err} comme l'équation d'un modèle linéaire mixte où $\beta$ représente tous les paramètres à effets fixes, avec sa matrice de design associée $X$, u les effets aléatoires et $\epsilon$ les résidus.
Dans l'optique de l'implémenter, nous avons calculé la formule explicite de l'approximation de Satterthwaite. 
A partir de \ref{eq:eq2err} on rappelle les valeurs suivantes:
\[
 Y \sim \mathcal{N}_n(X\beta, \sigma^2_{phy}K + \sigma^2_{err}I_n) \text{, }\theta=(\sigma^2_{phy}, \sigma^2_{err}) \text{ et } Var_\theta(Y) = V(\theta) = \sigma^2_{phy}K + \sigma^2_{err}I_n 
\]

De la documentation on obtient alors la covariance suivante:
\begin{equation}
    C(\theta) = (Cov(\beta_i , \beta_j))_{i,j} = (X^TV(\theta)^{-1}X)^{-1} = (X^T(\sigma^2_{phy}K + \sigma^2_{err}I_n)^{-1}X)^{-1}
\end{equation}

\begin{approximation}[F-statistique et approximation de Satterthwaite]
    \begin{align}
        &F_{approx}=\frac{||\hat{Y} - \bar{Y}||^2_{V_\lambda^{-1}}df_{approx}}{||Y - \hat{Y}||^2_{V_\lambda^{-1}}} \underset{\mathcal{H}_0}{\sim}\mathcal{F}\text{isher} (1, df_{approx})\\
        \text{Avec } &df_{approx} = \frac{2(f(\hat{\theta}))^2}{[\nabla f(\hat{\theta})]^T A[\nabla f(\hat{\theta})]} \\
        \text{où } &f(\theta) = l^TC(\theta)l \text{ et A matrice de variance-covariance de } \hat{\theta}=(\hat{\sigma}^2_{phy}, \hat{\sigma}^2_{err}) \notag
    \end{align}
\end{approximation}
\begin{proof}[Calcul explicite de l'approximation]
Toujours en suivant la documentation \cite{kuznetsovaLmerTestPackageTests2017} on part de l'expression pour les degrés de liberté $df$ et de l'approximation. Ce qui nous donne :
\begin{equation}
    df = \frac{2(l^T\hat{C}l)^2}{[Var(l^T\hat{C}l)]}=\frac{2(f(\hat{\theta}))^2}{[Var(f(\hat{\theta}))]}\approx \frac{2(f(\hat{\theta}))^2}{[\nabla f(\hat{\theta})]^T A[\nabla f(\hat{\theta})]}
\end{equation} 
\[\text{où} \quad \hat{C} = C(\hat\theta) \quad \text{et} \quad f(\theta) = l^TC(\theta)l\]
A partir de cette expression, on calcule $\nabla f(\theta)$ qu'on appliquera en $\hat{\theta}$ et $A$ la matrice de variance-covariance de $\hat{\theta}=(\hat{\sigma}^2_{phy}, \hat{\sigma}^2_{err})$ 
\newline
\newline
\textbf{Étape 1 :} Calcul du gradient
\newline
\newline
Nous voulons calculer les dérivées partielles $\partial_{\sigma^2_{phy}}f(\theta)$ et $\partial_{\sigma^2_{err}}f(\theta)$. Pour les premières étapes de calculs, on écrira seulement $\partial$ sans distinction car ce sont les mêmes expressions pour les 2 dérivées. 
On utilisera dans la suite les formules de \cite{petersenMatrixCookbook2012} pour les dérivées de matrice
\[
\partial f(\theta)=l^T\partial C(\theta)l
\]
\[
\partial C(\theta)=\partial (X^TV(\theta)^{-1}X)^{-1} = -C(\theta) \partial (X^TV(\theta)^{-1}X)C(\theta)
\]

\[
    \partial (X^TV(\theta)^{-1}X) = \partial (X^TV(\theta)^{-1})X + \cancel{X^TV(\theta)^{-1}\partial(X)} \quad (\partial_{\sigma^2_{phy}}(X)\text{ et } \partial_{\sigma^2_{err}}(X) \text{ sont nulles})
\]
\[\partial (X^TV(\theta)^{-1}) = \partial(X^T)V(\theta)^{-1} + X^T\partial(V(\theta)^{-1}) = \cancel{\partial(X)^TV(\theta)^{-1}} + X^T\partial(V(\theta)^{-1})
\]
\[\partial (V(\theta)^{-1}) = -V(\theta)^{-1}\partial(V(\theta))V(\theta)^{-1}
\]
\[\partial (V(\theta)) = \partial(\sigma^2_{phy}K + \sigma^2_{err}I_n)
\]
Ce qui donne :
\[\partial_{\sigma^2_{phy}}(V(\theta)) = K, \quad \text{et} \quad \partial_{\sigma^2_{err}}(V(\theta)) = I_n
\]
De là en remettant les formules explicite les unes dans les autres, on obtient :
\[[\nabla f(\hat{\theta})] = \begin{bmatrix} \partial_{\sigma^2_{phy}}f(\hat{\theta}) \\ \partial_{\sigma^2_{err}}f(\hat{\theta}) \end{bmatrix}=\begin{bmatrix} l^TC(\hat{\theta})X^TV(\hat{\theta})^{-1}KV(\hat{\theta})^{-1}XC(\hat{\theta})l \\ l^TC(\hat{\theta})X^TV(\hat{\theta})^{-1}I_nV(\hat{\theta})^{-1}XC(\hat{\theta})l\end{bmatrix}
\]
\newline
\newline
\textbf{Étape 2 :} Calcul de A
\newline
\newline
Par Cramer-Rao on sait que pour les estimateurs non biaisés, $Var(\hat{\theta})\geqslant I(\theta)^{-1}$. $I(\theta)^{-1}$ est l'inverse de l'information de Fisher.
Lorsque la fonction de vraissemblance est assez régulière comme ici, on a $I(\theta)= -\mathbb{E}\left[\frac{\partial^2}{\partial \theta^2} l(Y,\theta)\right]$
La matrice de variance-covariance $Var(\hat{\theta})$ peut alors être approximée par $I(\hat{\theta})^{-1}=A$.
\newline
Dans la plupart des cas il est plus simple d'estimer cette matrice par des méthodes numériques, pour autant une formule explicite de la Hessienne rend le calcul plus rapide et plus robuste quant aux erreurs numériques.  
\newline
On va donc abord calculer la log-vraissemblance du vecteur Y défini précédemment:

\begin{align*}
    l (\bf{Y}, \theta)&= \log (\frac{1}{(2\pi)^{n/2}|V(\theta)|^{1/2}} \exp\left( -\frac{1}{2}(Y - X\beta)^T V(\theta)^{-1} (Y - X\beta) \right)) \\
    &= - \frac{n}{2} \log(2\pi) -\frac{1}{2} \log(|V(\theta)|) - \frac{1}{2}(Y - X\beta)^T V(\theta)^{-1} (Y - X\beta) \\
\end{align*}
On calcule les dérivées premières de la log-vraissemblance
\begin{align*}
    \partial_{\sigma^2_{phy}} l  &= -\frac{1}{2} \partial_{\sigma^2_{phy}}(\log(|V(\theta)|)) - \frac{1}{2} \partial_{\sigma^2_{phy}}((Y - X\beta)^T V(\theta)^{-1} (Y - X\beta))\\
    &= -\frac{1}{2} Tr(V(\theta)^{-1}\partial_{\sigma^2_{phy}}(V(\theta))) - \frac{1}{2} (Y - X\beta)^T \partial_{\sigma^2_{phy}}(V(\theta)^{-1})(Y - X\beta)\\
    &= -\frac{1}{2} Tr(V(\theta)^{-1}K) + \frac{1}{2} (Y - X\beta)^T V(\theta)^{-1} K V(\theta)^{-1}(Y - X\beta)\\
\end{align*}
\begin{align*}
    \partial_{\sigma^2_{err}} l  &= -\frac{1}{2} \partial_{\sigma^2_{err}}(\log(|V(\theta)|)) - \frac{1}{2} \partial_{\sigma^2_{err}}((Y - X\beta)^T V(\theta)^{-1} (Y - X\beta))\\
    &= -\frac{1}{2} Tr(V(\theta)^{-1}\partial_{\sigma^2_{err}}(V(\theta))) - \frac{1}{2} (Y - X\beta)^T \partial_{\sigma^2_{err}}(V(\theta)^{-1})(Y - X\beta)\\
    &= -\frac{1}{2} Tr(V(\theta)^{-1}I_n) + \frac{1}{2} (Y - X\beta)^T V(\theta)^{-1} I_n V(\theta)^{-1}(Y - X\beta)\\
\end{align*}
Puis les dérivées secondes:
\begin{align*}
    \textbf{$\partial_{\sigma^2_{phy}\sigma^2_{phy}} $}l &= -\frac{1}{2} \partial_{\sigma^2_{phy}\sigma^2_{phy}} (Tr(V(\theta)^{-1}K)) + \frac{1}{2} \partial_{\sigma^2_{phy}\sigma^2_{phy}} \left( (Y - X\beta)^T V(\theta)^{-1} K V(\theta)^{-1}(Y - X\beta) \right)\\
    &= -\frac{1}{2} Tr(\partial_{\sigma^2_{phy}\sigma^2_{phy}} (V(\theta)^{-1})K) + \frac{1}{2}  (Y - X\beta)^T \partial_{\sigma^2_{phy}\sigma^2_{phy}}(V(\theta)^{-1} K V(\theta)^{-1})(Y - X\beta) \\
    &= \frac{1}{2} Tr(V(\theta)^{-1}KV(\theta)^{-1}K) -  (Y - X\beta)^T V(\theta)^{-1}KV(\theta)^{-1} K V(\theta)^{-1}(Y - X\beta) \\
\end{align*}
\[\text{car} \quad \partial ((Y - X\beta)^T V(\theta)^{-1} K V(\theta)^{-1}(Y - X\beta)) = (Y - X\beta)^T \partial ( V(\theta)^{-1} K V(\theta)^{-1})(Y - X\beta)\]
\[\text{et} \quad \partial \left( V(\theta)^{-1} K V(\theta)^{-1} \right) = -V(\theta)^{-1} \partial V(\theta) V(\theta)^{-1} K V(\theta)^{-1} - V(\theta)^{-1} K V(\theta)^{-1} \partial V(\theta) V(\theta)^{-1}\]
\[\text{ce qui donne} \quad\partial_{\sigma^2_{phy}\sigma^2_{phy}} \left( V(\theta)^{-1} K V(\theta)^{-1} \right) = -2V(\theta)^{-1} K V(\theta)^{-1} K V(\theta)^{-1}\]
\begin{align*}
    \textbf{$\partial_{\sigma^2_{err}\sigma^2_{phylo}}$}l& = \textbf{$\partial_{\sigma^2_{phy}\sigma^2_{err}}$}l \\
    &= -\frac{1}{2} Tr(\partial_{\sigma^2_{phy}\sigma^2_{err}} (V(\theta)^{-1})K) + \frac{1}{2}  (Y - X\beta)^T \partial_{\sigma^2_{phy}\sigma^2_{err}}(V(\theta)^{-1} K V(\theta)^{-1})(Y - X\beta) \\
    &= \frac{1}{2} Tr(V(\theta)^{-1}I_nV(\theta)^{-1}K) - \frac{1}{2} (Y - X\beta)^T (V(\theta)^{-1}V(\theta)^{-1}KV(\theta)^{-1} + \\
    & V(\theta)^{-1}KV(\theta)^{-1}V(\theta)^{-1})(Y - X\beta) \\
\end{align*}
\[\text{car} \quad \partial_{\sigma^2_{phy}\sigma^2_{err}} \left( V(\theta)^{-1} K V(\theta)^{-1} \right) = -(V(\theta)^{-1}V(\theta)^{-1}KV(\theta)^{-1} + V(\theta)^{-1}KV(\theta)^{-1}V(\theta)^{-1})\]
\begin{align*}
    \textbf{$\partial_{\sigma^2_{err}\sigma^2_{err}}$} l  &= -\frac{1}{2} \partial_{\sigma^2_{err}\sigma^2_{err}} (Tr(V(\theta)^{-1})) + \frac{1}{2} \partial_{\sigma^2_{err}\sigma^2_{err}} \left( (Y - X\beta)^T V(\theta)^{-1} V(\theta)^{-1}(Y - X\beta) \right)\\
    &=\frac{1}{2}Tr(V(\theta)^{-1}V(\theta)^{-1}) - (Y - X\beta)^T V(\theta)^{-1}V(\theta)^{-1}V(\theta)^{-1}(Y - X\beta)
\end{align*}

De là on obtient la Hessienne $\begin{pmatrix}
    \partial_{\sigma^2_{phy}\sigma^2_{phy}}l & \partial_{\sigma^2_{phy}\sigma^2_{err}}l \\
    \partial_{\sigma^2_{err}\sigma^2_{phy}}l & \partial_{\sigma^2_{err}\sigma^2_{err}}l \\
\end{pmatrix}$ puis A en inversant la matrice de Fisher liée, ce qui peut se faire par des méthodes numériques ou analytiques ici car cela signifie inverser une matrice 2 x 2 ce qui est facile.


\end{proof}

\subsection{REML}
Le REML, ou Maximum de Vraisemblance Restreint (Restricted Maximum Likelihood en anglais), est une méthode statistique utilisée dans l'estimation des paramètres de modèles linéaires mixtes (ou modèles à effets mixtes) et dans l'analyse de la variance (ANOVA).
Il s'agit d'une approche alternative à la méthode de maximum de vraisemblance (ML) standard, notamment lorsque l'on travaille avec des modèles à effets aléatoires.
\newline
\newline
L'une des formules pour la vraissemblance restreinte consite à regarder la vraissemblances des observations en intégrant sur les effets fixes, ici $\beta$ sur lesquels on a mis un prior impropre uniforme entre moins l'infini et plus l'infini.
C'est à dire sous l'hypothèse que les estimations des effets fixes sont éliminées (conditionnées)
\begin{equation*}
    L_{REML} (Y;\theta)=\int _{\beta \in \mathbb{R} ^2} \frac{1}{(2\pi)^{n/2}|V(\theta)|^{1/2}} \exp\left( -\frac{1}{2}(Y - X\beta)^T V(\theta)^{-1} (Y - X\beta) \right)d\beta
\end{equation*}
Lorsque dans un estimateur on a la division par $(n-p)$ au lieu de n, on fait sans le dire un REML dans le cas de la régression linéaire classique. 
\newline
Satterthewaite maximum likelihood ne fait pas l'astuce de diviser par $n-p$, au contraire de l'anova classique ou phylogénétique. Il est donc particulièrement mauvais. Le REML permet de mieux estimer la variance, et aussi le $\lambda$ (ratio des variances), ce qui est crucial dans la dérivation des degrés de libertés approchés.

\subsection{Méthode Likelihood Ratio Test}
La méthode du test de rapport de vraisemblance (LRT) est une technique statistique utilisée pour comparer deux modèles statistiques et déterminer s'ils diffèrent significativement en termes d'ajustement aux données.
Le rapport de vraisemblance est calculé comme le rapport des vraisemblances maximales de deux modèles (sous les deux hypothèses) emboîtés :
\begin{equation*}
    \text{LRT} = -2 \log\left(\frac{L(\theta_{\text{H}_0})}{L(\theta_{\text{H}_1})}\right)
    \end{equation*}
Sous l'hypothèse nulle que le modèle plus simple est correct, ce rapport suit approximativement une distribution du chi-deux avec un nombre de degrés de liberté égal à la différence dans le nombre de paramètres entre les deux modèles.
Ainsi, en comparant la valeur observée du rapport de vraisemblance à la distribution du chi-deux, on peut décider si l'ajout de paramètres dans le modèle conduit à une amélioration significative de l'ajustement aux données.
\newline
Il est à noter que cette méthode est plus coûteuse car il est nécessaire d'ajuster 2 modèles au lieu d'un seul dans les autres méthodes. 


\section{Simulations}
\label{chap:metho}
% Ou faire une partie à part entière avec 
% 1) ANOVA vs ANOVA phylo sans correction des degrés de liberté
%  b) avec une sous partie sur le REML

% 2) ANOVA phylo avec approximation de SAtterthwaite 
%  a) prez
%  a`) simulation et résultats 
%  b) instabilités numériques -> correction avec la Hessienne ? 
%  c) La hessienne analytique ? A voir si besoin d'une partie supplémentaire
% On importe le fichier

<<simulations-methodes, child='Rnw/simulations-methodes.Rnw'>>=
@

\section{Application aux données réelles}
\label{sec:data}
% Présentation des données utilisées.
<<donnees-reelles, child='Rnw/donnees-reelles.Rnw'>>=
@

\section{Conclusions sur le projet}
\label{sec:conclusion}
% Analyse critique des résultats, limites, perspectives.

Le travail portait sur l'élaboration de nouvelles méthodes statistiques 
spécifiquement conçues pour traiter les données phylogénétiques. Nous 
avons accompli cela en adaptant des approches classiques déjà existantes 
à ce contexte particulier. Pour y parvenir, nous avons approfondi notre 
compréhension de ces diverses méthodes en consultant la littérature spécialisée
et en échangeant de manière intensive avec Paul Bastide et Mélina Gallopin.
\newline
Ce projet nous a offert l'opportunité de renforcer nos connaissances théoriques 
en explorant diverses méthodes, tout en approfondissant notre compréhension de 
leurs relations et de la manière de formaliser un modèle.
Nous avons également effectué des calculs et programmé en R, utilisant à la fois
des packages existants et en développant nos propres fonctionnalités. De plus, 
nous avons analysé et interprété une variété de graphiques pour comparer les 
différentes méthodes.
Ce travail s'est déroulé en binôme, permettant ainsi d'exploiter les compétences
individuelles de chaque membre pour atteindre nos objectifs de manière 
collaborative.
\newline
Nos recherches ont abouti à des conclusions intéressantes concernant les 
méthodes d'ANOVA phylogénétique.
Nous avons observé que l'ANOVA phylogénétique surpasse l'ANOVA classique en 
termes de contrôle de l'erreur de type 1.
En revanche, l'utilisation isolée de la méthode de Satterthwaite présente des 
limitations, ce qui souligne l'importance d'intégrer REML dans nos analyses 
pour obtenir des résultats plus fiables.
\newline
Au cours de notre projet, plusieurs obstacles ont été identifiés. Tout d'abord, 
nous avons rencontré des difficultés avec la convergence de l'algorithme 
d'optimisation lors du calcul approximatif de la Hessienne dans l'approche de 
Satterthwaite. Comprendre le contexte biologique s'est avéré être un élément 
crucial pour le développement de méthodes statistiques appropriées. En outre, 
l'importance d'écrire un code propre, clair et réutilisable a été mise en 
évidence pour assurer l'efficacité et la durabilité de nos travaux.

De nombreux points restent à explorer, nous en proposons quelques-uns 
ci-dessous.\newline

Dans tout le rapport nous avons supposé que le processus de dérive était un
mouvement brownien mais nous aurions pu utiliser un autre processus comme le
processus d'Ornstein-Uhlenbeck. Ce dernier permet de modéliser des traits 
tendant vers une valeur, ce qui peut correspondre par exemple à un optimum 
écologique comme sur la figure~\ref{fig:OrnsteinUhlenbeck}.\newline

\begin{figure}[H]
    \centering
    \includegraphics[width=0.5\textwidth]{OrnsteinUhlenbeck3.png}
    \caption{Exemple de processus d'Ornstein-Uhlenbeck (tiré de Wikipédia)}
    \label{fig:OrnsteinUhlenbeck}
\end{figure}

Un autre point que l'on pourrait considérer, c'est que le LRT est un test
asymptotique, contrairement à l'ANOVA quequi est un test exact (sous les bonnes
hypothèses). Cela peut expliquer pourquoi leil est un peu moins bon sur 
les simulations.\newline
Par contre, il dépend moins des hypothèses (le ratio des vraisemblance converge
sous des hypothèses faible). Cela peut peut-être expliquer pourquoi il est plus
robuste dans le cas des données rééelles, qui n'ont pas de raison de suivre le 
bon modèle.

\newpage

% Bibliographie
\printbibliography
\nocite{*}

\appendix
\section{Application aux données réelles}

Comme nous l'avons remarqué dans la section~\ref{sec:data} l'application de la
méthode EVEmodel a produit des valeurs manquantes pour les gènes présentés dans
le tableau suivant.

\begin{table}[H]
    \centering
    <<'table_nas'>>=
    knitr::kable(evegenesNA, col.names = "Gènes ayant produits des NA", 
        align = "c", booktabs = TRUE, format = "latex", escape = TRUE)
    @
    \caption{Table des gènes pour lesquels la méthode \texttt{EVEmodel} a produit des NA}
    \label{tab:na-evemodel}
\end{table}

\end{document}