1. Moyennes

2. Variances

En ANOVA univariée, la variabilité est résumée par un scalaire (la variance).
En MANOVA, la variabilité est décrite par une matrice de covariance de taille \(D\times D\).

3. Décomposotion de la variance

L’idée est de décomposer la dispersion totale en une composante et une composante . En MANOVA, ces quantités sont des (Sum of Squares and Cross-Products).

1. Ecriture initiale

Pour chaque expérience \((i, j)\) (\(i\)-ième observation réalisée au niveau \(j\) de \(F\)), on écrit le vecteur aléatoire réponse \(Y_{ij}\) de \(\mathbb{R}^D\) sous la forme : \[ Y_{ij} = \beta_j + U_{ij} \]

Les trois éléments de cette écriture doivent être vus comme des vecteurs-lignes de \(\mathbb{R}^D\), comme précisé ci-dessous.

• Le vecteur \(\beta_j = (\beta_j^1 ...\beta_j^d ... \beta_j^D)\) est un paramètre à estimer ; il modélise la valeur de la réponse \(Y\) au niveau \(j\) de \(F\).

• Le terme \(U_{ij} = (U_ij^1 ... U_{ij}^D)\) est le vecteur aléatoire des erreurs. On suppose que les \(U_{ij}\) sont i.i.d., de loi \(\mathcal{N}(0_D,\Sigma)\), où \(\Sigma\) est une matrice symétrique et strictement définie-positive ; on doit également estimer \(\Sigma\). On notera que \(\Sigma\) ne dépend pas de \(j\), autrement dit on est toujours dans le cadre d’un modèle homoscédastique.

