Application pratique d'une Analyse de Variance

Avant Propos

Cette page WEB a pour but de donner quelques exemples concrêts de l'utilisation de l'Analyse de Variance (ANOVA) par le biais d'un logiciel libre de statistiques : R. Ce travail provient d'une étude réalisée conjointement par A.B. Dufour, D. Pontier et D. Chessel. Le modèle biologique choisit est le chat (Felis Catus), et la question principale est de savoir quel est le paramètre biologique qui influence la taille des testicules des chats domestiques.

I - La problématique biologique

Il est avant tout nécessaire de définir la notion de compétition spermatique hypothèse sousjacente à la théorie de la sélection sexuelle. D'après cette théorie, la taille des testicules des individus est d'autant plus importante qu'ils sont dans un milieu compétitif pour la reproduction. Ceci s'explique par le fait que l'augmentation de ces testicules permet la production d'une plus grande quantité de sperme ce qui multiplie naturellement les chances de fécondation d'une femelle lors de la reproduction. En effet les femelles des chats peuvent être fécondées par plusieurs mâles et donc avoir des chatons issus de mâles différents au sein d'une même portée.

L'Hypothèse de travail est donc qu'il existe un lien entre compétition spermatique et la taille des testicules.

A) Le contenu du fichier de données

Les variables quantitatives : Poids et Volume

Pour chaque individu capturé, deux variables morphométriques ont été relevées : le poids des testicules (en g) et leur volume (en cm3), ces deux variables sont considérées aléatoires. Leurs vecteurs respectifs sont poi et vol

Les variables qualitatives : Les populations de chats

Le fichier de données à proprement parlé

Voici ci-dessous le contenu exhaustif du fichier que nous allons traiter, ils vous est tout à fait possible de le sauvegarder afin de reproduire les commandes qui seront utilisées ici.

Si vous souhaitez enregistrer le fichier, faites un clic droit puis Enregistrer la cible du lien sous ici.

Haut de la page

II - Méthodes statistiques

A) Le cas des variables quantitatives

Etant donné que l'on suppose que chaque individu a été tiré de façon aléatoire, le poids et le volume sont considérés aléatoires, ce qui implique que le modèle est lui aussi aléatoire. Dans ce cas précis, on peut parler de mesure de corrélation entre le poids et le volume des testicules. En effet si le modèle avait été mixte (l'expérimentateur ayant fixé une des deux variables) le coefficient de détermination calculé entre le poids et le volume serait une mesure du pourcentage de variabilité de la variable aléatoire expliquée par le modèle, et non une mesure de corrélation entre ces deux variables puisque la distribution de la variable fixée influencerait ce coefficient.

B) Les facteurs et leur contrastes

Une variable dite qualitative, par exemple ici le vecteur pop et ses 5 modalités, se traduit par une matrice d'indicatrice de classe. Ceci permet de discriminer les différentes modalités et de fixer une référence c'est à dire un témoin. Il devient alors possible et même pertinent de créer des hypothèses précises et de les tester via cette matrice de contraintes.

C) Modèle linéaire et estimation de paramètres

Une manière géométrique de représenter les données et de se placer dans l'espace des variables, où chaque axe est un individu. Sur le schéma ci-dessous, on a par exemple 3 individus, Y est un vecteur observé, celui-ci se projette orthogonalement sur le vecteur Delta (1ère diagonale de constantes de l'espace) en un point qui représente l'estimation ponctuelle de la moyenne de la population. Il est possible de généraliser ce concept pour déterminer les variances et donc d'estimer les paramètres du modèle, à savoir la pente et l'ordonnée à l'origine de la droite de régression linéaire.

Explication Géométrique

Dans la pratique voici les formules que l'on utilise pour calculer ces paramètres :

D) Calculs des paramètres d'une régression linéaire (Javascript):

