2
$\begingroup$

I have a number of samples that I would like to test for normality. One of the samples exceeds 5,000 data points, the limit up to which the shapiro test accepts samples. This is the data:

c1 <- exp(rnorm(505))
c2 <- exp(rnorm(550))
c3 <- exp(rnorm(5500))

cluster.data <- c(c1, c2, c3)
cluster.factors <- c(rep("Cluster_1", length(c1)), 
                     rep("Cluster_2", length(c2)),
                     rep("Cluster_3", length(c3)))

# set up data for test:
cluster.df <- data.frame(cluster.data, cluster.factors)

To circumvent the 5,000 restriction, would it be statistically acceptable if I run the test on smallish subsamples of the data only? Here, for example, I draw a subsample of size 500 for all three variables:

tapply(cluster.df[,1], cluster.df[,2], function(x) shapiro.test(sample(x, 500)))

And the test returns sigificant results for all three:

$Cluster_1

    Shapiro-Wilk normality test

data:  sample(x, 500)
W = 0.59561, p-value < 2.2e-16


$Cluster_2

    Shapiro-Wilk normality test

data:  sample(x, 500)
W = 0.57891, p-value < 2.2e-16


$Cluster_3

    Shapiro-Wilk normality test

data:  sample(x, 500)
W = 0.67686, p-value < 2.2e-16
$\endgroup$
9
$\begingroup$

I have comments on five levels.

  1. This is a deficiency of a particular R function shapiro.test() and need not imply that that there aren't other ways to do it in R, on which I can't advise specifically. It may or may not be of practical relevance to you that no such limit applies to all software. For example, the Stata command swilk isn't limited in this way.

  2. I can't comment on why that particular function won't perform, but the larger question is why you are doing this kind of testing at all. A good reason not to care is generic: for sample sizes of that order, or even larger, such tests are arguably fairly useless as even minute deviations from normality will qualify as significant at conventional levels. More specifically: why is it important or interesting to test for normality? People often apply such tests to marginal distributions given a widespread myth that marginal normality is a requirement for very many procedures. Where normality is a relevant assumption, or ideal condition, it usually applies to distributions conditional on a structure of mean outcomes or responses.

  3. In response to your specific query of whether subsampling is acceptable, the serious reply in return is acceptable in what sense? A personal reply: as a reader, author and reviewer of statistical papers, and as a statistical journal editor, my reaction would be to suggest that such subsampling is at best awkward and at worst an avoidance of the main issue, which would be to find an implementation without such a limit, or more likely to think about the distribution in different terms.

  4. As often emphasised on CV, and elsewhere, the most helpful and informative way to check departure from normality is a normal quantile plot, often also called a normal probability plot, a normal scores plot, or a probit plot. Such a plot not only provides a visual assessment of degree of non-normality, it makes precise in what sense there are departures from the ideal shape. The lack of an associated P-value is not in practice much of a loss, although the procedure may be given some inferential impetus through confidence levels, simulations and so forth.

  5. Specifically, your examples consist of generating lognormal samples and then establishing that indeed they fail to qualify as normal with P-values $\ll 10^{-15}$. That has to seem puzzling, but be reassured that with larger samples your P-values will be, or should be, even more minute, subject to a machine level question of the minimum reportable P-value here. Conversely, it may well be that your real problem lies elsewhere and these examples are no more than incidental illlustrations.

$\endgroup$
  • 2
    $\begingroup$ +1. But maybe it would also help to point out the obvious: if one takes a single random subsample, tests it, and rejects the null, then the null can be rejected for the entire sample. Given the expectation of rejecting the null on a suitably large sample, this could be an effective strategy--FWIW. A second obvious point to make explicitly would be to ask what's the point of conducting a distributional test in the first place. $\endgroup$ – whuber Jan 24 at 14:23
  • 1
    $\begingroup$ @whuber Thanks for comments. I have tinkered with my answer to reflect your "second obvious point". $\endgroup$ – Nick Cox Jan 24 at 15:54
1
$\begingroup$

I think Nick Cox points out some of the difficulties with the approach.

A possible alternate recommendation would be to use another normality test. In classes I took we used a test based on skewness and kurtosis due to D'Agostino for larger samples. I implemented these tests in my lolcat statistical package. Consider:

#Install/load step
require(devtools)
install_github("burrm/lolcat")
require(lolcat)

set.seed(1)

#Normal distribution - no rejection
zz <- rnorm(5500)
skewness.test(zz)
kurtosis.test(zz)

# Log normal distribution - rejection on both skewness and kurtosis
zz1 <- exp(zz1)
skewness.test(zz1)
kurtosis.test(zz1)

Interestingly enough, even with a sample size of 5500, skewness/kurtosis would likely not reject with these tests. A log normal distribution would most likely reject, even at substantially lower sample sizes. As an example:

> set.seed(1)
> 
> #Normal distribution - no rejection
> zz <- rnorm(5500)
> skewness.test(zz)

    D'Agostino Skewness Normality Test

data:  input data
skewness = -0.035209, null hypothesis skewness = 0, p-value = 0.286
alternative hypothesis: true skewness is not equal to 0
95 percent confidence interval:
 -0.09992690  0.02950877
sample estimates:
   skewness           z      se.est     root.b1 
-0.03520907 -1.06683621  0.03301991 -0.03519946 

> kurtosis.test(zz)

    D'Agostino Kurtosis Normality Test

data:  input data
kurtosis = -0.052102, null hypothesis kurtosis = 0, p-value = 0.4362
alternative hypothesis: true kurtosis is not equal to 0
95 percent confidence interval:
 -0.18151406  0.07731029
sample estimates:
   kurtosis           z      se.est          b2 
-0.05210189 -0.77868046  0.06602783  2.94685476 

> 
> # Log normal distribution - rejection on both skewness and kurtosis
> zz1 <- exp(zz1)
> skewness.test(zz1)

    D'Agostino Skewness Normality Test

data:  input data
skewness = 5.2214, null hypothesis skewness = 0, p-value < 2.2e-16
alternative hypothesis: true skewness is not equal to 0
95 percent confidence interval:
 5.156675 5.286111
sample estimates:
   skewness           z      se.est     root.b1 
 5.22139319 63.31231869  0.03301991  5.21996907 

> kurtosis.test(zz1)

    D'Agostino Kurtosis Normality Test

data:  input data
kurtosis = 61.259, null hypothesis kurtosis = 0, p-value < 2.2e-16
alternative hypothesis: true kurtosis is not equal to 0
95 percent confidence interval:
 61.13006 61.38888
sample estimates:
   kurtosis           z      se.est          b2 
61.25946799 44.06817706  0.06602783 64.20270103 
$\endgroup$
  • $\begingroup$ There are good reasons to regard skewness and kurtosis as (sometimes) helpful measures of distribution shape but a poor basis for testing. That’s why the Shapiro-Wilk test and some others don’t use them. Specifically even if the parent is normal, sample skewness and kurtosis approach their asymptotic sampling distributions extraordinarily slowly. I don’t recall whether the D’Agostino test is smart about this, but many skewness and kurtosis tests are not. $\endgroup$ – Nick Cox Jan 26 at 10:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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