Mon premier article : Un design croisé non équilibré

Les analyses statistiques que j'ai réalisé pour mon premier article publié

Marie Vaugoyeau

8 minutes

En 2015, j’ai publié l’article : Is oxidative status influenced by dietary carotenoid and physical activity after moult in the great tit (Parus major)? sur ma première expérimentation.
Dans le but de tester les relations entre caroténoïdes alimentaires, activité physique et status oxidatif, j’ai capturé et gardé en volière des mésanges charbonnières.
Initialement, nous avions prévu de faire une expérience avec un design croisé équilibré : avec ou sans caroténoïdes dans l’alimentation et avec ou sans rémiges primaires enlevées (retirer des rémiges augmente l’activité physique des oiseaux).

Table 1: Design expérimental prévu
Caroténoïdes dans l’eau Eau pure
Activité physique augmentée 7F 7M 7F 7M
Activité physique normale 7F 7M 7F 7M

Avec F = femelles et M = mâles

A cause du temps, le nombre d’oiseaux capturés a été largement inférieur à celui prévu, nous n’avons donc fait que trois traitements sur les quatre car je n’ai capturé que 42 adultes (un est mort avant le début de l’expérience).

Table 2: Design expérimental réel
Caroténoïdes dans l’eau Eau pure
Activité physique augmentée 6F 8M 5F 8M
Activité physique normale 7F 7M

Ci-dessous, je détaille les analyses statistiques menées à l’époque.

Première étape : Vérifier la répétabilité des données

Pour certaines variables (comptage d’hétérophiles sur lymphocytes, couleur des plumes…), j’ai fais deux ou trois mesures. Afin de quantifier la répétabilité, j’ai utilisé un modèle linéaire qui m’a donnée la statistique de test F avec le degré de liberté et la p-value.

mod<-lm(Area~Identity, data=X)
anova(mod)
## Analysis of Variance Table
## 
## Response: Area
##            Df  Sum Sq Mean Sq F value   Pr(>F)   
## Identity    1  11.558 11.5584  10.346 0.001655 **
## Residuals 124 138.528  1.1172                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dans cet article j’ai aussi calculé manuellement le facteur de corrélation intraclass selon la méthode de Lessells and Boag 1987: r = (MSA - MSw)/(MSA + MSW(n-1)).

Une fois la répétabilité calculée, j’ai calculé et utilisé les moyennes pour le reste des analyses.

Deuxième étape : Modèles linéaires mixtes avec backward élimination

Pour tester les effets de l’activité physique et de la supplémentation en caroténoïdes sur le système antioxydant, la réponse au stress et la coloration du plumage, tous les modèles ont comme facteurs fixes la longueur du tarse, la masse et deux triples interactions entre la supplémentation en caroténoïdes ou l’activité physique avec le sexe et la période (au début ou à la fin de l’expérience).

Modèle linéaire mixte

Pour chaque variable, j’ai écrit une modèle linéaire mixte comme celui-ci :

library(nlme)
m1<-lme(Oxy ~ Tarsus+Mass+Exp*Sex*Carot+Exp*Sex*Plucking, random=list(Aviary=~1, Identity=~1), method="ML", data = Y, na.action=na.exclude)

Pour prendre en compte la variance biologique, l’identité de la volière et de l’individu étaient mis en facteurs aléatoires.

Distribution des résidus du modèle

Pour valider à posteriori l’utilisation d’un modèle linéaire, j’ai testé la distribution des résidus avec un test de Shapiro et un Q-Q plot.