Taper les valeurs de votre séquence explicatrice ici en les séparant par des virgules (comme dans l'exemple):

X :

Taper les valeurs de votre séquence à expliquer ici en les séparant par des virgules (comme dans l'exemple):

Y :

Intercept :     Pente :     Coeff. de corrélation :

E) Les modèles linéaires

Nous allons simplement évoquer quelques modèles de base de la régression linéaire.
  1. Le plus simple est le modèle constant que l'on écrit pour un chat i:

    Modèle Constant

    Ce modèle suppose que tous les individus ont des tailles de testicules identiques, c'est à dire de volume égal à la moyenne du vecteur volume. Ce modèle est bien évidemment beaucoup trop simple pour expliquer les données, nous devons le complexifier.

  2. Dans le second modèle, le poids influence directement le volume des testicules des chats, mais la provenance des individus (vecteur pop) n'est ici pas prise en compte, la pente et l'ordonnée à l'origine estimées sont identiques quelque soit l'individu.

    Modèle poids

  3. Le troisième modèle prend en compte à la fois le poids et l'appartenance à une population pour expliquer le volume des testicules. L'intercept varie pour chaque population, mais la pente est commune à chaque population. Ce modèle est additif, il ne suppose pas d'intéraction entre le poids la provenance géographique.

    Modèle poids et population

  4. Le quatrième modèle est identique au précédent mais il prend en compte l'intéraction entre le poids et la population. L'intercept et la pente varient pour chaque population, on peut le considérer comme le modèle parapluie.

    Modèle poids et population avec intéraction

Il est important de noter que le modèle linéaire suppose la normalité des résidus du modèle, ainsi que leur homoscedasticité.

Modèle poids et population avec intéraction

F) L'Analyse de Variance comme comparaison de modèles

Une analyse de variance à pour but de déterminer si la moyenne d'une variable d'intérêt, ici le volume des testicules, varie significativement entre k sous-groupes (les différentes populations), ce qui traduit un effet population. Voici comment se construit une analyse de variance :

Somme des carrés des écarts

Somme des carrés des écarts calculés

Pour tester par exemple l'effet de d'un facteur contrôlé (pop) à k modalités, on construit le tableau d'analyse de variance à un facteur suivant:

Tableau d'ANOVA

Ce tableau se généralise facilement à plusieurs facteurs, et son principe consiste à déterminer si une différence de moyenne d'une variable d'intérêt entre plusieurs groupes est significative ou non, en d'autre terme : y a-t-il un effet du facteur? Sous l'hypothèse nulle (absence d'effet du facteur) les carrés moyens estiment la variance de la population et leur rapport qui vaudrait 1 suit une loi de Fisher (loi dont la fonction de densité de probabilité est unilatérale). Sous l'hypothèse alternative, le carré moyen résiduel estime toujours la variance résiduelle.

Fonction de densité de probabilité d'une loi de Fisher

Si la statistique observée dépasse la statistique attendue sous l'hypothèse nulle (d'absence d'effet du facteur sur les moyennes) alors on rejette cette hypothèse nulle au profit de l'hypothèse complémentaire avec un risque de première espèce qui est exactement la pvalue.

Probabilité critique

Par exemple si l'on veut tester l'effet additif de la population (pop) sur le volume des testicules des chats en présence du facteur poids, on confronte en réalité deux modèles : le modèle poids et le modèle poids+population:

Modèle poids

Modèle poids+population

On teste le coefficient estimé b2 du second modèle à 0, si la probabilité critique est forte (>5%), alors ceci implique que l'on n'est pas dans les conditions qui nous permettent de rejetter l'hypothèse nulle, c'est à dire que b2 n'est pas significativement différent de 0, ce qui signifie que dans ce cas l'effet population n'est pas significatif, on choisit donc le modèle poids :

Modèle poids

Haut de la page

III - Applications sous R

A) Récupération des données et construction du data-frame

En premier lieu on va charger le fichier de données en mémoire, pour cela exécutez l'ordre :

data=read.table("donnees.txt",header=TRUE)
head(data)
summary(data)

Vous avez ainsi une idée assez un peu plus claire de comment est consruit le fichier, notamment grâce à la méthode summary qui donne la moyenne, la médiane, le minimum, le maximum de chaque vecteur quantitatif, et qui pour un vecteur qualitatif (ici pop) compte le nombre de répétition de chaque modalités du facteur.

Nous pouvons clairement émettre l'hypothèse que le volume des testicules augmente avec le poids des chats, avant de tester cette hypothèse, il conviendra de linéariser les variables afin de les rendre plus propices à l'utilisation de tests paramétriques comme l'ANOVA par exemple.

lnpoi=log(data$poi)
lnvol=log(data$vol)
pop=data$pop
new_data=cbind.data.frame(lnpoi,lnvol,pop)

B) Influence du poids sur le volume des testicules

On a reconstruit un data.frame avec les variables quantitatives linéarisées, on va donc tester la corrélation qu'il existe entre le poids et le volume (linéarisés).

cor.test(lnpoi,lnvol,data=new_data)

