分类:P值和样本大小的关系

来自Big Physics


问题背景

在有实验组和对照组的对比实验设计中,往往我们需要通过预实验来确定样本大小。而样本大小的确定通常可以通过分析预实验的数据,通过计算显著性和样本大小的关系,也就是p值和N的关系,画图,做延拓,来估计所需要的最小N值。

但是,是否存在着这样的情况,无论实验组和对照组的差别多么小,只要通过增加样本大小N,总是可以得到统计上显著的结果?如果有,那这个实验设计的思路,原则上有逻辑漏洞。

理论研究

假定我们的两组数据[math]\displaystyle{ \left(\left\{x^{a}_{i=1...N_{a}}\right\},\left\{x^{b}_{i=1...N_{b}}\right\}\right) }[/math]分别来自于两个正态分布[math]\displaystyle{ N(\mu^{a}, \sigma^{a}) }[/math][math]\displaystyle{ N(\mu^{b},\sigma^{b}) }[/math]。现在我们来从这两组数据估计出来这两个分布函数有差别的概率。其中一种定义有差别的方式是:从两个分布中随机各自取出来一个新的样本[math]\displaystyle{ \left(x^{a}_{*},x^{b}_{*}\right) }[/math],其中一个大于另一个的概率,[math]\displaystyle{ P\left(x^{a}_{*}\gt x^{b}_{*}\right)=P\left(x^{a}_{*}-x^{b}_{*}\gt 0\right)\gt p_{*} }[/math]

我们先在参数[math]\displaystyle{ \mu^{a}, \sigma^{a}, \mu^{b},\sigma^{b} }[/math]全都已知的条件下把这个公式推导出来,然后,再来解决如何从样本数据[math]\displaystyle{ \left(\left\{x^{a}_{i=1...N_{a}}\right\},\left\{x^{b}_{i=1...N_{b}}\right\}\right) }[/math]中估计出来这些参数的问题,以及估计出来以后如何更好地给出来[math]\displaystyle{ P\left(x^{a}_{*}-x^{b}_{*}\gt 0\right) }[/math]的问题。注意,后者由于估计出来的量[math]\displaystyle{ \mu^{a}, \sigma^{a}, \mu^{b},\sigma^{b} }[/math]都有误差,因此,会比直接代入以[math]\displaystyle{ \mu^{a}, \sigma^{a}, \mu^{b},\sigma^{b} }[/math]表达的公式要复杂一些。

概率论问题:从[math]\displaystyle{ \mu^{a}, \sigma^{a}, \mu^{b},\sigma^{b} }[/math]得到[math]\displaystyle{ P\left(x^{a}_{*}-x^{b}_{*}\gt 0\right) }[/math]

[math]\displaystyle{ X=x^{a}-x^{b} }[/math]的分布函数是[math]\displaystyle{ N(\mu^{a}-\mu^{b},\sqrt{\left(\sigma^{a}\right)^{2}+\left(\sigma^{b}\right)^{2}}) }[/math]。也就是说,只要这四个参数满足下面的条件,两个分布函数就有区别,

[math]\displaystyle{ \int_{0}^{\infty} d\xi \frac{1}{\sqrt{2\pi}}e^{-\frac{\left(\xi-\left(\mu^{a}-\mu^{b}\right)\right)^{2}}{2\left(\sigma^{a}\right)^{2}+2\left(\sigma^{b}\right)^{2}}}\geq p_{*} }[/math].

这个概率实际上,由两个分布函数的均值和方差四个参数一起决定的。就算两个均值相差很大,但是方差也比较大的时候,两个分布函数重叠的地方特别多,则仍然这个概率会很小,不一定能够达到大于阈值[math]\displaystyle{ p_{*} }[/math]。用通常的p值的语言,可以重新表述为

[math]\displaystyle{ p=1-\int_{0}^{\infty} d\xi \frac{1}{\sqrt{2\pi}}e^{-\frac{\left(\xi-\left(\mu^{a}-\mu^{b}\right)\right)^{2}}{2\left(\sigma^{a}\right)^{2}+2\left(\sigma^{b}\right)^{2}}}=\int_{-\infty}^{0} d\xi \frac{1}{\sqrt{2\pi}}e^{-\frac{\left(\xi-\left(\mu^{a}-\mu^{b}\right)\right)^{2}}{2\left(\sigma^{a}\right)^{2}+2\left(\sigma^{b}\right)^{2}}} }[/math].

