选择样本大小

此示例说明如何确定执行统计测试所需的样本数或观察数。它说明了一个简单问题的样本量计算，然后展示了如何使用 sampsizepwr 函数为两个更实际的问题计算功效和样本量。最后，它说明了如何使用统计库函数来计算 sampsizepwr 函数不支持的测试所需的样本量。

单边测试具有已知标准差的正态平均值

只是为了介绍一些概念，让我们考虑一个不切实际的简单示例，我们想要测试均值并且我们知道标准差。我们的数据是连续的，可以用正态分布建模。我们需要确定样本量 N，以便我们可以区分平均值 100 和平均值 110。我们知道标准差为 20。

当我们进行统计检验时，我们通常会针对备择假设检验原假设。我们找到一个检验统计量 T，并查看其在零假设下的分布。如果我们观察到一个不寻常的值，比如如果零假设为真，则发生概率小于 5% 的值，那么我们将拒绝零假设而支持替代假设。 （5% 的概率称为检验的显着性水平。）如果该值没有异常，我们不会拒绝原假设。

在这种情况下，检验统计量 T 是样本均值。在原假设下，它的均值为 100，标准差为 20/sqrt(N)。首先让我们看一个固定的样本量，假设 N=16。如果 T 在阴影区域，即其分布的上尾，我们将拒绝零假设。这使得这是一个单边测试，因为如果 T 位于下尾，我们将不会拒绝。此阴影区域的截止值比平均值高 1.6 个标准差。

using TyPlot
using TyMath
using Printf
using TyStatistics
figure()
mu0 = 100
sig = 20
N = 16
alpha = 0.05
conf = 1-alpha
cutoff = norminv(conf, mu0, sig/sqrt(N))
x = [LinRange(90,cutoff,100); LinRange(cutoff,127,100)]
y = normpdf.(x,mu0,sig/sqrt(N))
h1 = plot(x,y)
xhi = [cutoff; x[x.>=cutoff]]
yhi = [0; y[x.>=cutoff]]
patch(xhi,yhi,"b")
title("Distribution of sample mean, N=16")
xlabel("Sample mean")
ylabel("Density")
text(96, 0.01, @sprintf("Reject if mean>%.4g\nProb = 0.05", cutoff), color="b")

这就是 T 在原假设下的表现，但在替代假设下呢？我们的替代分布的平均值为 110，如红色曲线所示。

hold("on")
mu1 = 110
y2 = normpdf.(x,mu1,sig/sqrt(N))
h2 = line(x,y2,color="r")
yhi = [0; y2[x.>=cutoff]]
patch(xhi,yhi,"r",facealpha=0.25)
P = 1 .- normcdf.(cutoff,mu1,sig/sqrt(N))
text(115,.06,@sprintf("Reject if mean>%.4g\nProb = %.2g", cutoff,P),color=[1,0,0])
legend(["Null hypothesis","Alternative hypothesis"])
hold("off")

如果备选方案为真，则更有可能拒绝原假设。这正是我们想要的。如果我们查看累积分布函数 (cdf) 而不是密度 (pdf)，则更容易可视化。我们可以直接从这张图中读取概率，而不必计算面积。

figure()
ynull = normcdf(x,mu0,sig/sqrt(N))
yalt = normcdf(x,mu1,sig/sqrt(N))
h12 = plot(x,ynull,"b-",x,yalt,"r-")
hold("on")
zval = norminv(conf)
cutoff = mu0 + zval * sig / sqrt(N)
line([90, cutoff, cutoff], [conf, conf, 0], linestyle = ":")
msg = @sprintf("Cutoff = 100 + %.2g * 20 / √n",zval)
text(cutoff,.15,msg,color="b")
text(minimum(x),conf,@sprintf("   %g%% test",100*alpha),color = "b",verticalalignment="top")
palt = normcdf(cutoff,mu1,sig/sqrt(N))
line([90, cutoff, cutoff], [conf, conf, 0], linestyle=":")
line([90,cutoff],[palt,palt],color="r",linestyle=":")
text(91,palt+.02,@sprintf(" Power is 1-%.2g = %.2g",palt,1-palt),color=[1,0,0])
legend(["Null hypothesis","Alternative hypothesis"])
hold("off")

该图显示了当 N=16 时，对于两个不同的 mu 值获得显着统计量（拒绝原假设）的概率。

