T-tests for power-law distributed data?

Question

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?

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

Coefficients:
            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.test(values~ind,data=pdata)

    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 
             0.0668152 

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