12
$\begingroup$

Is it possible to do a power analysis for a 2-sided Kolmogorov Smirnov test in R?

I am testing whether two empirical distributions differ using ks.test(), and am looking to add a power analysis.

I haven't been able to find any built-in power analyses for KS tests in R. Any suggestions?

Edit: These are randomly generated distributions that closely approximate my data (with true sample sizes and estimated decay rates for exponential distributions)

set.seed(100)
x <- rexp(64, rate=0.34)
y <- rexp(54,rate=0.37)

#K-S test: Do x and y come from same distribution?
ks.test(x,y)

These data are measures of body size in two different groups. I want to show that the two groups have essentially the same distribution, but was asked by a collaborator if I had the power to say that based on sample sizes. I've randomly drawn from an exponential distribution here, but these are close to the real data.

So far, I've said that there is no significant difference in these distributions based on the two-sided K-S test. I've also plotted the two distributions. How can I show that I have the power to make such a statement, given the sample sizes and decay rates for x and y?

$\endgroup$
7
  • 5
    $\begingroup$ Power will depend on many things which is why there wouldn't be a built in for the two sample test. You can simulate for given situations. So: Power given what assumptions about the situation? Against what alternative or sequence of alternatives? For example, you could compute (simulate) a power curve for exponential distributed data, against a set of scale shift alternatives. Or you could compute power for a normal against a location-shift. Or you could compute power in a Weibull as you vary the shape parameter. Do you have some additional details? $\endgroup$
    – Glen_b
    Mar 24, 2015 at 22:42
  • $\begingroup$ To actually compute power you'd need sample sizes as well. If you're trying to identify sample size given specified power against a particular alternative, that could be done via root-finding, but often you can find the point by simple approaches (trying a couple of sample sizes is usually sufficient to get very close). $\endgroup$
    – Glen_b
    Mar 24, 2015 at 22:51
  • $\begingroup$ What variable is being measured? Are these times? $\endgroup$
    – Glen_b
    Mar 25, 2015 at 22:22
  • $\begingroup$ @Glen_b These are not times. They are measures of body size in two different groups. I want to show that the two groups have essentially the same distribution, but was asked if I had the power to say that based on sample sizes. $\endgroup$
    – Sarah
    Mar 26, 2015 at 18:24
  • 1
    $\begingroup$ Ah! That's two pieces of useful context that may help to have in your question. So the idea is if you show that power to identify some notionally modest differences was reasonable, one could take failure to reject as an indication that the difference was small. Yes, a prior power analysis can help with making that argument. After the fact I'd probably focus more on something like a estimate (and confidence interval perhaps) of the change in scale as an indication that the difference was actually small in size, as well as a plot of the two sample cdfs. $\endgroup$
    – Glen_b
    Mar 27, 2015 at 0:45

2 Answers 2

18
$\begingroup$

Finding power against exponential scale-shift alternatives is reasonably straightforward.

However, I don't know that you should be using values calculated from your data to work out what the power might have been. That sort of post hoc power calculation tends to result in counter-intuitive (and often misleading) conclusions.

Power, like significance level, is a population-level phenomenon and one you deal with before the fact; you'd use a priori understanding (including theory, reasoning or any previous studies) to decide on a reasonable set of alternatives to consider, and a desirable effect size

You can also consider a variety of other alternatives (e.g. you could embed the exponential inside a gamma family to consider the impact of more or less skew cases).

The usual questions one might try to answer by a power analysis are:

  1. what's the power, for a given sample size, at some effect size or set of effect sizes*?

  2. given a sample size and power, how large an effect is detectable?

  3. Given a desired power for a particular effect size, what sample size would be required?

*(where here 'effect size' is intended generically, and might be for example, a particular ratio of means, or difference of means, not necessarily standardized).

Clearly you already have a sample size, so you're not in case (3). You might reasonably consider case (2) or case (1).

I'd suggest case (1) (which also gives a way to deal with case (2)).

To illustrate an approach to case (1) and see how it relates to case (2), let's consider a specific example, with:

  • scale shift alternatives

  • exponential populations

  • sample sizes in the two samples of 64 and 54

Because the sample sizes are different, we have to consider the case where the relative spread in one of the samples is both smaller and larger than 1 (if they were the same size, symmetry considerations make it possible to consider just one side). However, because they're quite close to the same size, the effect is very small. In any case, fix the parameter for one of the samples and vary the other.

So what one does is:

Beforehand:

choose a set of scale multipliers representing different alternatives
select an nsim (say 1000)
set mu1=1

To do the calculations:

for each possible scale multiplier, kappa 
  repeat nsim times
    generate a sample of size n1 from Exp(mu1) and n2 from Exp(kappa*mu1)
    perform the test
  compute the rejection rate across nsim tests at this kappa

In R, I did this:

alpha = 0.05
n1 = 54
n2 = 64
nsim = 10000
s = c(1.1,1.2,1.5,2,2.5,3) # set up grid for kappa
s = c(1/rev(s),1,s)        #  also below and at 1
rr = array(NA,length(s))   # to hold rejection rates

for(i in seq_along(s)) rr[i]=mean(replicate(nsim,
                                    ks.test(rexp(n1,1),rexp(n2,s[i]))$p.value)<alpha
                                 )

plot(rr~s,log="x",ylim=c(0,1),type="n") #set up plot
points(rr~rev(s),col=3) # plot the reversed case to show the (tiny) asymmetry+noise
points(rr~s,col=1) # plot the "real" case last 
abline(h=alpha,col=8,lty=2) # draw in alpha

which gives the following power "curve"

enter image description here

The x-axis is on a log-scale, the y-axis is the rejection rate.

It's hard to tell here, but the black points are slightly higher on the left than the right (that is, there's fractionally more power when the larger sample has the smaller scale).

Using the inverse normal cdf as a transformation of the rejection rate, we can make the relationship between the transformed rejection rate and log kappa (kappa is s in the plot, but the x-axis is log-scaled) very nearly linear (except near 0), and the number of simulations was high enough that the noise is very low -- we can just about ignore it for present purposes.

So we can just use linear interpolation. Shown below are approximate effect sizes for 50% and 80% power at your sample sizes:

enter image description here

The effect sizes on the other side (larger group has smaller scale) are only slightly shifted from that (can pick up a fractionally smaller effect size), but it makes little difference, so I won't labor the point.

So the test will pick up a substantial difference (from a ratio of scales of 1), but not a small one.


Now for some comments: I don't think hypothesis tests are particularly relevant to the underlying question of interest (are they quite similar?), and consequently these power calculations don't tell us anything directly relevant to that question.

I think you address that more useful question by prespecifying what you think "essentially the same" actually means, operationally. That - pursued rationally to a statistical activity - should lead to meaningful analysis of the data.

$\endgroup$
1
  • $\begingroup$ Thank you so much! This is really helpful, much appreciated. $\endgroup$
    – Sarah
    Apr 1, 2015 at 15:01
1
$\begingroup$

Since Kolmogorov-Smirnov is non-parametric so by definition there can be no applicable power analysis. To have some sort of estimation, you need to assume a background model (and thus divert from the non-parametric world...) and use it to compute one of the following: sample-size, MDE, or power (i.e., you fix/choose two and compute the third).

$\endgroup$

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.