也可以类似地讨论[math]\displaystyle{ X=\left|x^{a}-x^{b}\right| }[/math]的分布函数,然后把统计检验构建在这个分布函数之上,不过这里暂时不展开讨论。

Monte Carlo问题:从[math]\displaystyle{ \mu^{a}, \sigma^{a}, \mu^{b},\sigma^{b} }[/math]到样本数据

假设我们已知一个正态分布[math]\displaystyle{ N\left(\mu, \sigma\right) }[/math]我们来看看从中得到的大小为[math]\displaystyle{ N }[/math]的样本数据[math]\displaystyle{ \left\{x_{i=1...N}\right\} }[/math]的统计性质。

首先,

[math]\displaystyle{ \left\langle\bar{x}\right\rangle=\frac{1}{N}\sum_{i=1}^{N}\left\langle x_{i} \right\rangle=\mu }[/math]

也就是样本均值的均值等于原始分布函数的均值。注意,这里的[math]\displaystyle{ \left\langle \cdot \right\rangle }[/math]是系综平均,也就是如果生成这样的样本很多很多次例如L次(在L趋向无穷大的极限下),每次N个,然后对每一次得到的结果先求出来单次的量,接着做多次平均。

其次,

[math]\displaystyle{ \left\langle s^{2}_p \right\rangle=\left\langle \frac{1}{N}\sum_{i=1}^{N}\left(x_{i}-\bar{x}\right)^2 \right\rangle=\sigma^{2} }[/math]

也就是样本方差的均值等于原始分布函数的方差。

最后,我们计算一下样本均值的方差,也就是

[math]\displaystyle{ \sigma_{N}=\left\langle\bar{x}^2-\left\langle\bar{x}\right\rangle^2\right\rangle=\frac{1}{N}\sigma^2 }[/math]

也就是样本均值的方差是样本方差的[math]\displaystyle{ \frac{1}{N}\sigma^2 }[/math]

经过这个计算,我们一定要看到样本方差和样本均值的方差之间的差异和联系。前者是对整个样本来说的,后者是把每一次抽样以后的均值当作随机变量,考察这个随机变量的涨落得到的。这个涨落,在样本大小N越大的时候,越小。

统计学问题:从样本数据到[math]\displaystyle{ \mu^{a}, \sigma^{a}, \mu^{b},\sigma^{b} }[/math]

从样本数据[math]\displaystyle{ \left(\left\{x^{a}_{i=1...N_{a}}\right\},\left\{x^{b}_{i=1...N_{b}}\right\}\right) }[/math],通过极大似然估计(或者别的估计方法),我们可以得到参数值[math]\displaystyle{ \mu^{a}, \sigma^{a}, \mu^{b},\sigma^{b} }[/math]以及每个参数的估计误差,例如(先暂时写成一个[math]\displaystyle{ \sigma }[/math]的形式,如果需要转化成95%置信区间,则需要做相应转换),

[math]\displaystyle{ \mu^{a}=\bar{x}^{a}\pm \sigma_{N_{a}} = \bar{x}^{a}\pm \frac{1}{\sqrt{N_{a}}}s^{a}_{p} }[/math]

其中[math]\displaystyle{ \bar{x}^{a}=\frac{1}{N_{a}}\sum_{i=1}^{N_{a}}x^{a}_{i} }[/math][math]\displaystyle{ s^{a,2}_p=\frac{1}{N_{a}}\sum_{i=1}^{N_{a}}\left(x^{a}_{i}-\bar{x}\right)^2 }[/math],也被称为样本方差。注意这里的几个公式和上一节的系综平均公式很像,但是,它们意义很不相同。在这里,我们不需要引入系综平均。除非我们换一个思路,用Bootstrap方法来解决这个误差估计的问题,则需要从实际样本中产生多次抽样数据,于是刚好回到系综平均。注意,这里有个非常微妙的地方,参数[math]\displaystyle{ \mu^{a} }[/math]的估计误差不是直接就是样本方差[math]\displaystyle{ s^{a,2}_p }[/math],而是[math]\displaystyle{ \sigma_{N_{a}}=\frac{1}{\sqrt{N_{a}}}s^{a,2}_p }[/math],其额外还增加了一个[math]\displaystyle{ \frac{1}{\sqrt{N_{a}}} }[/math]的系数。这也是合理的,样本数量越多,对参数的估计就更准确。

