# Dunnett's test in R returning different values each time

I am using the R 'multcomp' library (http://cran.r-project.org/web/packages/multcomp/) to calculate Dunnett's test. I am using the script below:

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")))


Now if I run this script through the R Console multiple times I get very slightly different results each time. Here's one example:

         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)


And here's another:

         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)


As you can see, the above two results differ very slightly, but it's enough to move the final group (F) from two stars to three stars, which I find worrying.

I have several questions regarding this:

1. Why is this happening?! Surely if you put the same data in each time you should get the same data out.
2. Is there some kind of random number being used somewhere in the Dunnett's calculation?
3. Is this slight variation each time actually a problem?

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)


Results:

     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)


Run again (without setting the seed):

summary(Dunnet)


Different results:

     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)


Run again (with a set seed):

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


Same results:

     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)


By setting the seed before each run, you get consistent results. Therefore it appears that a random number is being used in the calculation of the p-values.

Do I think this slight variation is a problem? I don't really like it, but I'd live with it. Using a set seed will make your results reproducible. I'd recommend not thinking about p-values in terms of how many stars are next to them – rather choose an $alpha$ that is meaningful and useful for you. I try not to get caught up in what is happening 5 or 6 decimal points out unless the project I am working on really requires that level of precision. In this case I think most people would agree that even if the p-value calculation changes slightly the interpretation of the results would be the same.

• Thanks very much for your answer. I think you are right about not thinking of how many stars are there - people should be looking at the P-value anyway. I think I will have to set the seed to a known value though, because in order to validate my program the resuls must be reproducible exactly. Just one more question - do you know why the random seed is used? Jan 23, 2014 at 14:57
• See the answer written by @Aniko which provides a more detailed explanation. Notice I used today's date as the seed. Jan 23, 2014 at 15:43

You are correct, there is a random number generation involved, and it makes the calculations vary from run-to-run. The culprit is actually not Dunnett's procedure, but the multivariate t distribution required for the single-step adjustment.

The following code shows an example calculating $P(X<0)$ with a 5-dimensional vector $X$ having multivariate $T_5$ distribution with exchangeable correlation:

> 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"


If this is of concern, just call set.seed with any argument before the calculation to make it exactly reproducible.

By the way, there is an acknowledgement and quantification of the error in the output of glht:

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