• Les vecteurs aléatoires \(Y_{ij}\) sont donc indépendants, de loi \(\mathcal{N}_D(\beta_j,\Sigma)\). Finalement, il y a \(J \times D\) paramètres de moyenne \(\beta_j^d\) à estimer, ainsi que \(\frac{D(D+1)}{2}\) paramètres de variance \((\Sigma)_d^{d^{\prime}} (1 \leq d \leq D ; 1 \leq d' \leq D)\). Comme on dispose de \(nD\) observations, on doit veiller à ce que la taille \(n\) de l’échantillon utilisé vérifie : \(n \geq J + \frac{D+1}{2}\)

2. Écriture matricielle

L’ensemble des \(nD\) observations réalisées peut se mettre sous la forme matricielle suivante : \[Y = X\beta + U\]

Dans l’écriture ci-dessus, \(X\) et \(\beta\) sont des matrices réelles (non aléatoires) de dimensions respectives \(n \times J\) et \(J \times D\).

Comme dans le cas unidimensionnel, les colonnes de la matrice d’incidence \(X\) sont les indicatrices \(Z^j\) des niveaux du facteur \(F\), de sorte que \(X\) ne comporte que des 0 et des 1.

Les termes \(Y\) et \(U\) sont des matrices aléatoires de dimension \(n \times D\). Elles sont gaussiennes et vérifient \(\mathbb{E}(U) = 0_{n \times D}; \quad \mathbb{E}(Y) = X\beta; \quad \text{Var}(Y) = I_n \otimes \Sigma\) où \(I_n\) est matrice identité d’ordre \(n\) et \(\otimes\) le produit matriciel direct, ou produit de Kronecker.

Afin d’estimer les paramètres, deux méthodes peuvent être utilisées :

• Paramétrage centré

Ce paramétrage consiste à décomposer chaque vecteur-ligne \(\beta_j\) sous la forme :

\[ \beta_j=\mu+\alpha_j \text{ avec } \mu=\frac{1}{J}\sum_{j=1}^J\beta_j\text{, } \alpha_j=\beta_j-\mu \text{ et } \sum_{j=1}^J\alpha_j=0_D \]

• Paramétrage SAS

Pour ce paramétrage, on pose \(m = \beta_J\) et \(a_j = \beta_j - \beta_J\) (de sorte que, encore une fois, \(a_J = 0_D\)). Les paramètres \(m\) et \(a_j\) sont également des vecteurs de \(\mathbb{R}^D\).

1. Vraisemblance et log-vraisemblance

la vraisemblance de l’échantillon des \(y_{ij}\) est

\[ \begin{aligned} L(y_{ij},\beta,\Sigma)&=\prod_{j=1}^J\prod_{i=1}^{n_j} \frac{1}{(2\pi)^{D/2}(\det\Sigma)^{1/2}} \exp\left[- \frac12(y_{ij}-\beta_j)\Sigma^{-1}(y_{ij}- \beta_j)^\prime\right]\\ &=C_1(\det \Sigma)^{- n/2} \exp\left[- \frac12\sum_{j=1}^{J}\sum_{i=1}^{n_j}(y_{ij}-X_j\beta)\Sigma^{-1}(y_{ij}- \beta_j)^\prime\right] \end{aligned} \]

Où \(C_1\) est une constante et \(X_j\) un vecteur à \(J\) éléments, comportant un comportant un 1 en \(j\)-ième colonne et des 0 partout ailleurs (en fait, \(X_j\) est n’importe laquelle des lignes de la matrice \(X\) correspondant aux observations du niveau \(j\) de \(F\)).

La log-vraisemblance s’écrit

\[ \begin{aligned} l(y_{ij},\beta,\Sigma)&=\log[L(y_{ij},\beta,\Sigma)]\\ &=C_2-\frac{n}{2}\log(\det\Sigma)-\frac12\|Y-X\beta\|^2_{I_n,\Sigma^{-1}} \end{aligned} \]

Où \(C_2=\log(C_1)\) et la norme d’une matrice \(A\) (de dimension \(n \times p\) est définie par

\[ \|A\|^2_{N,P}=\text{tr}(APA'N) \]

2. Estimation du maximum de vraisemblance

En maximisant la log-vraisemblance, on obtient \(\hat{B}\) l’estimateur de \(\beta\) et qui est défini par \(\hat{B}=(X^\prime X)^{-1}X^{\prime}Y\)

Pour un niveau \(j\) donné, et pour toute observation \(i\) faite à ce niveau, la valeur prédite correspondante est : \(\hat{y}_{ij}=\hat{\beta}_j=\bar{y}_{\cdot j}\) (vecteur de \(\mathbb{R}^D\)). On note \(\hat{Y}\) la matrice aléatoire \(n\times D\) de l’ensemble des valeurs prédites.

Les résidus sont donnés par \(\hat{u}_{ij}=y_{ij}-\hat{y}_{ij}=y_{ij}-\bar{y}_{\cdot j}\) (vecteur de \(\mathbb{R}^D\)). On note \(\hat{U}\) la matrice aléatoire \(n\times D\) des résidus ainsi définis.

La matrice de covariance est estimée par \(\hat{\Sigma}=\frac1{n-J}\hat{U}^\prime \hat{U}\) (matrice \(D\times D\))

3. Propirétés des estimateurs du maximum de vraisemblance

• Les matrices \(\hat{B}\) , de dimension \(J \times D\), \(\hat{Y}\), de dimension \(n \times D\), et \(\hat{U}\) , de dimension \(n \times D\), sont des matrices aléatoires gaussiennes, d’espérances respectives \(\beta \text{, } X\beta \text{ et } 0n\times D\).

• La matrice aléatoire \(\hat{U}\) est indépendante des matrices aléatoires \(\hat{B}\) et \(\hat{Y}\).

• Enfin, \((n-j)\hat{\Sigma}=\hat{U}^\prime\hat{U}\) est une matrice aléatoire distribuée selon une loi de Wishart (une sorte de généralisation multidimentionnelle de la loi du khi-deux) de dimension \(D\), à \(n-J\) degrés de liberté et de matrice associée \(\sigma\).

1. Le test de Wilks

Il s’agit du test le plus courant dans le contexte de la MANOVA qui est, en fait, une adaptation du test du rapport des vraisemblances. Notons \(\theta\) le vecteur de tous les paramètres du modèle, \(\hat\theta\) son estimation par le maximum de vraisemblance et \(\hat{\theta}_0\) son estimation par le maximum de vraisemblance sous la contrainte définie par \(\text{H}_0\).

La statistique du test du rapport des vraisemblances est \(\frac{L(y,\hat{\theta}_0)}{L(y,\hat{\theta})}\), dont on peut vérifier qu’elle vaut : \([\frac{\det(E)}{\det(H+E)}]^{n/2}\). Le test de Wilks consiste à considérer la puissance \(2/n\) de cette quantité, autrement dit sa statistique est définie par : \[ \Lambda=\frac{\det(E)}{\det(H+E)}=\prod_{k=1}^{s}\frac{1}{1+\lambda_k} \]

Où les \(\lambda_k\) sont les valeurs propres de la matrice \(HE^{-1}\) et \(s =\inf(D, J - 1)\) est le nombre de valeurs propres non nulles de cette matrice.

On peut considerer 3 cas particuliers et généraliser pour les autres cas.

2. Autres tests

Dans la littérature statistique, on trouve d’autres tests permettant d’éprouver la même hypothèse nulle.

1. Effectifs et dimension

On note \(n_{rs}\) l’effectif de la cellule \((r,s)\), et l’effectif total : \[ n=\sum_{r=1}^{a}\sum_{s=1}^{b} n_{rs}. \] On définit également les effectifs marginaux : \[ n_{r\cdot}=\sum_{s=1}^{b}n_{rs},\qquad n_{\cdot s}=\sum_{r=1}^{a}n_{rs}. \]

Chaque individu est observé sur \(D\) variables quantitatives (profil multivarié). La réponse est donc de dimension \(D\).

On reprend notre matrice d’observations \(\mathbf{Y}\) de taille \(n\times D\) : \[ \mathbf{Y}= \begin{pmatrix} y_{11} & y_{12} & \cdots & y_{1D}\\ y_{21} & y_{22} & \cdots & y_{2D}\\ \vdots & \vdots & \ddots & \vdots\\ y_{n1} & y_{n2} & \cdots & y_{nD} \end{pmatrix} \in \mathbb{R}^{n\times D}. \]

L’appartenance à une cellule \((r,s)\) est une information qualitative associée à chaque individu : elle permet de définir les moyennes de cellule et les moyennes marginales présentées ci-dessous.

2. Moyennes

3. Variabilité multivariée

Comme dans toute MANOVA, la variabilité conjointe entre les \(D\) variables est décrite par une matrice de covariance \(\boldsymbol{\Sigma}\in\mathbb{R}^{D\times D}\).

4. Décomposition de la variabilité (SSCP)

Ces quantités sont des matrices SSCP de taille \(D\times D\).

1. Ecriture scalaire du modèle

Le modèle MANOVA à deux facteurs sans interaction s’écrit, pour chaque observation : \[ Y_{ijk} = \mu + \alpha_i + \beta_j + U_{ijk} \text{, où } U_{ijk} \overset{\text{i.i.d}}{\sim} \mathcal{N}_D(0_D, \Sigma) \text{ (1.1)} \]

Avec :

• \(\mu = (\mu^1,..., \mu^D)^\prime \in \mathbb{R}\) D : vecteur moyenne générale ;

• \(\alpha_i \in \mathbb{R}^D\) : effet principal du niveau \(i\) du facteur \(A\) ;

• \(\beta_j \in \mathbb{R}^D\) : effet principal du niveau \(j\) du facteur \(B\) ;

• \(U_{ijk}\) : vecteur d’erreur aléatoire, i.i.d. \(\mathcal{N}_D(0_D, \Sigma)\).

La matrice \(\Sigma (D \times D)\) est symétrique définie positive.

Le modèle (1.1) est sur-paramétré. Pour assurer l’unicité des estimateurs, on impose les contraintes de somme nulle (vecteur nul de \(\mathbb{R}^D\) :

\[ \sum_{i=1}^{a} \alpha_i=0_D \text{ } \sum_{j=1}^{b} \beta_j=0_D \text{ (1.2)} \] Ainsi, le vecteur moyenne théorique de la cellule \((i, j)\) est \(\mu_{ij} = \mathbb{E}[Y_{ijk}] = \mu + \alpha_i + \beta_j\).

2. Ecriture matricielle du modèle

L’ensemble des \(n\) observations et des paramètres se représentent par des matrices.

• Matrice des observations : \(Y \in \mathbb{R}^{n \times D}\)

• Matrice d’incidence : \(X \in \mathbb{R}^{n \times p}\), où \(p = 1 + (a − 1) + (b − 1) = a + b − 1\) est le nombre de colonnes nécessaires pour coder le modèle (par exemple, une colonne pour la moyenne générale, \(a − 1\) colonnes pour \(A\), \(b − 1\) colonnes pour \(B\)).

• Matrice des paramètres : \(\beta \in \mathbb{R}^{p \times D}\) (chaque colonne correspond à une variable réponse).

• Matrice des erreurs : \(U \in \mathbb{R}^{n \times D}\), avec \(\mathbb{E}[U] = 0_{n \times D}\) et \(Var(vec(U))= \Sigma \otimes I_n\), où \(\otimes\) désigne le produit de Kronecker.

Le modèle s’écrit alors sous la forme compacte : \[ Y = X\beta + U \text{, } \mathbb{E}[Y] = X\beta \text{, } Var(vec(Y))= \Sigma \otimes I_n \text{ (1.3)} \]

Remarque : Une autre paramétrisation possible est la paramétrisation cellmeans où contient directement les \(ab\) moyennes \(\mu_{ij}\) . Dans ce cas, \(X\) est la matrice d’indicatrices des cellules (dimension \(n \times ab\)), et les contraintes (1.2) servent ensuite à définir les effets.

1. Estimateurs des paramètres

Pour chaque variable réponse \(d = 1,..., D\), le modèle univarié \(y_d = X\beta_d + u_d\) (où \(y_d\) est la \(d\)-ième colonne de \(Y\)) admet l’estimateur des moindres carrés ordinaires : \(\hat{\beta}=(X^\prime X)^{-1}X^\prime y_d\) où \((X^\prime X)^{-1}\) désigne une inverse généralisée. En empilant les colonnes, on obtient la matrice des estimateurs.

• \(\hat{\beta} \in \mathbb{R}^{p \times D}\)

• \(\hat{\beta}\) est sans biais : \(\mathbb{E}[\hat{\beta}]=\beta\)

Les valeurs ajustées sont \(\hat{Y}=X\hat{\beta}=P_X Y\) où \(P_X = X(X^\prime X)^{-1}X^\prime\) est le projecteur orthogonal (par rapport au produit scalaire usuel) sur l’espace engendré par les colonnes de \(X\). La matrice des résidus est \(\hat{U} = Y − \hat{Y} = (I_n − P_X)Y\).

L’estimateur sans biais de \(\Sigma\) est \(\hat{\Sigma}=\frac{1}{n-p} \hat{U}^{\prime} \hat{U}=\frac{1}{n-p}E\) où \(p = \text{rang}(X)\) et \(E = \hat{U}^{\prime} \hat{U}\) est la matrice de dispersion résiduelle.

2. Propriétés des estimateurs

• Les estimateurs \(\hat{\beta}\) et \(\hat{U}\) (donc \(\hat{\Sigma}\)) sont indépendants.

• Sous les hypothèses du modèle,

  • vec\((\hat{\beta})\sim \mathcal{N}_{pD}(\text{vec}(\beta), \Sigma \otimes (X^{\prime}X)^{-1})\)

  • \(E =\hat{U}^{\prime}\hat{U} \sim \mathcal{W}_D(n − p, \Sigma)\) (loi de Wishart)

  • \(\hat{\beta}\) et \(E\) sont indépendantes.

1. Formulation des hypothèses

On s’intéresse aux deux hypothèses nulles classiques :

\[ \text{H}_0^A \text{: } \alpha_1 =...= \alpha_a = 0_D \text{, } \text{H}_0^B \text{: } \beta_1 =...= \beta_b = 0_D \]

2. Statistiques de test

Pour tester \(\text{H}_0^A\), on utilise la matrice \(H_A\) définie précédemment. Sous \(\text{H}_0^A\), \(H_A\) mesure la dispersion due à \(A\) et doit être “petite” par rapport à \(E\). La matrice \(E\) est la même pour tous les tests. On note \(\nu_H\) le degré de liberté de l’hypothèse (\(\nu_H = a−1\) pour \(A, b − 1\) pour \(B\)) et \(\nu_E = n − (a + b − 1)\).

Soit \(s = \min(D, \nu_H)\). Soient \(\lambda_1 \geq \lambda_2 \geq ... \geq \lambda_s > 0\) les valeurs propres de \(E^{-1}H\) (ou les racines de \(|H − \lambda E| = 0\)). Les quatre statistiques classiques sont : \[ \begin{aligned} \Lambda_{\text{Wilks}} &= \prod_{i=1}^{s} \frac{1}{1+\lambda_i} =\frac{\text{det}(E)}{\text{det}(H+E)}\\ \Lambda_{\text{Pillai}} &= \sum_{i=1}^{s} \frac{\lambda_i}{1+\lambda_i} =\text{tr}[H(H+E)^{-1}]\\ U_{\text{Hotelling-Lawley}} &= \sum_{i=1}^{s} \lambda_i =\text{tr}(E^{-1}H)\\ \Theta_{\text{Roy}} &= \frac{\lambda_1}{1+\lambda_1}\\ \end{aligned} \]

Ces statistiques sont des fonctions monotones des \(\lambda_i\). Par exemple, \(\Lambda\) est d’autant plus petit que les \(\lambda_i\) sont grands, ce qui conduit au rejet de \(\text{H}_0\).

3. Distributions sous H0

Sous l’hypothèse nulle, on a les résultats suivants : \(E \sim \mathcal{W}_D (\nu_E, \Sigma)\), \(H \sim \mathcal{W}_D (\nu_H, \Sigma)\) et \(E\) et \(H\) sont indépendantes.

• Si \(\nu_H = 1\), alors :

\[ \frac{1- \Lambda}{\Lambda} . \frac{\nu_E -D+1}{D} \sim F_{D,\nu_E -D+1} \]

• si \(\nu_H=2\), alors :

\[ \frac{1- \sqrt{\Lambda}}{\sqrt{\Lambda}} . \frac{\nu_E -D+1}{D} \sim F_{2D,2(\nu_E -D+1)} \] * Si \(D=2\), alors :

\[ \frac{1- \sqrt{\Lambda}}{\sqrt{\Lambda}} . \frac{\nu_E -1}{\nu_H} \sim F_{2\nu_H,2(\nu_E -1)} \]

Pour les cas généraux, on utilise l’approximation de Rao pour Wilks : \[ \frac{1-\Lambda^{1/t}}{\Lambda^{1/t}}.\frac{f_t-g}{D_{\nu_H}}\overset{\text{approx}}{\sim} F_{D_{\nu_H},f_t-g} \] Avec \(f=\nu_E+\nu_H-\frac{\nu_H+D+1}{2}\), \(g=\frac{D_{\nu_H}}{2}-1\), \(t=\sqrt{\frac{D^2\nu^2_H-4}{D^2+\nu^2_H-5}}\). De plus, asymptotiquement (\(n \to +\infty\), \(\nu_H\) fixé) : \(\nu_E \log \Lambda \overset{\mathcal{L}}{\to} \chi^2_{D\nu_H}\),\(\nu_E V \overset{\mathcal{L}}{\to} \chi^2_{D\nu_H}\), \(\nu_E U \overset{\mathcal{L}}{\to} \chi^2_{D\nu_H}\).

1. Ecriture scalaire du modèle

Le modèle avec interaction s’écrit : \[ Y_{ijk} = \mu + \alpha_i + \beta_j + \gamma_{ij} + U_{ijk} \text{, } U_{ijk} \overset{\text{i.i.d.}}{\sim} \mathcal{N}_D(O_D, \Sigma) \] Où \(\gamma_{ij} \in \mathbb{R}^D\) est l’effet d’interaction entre le niveau \(i\) de \(A\) et le niveau \(j\) de \(B\).

Pour assurer l’unicité des estimateurs, on impose les contraintes de somme nulle :

\[ \sum_{i=1}^{a} \alpha_i=0_D \text{, } \sum_{j=1}^{b} \beta_j=0_D \text{, } \sum_{i=1}^{a} \gamma_{ij}=O_D \text{ }(\forall j) \text{, } \sum_{j=1}^{b} \gamma_{ij}=O_D \text{ } (\forall j) \] Ainsi, le vecteur moyenne théorique de la cellule \((i,j)\) est \(\mu_{ij} = \mathbb{E}[Y_{ijk}] = \mu+ \alpha_i+\beta_j + \gamma_{ij}\) .

2. Ecriture matricielle du modèle

On peut toujours écrire \(Y=X\beta+U\). Ici, on utilise la paramétrisation cell-means où \(\beta\) est la matrice \(ab \times D\) contenant les \(\mu_{ij}\) , et \(X\) est la matrice d’indicatrices des cellules (dimension \(n \times ab\)). Les contraintes servent ensuite à définir les effets à partir des \(\mu_{ij}\).

1. Estimateurs des paramètres

L’estimateur des moindres carrés est toujours \(\hat{\beta}=(X^{\prime}X)^{-1}X^{\prime}Y\). Dans la paramétrisation cell-means, \(X^{\prime}X\) est diagonale par blocs si les observations sont ordonnées par cellule, et son inverse est simple. En pratique, on calcule directement les moyennes des cellules : \(\hat{\mu}_{ij}=\bar{y}_{ij}\).

La matrice résiduelle est \(\hat{U}=Y-X\hat{\beta}\). L’estimateur sans biais de \(\Sigma\) est \(\hat{\Sigma}\frac{1}{n-ab}\hat{U}^{\prime}\hat{U}=\frac{1}{n-ab}\), car \(\text{rang}(X) = ab\) (si toutes les cellules sont observées).

2. Propriétés des estimateurs

Comme précédemment, \(\hat{\beta}\) et \(\hat{U}\) sont indépendants, \(\hat{\beta}\) est sans biais, et \(E \sim \mathcal{W}_D(n-ab, \Sigma)\).

1. Formulation des hypothèses

On teste maintenant trois hypothèses :

• \(H_0^{AB}\) : \(\gamma_{ij}=0_D\) pour tout \(i\), \(j\) (pas d’interaction) ;

• \(H_0^{A}\) : \(\alpha_{i}=0_D\) pour tout \(i\) (pas d’effet principal de \(A\)) ;

• \(H_0^{A}\) : \(\alpha_{i}=0_D\) pour tout \(i\) (pas d’effet principal de \(A\)) ;

• \(H_0^{B}\) : \(\beta_{j}=0_D\) pour tout \(j\).

Pour \(H_0^{AB}\), on utilise \(H_{AB}\) avec \(\nu_H = (a − 1)(b − 1)\). Pour \(H_0^{A}\), on utilise \(H_A\) avec \(\nu_H = a − 1\), et pour \(H_0^{B}\), \(H_B\) avec \(\nu_H = b − 1\). Dans tous les cas, \(\nu_E = n − ab\).

2. Statistiques de test

Les statistiques de test sont les mêmes que dans la partie I, avec les degrés de liberté appropriés. Sous chaque hypothèse nulle, on a \(H \sim \mathcal{W}_D(\nu_H, \Sigma)\) et \(E \sim \mathcal{W}_D(\nu_E, \Sigma)\), indépendantes. Les transformations exactes en \(F\) et les approximations décrites précédemment s’appliquent en remplaçant \(\nu_H\) et \(\nu_E\) par leurs valeurs.

1.3.1 valeurs manquantes

# Compte NA par variable (robuste)
na_count <- colSums(is.na(df))
na_count

##     Zone  Especes  Hauteur Diametre 
##        0        0        0        0

On voit qu’uncune variable ne présnete de valeurs manquantes. On peut donc passer à l’analyse des valeurs abérrantes

1.3.2 valeurs abérrantes

# Top 8 espèces
k <- 8
top_especes <- df %>%
  count(Especes, sort = TRUE) %>%
  slice_head(n = k) %>%
  pull(Especes)

df_top <- df %>% filter(Especes %in% top_especes)

# Graphiques
p1 <- ggplot(df, aes(x = Zone, y = Hauteur, fill = Zone)) +
  geom_boxplot(alpha = 0.8) +
  labs(title = "Hauteur par Zone", x = "Zone", y = "Hauteur") +
  theme_minimal() + guides(fill = "none")

p2 <- ggplot(df, aes(x = Zone, y = Diametre, fill = Zone)) +
  geom_boxplot(alpha = 0.8) +
  labs(title = "Diamètre par Zone", x = "Zone", y = "Diamètre") +
  theme_minimal() + guides(fill = "none")

p3 <- ggplot(df_top, aes(x = reorder(Especes, Hauteur, FUN = median), y = Hauteur, fill = Especes)) +
  geom_boxplot(alpha = 0.8) +
  labs(title = "Hauteur par Espèces (Top 8)", x = "Espèces", y = "Hauteur") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  guides(fill = "none")

p4 <- ggplot(df_top, aes(x = reorder(Especes, Diametre, FUN = median), y = Diametre, fill = Especes)) +
  geom_boxplot(alpha = 0.8) +
  labs(title = "Diamètre par Espèces (Top 8)", x = "Espèces", y = "Diamètre") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  guides(fill = "none")

# Carte de graphiques 2x2
gridExtra::grid.arrange(p1, p2, p3, p4, ncol = 2)

L’analyse des valeurs aberrantes par groupe (zone/espèce) montre que la présence d’observations extrêmes n’est pas uniforme : au niveau des , on observe des points atypiques surtout pour la (notamment en hauteur, avec des valeurs très élevées) et, dans une moindre mesure, pour la (quelques diamètres élevés), tandis que la présente globalement une dispersion plus large, en particulier pour le diamètre. Au niveau des (Top 8), certaines modalités se distinguent par une variabilité plus marquée et/ou des valeurs élevées : (et dans une moindre mesure et ) apparaît associée à des hauteurs et diamètres plus importants, avec des distributions plus étalées, alors que d’autres espèces (ex. , ) montrent des profils plus concentrés.

Ces observations indiquent que la base contient des points extrêmes potentiellement influents, mais qui peuvent aussi refléter une hétérogénéité réelle entre groupes (et non nécessairement des erreurs). Avant de poursuivre l’analyse multivariée, il est donc essentiel de vérifier la qualité des données et de s’assurer que les valeurs manquantes ne biaisent pas les comparaisons entre groupes.

df_wins <- df %>%
  group_by(Zone) %>%
  mutate(
    # Bornes Tukey Hauteur
    Q1_h  = quantile(Hauteur, 0.25, na.rm = TRUE),
    Q3_h  = quantile(Hauteur, 0.75, na.rm = TRUE),
    IQR_h = Q3_h - Q1_h,
    low_h = Q1_h - 1.5 * IQR_h,
    up_h  = Q3_h + 1.5 * IQR_h,
    
    # Bornes Tukey Diametre
    Q1_d  = quantile(Diametre, 0.25, na.rm = TRUE),
    Q3_d  = quantile(Diametre, 0.75, na.rm = TRUE),
    IQR_d = Q3_d - Q1_d,
    low_d = Q1_d - 1.5 * IQR_d,
    up_d  = Q3_d + 1.5 * IQR_d,
    
    # Winsorisation (cap aux bornes)
    Hauteur  = pmin(pmax(Hauteur,  low_h), up_h),
    Diametre = pmin(pmax(Diametre, low_d), up_d)
  ) %>%
  ungroup() %>%
  select(-Q1_h,-Q3_h,-IQR_h,-low_h,-up_h,
         -Q1_d,-Q3_d,-IQR_d,-low_d,-up_d)

df_wins <- df_wins %>%
  group_by(Especes) %>%
  mutate(
    # Bornes Tukey Hauteur
    Q1_h  = quantile(Hauteur, 0.25, na.rm = TRUE),
    Q3_h  = quantile(Hauteur, 0.75, na.rm = TRUE),
    IQR_h = Q3_h - Q1_h,
    low_h = Q1_h - 1.5 * IQR_h,
    up_h  = Q3_h + 1.5 * IQR_h,
    
    # Bornes Tukey Diametre
    Q1_d  = quantile(Diametre, 0.25, na.rm = TRUE),
    Q3_d  = quantile(Diametre, 0.75, na.rm = TRUE),
    IQR_d = Q3_d - Q1_d,
    low_d = Q1_d - 1.5 * IQR_d,
    up_d  = Q3_d + 1.5 * IQR_d,
    
    # Winsorisation (cap aux bornes)
    Hauteur  = pmin(pmax(Hauteur,  low_h), up_h),
    Diametre = pmin(pmax(Diametre, low_d), up_d)
  ) %>%
  ungroup() %>%
  select(-Q1_h,-Q3_h,-IQR_h,-low_h,-up_h,
         -Q1_d,-Q3_d,-IQR_d,-low_d,-up_d)


df <- df_wins

res <- df %>%
  summarise(
    n = n(),
    hauteur_min = min(Hauteur, na.rm = TRUE),
    hauteur_max = max(Hauteur, na.rm = TRUE),
    hauteur_moy = mean(Hauteur, na.rm = TRUE),
    hauteur_sd  = sd(Hauteur, na.rm = TRUE),
    diam_min    = min(Diametre, na.rm = TRUE),
    diam_max    = max(Diametre, na.rm = TRUE),
    diam_moy    = mean(Diametre, na.rm = TRUE),
    diam_sd     = sd(Diametre, na.rm = TRUE)
  )

knitr::kable(res, digits = 2, caption = "Résumé global des variables quantitatives")

tab_zone <- df %>%
  count(Zone) %>%
  mutate(pourcentage = round(100 * n / sum(n), 2)) %>%
  arrange(desc(n))

knitr::kable(tab_zone, caption = "Répartition des observations par Zone (effectifs et pourcentages)")

ggplot(tab_zone, aes(x = reorder(Zone, -n), y = n, fill = Zone)) +
  geom_col() +
  labs(x = "Zone", y = "Effectif", title = "Effectifs par Zone") +
  theme_minimal() +
  guides(fill = "none")

tab_especes <- df %>%
  count(Especes) %>%
  mutate(pourcentage = round(100 * n / sum(n), 2)) %>%
  arrange(desc(n))

knitr::kable(tab_especes, caption = "Répartition des observations par Espèces")

ggplot(tab_especes, aes(x = reorder(Especes, -n), y = n, fill = Especes)) +
  geom_col() +
  labs(x = "Espèces", y = "Effectif", title = "Effectifs par Espèces") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  guides(fill = "none")

0.1 Vérification graphique de l’indépendance des observations

Une matrice de dispersion des variables Hauteur et Diametre a été réalisée afin d’examiner visuellement la structure des données.

Le graphique montre un nuage de points continu et homogène, sans structure particulière ni regroupement artificiel.

La relation positive observée entre la hauteur et le diamètre est cohérente avec les caractéristiques biologiques des arbres et ne constitue pas une violation de l’indépendance.

Aucune structure anormale ou dépendance systématique entre les observations n’est observée.

Nous pouvons donc considérer que l’hypothèse d’indépendance des observations est raisonnablement vérifiée.

0.2 Test de normalité multivariée

Nous utilisons le test de Shapiro-Wilk multivarié.

mshapiro.test(t(Base[, c("Hauteur", "Diametre")]))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.97475, p-value = 0.1896

Interprétation du test de normalité multivariée

Le test de Shapiro-Wilk multivarié a été appliqué aux variables dépendantes Hauteur et Diametre.

Les résultats obtenus sont :

• Statistique W = 0.97475
• p-value = 0.1896

La p-value est supérieure au seuil de signification de 5 % (α = 0.05).

Nous ne rejetons donc pas l’hypothèse nulle de normalité multivariée.

Conclusion

Les variables Hauteur et Diametre suivent une distribution compatible avec la normalité multivariée.

Cette condition est nécessaire à l’application de la MANOVA et est donc vérifiée.

Nous vérifions également la normalité multivariée dans chaque zone écologique.

by(Base[, c("Hauteur", "Diametre")], Base$Zone,
   function(x) mshapiro.test(t(x)))

## Base$Zone: zone1
## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.95841, p-value = 0.4579
## 
## ------------------------------------------------------------ 
## Base$Zone: zone2
## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.8392, p-value = 0.005749
## 
## ------------------------------------------------------------ 
## Base$Zone: zone3
## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.83674, p-value = 0.0006365

Interprétation de la normalité multivariée dans chaque zone

Le test de Shapiro-Wilk multivarié a été appliqué séparément dans chaque zone écologique.

Les résultats obtenus sont les suivants :

• Zone 1 :
  • W = 0.95841
  • p-value = 0.4579
• Zone 2 :
  • W = 0.8392
  • p-value = 0.005749
• Zone 3 :
  • W = 0.83674
  • p-value = 0.0006365

Décision statistique

Au seuil de signification de 5 % :

• Pour la zone 1, la p-value est supérieure à 0.05.

  Nous ne rejetons pas l’hypothèse de normalité multivariée.

• Pour les zones 2 et 3, les p-values sont inférieures à 0.05.

  Nous rejetons l’hypothèse de normalité multivariée.

Conclusion

L’hypothèse de normalité multivariée n’est pas entièrement vérifiée dans toutes les zones écologiques.

Cela indique que la distribution conjointe de la hauteur et du diamètre n’est pas normale dans certaines zones.

Implication pour la MANOVA

La violation de l’hypothèse de normalité peut affecter la validité de la MANOVA.

Cependant, la MANOVA est relativement robuste à des violations modérées de la normalité, en particulier lorsque la taille de l’échantillon est suffisante.

Vérification de la normalité multivariée dans chaque espèce

Nous souhaiterions vérifier également la normalité multivariée dans chaque groupe défini par la variable Especes. Seulement,on ne peut pas utiliser le test de Shapiro-Wilk multivarié pour les groupes définis par la variable Especes, car certains groupes ont des effectifs très faibles (n=1), ce qui rend le test inapplicable.

0.3 Test d’homogénéité des matrices de covariance

library(biotools)

boxM(Base[, c("Hauteur", "Diametre")], Base$Zone)

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Base[, c("Hauteur", "Diametre")]
## Chi-Sq (approx.) = 17.228, df = 6, p-value = 0.008482

Interprétation du test d’homogénéité des matrices de covariance

Le test de Box permet de vérifier si les matrices de covariance des variables dépendantes sont égales entre les différentes zones écologiques.

Les résultats obtenus sont :

• Chi-Square = 17.228
• degrés de liberté = 6
• p-value = 0.008482

La p-value est inférieure au seuil de signification de 5 % (α = 0.05).

Nous rejetons donc l’hypothèse nulle d’égalité des matrices de covariance.

Conclusion

L’hypothèse d’homogénéité des matrices de covariance n’est pas vérifiée.

Cela signifie que la structure de variabilité conjointe de la hauteur et du diamètre diffère selon les zones écologiques.

0.4 Vérification de la corrélation entre les variables dépendantes

La MANOVA repose sur l’analyse simultanée de plusieurs variables dépendantes.

Pour que la MANOVA soit pertinente, ces variables doivent :

• être corrélées (sinon plusieurs ANOVA séparées seraient suffisantes),
• mais ne pas être excessivement corrélées.

En effet, si la corrélation entre les variables dépendantes est trop élevée (par exemple supérieure à 0.9 en valeur absolue), cela indique une situation de multicolinéarité.

Une multicolinéarité forte signifie que les variables mesurent pratiquement la même information, ce qui peut rendre les estimations instables et diminuer la validité de l’analyse multivariée.

Nous vérifions donc la corrélation entre Hauteur et Diametre.

cor(Base[, c("Hauteur", "Diametre")])

##            Hauteur  Diametre
## Hauteur  1.0000000 0.8662996
## Diametre 0.8662996 1.0000000

Interprétation de la corrélation entre les variables dépendantes

La matrice de corrélation indique une corrélation de :

r = 0.8663

entre les variables Hauteur et Diametre.

Cette corrélation est :

• positive,
• relativement élevée,
• mais inférieure au seuil critique de 0.9 généralement utilisé pour détecter une multicolinéarité problématique.

Cette valeur indique que les deux variables sont fortement liées, ce qui est biologiquement cohérent, puisque les arbres plus hauts tendent à avoir un diamètre plus important.

Cependant, la corrélation n’est pas suffisamment élevée pour suggérer que les deux variables mesurent exactement la même information.

Conclusion

Il n’existe pas de multicolinéarité excessive entre les variables dépendantes.

Les variables Hauteur et Diametre sont suffisamment corrélées pour justifier l’utilisation d’une MANOVA, tout en conservant une information propre à chacune.

La condition relative à l’absence de multicolinéarité problématique est donc vérifiée.

Interprétation du test MANOVA (statistique de Pillai)

Les résultats du test MANOVA basé sur la statistique de Pillai sont les suivants :

• Pillai’s Trace = 0.13827
• F approximatif = 2.3767
• degrés de liberté du numérateur = 4
• degrés de liberté du dénominateur = 128
• p-value = 0.05535

La statistique de Pillai mesure la proportion de la variance totale des variables dépendantes expliquée par le facteur Zone.

La p-value obtenue est de 0.05535.

Au seuil de signification de 5 % (α = 0.05), cette p-value est légèrement supérieure au seuil critique.

Nous ne rejetons donc pas l’hypothèse nulle d’égalité des vecteurs moyens entre les zones écologiques.

Cela signifie qu’il n’existe pas de différence statistiquement significative entre les zones écologiques en ce qui concerne la combinaison des variables Hauteur et Diametre, au seuil de 5 %.

Cependant, la proximité de la p-value avec le seuil de signification suggère l’existence d’un effet marginal de la zone écologique.

Afin d’examiner séparément l’effet de la zone écologique sur chacune des variables, nous réalisons des analyses de variance univariées complémentaires. ## 1.5 Analyses univariées complémentaires

Les analyses de variance univariées permettent d’examiner séparément l’effet de la variable Zone sur chacune des variables dépendantes.

Ces analyses complètent la MANOVA en identifiant les variables qui contribuent aux différences observées entre les groupes.

summary.aov(modele_manova_zone)

##  Response Hauteur :
##             Df  Sum Sq Mean Sq F value Pr(>F)
## Zone         2   5.972  2.9861  0.7311 0.4854
## Residuals   64 261.413  4.0846               
## 
##  Response Diametre :
##             Df Sum Sq Mean Sq F value  Pr(>F)  
## Zone         2  29.98 14.9904  2.5418 0.08665 .
## Residuals   64 377.45  5.8976                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interprétation des analyses univariées

Les analyses de variance univariées permettent d’examiner séparément l’effet de la variable Zone sur chacune des variables dépendantes.

Synthèse des résultats univariés

Les analyses univariées confirment les résultats de la MANOVA : la zone écologique n’a pas d’effet statistiquement significatif sur la hauteur ni sur le diamètre des arbres.

Interprétation du test MANOVA (statistique de Pillai)

Le test MANOVA basé sur la statistique de Pillai permet d’évaluer l’effet des facteurs Zone et Especes sur la combinaison des variables dépendantes Hauteur et Diametre.

Les résultats obtenus sont les suivants :

Effet de la zone écologique (Zone)

• Pillai’s Trace = 0.26578
• F approximatif = 3.9847
• degrés de liberté du numérateur = 4
• degrés de liberté du dénominateur = 104
• p-value = 0.004765

La p-value est inférieure au seuil de signification de 5 % (α = 0.05).

Nous rejetons donc l’hypothèse nulle d’égalité des vecteurs moyens entre les zones écologiques.

Cela signifie que la zone écologique a un effet statistiquement significatif sur la combinaison des variables Hauteur et Diametre.

Effet du type d’espèce (Especes)

• Pillai’s Trace = 0.99274
• F approximatif = 4.2709
• degrés de liberté du numérateur = 24
• degrés de liberté du dénominateur = 104
• p-value = 1.129 × 10⁻⁷

La p-value est extrêmement inférieure au seuil de signification de 5 %.

Nous rejetons donc fortement l’hypothèse nulle d’égalité des vecteurs moyens entre les espèces.

Cela indique que le type d’espèce a un effet hautement significatif sur les dimensions des arbres.

Conclusion du test MANOVA

Les résultats montrent que les facteurs Zone et Especes ont tous deux un effet statistiquement significatif sur la combinaison des variables Hauteur et Diametre.

Le facteur Especes présente un effet particulièrement important, comme l’indique la valeur élevée de la statistique de Pillai.

Afin d’identifier précisément quelles variables sont affectées par ces facteurs, nous réalisons des analyses de variance univariées complémentaires.

Interprétation des analyses univariées

Les analyses de variance univariées permettent d’examiner séparément l’effet des facteurs Zone et Especes sur chacune des variables dépendantes.

Effet sur la hauteur des arbres (Hauteur)

Les résultats obtenus sont les suivants :

• Effet de la Zone :
  • F = 2.0732
  • p-value = 0.1361

La p-value est supérieure au seuil de signification de 5 %.

Nous ne rejetons donc pas l’hypothèse nulle.

La zone écologique n’a pas d’effet statistiquement significatif sur la hauteur des arbres.

• Effet des Especes :
  • F = 19.8189
  • p-value = 2.625 × 10⁻¹⁵

La p-value est extrêmement inférieure au seuil de signification de 5 %.

Nous rejetons donc l’hypothèse nulle.

Le type d’espèce a un effet hautement significatif sur la hauteur des arbres.

Effet sur le diamètre des arbres (Diametre)

Les résultats obtenus sont les suivants :

• Effet de la Zone :
  • F = 6.3732
  • p-value = 0.003346

La p-value est inférieure au seuil de signification de 5 %.

Nous rejetons donc l’hypothèse nulle.

La zone écologique a un effet statistiquement significatif sur le diamètre des arbres.

• Effet des Especes :
  • F = 20.6280
  • p-value = 1.153 × 10⁻¹⁵

La p-value est extrêmement inférieure au seuil de signification de 5 %.

Nous rejetons donc l’hypothèse nulle.

Le type d’espèce a un effet hautement significatif sur le diamètre des arbres.

Synthèse des analyses univariées

Les résultats montrent que :

• le facteur Especes a un effet hautement significatif sur la hauteur et le diamètre,
• le facteur Zone a un effet significatif sur le diamètre, mais pas sur la hauteur.

Interprétation du modèle MANOVA avec interaction

Le modèle MANOVA avec interaction a été estimé afin d’évaluer les effets des facteurs Zone, Especes, et de leur interaction sur les variables dépendantes Hauteur et Diametre.

Les degrés de liberté sont répartis comme suit :

• Zone : 2 degrés de liberté,
• Especes : 12 degrés de liberté,
• Interaction Zone × Especes : 5 degrés de liberté,
• Résidus : 47 degrés de liberté.

Les erreurs standards résiduelles sont de :

• 1.031 pour la variable Hauteur,
• 1.275 pour la variable Diametre.

Ces valeurs sont similaires à celles obtenues dans le modèle sans interaction, ce qui indique que l’ajout du terme d’interaction n’améliore pas fortement l’explication de la variabilité.

Le message suivant apparaît :

“19 out of 39 effects not estimable”

Cela signifie que certains effets d’interaction ne peuvent pas être estimés.

Cette situation se produit lorsque certaines combinaisons Zone × Especes ne sont pas présentes dans les données ou ne possèdent pas suffisamment d’observations.

Cela correspond à un plan expérimental déséquilibré ou incomplet.

Dans ce cas, le modèle peut toujours être estimé, mais l’interprétation de certains effets d’interaction doit être faite avec prudence.

Nous procédons maintenant au test MANOVA afin d’évaluer la significativité statistique des effets de Zone, Especes, et de leur interaction.

Interprétation du test MANOVA avec interaction (statistique de Pillai)

Le test MANOVA basé sur la statistique de Pillai permet d’évaluer l’effet des facteurs Zone, Especes, et de leur interaction sur la combinaison des variables Hauteur et Diametre.

Les résultats obtenus sont les suivants :

Effet principal de la zone écologique (Zone)

• Pillai’s Trace = 0.27860
• F approximatif = 3.8033
• degrés de liberté du numérateur = 4
• degrés de liberté du dénominateur = 94
• p-value = 0.006544

La p-value est inférieure au seuil de signification de 5 %.

Nous rejetons donc l’hypothèse nulle.

La zone écologique a un effet statistiquement significatif sur la combinaison des variables Hauteur et Diametre.

Effet principal du type d’espèce (Especes)

• Pillai’s Trace = 1.00598
• F approximatif = 3.9638
• degrés de liberté du numérateur = 24
• degrés de liberté du dénominateur = 94
• p-value = 8.156 × 10⁻⁷

La p-value est extrêmement inférieure au seuil de signification de 5 %.

Nous rejetons donc fortement l’hypothèse nulle.

Le type d’espèce a un effet hautement significatif sur les dimensions des arbres.

Effet d’interaction Zone × Especes

• Pillai’s Trace = 0.11312
• F approximatif = 0.5635
• degrés de liberté du numérateur = 10
• degrés de liberté du dénominateur = 94
• p-value = 0.839635

La p-value est largement supérieure au seuil de signification de 5 %.

Nous ne rejetons donc pas l’hypothèse nulle.

Cela indique qu’il n’existe pas d’interaction statistiquement significative entre la zone écologique et le type d’espèce.

Autrement dit, l’effet du type d’espèce sur les dimensions des arbres est similaire dans toutes les zones écologiques.

a. Projecteurs orthogonaux

Pour toute matrice \(X\) de plein rang, \(P_X = X(X^{\prime}X)^{-1}X^{\prime}\) est le projecteur orthogonal sur l’espace des colonnes de \(X\). Il vérifie \(P^2_X = P_X\), \(P^{\prime}_X = P_X\) et \(I_n-P_X\) est le projecteur sur l’orthogonal.

b. Produit de Kronecker

Pour des matrices \(A (m\times n)\) et \(B (p\times q)\), le produit de Kronecker \(A\otimes B\) est la matrice \(mp \times nq\) définie par blocs \((A \otimes B)_{ik,jl} = a_{ij} b_{kl}\).

Propriétés : \[ (A \otimes B)(C \otimes D) = AC \otimes BD, vec(ABC) = (C^{\prime} \otimes A) vec(B) \]

a. Loi de Wishart

Soit \(Z\) une matrice \(m \times D\) dont les lignes sont i.i.d. \(\mathcal{N}_D(0, \Sigma)\). Alors \(W = Z^{\prime}Z\) suit une loi de Wishart à \(m\) degrés de liberté et matrice d’échelle \(\Sigma\), notée \(\mathcal{W}_D(m, \Sigma)\).

b. Théorème de Cochran multivarié

Soit \(Y\) une matrice \(n \times D\) dont les lignes sont i.i.d. \(\mathcal{N}_D(0,\Sigma)\). Soient \(Q_1, ..., Q_k\) des matrices symétriques idempotentes de rangs \(r_i\) telles que \(\sum_{i=1}^{k}Q_i=I_n\) et \(Q_iQ_j=0\) pour \(i \neq j\). Alors les formes quadratiques \(Y^{\prime}Q_iY\) sont indépendantes et suivent des lois de Wishart \(\mathcal{W}_D(r_i,\Sigma\)).