幂函数定义为当备选方案为真时拒绝零假设的概率。这取决于备选方案的价值和样本量。我们将绘制作为 N 函数的功效（即 1 减去 cdf），将备选方案固定为 110。我们将选择 N 以获得 80% 的功效。该图显示我们需要大约 N=25。

figure()
DesiredPower = 0.80
Nvec = 1:30
cutoff = @. mu0 + norminv(conf)*sig/sqrt(Nvec)
power = @. 1 - normcdf(cutoff, mu1, sig/sqrt(Nvec))
plot(Nvec,power,"bo-",[0,30],[DesiredPower,DesiredPower],"k:")
xlabel("N = sample size")
ylabel("Power")
title("Power function for the alternative:  μ = 110")

在这个过于简单的示例中，有一个公式可以直接计算所需的值以获得 80% 的幂：

mudiff = (mu1 - mu0) / sig
N80 = ceil(((norminv(1-DesiredPower)-norminv(conf)) / mudiff)^2)

N80 = 25.0

为了验证这是否有效，让我们进行蒙特卡洛模拟并在均值为 100 的原假设和均值为 110 的备择假设下生成 400 个大小为 25 的样本。如果我们测试每个样本以查看其均值是否为 100，我们应该预计第一组的大约 5% 和第二组的 80% 是重要的。

figure()
rng = MT19937ar(5489)
nsamples = 400
samplenum = 1:nsamples
N = 25
h0 = Vector{Float64}(undef,nsamples)
h1 = similar(h0)
for j = 1:nsamples
    Z0 = normrnd(rng,mu0,sig,N)
    h0[j] = ztest(Z0,mu0,sig,alpha=alpha,tail="right")[1]
    Z1 = normrnd(rng,mu1,sig,N)
    h1[j] = ztest(Z1,mu0,sig,alpha=alpha,tail="right")[1]
end
p0 = cumsum(h0) ./ samplenum
p1 = cumsum(h1) ./ samplenum
plot(samplenum,p0,"b-",samplenum,p1,"r-")
xlabel("Sample number")
ylabel("Proportion significant")
title("Verification of power computatio")
legend(["Null hypothesis","Alternative hypothesis"],loc="east")

双侧测试标准差未知的正态平均值

现在假设我们不知道标准偏差，并且我们想要进行双侧检验，即拒绝零假设样本均值过高或过低的检验。

检验统计量是t统计量，它是样本均值与被检验均值之差除以均值的标准误差。在原假设下，这具有自由度为 N-1 的学生 t 分布。在备择假设下，分布是非中心 t 分布，其非中心参数等于真实均值与被检验均值之间的归一化差值。

对于这个双侧检验，我们必须将原假设下 5% 的错误概率平均分配给两条尾部，如果检验统计量在任一方向过于极端，则拒绝。我们还必须考虑任何替代方案下的两条尾巴。

figure()
N = 16
df = N-1
alpha = 0.05
conf = 1-alpha
cutoff1 = tinv(alpha/2,df)
cutoff2 = tinv(1-alpha/2,df)
x = [LinRange(-5,cutoff1,100); LinRange(cutoff1,cutoff2,100);LinRange(cutoff2,5,100)]
y = tpdf(x,df)
h1 = plot(x,y)
xlo = [x[x.<=cutoff1];cutoff1]
ylo = [y[x.<=cutoff1];0]
xhi = [cutoff2;x[x.>=cutoff2]]
yhi = [0; y[x.>=cutoff2]]
patch(xlo,ylo,"b")
patch(xhi,yhi,"b")
title("Distribution of t statistic, N=16")
xlabel("t")
ylabel("Density")
text(2.5,.05,@sprintf("Reject if t>%.4g\nProb = 0.025",cutoff2),color="b")
text(-4.5,.05,@sprintf("Reject if t<%.4g\nProb = 0.025",cutoff1),color="b")

我们可以将其视为 mu 的函数，而不是仅在原假设和 mu 的单个替代值下检查幂函数。随着 mu 在任一方向上远离原假设值，功效会增加。我们可以使用 sampsizepwr 函数来计算功率。对于功效计算，我们需要为标准差指定一个值，我们怀疑该值大约为 20。这是样本大小 N=16 的功效函数的图片。

figure()
N = 16
x = LinRange(90,127,100)
power = sampsizepwr("t",[100,20],x,nothing,N)
plot(x,power)
xlabel("True mean")
ylabel("Power")
title("Power function for N=16")