Le coefficient de corrélation est estimé à 0.467, la probabilité critique étant très significative (pv=1.718e-09) on rejette l'hypotèse nulle d'absence de corrélation entre ces deux vecteurs, ce qui implique que le volume est positivement corrélé au poids. Les gros chats ont de plus gros testicules que les chats plus petits. Notez que l'on peut retrouver cette estimation du coefficient de corrélation via la formule de son estimation ponctuelle :

cov(lnpoi,lnvol)/sqrt(var(lnpoi)*var(lnvol))

On peut réaliser la régression linéaire du modèle volume en fonction du poids et tracer la représentation graphique correspondante. On estime par ailleurs la pente et l'intercept du modèle et on trace la droite de régression linéaire sur le nuage de point.

lm1=lm(lnvol~lnpoi,data=new_data)
coefficients(lm1)
pente=cov(lnpoi,lnvol)/var(lnpoi)
intercept=mean(lnvol)-pente*mean(lnpoi)
plot(new_data$lnpoi,new_data$lnvol,main="log(volume) en fonction de log(poids)",xlab="log(poids)",ylab="log(volume)")
abline(lm1,col="blue")
abline(intercept,pente,col="red",lty=2,lwd=2)

Régression Linéaire

Les formules permettent là encore de réaliser des estimations ponctuelles des paramètres du modèle, on vérifie que l'on retrouve les mêmes estimations que le logiciel R, en effet les deux droites de régression linéaire sont parfaitement confondues.

anova(lm1)

L'utilisation de la méthode anova permet ici de réaliser la régression linéaire du modèle, l'effet du poids sur le volume est très significatif, ce qui est exactement la même conclusion que ce qui précède.

C) Effet population et vérification des pré-supposées

Nous allons maintenant tester s'il existe un effet population sur le volume des testicules, en d'autres termes, on va vérifier si la moyenne du volume des testicules est la même pour toute les population, pour celà on construit le modèle linéaire approprié et on teste l'effet population:

lm2=lm(lnvol~pop,data=new_data)
anova(lm2)

L'effet population est très significatif, on a donc clairement une différence moyenne du volume des testicules de chats suivant leur situation géographique, toutefois avant d'aller plus loin, nous devons vérifier que les conditions pré-requises à l'utilisation d'une analyse de variance sont réunies, il nous faut vérifier l'indépendance des résidus (absence d'autocorrélation), leur homoscedasticité entre chaque population et enfin nous devons vérifier qu'ils sont distribués normalement. Notez qu'un test doit toujours être accompagné d'un graphique, bien souvent l'aspec visuel renseigne sur les tendances générales et permet de mieux cibler quelle analyse on va mettre en place. On peut représenter les résidus d'un modèle de nombreuses manières, en voici quelques unes:

par(mfrow=c(2,2))
plot(lm2,ask=F)

Plot(lm)

Cette série de graphique permet de visuellement apprécier que ces conditions sont réunies, en haut à gauche : les résidus du modèle sont confrontés aux valeurs ajustées, cela permet de savoir s'il existe une organisation des résidus. En moyenne ils valent 0, et on ne note pas d'organisation particulière du nuage de points. Les résidus ne sont pas auto-corrélés. En haut à droite : la normalité des résidus semble avérée (graphe quantile/quantile). En bas à gauche : permet de voir si la variance des résidus est constante (homoscedasticité), on ne note pas de modification particulière de la morphologie du nuage de point, l'homoscedasticité des résidus est admise. En bas à droite : le graphe de Cook détermine le "poids" de chaque point dans la régression linéaire.

On trace ensuite un graphe quantile-quantile (droite de Henry) qui confronte les résidus empiriques du modèle aux résidus théoriques de la loi normale, la normalité est difficilement contestable. De même on trace l'histogramme des résidus du modèle, la courbe représente la fontion de densité de probabilité de la distribution, l'ajustement est plutôt bon. Enfin le dernier graphique représente des boîtes à moustache (boxplot) horizontales, des résidus pour chaque population, là encore l'homoscedasticité semble avérée.

par(mfrow=c(1,2))
qqnorm(lm2$residuals)
qqline(lm2$residuals)
hist(lm2$residuals,nclass=15,proba=TRUE,col="lightblue",xlab="Résidus",main="Histogramme des résidus")
x1 <- seq(-1.5, 1, length = 100)
lines(x1,dnorm(x1,mean(lm2$residuals),sd(lm2$residuals)),col="red")

Droite de Henry et Histogramme des résidus

library(lattice)
bwplot(pop ~ residuals(lm2), data=new_data, xlab="Résidus")