shapiro.test(residuals(m1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m1)
## W = 0.98678, p-value = 0.5759
qqnorm(residuals(m1));qqline(residuals(m1))

Si les résidus ne suivaient pas une loi normale, la variable était log-transformée puis un modèle linéaire mixte était refait.

Sélection de modèle avec élimination

Après le premier modèle, j’en ai réalisé un autre en retirant le facteur fixe le moins significatif, puis j’ai utilisé le critère d’information d’Akaike pour choisir le meilleur modèle.

anova(m1)
##                  numDF denDF  F-value p-value
## (Intercept)          1    33 3894.578  <.0001
## Tarsus               1    25    9.314  0.0053
## Mass                 1    33    0.414  0.5243
## Exp                  1    33   37.679  <.0001
## Sex                  1    25    5.208  0.0313
## Carot                1     9    0.519  0.4897
## Plucking             1     9    1.001  0.3432
## Exp:Sex              1    33    0.888  0.3529
## Exp:Carot            1    33    2.560  0.1191
## Sex:Carot            1    25    0.017  0.8986
## Exp:Plucking         1    33    1.918  0.1753
## Sex:Plucking         1    25    2.002  0.1694
## Exp:Sex:Carot        1    33    0.113  0.7391
## Exp:Sex:Plucking     1    33    0.782  0.3831
summary(m1)
## Linear mixed-effects model fit by maximum likelihood
##  Data: Y 
##        AIC      BIC    logLik
##   772.4558 813.1614 -369.2279
## 
## Random effects:
##  Formula: ~1 | Aviary
##         (Intercept)
## StdDev:     7.01982
## 
##  Formula: ~1 | Identity %in% Aviary
##         (Intercept) Residual
## StdDev:    7.071078 21.21047
## 
## Fixed effects: Oxy ~ Tarsus + Mass + Exp * Sex * Carot + Exp * Sex * Plucking 
##                               Value Std.Error DF   t-value p-value
## (Intercept)               289.61059 143.81205 33  2.013813  0.0522
## Tarsus                     -3.63000   6.48630 25 -0.559640  0.5807
## Mass                       -0.09082   4.43841 33 -0.020462  0.9838
## ExpStart                   44.13499  13.03342 33  3.386293  0.0018
## SexM                      -14.94850  15.46285 25 -0.966736  0.3429
## CarotYes                   -1.46163  16.34930  9 -0.089400  0.9307
## PluckingYes                23.53261  15.67578  9  1.501208  0.1675
## ExpStart:SexM               5.54477  18.34549 33  0.302241  0.7644
## ExpStart:CarotYes         -11.31881  19.99812 33 -0.565994  0.5752
## SexM:CarotYes              10.80621  20.00440 25  0.540192  0.5938
## ExpStart:PluckingYes       -5.26320  19.44272 33 -0.270703  0.7883
## SexM:PluckingYes           -9.06389  20.22299 25 -0.448197  0.6579
## ExpStart:SexM:CarotYes      4.51084  25.90055 33  0.174160  0.8628
## ExpStart:SexM:PluckingYes -23.19101  26.23251 33 -0.884056  0.3831
##  Correlation: 
##                           (Intr) Tarsus Mass   ExpStr SexM   CartYs PlcknY
## Tarsus                    -0.873                                          
## Mass                      -0.247 -0.251                                   
## ExpStart                   0.031  0.073 -0.292                            
## SexM                       0.283 -0.133 -0.381  0.497                     
## CarotYes                   0.042 -0.144  0.205 -0.060 -0.056              
## PluckingYes               -0.168  0.147 -0.050  0.395  0.353 -0.558       
## ExpStart:SexM             -0.005 -0.056  0.182 -0.703 -0.662  0.038 -0.279
## ExpStart:CarotYes          0.013  0.013 -0.052  0.015  0.020 -0.621  0.350
## SexM:CarotYes              0.064  0.047 -0.222  0.065  0.087 -0.732  0.402
## ExpStart:PluckingYes      -0.001 -0.029  0.116 -0.647 -0.303  0.366 -0.618
## SexM:PluckingYes           0.075 -0.131  0.175 -0.346 -0.656  0.417 -0.690
## ExpStart:SexM:CarotYes    -0.002 -0.002  0.008 -0.002 -0.003  0.473 -0.268
## ExpStart:SexM:PluckingYes  0.001  0.036 -0.115  0.488  0.458 -0.278  0.459
##                           ExS:SM ExS:CY SxM:CY ExS:PY SxM:PY ES:SM:C
## Tarsus                                                              
## Mass                                                                
## ExpStart                                                            
## SexM                                                                
## CarotYes                                                            
## PluckingYes                                                         
## ExpStart:SexM                                                       
## ExpStart:CarotYes         -0.009                                    
## SexM:CarotYes             -0.040  0.510                             
## ExpStart:PluckingYes       0.457 -0.566 -0.305                      
## SexM:PluckingYes           0.488 -0.278 -0.537  0.495               
## ExpStart:SexM:CarotYes     0.001 -0.770 -0.649  0.433  0.339        
## ExpStart:SexM:PluckingYes -0.697  0.421  0.362 -0.745 -0.672 -0.521 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.67403721 -0.67266044  0.03448938  0.58915181  2.07277186 
## 
## Number of Observations: 81
## Number of Groups: 
##               Aviary Identity %in% Aviary 
##                   12                   41
m2<-lme(Oxy ~ Tarsus+Mass+Exp*Carot+Sex*Carot+Exp*Sex*Plucking, random=list(Aviary=~1, Identity=~1), method="ML", data = Y, na.action=na.exclude)
anova(m1,m2)
##    Model df      AIC      BIC    logLik   Test    L.Ratio p-value
## m1     1 17 772.4558 813.1614 -369.2279                          
## m2     2 16 770.4925 808.8036 -369.2462 1 vs 2 0.03665401  0.8482

Je garde le modèle avec l’AIC le plus faible ou si la différence entre les AIC est plus petite que 2, le plus parcimonieux.

Modèle final

Les modèles finaux sont des modèles mixtes linéaire REML réalisés avec la même fonction lme et la fonction Anova du package car pour avoir les erreur de type 3.

mf<-lme(Oxy ~ Sex+Exp*Plucking, random=list(Aviary=~1, Identity=~1), method="REML", data = Y, na.action=na.exclude)
shapiro.test(residuals(mf))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mf)
## W = 0.99136, p-value = 0.8646
qqnorm(residuals(mf));qqline(residuals(mf))