当平均值为 110 时，我们想要 80% 的功效。根据此图，我们的功效小于 50%，样本量为 N=16。多大的样本量可以提供我们想要的功效？

N = sampsizepwr("t",[100,20],110,0.8)

N = 34

我们需要大约 34 的样本量。与前面的示例相比，我们需要进行九次额外的观察以允许进行双侧测试并弥补不知道真实标准差的情况。

我们可以绘制不同 N 值的幂函数图。

figure()
Nvec = 2:40;
power = sampsizepwr("t",[100,20],110,nothing,Nvec)
plot(Nvec,power,"bo-",[0 40],[DesiredPower DesiredPower],"k:")
xlabel("N = sample size")
ylabel("Power")
title("Power function for the alternative: μ = 110")

我们可以进行类似于之前的模拟，以验证我们是否获得了所需的功率。

figure()
nsamples = 400
samplenum = 1:nsamples
N = 34
h0 = Vector{Float64}(undef,nsamples)
h1 = similar(h0)
for j = 1:nsamples
    Z0 = normrnd(rng,mu0,sig,N)
    h0[j] = ttest(Z0,mu0,alpha=alpha)[1]
    Z1 = normrnd(rng,mu1,sig,N)
    h1[j] = ttest(Z1,mu0,alpha=alpha)[1]
end
p0 = cumsum(h0) ./ samplenum
p1 = cumsum(h1) ./ samplenum
plot(samplenum,p0,"b-",samplenum,p1,"r-")
xlabel("Sample number")
ylabel("Proportion significant")
title("Verification of power computation")
legend(["Null hypothesis","Alternative hypothesis"],loc="east")

假设我们可以负担得起 50 的样本量。大概我们检测替代值 mu=110 的能力将大于 80%。如果我们将功率保持在 80%，我们可以检测到什么替代方案？

mu1 = sampsizepwr("t",[100,20],nothing,0.8,50)

mu1 = 108.08366004271574

测试比例

现在让我们转向确定区分两个比例所需的样本量的问题。假设我们正在对一个人口进行抽样，其中大约 30% 的人支持某个候选人，并且我们希望对足够多的人进行抽样，以便我们可以将这个值与 33% 区分开来。

这里的思路和之前一样。这里我们可以使用样本计数作为我们的测试统计量。此计数服从二项分布。对于任何样本大小 N，我们可以计算拒绝零假设 P=0.30 的截止值。例如，对于 N=100，如果样本计数大于按如下方式计算的截止值，我们将拒绝原假设：

N = 100
alpha = 0.05
p0 = 0.30
p1 = 0.33
cutoff = binoinv(1-alpha, N, p0)

cutoff = 38.0

二项式分布的复杂性在于它是离散的。超过截止值的概率不正好是 5%：

1 - binocdf(cutoff, N, p0)

0.03397899826893125

再一次，让我们针对一系列样本量计算备选方案 P=0.33 的功效。我们将使用单侧（右尾）测试，因为我们只对大于 30% 的替代值感兴趣。

figure()
Nvec = 50:50:2000
power = sampsizepwr("p",p0,p1,nothing,Nvec,tail="right")
plot(Nvec,power,"bo-",[0 2000],[DesiredPower DesiredPower],"k:")
xlabel("N = sample size")
ylabel("Power")
title("Power function for the alternative: p = 0.33")

我们可以使用 sampsizepwr 函数来请求 80% 的功效所需的样本量。

approxN = sampsizepwr("p",p0,p1,0.80,nothing,tail="right")

┌ Warning: 值 N>=200 是近似值。将幂绘制为 N 的函数可以发现具有所需幂
的更小的 N 值。
└ @ TyStatistics 
1500

一条警告消息告诉我们答案只是近似值。如果我们查看不同样本量的幂函数，我们可以看到该函数通常是递增的，但不规则，因为二项分布是离散的。让我们看看在 1470 到 1480 的样本量范围内拒绝 p=0.30 和 p=0.33 的原假设的概率。

figure()
subplot(3,1,1)
Nvec = 1470:1480
power = sampsizepwr("p",p0,p1,nothing,Nvec,tail="right")
plot(Nvec,power,"ro-",[minimum(Nvec),maximum(Nvec)],[DesiredPower,DesiredPower],"k:")
ylabel(@sprintf("Prob[T>cutoff]\nif p=0.33"))
h_gca = gca()
h_gca.set_xticklabels("")
ylim([.78, .82])

