I have empirical data consisting of a variable that is arranged in 4 groups. In each of the groups the variable roughly follows a power-law distribution (Pareto distribution). I want to test for differences in means across the group, as I would with a t-test of the mean for normally distributed data.

What is the best way of going about doing this?

  • $\begingroup$ When you say it follows 'a power-law distribution', do you mean a Pareto distribution? or something else? $\endgroup$
    – Glen_b
    Feb 13, 2014 at 20:37
  • $\begingroup$ @Glen_b - yes, a Pareto distribution. $\endgroup$
    – histelheim
    Feb 13, 2014 at 20:43
  • 1
    $\begingroup$ With a Pareto distributions, if the shape parameter is $\alpha$ (so the exponent in the power law is $-(\alpha+1)$ ) with $0 \lt \alpha \le 1$, then the mean will be infinite. You can have different $\alpha$s each giving an infinite mean. $\endgroup$
    – Henry
    Feb 14, 2014 at 7:45

1 Answer 1


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
AIC: 135.73

(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:
> wilcox.test(log(log(values))~ind,data=pdata,conf.int=TRUE)

    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:
 -0.8001905  0.8935245
sample estimates:
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)

  • 1
    $\begingroup$ +1, I think this is good, but I think it would be even better if you could add a simple R based example for each option, if that's not too inconvenient. $\endgroup$ Feb 13, 2014 at 22:11
  • 1
    $\begingroup$ @gung I agree; I've been working on doing that for the first two already. $\endgroup$
    – Glen_b
    Feb 13, 2014 at 22:13
  • $\begingroup$ Empirically, I can see that the log of the log makes my data normal. However, what is the substantive interpretation of "double-logged" data (i.e. log(log(variable)))? I understand that taking the logarithm once enables you to make relative as opposed to absolute comparisons between, say, averages, but what does the log of a logged value "mean"? $\endgroup$
    – histelheim
    Feb 18, 2014 at 19:07
  • 1
    $\begingroup$ log-log doesn't make the a Pareto-distribution normal - the result is definitely left skew. It makes it constant variance with mean the log of the index. Sorry, I have no good interpretation of the log-log, it was just presented as a way of doing inference on the parameter. The first log reduces the shape parameter (the index of the Pareto) to a scale parameter and the second one makes it (or rather, its log) a location parameter. My own instinct would be to work on the log scale with the exponential distribution. $\endgroup$
    – Glen_b
    Feb 18, 2014 at 21:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.