library(car)
Anova(mf, type="III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: Oxy
##                 Chisq Df Pr(>Chisq)    
## (Intercept)  711.1752  1  < 2.2e-16 ***
## Sex           17.7193  1  2.560e-05 ***
## Exp           31.2056  1  2.321e-08 ***
## Plucking       5.6370  1    0.01759 *  
## Exp:Plucking   4.9661  1    0.02585 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mf)
## Linear mixed-effects model fit by REML
##  Data: Y 
##        AIC      BIC    logLik
##   740.8622 759.6127 -362.4311
## 
## Random effects:
##  Formula: ~1 | Aviary
##         (Intercept)
## StdDev:    7.375749
## 
##  Formula: ~1 | Identity %in% Aviary
##         (Intercept) Residual
## StdDev:    8.367793 22.21293
## 
## Fixed effects: Oxy ~ Sex + Exp * Plucking 
##                          Value Std.Error DF   t-value p-value
## (Intercept)          211.16883  7.918474 39 26.667869  0.0000
## SexM                 -24.19386  5.747544 28 -4.209425  0.0002
## ExpStart              46.90000  8.395700 39  5.586193  0.0000
## PluckingYes           21.58673  9.092049 10  2.374242  0.0390
## ExpStart:PluckingYes -23.05556 10.345869 39 -2.228479  0.0317
##  Correlation: 
##                      (Intr) SexM   ExpStr PlcknY
## SexM                 -0.363                     
## ExpStart             -0.530  0.000              
## PluckingYes          -0.736 -0.054  0.462       
## ExpStart:PluckingYes  0.430  0.000 -0.812 -0.569
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.64507717 -0.62459656  0.07843944  0.59250358  2.26879382 
## 
## Number of Observations: 82
## Number of Groups: 
##               Aviary Identity %in% Aviary 
##                   12                   41

Effets aléatoires

Pour les modèles finaux, les facteurs aléatoires ne sont pas significatifs. Je l’ai testé comme ça :

mfR<-lme(Oxy ~ Sex+Exp*Plucking, random=list(Aviary=~1, Identity=~1), method="REML", data = Y, na.action=na.exclude)
mfNR<-gls(Oxy ~ Sex+Exp*Plucking, method="REML", data = Y, na.action=na.exclude)
anova(mfR,mfNR)
##      Model df      AIC      BIC    logLik   Test L.Ratio p-value
## mfR      1  8 740.8622 759.6127 -362.4311                       
## mfNR     2  6 739.2564 753.3192 -363.6282 1 vs 2 2.39415  0.3021

Troisième et dernière étape : Expliquer les résultats trouvés

La dernière étape et pas des moindre est d’analyser et expliquer les résultats pour écrire un article !

Si vous avez lu jusqu’ici, merci !
Je reste à votre disposition pour toutes questions ou tous commentaires. Si vous voulez le pdf de l’article, le texte intégral est disponible ici.