subplot(3,1,2)
alf = sampsizepwr("p",p0,p0,nothing,Nvec,tail="right")
plot(Nvec,alf,"bo-",[minimum(Nvec),maximum(Nvec)],[alpha,alpha],"k:")
ylabel(@sprintf("Prob[T>cutoff]\nif p=0.30"))
h_gca = gca()
h_gca.set_xticklabels("")
ylim([.04, .06])

subplot(3,1,3)
cutoff = binoinv.(1-alpha, Nvec, p0)
plot(Nvec,cutoff,"go-")
xlabel("N = sample size")
ylabel("Cutoff")

该图显示幂函数曲线（顶部图）不仅不规则，而且在某些样本量下还会减少。这些样本大小需要增加截止值（底部图）以保持显着性水平（中间图）不大于 5%。我们可以在此范围内找到较小的样本量，以提供 80% 的所需功效：

minimum(Nvec[power.>=0.80])

1478

测试相关性

在我们目前考虑的示例中，我们能够计算出检验统计量达到特定显着性水平的截止值，并计算在替代假设下超过该截止值的概率。对于我们的最后一个例子，让我们考虑一个不太容易的问题。

想象一下，我们可以从两个变量 X 和 Y 中抽取样本，并且我们想知道我们需要多少样本量来测试它们是否不相关，而不是它们的相关性高达 0.4 的备选方案。虽然可以通过将样本相关性转换为 t 分布来计算出样本相关性的分布，但让我们使用一种即使在无法计算出检验统计量分布的问题中也可以使用的方法。

对于给定的样本量，我们可以使用蒙特卡罗模拟来确定相关性检验的近似截止值。让我们进行一次大型模拟运行，这样我们就可以准确地得到这个值。我们将从 25 的样本量开始。

nsamples = 10000
N = 25
alpha = 0.05
conf = 1-alpha
r = zeros(nsamples)
for j = 1:nsamples
    xy = normrnd(rng,0,1,N,2)
    r[j:j], = corr(xy[:,1],xy[:,2])
end
cutoff = ty_quantile(r,conf)

0.33718635817903886

然后我们可以在备择假设下生成样本，并估计检验的功效。

nsamples = 1000
mu = [0; 0]
sig = [1 0.4; 0.4 1]
r = zeros(nsamples)
for j = 1:nsamples
    xy, = mvnrnd(rng,mu,sig,N)
    r[j:j], = corr(xy[:,1],xy[:,2])
end
power,powerci = binofit(count(>(cutoff),r),nsamples)

power = 0.647

powerci = 1×2 Matrix{Float64}:
 0.616481  0.67665

我们估计功效为 65%，并且我们有 95% 的置信度认为真实值介于 62% 和 68% 之间。要获得 80% 的功效，我们需要更大的样本量。我们可能会尝试将 N 增加到 50，估计此样本大小的截止值，然后重复功效模拟。

nsamples = 10000
N = 50
alpha = 0.05
conf = 1-alpha
r = zeros(nsamples)
for j = 1:nsamples
    xy = normrnd(rng,0,1,N,2)
    r[j:j], = corr(xy[:,1],xy[:,2])
end
cutoff = ty_quantile(r,conf)

cutoff = 0.2314634086611933

nsamples = 1000;
mu = [0; 0];
sig = [1 0.4; 0.4 1];
r = zeros(nsamples);
for j = 1:nsamples
    xy, = mvnrnd(rng,mu,sig,N)
    r[j:j], = corr(xy[:,1],xy[:,2])
end
power, powerci = binofit(count(>(cutoff), r), nsamples)

power = 0.899

powerci = 1×2 Matrix{Float64}:
 0.878633  0.91698

此样本量提供的功效优于我们 80% 的目标。我们可以继续以这种方式进行试验，尝试找到满足我们要求的小于 50 的样本量。

结论

统计工具箱中的概率函数可用于确定在假设检验中达到所需功效水平所需的样本量。在某些问题中，样本量可以直接计算；在其他情况下，有必要搜索一系列样本大小，直到找到正确的值。随机数生成器可以帮助验证是否满足所需功效，也可用于研究特定测试在替代条件下的功效。