对于参数[math]\displaystyle{ \sigma^{a} }[/math]我们有估计方式

[math]\displaystyle{ \sigma^{a,2}=\frac{1}{N_{a}}\sum_{i=1}^{N_{a}}\left(x^{a}_{i}-\bar{x}\right)^2=s^{a,2}_p }[/math].

这个参数的估计误差的公式就不再推导了,但是,类似可以得到。

对什么东西做假设检验?对比两个均值,还是两个分布

有了上面的理解,现在我们可以来回答统计检验的问题了。核心问题是,在一个统计检验里面,我们是希望对比两组样本背后的分布函数的均值,还是说对比两组样本背后的分布函数?

如果我们对比均值,则相当于我们在比较,[math]\displaystyle{ \mu^{a} }[/math][math]\displaystyle{ \mu^{b} }[/math]。于是,大概我们需要计算下面的量(假设a样本均值更大)然后看看这个量大于零的概率,

[math]\displaystyle{ \mu^{a}-\mu^{b}=\left(\bar{x}^{a}- \frac{1}{\sqrt{N_{a}}}s^{a}_{p}\right)-\left(\bar{x}^{b}+ \frac{1}{\sqrt{N_{b}}}s^{b}_{p}\right) }[/math],

[math]\displaystyle{ =\left(\bar{x}^{a}- \bar{x}^{b} \right)-\left(\frac{1}{\sqrt{N_{a}}}s^{a}_{p}+ \frac{1}{\sqrt{N_{b}}}s^{b}_{p}\right) }[/math].

于是,第一项和第二项的比例,决定了这个量大于零的概率,也就是说我们需要计算

[math]\displaystyle{ t=\frac{\left(\bar{x}^{a}- \bar{x}^{b} \right)}{\left(\frac{1}{\sqrt{N_{a}}}s^{a}_{p}+ \frac{1}{\sqrt{N_{b}}}s^{b}_{p}\right)} }[/math].

这正好就是传统的t检验的内容(除了样本方差计算的时候需要修正一下自由度这些细节)。

如果我们对比的是两个分布函数,则其实我们需要同时考虑[math]\displaystyle{ \mu^{a}, \sigma^{a}, \mu^{b},\sigma^{b} }[/math]。也就是把这四个参数都估计出来以后代入[math]\displaystyle{ P\left(x^{a}_{*}-x^{b}_{*}\gt 0\right) }[/math]的公式,计算出来分布函数有差异的概率。而由于[math]\displaystyle{ \sigma^{a},\sigma^{b} }[/math]不随着样本数量的改变而减小,因此,不存在任何不等于零的差别只要增加样本都会显著这回事。

和传统统计检验的联系

KS检验似乎是对比的分布函数,好像不会出现随着。ks检验实际上对比的是两个分布函数之间的差。因此,从逻辑上说,ks检验先从实际数据中拟合出来分布函数,然后,计算两个分布函数的最大区别。因此,增加样本量只能使得分布函数估计的更加准确,但是不会增加那个最大区别。这个结论需要做数值检验。

当然,也有可能,其实KS检验,也经过了人工调整,使得它其实检验的是分布函数的均值。因此,这需要数值实验来检验。

从正态分布的双样本t检验的角度来看,由于其检验目标就是两个分布的均值是否有区别,因此,统计显著性会随着样本数量的增加而提高。这个结论需要做数值检验。

理论分析得到的解决方式

明确指出来,到底是对比的两个分布函数的均值,还是对比的两个分布函数。然后,如果检验的是均值的差异,用传统的统计检验,并且必然存在着样本量和统计显著性的平庸关系。如果检验的是分布函数的差异,则我们提出新(是不是新的,需要做文献调研和请教统计学专家)的检验方法,用我们的方法,并且,显著性和样本量之间不存在着这个增加样本量必然导致显著性增加的平庸关系。

另外,对比均值的检验,如果看作是分布函数,则正好是把大小为[math]\displaystyle{ N }[/math]的样本当作一个整体看,得到的分布函数(也就是看大小为[math]\displaystyle{ N }[/math]的样本均值的分布函数)和对这个函数的统计检验结果[1]

数值实验

