3
$\begingroup$

When simulating the Shapiro test for normal distributed data in R I encountered that the p-values for rejecting the Null-Hypothesis often increase in the range of ~30 to ~200 observations, than decrease again and than somehow fluctuate. How can I explain this behavior? Intuitively I thought that with a larger number of observations, the p-values would increase as well and that I would get low p-values especially if I have a small number of observations.

test <- vector(length = length(10:5000))
for(i in 10:5000){
  set.seed(500)
  test[i] <- shapiro.test(rnorm(i))[[2]]}
plot(test)

enter image description here

If you change the seed the results look different, however the general pattern is often similar.

$\endgroup$
3
  • 1
    $\begingroup$ This isn't telling you anything useful about the Shapiro test or anything else. Remove the seed setting. If you don't, you induce extremely strong relationships among successive results and all you get is a variation of a random walk. $\endgroup$
    – whuber
    Commented Feb 21, 2018 at 15:31
  • $\begingroup$ If I remove the seed setting than there is no pattern at all any more with p-values between 0-1 for any level of observations... shouldn't be the general pattern an increase in robustness and hence higher p-values? Is this a specific problem of the rnorm function? $\endgroup$
    – joaoal
    Commented Feb 22, 2018 at 8:16
  • 4
    $\begingroup$ No, the p-value is defined in such a way that when it is correctly and accurately calculated, when the distribution of the test statistic is continuous, and the null hypothesis holds, it will be randomly and uniformly distributed between $0$ and $1$. All of those circumstances apply to your experiment. $\endgroup$
    – whuber
    Commented Feb 22, 2018 at 14:32

1 Answer 1

5
$\begingroup$

This question is resolved by applying software engineering principles.


The plot is a very roundabout and computationally expensive way to construct a random walk.

To see what this means, let's redo the code to make it faster and clearer:

set.seed(500)
x <- rnorm(5000)
test <- c(rep(0, 9), sapply(10:5000, \(n) shapiro.test(x[1:n])$p.value))

You can plot test to verify the results are identical. Notice that the data vector x never changes. The sapply function loops over the indices $n = 10, 11, \ldots, 5000$ and applies the function shapiro.test to the first $n$ entries in the constant vector x.

The same qualitative behavior arises with most other functions. Consider this variant where shapiro.test is replaced by a t test:

test <- c(rep(0, 9), sapply(10:5000, \(n) t.test(x[1:n])$p.value))
plot(test, type = "l")

When applied to the very same vector x, the plot is qualitatively like that in the question:

enter image description here

Abstractly, let $f$ represent the function shapiro.test or t.test or whatever. Its values are p-values in the range $[0,1].$ The plot maps each $n$ to $f(x_1,x_2,\ldots, x_n).$

This abstract description helps clarify how the points on the plot evolve from left to right: $f(x_1,\ldots, x_n)$ is followed by an updated version $f(x_1,\ldots, x_n, x_{n+1}),$ where $x_{n+1}$ is a random value independent of the first $n$ values in $x.$

Let's express this conventionally. The vector x is the computer's realization of the mathematical sequence $X_1, X_2,\ldots, X_n,\ldots$ of independent random variables (all with standard Normal distributions in this example). The sequence being plotted is

$$S = f(X_1), f(X_1,X_2), \ldots, f(X_1,X_2,\ldots,X_n), \ldots$$

It generalizes the archetypical random walk where $S_n = X_1 + X_2 + \cdots + X_n.$ Unless $f$ doesn't depend at all on $x_1,\ldots, x_n,$ the updated value must be related to the previous value: successive points on the plot are not independent. That's what makes the plot wander around rather than producing a bunch of scattered points.


The intention behind this code might have been to study what happens to the distributions of the Shapiro-Wilk p-values as they might depend on the sample size. If so, the code is erroneous.

A good approach to implementing the intention is to encapsulate the problem by putting the data generation and testing into a function. Below, the function simulate generates n.sim samples of size n using a common random number generator distribution, applies the function $f$ (passed as statistic), and returns n p-values. By wrapping all this in a call to sapply it allows the caller to pass a vector of various sample sizes in n, returning an array whose columns correspond to those sample sizes. All this can be done in one line.

simulate <- function(n, statistic, distribution, n.sim = 5000, ...) 
  sapply(n, \(n) apply(matrix(distribution(n * n.sim, ...), n), 2, statistic))

(The ellipsis ... permits the caller to pass named arguments to distribution to vary its parameters if desired.)

Let's use this function to study the Shapiro-Wilk p-values for iid Normal samples (the null distribution). This one uses sample sizes of $10,$ $20,$ $50,$ and $100.$ Because the main computation was encapsulated, this requires only one line, making it easy to verify the correctness of the implementation (subject to verifying that simulate works as intended).

n <- c(10, 20, 50, 100) # Sample sizes
set.seed(17)
Sim <- simulate(n, function(x) shapiro.test(x)$p.value, rnorm)

Here is what plots of the columns look like.

for (i in seq_along(n)) 
  plot(Sim[, i], col = "#00000020", ylab = "p value", main = bquote(n==.(n[i])))

enter image description here

The p-values are (literally) all over the place, exactly as one would hope when the data are generated from the null distribution for the test.

Maybe it's better to examine histograms of these p-values, now that we have seen there's no sequential pattern to them.

for (i in seq_along(n)) {
  hist(Sim[, i], freq = FALSE, xlim = 0:1, xlab = "p value", main = bquote(n==.(n[i])))

enter image description here

They look just like samples of a uniform distribution. Everything works as it should, regardless of the sample size.

The extreme abstraction of simulate pays off here. By making one tiny change we can study how well the Shapiro-Wilk test performs under a specific alternative hypothesis. Consider, for instance, what happens when the data are drawn from a Gamma$(5)$ distribution:

Sim <- simulate(n, function(x) shapiro.test(x)$p.value, rgamma, shape = 5)

The p-value distributions now look like these:

enter image description here

As the sample size grows, a small p-value becomes more and more likely. The test grows more powerful with increasing 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.