Actually with 4 groups, you'd normally compare means using one ANOVA, not six t-tests (though t-tests would still come in with multiple comparisons or planned contrasts).
If one assumes Pareto distributions, then there are a number of possible approaches. I'll mention only a few, starting with one I think is perhaps easiest and also some nonparametric tests that would work. I'll assume that the left boundary ($x_m$, the scale parameter) is known; if that's not the case, things are more complicated (especially if they're unknown and not necessarily equal).
(1) comparison of the power parameter (shape parameter, $\alpha$) implies a comparison of population means (that is, equal $\alpha$ implies equal mean, different $\alpha$ implies different mean). Take $z = \log(y/x_m)$ and compare scale parameters of the resulting exponential distributions via a generalized linear model; it's straightforward to do ANOVA-type comparisons as well as pairwise comparisons.
(1a) A quick way to deal with unknown-but-equal $x_m$ is to take logs and subtract the smallest value (of the whole set) from all samples (losing that value from the data in the process). You can then proceed as above.
Here's an R example (with $x_m=1$). Similar fits can be done in pretty much any decent statistics package:
# create some (Pareto) data:
y1 <- c(814.660, 1.47520, 1.28029, 2.08808, 13.5882, 25.1290, 10.7137,
10.3032, 13.9075, 1556.73, 1.73512, 1783.04, 2.10658, 56.7400,
1.34085, 4.01592, 1.19537, 2.23376, 22.5796, 12.3961)
y2 <- c(332.949, 13.0680, 1.19512, 9.19466, 1.10640, 11.5778, 4.69242, 2.50173,
1.51986, 184.397, 2.61102, 17.86237, 6.01949, 76.9210, 3.66999)
pdata <- stack(list(y1=y1,y2=y2))
pcompfit <- glm(log(values)~ind,family=Gamma(link=log),data=pdata)
summary(pcompfit, dispersion=1) # dispersion = 1 for exponential
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9106 0.2236 4.072 4.66e-05 ***
indy2 -0.1311 0.3416 -0.384 0.701
(Dispersion parameter for Gamma family taken to be 1)
Null deviance: 34.674 on 34 degrees of freedom
Residual deviance: 34.528 on 33 degrees of freedom
(Some non-essential lines of output removed)
The test has (correctly) picked up that the parameter was smaller for the second sample, but for such small samples the effect size is too small to tell from random variation.
Exponential-distributed data might also be compared via
survreg in R.
(2) With the Pareto assumption (and $x_m\geq 1$), the log of the log (which is monotonic) reduces to a comparison of location with a shift alternative. You could therefore meaningfully compare the original data using Kruskal-Wallis, and it would be fairly easy to interpret the results under such a transformed shift-alternative.
# example - same data as before
# the 'anova-like' test (you don't need to transform to test for equality:
Kruskal-Wallis rank sum test
data: values by ind
Kruskal-Wallis chi-squared = 0.0544, df = 1, p-value = 0.8155
If you work on the loglog scale you can even produce a confidence interval for the log of the ratio of parameters. [We're not used to thinking of these tests as being a comparison of means, but with the specific distributional assumption and common $x_m$, it is a comparison of population means (as well as medians and so on).]
Here that's shown with a Mann-Whitney, for which interval construction is easy --
# the estimate of shift in log-parameter:
Wilcoxon rank sum test
data: log(log(values)) by ind
W = 157, p-value = 0.8307
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
difference in location
[However, something more appropriate like Dunn's test should be used for post hoc comparisons with the Kruskal-Wallis; though often suggested (e.g. at the bottom of this section, the Mann-Whitney is not quite the most appropriate choice for that purpose. One reason for that is discussed here.]
Note that this comparison is done 'the opposite way around' from the GLM ("first-second", not "second-first" as in the GLM). In any case it's trivial to flip the sign of the estimate and the end of the confidence interval, but it's important to know that it happens.
[The small difference in p-values is because the K-W is using a chi-square approximation, not the exact distribution of the test statistic. In large samples the approximation is very good.]
In large samples the estimate of the shift parameter should be more consistent between the two tests.
(3) you can do a likelihood ratio test. (This can deal with unknown $x_m$ without difficulty, but you may need to rely on the asymptotic chi-square result)