是否可以使用 shapiro.test 通过将测试应用于子样本来测试大于 5,000 个数据点的样本的正态性?

机器算法验证 正态分布 采样 夏皮罗威尔克测试
2022-03-14 20:27:55

我有一些我想测试正常性的样本。其中一个样本超过 5,000 个数据点,这是 shapiro 测试接受样本的上限。这是数据:

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)

为了规避 5,000 的限制,如果我只对较小的数据子样本进行测试,在统计上是否可以接受?例如,在这里,我为所有三个变量绘制了一个大小为 500 的子样本:

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

并且测试返回了所有三个显着的结果:

$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
4个回答

我有五个层面的意见。

  1. 根据这个证据,这是对特定 R 函数的限制, shapiro.test()并不意味着在 R 中没有其他方法可以做到这一点,对此我无法具体建议。对于所有软件都没有此类限制,这可能与您实际相关,也可能不相关。例如,Stata 命令swilk并没有以这种方式受到限制,但是手册和命令输出警告说,对于大约 5000 以上的样本量,P 值计算不能得到太多信任。(编辑:本段于 2021 年 1 月 26 日编辑根据@Ben Bolker 的评论和单独的答案。)

  2. 我无法评论为什么该特定功能无法执行,但更大的问题是您为什么要进行这种测试。不关心的一个很好的理由是通用的:对于该顺序甚至更大的样本量,这样的测试可以说是相当无用的,因为即使是与正态性的微小偏差也将在传统水平上被视为显着。更具体地说:为什么测试正态性很重要或很有趣?人们经常将此类检验应用于边际分布,因为人们普遍认为边际正态性是许多程序的要求。在正态性是相关假设或理想条件的情况下,它通常适用于以平均结果或响应结构为条件的分布。

  3. 针对您对二次抽样是否可以接受的具体疑问,严肃的回复在什么意义上可以接受?个人回复:作为统计论文的读者、作者和审稿人,以及统计期刊的编辑,我的反应是建议这样的二次抽样充其量是尴尬的,最坏的情况是回避主要问题,即发现没有这种限制的实现,或者更有可能以不同的方式考虑分布。

  4. 正如在 CV 和其他地方经常强调的那样,检查偏离正态性的最有用和信息量最大的方法是正态分位数图,通常也称为正态概率图、正态分数图或概率图。这样的图不仅提供了对非正态性程度的视觉评估,而且在何种意义上精确地说明了与理想形状的偏差。缺乏相关的 P 值实际上并没有太大的损失,尽管可以通过置信水平、模拟等为过程提供一些推论动力。(编辑 2021 年 1 月 26 日:还有其他术语是高斯百分位图和高斯概率图。)

  5. 具体来说,您的示例包括生成对数正态样本,然后确定它们确实不符合P 值的正常条件。这似乎令人费解,但请放心,对于更大的样本,您的 P 值将或应该更小,受制于机器级别的最小可报告 P 值问题。相反,您的真正问题很可能出在其他地方,而这些示例只不过是偶然的说明。1015

一个小的历史记录/更正:与此处和其他地方的其他答案中所说的(或可能从中推断出的)相反,R 的 Shapiro-Wilk 测试对 <=5000 观察的限制不是

  • R 实现中的一个意外限制
  • 为保护用户免于执行有问题的测试而施加的故意限制(如此处可能建议的那样

出现限制是因为 R 拒绝提供p-未验证原始函数的范围的值。相比之下,Royston (1995) 的原始实现和 Stata 的swilk函数确实提供了 p 值,但给出了错误代码/警告,说明p-值可能不可靠。

p- 对应于给定的值W统计数据很难计算:文献中有一系列论文(参见下面的参考文献)使用复杂的数学技术来提出在给定范围内计算有效且足够准确的近似值n提供对 Shapiro-Wilk 统计量 p 值的可靠估计。罗伊斯顿 (1995) 说:

所有计算均针对大于 5000 的样本进行,但IFAULT返回为 2。虽然W将被正确计算,其准确性P-价值无法保证。

换句话说,这超出了 Royston 和其他作者精心构建的有效函数的范围,这些函数可以很好地近似于p- 对应于给定值的值W.

我怀疑夏皮罗威尔克的实现p现代统计软件包中的值都基于 Royston (1995) 中描述的 Fortran 代码。如果你想计算可靠的 Shapiro-Wilkp- 样本的值n>5000(忽略此处和其他地方给出的所有关于为什么在非常大的数据集上进行正态性测试通常只是愚蠢的建议),您将不得不回到 Royston 1992 和 Verrill 和 Johnson 1988 的论文并重新执行/扩展这些方法以较大的值n- 不是胆小的人的项目。


  • 罗伊斯顿,帕特里克。“逼近非正态性的 Shapiro-Wilk W 检验。” 统计与计算 2,没有。3(1992 年 9 月 1 日):117–19。https://doi.org/10.1007/BF01891203
  • ———。“备注 AS R94:关于算法 AS 181 的备注:正态性 W 检验。” 皇家统计学会杂志。系列 C(应用统计)44,没有。4 (1995): 547–51。https://doi.org/10.2307/2986146
  • 维里尔、史蒂夫和理查德 A. 约翰逊。“用于检验正态性的删失数据相关统计的表和大样本分布理论。” 美国统计协会杂志 83,没有。404(1988):1192-97。https://doi.org/10.2307/2290156

我认为尼克考克斯指出了这种方法的一些困难。

一个可能的替代建议是使用另一个正态性检验。在我参加的课程中,我们对更大的样本使用了基于 D'Agostino 的偏度和峰度的测试。我在我的 lolcat 统计包中实现了这些测试。考虑:

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

有趣的是,即使样本大小为 5500,偏度/峰度也可能不会拒绝这些测试。对数正态分布很可能会拒绝,即使样本量大大降低。举个例子:

> 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 

我为恢复旧线程而道歉,但我在搜索中遇到了这个问题,我想添加我的输入以防其他人有同样的问题。尼克考克斯为这个问题提供了一些很好的意见,但我想从不同的角度提出一个答案。在诸如公司之类的标准化环境中,通常必须满足某些约束和标准,例如 p 值。虽然我绝对同意应该完成诸如正态图之类的图形评估,但这些类型的分析很难制定政策。建立量化约束在受监管的环境中很重要,而 p 值是对此的合理解决方案。

抽取一部分样本确实是一个合理的解决方案。可以这样想:您有一个制造系统,其中不断生产零件,比如说乐高积木。数百万(数十亿?)被创造出来,这就是人口。您选择一个桶来进行评估,该桶包含 100,000 个乐高积木。这远远超过你的需要!你拿起一个勺子,拉出大约 100 个乐高积木。看看你做了什么?你拿了一个样品,桶,然后你重新取样,勺子。现在在这个例子中,我在随机化组件方面做得很糟糕,但我认为它仍然很好地说明了我们每天都在做这种事情,无论是使用物理组件还是对数据列表进行采样。

总而言之,对您的样本进行抽样绝对是可以接受的,但是您要确保获得足够的随机样本,不仅代表较大的样本数据点,而且代表总体上的总体。您使用 500 个数据点的示例仍然是需要分析的大量样本。