9
$\begingroup$

I want to assess the accuracy of normality tests over different sample sizes in R (I realize that normality tests may be misleading). For example, to look at the Shapiro-Wilk test, I'm conducting the following simulation (as well as plotting the results) and would expect that as the sample size increases the probability of rejecting the null decreases:

n <- 1000
pvalue_mat <- matrix(NA, ncol = 1, nrow = n)

for(i in 10:n){
    x1 <- rnorm(i, mean = 0, sd = 1)
    pvalue_mat[i,] <- shapiro.test(x1)$p.value
}   

plot(pvalue_mat)

My thought would be that as the sample size grows there should be a lower rejection rate, however it seems pretty uniform. I think I am misunderstanding this - any and all thoughts welcome.

$\endgroup$
1

3 Answers 3

8
$\begingroup$

You are simulating under the null hypothesis (normal distribution), therefore the rejection rate will tend to the significance level as expected. To assess the power, you need to simulate under any non-normal distribution. There are infinite possibilities/scenarios (e.g. gamma distributions with increasing skewness, t-distribution with decreasing df etc.) to choose from, depending on the scope of your study.

$\endgroup$
2
  • $\begingroup$ Thanks for the response. When I simulate over non-normal distributions, I observe a convex pattern with respect to the origin - i.e. as the sample size gets bigger for any non-normal distribution, the probability of rejecting the null of normality increases. However, I don't understand why it isn't the inverse when drawing from a normal distribution: why doesn't the probability of rejecting the null decrease as the sample size gets bigger? Thanks $\endgroup$
    – user94759
    Oct 29, 2013 at 19:28
  • 3
    $\begingroup$ Because the probability of commiting such type-1-error is by definition equal to the significance level, which is constant. Or put differently, p values are uniformly distributed under the null. Btw, it is advisable to do many simulations per setting, including the choice of n, not just one as in your code. $\endgroup$
    – Michael M
    Oct 29, 2013 at 20:11
9
$\begingroup$

Understanding power analysis of statistical hypothesis tests can be enhanced by carrying some out and looking closely at the results.


By design, a test of size $\alpha$ is intended to reject the null hypothesis with a chance of at least $\alpha$ when the null is true (its expected false positive rate). When we have the ability (or luxury) of choosing among alternative procedures with this property we would prefer those that (a) actually come close to the nominal false positive rate and (b) have relatively higher chances of rejecting the null hypothesis when it is not true.

The second criterion requires us to stipulate in what way(s) and by how much the null fails to be true. In textbook cases this is easy, because the alternatives are limited in scope and clearly specified. With distribution tests like the Shapiro-Wilk, the alternative are much more vague: they are "non-normal." When choosing among distribution tests, then, the analyst is likely to have to conduct their own one-off power study to assess how well the tests work against more specific alternative hypotheses that are of concern in the problem at hand.

An example motivated by Michael Mayer's answer posits that the alternative distribution may have qualities similar to those of the family of Student t distributions. This family, parameterized by a number $\nu\ge 1$ (as well as by location and scale) includes in the limit of large $\nu$ the Normal distributions.

In either situation--whether evaluating the actual test size or its power--we must generate independent samples from a specified distribution, run the test on each sample, and find the rate at which it rejects the null hypothesis. However, there is more information available in any test result: its P-value. By retaining the set of P-values produced during such a simulation, we can later assess the rate at which the test would reject the null for any value of $\alpha$ we might care about. The heart of the power analysis, then, is a subroutine that generates this P-value distribution (either by simulation, as just described, or--occasionally--with a theoretical formula). Here is an example coded in R. Its arguments include

  • rdist, the name of a function to produce a random sample from some distribution

  • n, the size of samples to request of rdist

  • n.iter, the number of such samples to obtain

  • ..., any optional parameters to be passed on to rdist (such as the degrees of freedom $\nu$).

The remaining parameters control the display of the results; they are included mainly as a convenience for generating the figures in this answer.

sim <- function(rdist, n, n.iter, prefix="",
                breaks=seq(0, 1, length.out=20), alpha=0.05,
                plot=TRUE, ...) {

  # The simulated P-values.
  # NB: The optional arguments "..." are passed to `rdist` to specify
  #     its parameters (if any).
  x <- apply(matrix(rdist(n*n.iter, ...), ncol=n.iter), 2, 
             function(y) shapiro.test(y)$p.value)

  # The histogram of P-values, if requested.
  if (plot) {
    power <- mean(x <= alpha)
    round.n <- 1+ceiling(log(1 + n.iter * power * (1-power), base=10) / 2)
    hist(x[x <= max(breaks)], xlab=paste("P value (n=", n, ")", sep=""), 
         breaks=breaks, 
         main=paste(prefix, "(power=", format(power, digits=round.n), ")", sep=""))
    # Specially color the "significant" part of the histogram
    hist(x[x <= alpha], breaks=breaks, col="#e0404080", add=TRUE)
  }

  # Return the array of P-values for any further processing.
  return(x)
}