BwPlot

Pour rassurer les plus sceptiques, on réalise un test de normalité sur les résidus (Shapiro) ainsi qu'un test d'égalité des variances (Bartlett), les conclusions sont les mêmes.

shapiro.test(lm2$residuals)
bartlett.test(lm2$residuals,pop,data=new_data)

D) Détermination du meilleur modèle

Nous allons maintenant tenter de déterminer quel est le modèle le plus adapté à ces données en intégrant chaque facteur dans un modèle linéaire. Par ailleurs nous verrons que réaliser une analyse de variance sur le modèle le plus complexe (modèle parapluie) revient exactement à réaliser un test de linéarité entre le modèle dit additif (poi+pop)et le modèle interactif (poi*pop).

lm4=lm(lnvol~lnpoi*pop,data=new_data)
anova(lm4)

L'interaction entre le facteur poids (linéarisé) et le facteur population n'est pas significative, ce modèle est donc inutilement sur-paramétré par rapport au modèle additif qui est donc plus simple et tout aussi informatif. On vérifie que l'on retrouve la même statistique F et donc la même probabilité critique avec un test de linéarité:

lm3=lm(lnvol~lnpoi+pop,data=new_data)
anova(lm3,lm4)

On peut très facilement retrouver la statistique F et donc la probabilité critique associée au test de linéarité par un calcul algébrique. Ici Dm1 et Dm2 sont les déviances des 2 modèles, df correspond au calcul p-q (différence de degrés de liberté de chacun des modèles), on applique alors la formule suivante:

Test de linéarité

D1=sum(residuals(lm4)^2)
D2=sum(residuals(lm3)^2)
df=summary(lm4)$fstatistic[[2]]-summary(lm3)$fstatistic[[2]]
F=((D2-D1)/df)/(D1/lm4$df.residual)
1-pf(F,df,lm4$df.residual)

E) Graphique du modèle sélectionné

On va réaliser le graphique correspondant au modèle additif, chaque droite aura donc la même pente, mais l'intercept va différer selon la population; on définit ces deux paramètres de régression multiple grâce au modèle linéaire additif lm3. Attention, ce script n'a pas la prétention d'être optimisé, il a pour but principal de faciliter la compréhension de la consruction d'un graphique et de la manipulation des estimations ponctuelles des paramètres du modèle.

coeff=coefficients(lm3)
plot(lnpoi[pop=="BAC"],lnvol[pop=="BAC"],xlim=c(range(lnpoi)[1],range(lnpoi)[2]),ylim=c(range(lnvol)[1],
range(lnvol)[2]),col="blue",main="Modèle Additif",xlab="log(poids)",ylab="log(volume)")
points(lnpoi[pop=="CxR"], lnvol[pop=="CxR"],col="red")
points(lnpoi[pop=="KER"], lnvol[pop=="KER"],col="green")
points(lnpoi[pop=="ROM"], lnvol[pop=="ROM"],col="black")
points(lnpoi[pop=="SJT"], lnvol[pop=="SJT"],col="orange")
abline(coeff[[1]],coeff[[2]],col="blue",lwd=2)
abline(coeff[[1]]+coeff[[3]],coeff[[2]],col="red",lwd=2)
abline(coeff[[1]]+coeff[[4]],coeff[[2]],col="green",lwd=2)
abline(coeff[[1]]+coeff[[5]],coeff[[2]],col="black",lwd=2)
abline(coeff[[1]]+coeff[[6]],coeff[[2]],col="orange",lty=2,lwd=2)
legend(7.6,8,col=c("blue","red","green","black","orange"),pch=1, legend=levels(pop),lty=c(rep(1,4),2))

Modèle additif

F) Utilisation d'une matrice de contrastes

Nous avons vu que globalement la population dont sont originaires les chats influence le volume de leurs testicules, mais cela ne nous permet de pas de répondre implicitement à notre question biologique. En effet comment la provenance géographique influence-t-elle vraiment ce volume? Par exemple il serait intéressant de comparer le volume des testicules de chats de deux populations à système d'appariement polygame, ou de comparer l'effet insulaire par rapport au continent ou même de voir si des regroupements sont possibles, bref on souhaite avoir plus de flexibilité dans la manière de tester nos hypothèses. L'ANOVA ne permet pas cela puisqu'elle teste l'effet global d'un facteur, en revanche un test sur chaque coefficients résoud ce problème et permet de tester absolument toutes les combinaisons de regroupement et de confrontation de modalités d'un même facteur qualitatif, c'est la méthode des contrastes. On va en premier lieu prendre connaissance de la construction de la matrice des contrastes du facteur population par le logiciel.