做一下数值实验,看看针对对比均值的统计检验怎么做(传统上就是),针对分布函数的统计检验怎么做(可能也已经有了,也可以按照上面的理论分析自己发明),结果是否根理论分析一致——对于均值的检验统计显著性随着样本增加而增加,对于分布函数检验统计显著性没有和样本直接相关(由于增加样本会使得分布函数估计更准确,因此,显著性也会和样本数量有关系,但是,不应该正好是[math]\displaystyle{ \frac{1}{\sqrt{N}} }[/math]的关系)。

使用MATLAB开展数值检验

文件:Figure PN.tiff

   clc; clear;
   ES = [0, 0.01, 0.2, 0.3]; % set 4 effect size
   for ef = 1:length(ES) % simulate for each effect size
       EffectSize = ES(ef);
       
       for n = 1:1000 % simulate starting from n = 1 to 10000
           G1 = rand(n,1) + EffectSize; % sample from the normal distribution and add effect size
           G2 = rand(n,1);
           [~,ttest2_p,~,~] = ttest2(G1,G2); % two sample t test
           [~,kstest2_p,~] = kstest2(G1,G2); % two sample ks test
           data0(n,:) = [n, ttest2_p,kstest2_p]; % collect data
       end
       data = log(data0); % log n and log ps 
       N = data(:,1); 
      ttest2_P = data(:,2);
       kstest2_P = data(:,3);
   
       % plot p-n for t test
       subplot(length(ES),2,1+2*(ef-1))
       plot(N, ttest2_P)
       title(['Effect size = ',num2str(EffectSize),'; Two sample t test'])
       ylabel('Significance (p)')
       xlabel('Sample size (n)')
       set(gca,'XTick',0:max(N)/4:max(N))
       set(gca,'XTicklabel',round(exp(0:max(N)/4:max(N))))
       Ytick = [min(ttest2_P(~isinf(ttest2_P))),log(0.05)];
       set(gca,'YTick',Ytick)
       set(gca,'YTicklabel',round(exp(Ytick),3))
       % plot p-n for ks test
       subplot(length(ES),2,2+2*(ef-1))
       plot(N, kstest2_P)
       title(['Effect size = ',num2str(EffectSize),'; Two sample ks test'])
       ylabel('Significance (p)')
       xlabel('Sample size (n)')
       set(gca,'XTick',0:max(N)/4:max(N))
       set(gca,'XTicklabel',round(exp(0:max(N)/4:max(N))))
       Ytick = [min(kstest2_P(~isinf(kstest2_P))),log(0.05)];
       set(gca,'YTick',Ytick)
       set(gca,'YTicklabel',round(exp(Ytick),3))
   end
   saveas(gcf,'Figure_PN.tif')

下一步工作

  1. 文献调研,看一下,这个问题其他人是否已经回答过,是否已经有好的答案。
  2. 数值实验
    1. 选择一个分布函数,例如正态分布,调整两者的均值和方差(可以暂时不调整),使得两者分别处于有显著区别和没显著区别的状态(例如之前的[math]\displaystyle{ P\left(x^{a}_{*}-x^{b}_{*}\gt 0\right)\gt 0.95 }[/math]
    2. 从每个分布函数中产生样本,给定样本数N
    3. 完成针对均值的检验
      1. 对这两组样本做t检验和KS统计检验,得到p值,看是否统计显著
      2. 画p-N曲线,看是否随着N增加,p减小
    4. 完成针对分布函数的检验
      1. 对这两组样本做t检验和KS统计检验,得到p值,看是否统计显著
      2. 画p-N曲线,看是否随着N增加,p减小
  3. 选择一组足够大的实际数据,有实验组和对照组的,分开具有显著不同的和没有显著不同的,做bootstrap抽样,大小为N
    1. 完成针对均值的检验
      1. 对这两组样本做t检验和KS统计检验,得到p值,看是否统计显著
      2. 画p-N曲线,看是否随着N增加,p减小
    2. 完成针对分布函数的检验
      1. 对这两组样本做t检验和KS统计检验,得到p值,看是否统计显著
      2. 画p-N曲线,看是否随着N增加,p减小

参考文献

  1. Zhesi Shen, Liying Yang, Zengru Di, Jinshan Wu. Large enough sample size to rank two groups of data reliably according to their means. Scientometrics 118: 653-671 (2019). https://doi.org/10.1007/s11192-018-2995-0

本分类目前不含有任何页面或媒体文件。