You can see the computation actually takes just one line; the rest of the code plots the histogram. To illustrate, let's use it to compute the expected false positive rates. "Rates" is in the plural because the properties of a test usually vary with the sample size. Since it is well-known that distributional tests have high power against qualitatively small alternatives when sample sizes are large, this study focuses on a range of small sample sizes where such tests of often applied in practice: typically about $5$ to $100.$ To save computation time, I report only on values of $n$ from $5$ to $20.$

n.iter <- 10^5                 # Number of samples to generate
n.spec <- c(5, 10, 20)         # Sample sizes to study
par(mfrow=c(1,length(n.spec))) # Organize subsequent plots into a tableau
system.time(
  invisible(sapply(n.spec, function(n) sim(rnorm, n, n.iter, prefix="DF = Inf ")))
)

After specifying the parameters, this code also is just one line. It yields the following output:

Histograms for the null

This is the expected appearance: the histograms show nearly uniform distributions of P-values across the full range from $0$ to $1$. With the nominal size set at $\alpha=0.05,$ the simulations report between $.0481$ and $0.0499$ of the P-values were actually less than that threshold: these are the results highlighted in red. The closeness of these frequencies to the nominal value attests that the Shapiro-Wilk test does perform as advertised.

(There does seem to be a tendency towards an unusually high frequency of P-values near $1$. This is of little concern, because in almost all applications the only P-values one looks at are $0.2$ or less.)

Let's turn now to assessing the power. The full range of values of $\nu$ for the Student t distribution can adequately be studied by assessing a few instances from around $\nu=100$ down to $\nu=1$. How do I know that? I performed some preliminary runs using very small numbers of iterations (from $100$ to $1000$), which takes no time at all. The code now requires a double loop (and in more complex situations we often need triple or quadruple loops to accommodate all the aspects we need to vary): one to study how the power varies with the sample size and another to study how it varies with the degrees of freedom. Once again, though, everything is done in just one line of code (the third and final):

df.spec <- c(64, 16, 4, 2, 1)
par(mfrow=c(length(n.spec), length(df.spec)))
for (n in n.spec) 
  for (df in df.spec)
    tmp <- sim(rt, n, n.iter, prefix=paste("DF =", df, ""), df=df)

Histograms for the alternatives

A little study of this tableau provides good intuition about power. I would like to draw attention to its most salient and useful aspects:

  • As the degrees of freedom reduce from $\nu=64$ on the left to $\nu=1$ on the right, more and more of the P-values are small, showing that the power to discriminate these distributions from a Normal distribution increases. (The power is quantified in each plot title: it equals the proportion of the histogram's area that is red.)

  • As the sample size increase from $n=5$ on the top row to $n=20$ on the bottom, the power also increases.

  • Notice how as the alternative distribution differs more from the null distribution and the sample size increases, the P-values start collecting to the left, but there is still a "tail" of them stretching all the way to $1$. This is characteristic of power studies. It shows that testing is a gamble: even when the null hypothesis is flagrantly violated and even when our sample size is reasonably large, our formal test may fail to produce a significant result.

  • Even in the extreme case at the bottom right, where a sample of $20$ is drawn from a Student t distribution with $1$ degree of freedom (a Cauchy distribution), the power is not $1$: there is a $100 - 86.57 = 13\%$ chance that a sample of $20$ iid Cauchy variates will not be considered significantly different from Normal at a level of $5\%$ (that is, with $95\%$ confidence).

  • We could assess the power at any value of $\alpha$ we choose by coloring more or fewer of the bars on these histograms. For instance, to evaluate the power at $\alpha=0.10$, color in the left two bars on each histogram and estimate its area as a fraction of the total.

    (This won't work too well for values of $\alpha$ smaller than $0.05$ with this figure. In practice, one would limit the histograms to P-values only in the range that would be used, perhaps from $0$ to $20\%$, and show them in enough detail to enable visual assessment of power down to $\alpha=0.01$ or even $\alpha=0.005$. (That is what the breaks option to sim is for.) Post-processing of the simulation results can provide even more detail.)


It is amusing that so much can be gleaned from what, in effect, amounts to three lines of code: one to simulate i.i.d. samples from a specified distribution, one to apply that to an array of null distributions, and the third to apply it to an array of alternative distributions. These are the three steps that go into any power analysis: the rest is just summarizing and interpreting the results.

$\endgroup$
1
$\begingroup$

(More than a comment, perhaps not a complete answer)

[I] would expect that as the sample size increases the probability of rejecting the null decreases

Leaving aside considerations of biased tests (which are not uncommon in goodness of fit, so it's worth a mention), there are three situations relating to rejection rate one might want to consider:

1) the rejection rate when simulating from the null (as you seem to be doing in your question)

Here, the rejection rate should be at or near the significance level, so, no, if you hold the significance level constant, the rejection rate doesn't decrease as n increases, but stays at/near $\alpha$.

2) the rejection rate when simulating from some alternative

Here the rejection rate should increase as n increases.

3) the rejection rate for some collection of real data

Practically, the null is never actually true, and real data will have some mixture of amounts of non-normality (as measured by the test statistic). If the degree of non-normality is not related to sample size, the rejection rate should increase as n increases.

So in fact, in none of these situations should we see the rejection rate decrease with sample size.

$\endgroup$

Your Answer

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

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