head(model.matrix(~new_data$pop))
attributes(model.matrix(~new_data$pop))$assign

Etant donné que la modalité "BAC" est la première de l'ordre alphabétique, elle sera choisie comme témoin, c'est à dire que les tests sur les coefficients vont comparer chaque population à Barisey-la-Côte. Comme ce n'est pas ce que l'on cherche à tester, on va devoir redéfinir la matrice des contrastes. Dans ce type de plan expérimental purement prospectif, il n'existe tout simplement pas de témoin, on va alors construire un vecteur pour chaque hypothèse biologique à tester:

  1. Le volume des testicules est-il plus élevé en milieu urbain?
    Le premier contraste va comparer les population urbaines aux autres (CxR,ROM VS KER,SJT,BAC)
  2. La monogamie de la population insulaire diffère-t-elle de la polyginie des populations rurales?
    Le second contraste opposera donc BAC,SJT VS KER .
  3. Les populations rurales diffèrent-elles?
    Ce contraste confrontera CxR VS ROM .
  4. Les populations urbaines diffèrent-elles significativement l'une de l'autre?
    Ce dernier contraste sera SJT VS BAC

On obtient alors la matrice de contraste suivante, chacune des hypothèses se 'lit' en colonne:

Matrice de contrastes

contrastes <- matrix(c(-1, 1, -1, 1, -1, 1, 0, -1, 0, 1, -1, 0, 0, 0, 1, 0, 1, 0, -1, 0), ncol = 4)
dimnames(contrastes) <- list(levels(pop), c("Urb", "Mono_Poly", "WRur", "WUrb"))
contrasts(pop) <- contrastes
summary(lm(lnvol ~ lnpoi + pop))

Les tests sur les coefficients permettent donc de conclure facilement sur chacune des hypothèses:

  1. Il existe une différence significative du volume des testicules des chats entre les populations de faible et de forte densité
  2. Le volume est différent entre les deux systèmes d'appariemment (polygames VS polygines).
  3. Le volume n'est pas significativement différent entre les deux populations rurales.
  4. Les chats ont des volumes de testicules significativement différents entre les deux populations urbaines.

Haut de la page

IV - Interprétation biologique

Pour voir de quelle manière se différencient le poids et le volume des testicules à travers les populations, on réalise deux boxplot:

par(mfrow=c(1,2))
boxplot(lnpoi~pop,data=new_data,xlab="populations",ylab="log(poids)")
boxplot(lnvol~pop,data=new_data,xlab="populations",ylab="log(volume)")

Boxplot

Le poids des testicules varie, comme on pouvait s'y attendre, très peu entre les populations. Par contre ce qui frappe, c'est que ce sont les chats de Barisey-la-Côte et de Saint-Just-Chaleyssin (les deux populations rurales) qui semblent avoir le plus gros volumes de testicules, ce qui réfute totalement notre hypothèse biologique, en effet les chats urbains sencés avoir la compétition spermatique la plus forte sont loins derrière en terme de volume de testicules. Il semblerait qu'on ait omis une réalité biologique majeure : les testicules ne produisent pas seulement des spermatozoïdes, ils sont aussi à l'origine de la sécretion d'hormone mâles (testostérone) qui se libère notamment lors de gros efforts physiques comme... les combats par exemple! En effet, les chats vivant dans un milieu rural ayant moins d'opportunité de reproduction que les chats urbains, ceux-ci monopolisent les femelles et cela passe nécessairement par des affrontements entre mâles, d'où la nécessité de pouvoir produire suffisamment d'hormones mâles pour entretenir cette agressivité, ce qui implique donc des volumes de testicules supérieurs à la moyenne. Toutefois l'histoire ne dit pas si ces chats ruraux sont plus sujets aux crises cardiaques que leurs congénères des villes!

Haut de la page

Références et contacts

Vous trouverez la fiche originale au format pdf ici. Pour toute question concernant cette page HTML, n'hésitez pas me contacter, pour des questions relatives à la fiche originale, vous pouvez de même contacter les auteurs.

Haut de la page

Donnez votre avis :

Nom      
Prénom

Quelles améliorations suggérez-vous pour cette page?

Utilisez-vous souvent la statistique? Oui   Non  
Que pensez-vous de cette page? Très Bien   Bien   Moyen   Mauvais