Le test de Dunnett dans R renvoyant des valeurs différentes à chaque fois

J'utilise la bibliothèque R 'multcomp' ( http://cran.r-project.org/web/packages/multcomp/ ) pour calculer le test de Dunnett. J'utilise le script ci-dessous:

Group <- factor(c("A","A","B","B","B","C","C","C","D","D","D","E","E","F","F","F"))
Value <- c(5,5.09901951359278,4.69041575982343,4.58257569495584,4.79583152331272,5,5.09901951359278,4.24264068711928,5.09901951359278,5.19615242270663,4.58257569495584,6.16441400296898,6.85565460040104,7.68114574786861,7.07106781186548,6.48074069840786)
data <- data.frame(Group, Value)
aov <- aov(Value ~ Group, data)
summary(glht(aov, linfct=mcp(Group="Dunnett")))

Maintenant, si j'exécute ce script plusieurs fois dans la console R, j'obtiens à chaque fois des résultats très légèrement différents. Voici un exemple:

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0 -0.35990    0.37009  -0.972  0.76545   
C - A == 0 -0.26896    0.37009  -0.727  0.90019   
D - A == 0 -0.09026    0.37009  -0.244  0.99894   
E - A == 0  1.46052    0.40541   3.603  0.01710 * 
F - A == 0  2.02814    0.37009   5.480  0.00104 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Et voici un autre:

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)    
B - A == 0 -0.35990    0.37009  -0.972   0.7654    
C - A == 0 -0.26896    0.37009  -0.727   0.9001    
D - A == 0 -0.09026    0.37009  -0.244   0.9989    
E - A == 0  1.46052    0.40541   3.603   0.0173 *  
F - A == 0  2.02814    0.37009   5.480   <0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Comme vous pouvez le voir, les deux résultats ci-dessus diffèrent très légèrement, mais il suffit de déplacer le groupe final (F) de deux étoiles à trois étoiles, ce que je trouve inquiétant.

J'ai plusieurs questions à ce sujet:

1. Pourquoi cela arrive-t-il?! Sûrement, si vous mettez les mêmes données à chaque fois, vous devriez obtenir les mêmes données.
2. Y a-t-il une sorte de nombre aléatoire utilisé quelque part dans le calcul de Dunnett?
3. Cette légère variation à chaque fois est-elle réellement un problème?

Je réponds ensemble à vos deux premières questions par l'exemple.

library(multcomp)

Group <- factor(c("A","A","B","B","B","C","C","C","D","D","D","E","E","F","F","F"))
Value <- c(5,5.09901951359278,4.69041575982343,4.58257569495584,4.79583152331272,5,5.09901951359278,4.24264068711928,5.09901951359278,5.19615242270663,4.58257569495584,6.16441400296898,6.85565460040104,7.68114574786861,7.07106781186548,6.48074069840786)
data <- data.frame(Group, Value)

fit <- aov(Value ~ Group, data)

set.seed(20140123)
Dunnet <- glht(fit, linfct=mcp(Group="Dunnett"))
summary(Dunnet)

Résultats:

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0 -0.35990    0.37009  -0.972  0.76536   
C - A == 0 -0.26896    0.37009  -0.727  0.90012   
D - A == 0 -0.09026    0.37009  -0.244  0.99895   
E - A == 0  1.46052    0.40541   3.603  0.01794 * 
F - A == 0  2.02814    0.37009   5.480  0.00112 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Exécutez à nouveau (sans définir la valeur de départ):

summary(Dunnet)

Différents résultats:

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0 -0.35990    0.37009  -0.972  0.76535   
C - A == 0 -0.26896    0.37009  -0.727  0.90020   
D - A == 0 -0.09026    0.37009  -0.244  0.99895   
E - A == 0  1.46052    0.40541   3.603  0.01767 * 
F - A == 0  2.02814    0.37009   5.480  0.00105 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Exécutez à nouveau (avec une graine définie):

set.seed(20140123)
Dunnet <- glht(fit, linfct=mcp(Group="Dunnett"))
summary(Dunnet)

Mêmes résultats:

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = Value ~ Group, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0 -0.35990    0.37009  -0.972  0.76536   
C - A == 0 -0.26896    0.37009  -0.727  0.90012   
D - A == 0 -0.09026    0.37009  -0.244  0.99895   
E - A == 0  1.46052    0.40541   3.603  0.01794 * 
F - A == 0  2.02814    0.37009   5.480  0.00112 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

En définissant la graine avant chaque exécution, vous obtenez des résultats cohérents. Il apparaît donc qu'un nombre aléatoire est utilisé dans le calcul des valeurs de p.

$alpha$

Vous avez raison, une génération de nombres aléatoires est impliquée et les calculs varient d'une exécution à l'autre. Le coupable n'est en fait pas la procédure de Dunnett, mais la distribution t à plusieurs variables requise pour l'ajustement en une seule étape.

Le code suivant montre un exemple de calcul $P(X<0)$ avec un vecteur à 5 dimensions $X$ avoir plusieurs variables $T_5$ distribution avec corrélation échangeable:

> library(mvtnorm)
> cr2 <- matrix(rep(0.3, 25), nr=5); diag(cr2) <- 1
> cr2
     [,1] [,2] [,3] [,4] [,5]
[1,]  1.0  0.3  0.3  0.3  0.3
[2,]  0.3  1.0  0.3  0.3  0.3
[3,]  0.3  0.3  1.0  0.3  0.3
[4,]  0.3  0.3  0.3  1.0  0.3
[5,]  0.3  0.3  0.3  0.3  1.0
> b <- pmvt(lower=rep(-Inf,5), upper=rep(0,5), delta=rep(0,5), df=5, corr=cr2)
> a <- pmvt(lower=rep(-Inf,5), upper=rep(0,5), delta=rep(0,5), df=5, corr=cr2)
> all.equal(a,b)
[1] "Attributes: < Component 1: Mean relative difference: 0.1527122 >"
[2] "Mean relative difference: 0.0003698006"     

Si cela vous inquiète, appelez simplement set.seedavec n'importe quel argument avant le calcul pour le rendre exactement reproductible.

Soit dit en passant, il y a un accusé de réception et une quantification de l'erreur dans la sortie de glht:

> ss <- summary(glht(aov, linfct=mcp(Group="Dunnett")))
> attr(ss$test$pvalues, "error")
[1] 0.0006597562