Se for balanceado, use a função agricolae::HSD.test() que solta essa estimativa na saída da função, veja

> x <- gl(5,5)
> y <- rnorm(x)
> agricolae::HSD.test(y, x, DFerror=20, MSError=1)

Study:

HSD Test for y

Mean Square Error:  1

x,  means

           y   std.err replication
1 -0.3230810 0.1327155           5
2  0.3595503 0.5492514           5
3 -0.1703578 0.3377758           5
4 -0.1483504 0.5129000           5
5 -0.1101464 0.2727455           5

alpha: 0.05 ; Df Error: 20
Critical Value of Studentized Range: 4.231857

Honestly Significant Difference: 1.892544 <------------------------------------------ DMS

Means with the same letter are not significantly different.

Groups, Treatments and means
a      2      0.3595503
a      5      -0.1101464
a      4      -0.1483504
a      3      -0.1703578
a      1      -0.323081
>

==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: walmes@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218