Тестирование нормального среднего значения с известным стандартным отклонением, одностороннее

Просто чтобы ввести некоторые понятия, давайте рассмотрим нереально простой пример, где мы хотим проверить среднее и мы знаем стандартное отклонение. Наши данные непрерывны и могут быть смоделированы с нормальным распределением. Нам нужно определить размер выборки N, чтобы мы могли различать среднее значение 100 и среднее значение 110. Мы знаем, что стандартное отклонение равно 20.

Когда мы проводим статистический тест, мы обычно проверяем нулевую гипотезу против альтернативной гипотезы. Мы находим тестовую статистику T и смотрим на ее распределение под нулевой гипотезой. Если мы наблюдаем необычное значение, скажем, которое имеет менее 5% шансов на то, чтобы это произошло, если нулевая гипотеза верна, то мы отвергаем нулевую гипотезу в пользу альтернативы. (Вероятность 5% известна как уровень значимости теста.) Если значение не является необычным, мы не отвергаем нулевую гипотезу.

В этом случае проверочная статистика Т представляет собой среднее значение выборки. При нулевой гипотезе он имеет среднее значение 100 и стандартное отклонение 20/sqrt (N). Сначала рассмотрим фиксированный размер выборки, скажем, N = 16. Мы отвергнем нулевую гипотезу, если Т находится в затененной области, которая является верхним хвостом его распределения. Это делает это односторонним тестом, так как мы не отвергнем, если Т находится в нижнем хвосте. Отсечение для этой затененной области составляет 1,6 стандартных отклонений от среднего значения.

rng(0,'twister');
mu0 = 100;
sig = 20;
N = 16;
alpha = 0.05;
conf = 1-alpha;
cutoff = norminv(conf, mu0, sig/sqrt(N));
x = [linspace(90,cutoff), linspace(cutoff,127)];
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,.01,sprintf('Reject if mean>%.4g\nProb = 0.05',cutoff),'Color','b');

Вот как Т ведет себя при нулевой гипотезе, а как же при альтернативе? Наше альтернативное распределение имеет среднее значение 110, как представлено красной кривой.

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 T>%.4g\nProb = %.2g',cutoff,P),'Color',[1 0 0]);
legend([h1 h2],'Null hypothesis','Alternative hypothesis');

Существует большая вероятность отклонения нулевой гипотезы, если альтернатива верна. Это то, чего мы хотим. Проще визуализировать, если посмотреть на кумулятивную функцию распределения (cdf) вместо плотности (pdf). Мы можем считывать вероятности непосредственно из этого графа, вместо того, чтобы вычислять площади.

ynull = normcdf(x,mu0,sig/sqrt(N));
yalt = normcdf(x,mu1,sig/sqrt(N));
h12 = plot(x,ynull,'b-',x,yalt,'r-');
zval = norminv(conf);
cutoff = mu0 + zval * sig / sqrt(N);
line([90,cutoff,cutoff],[conf, conf, 0],'LineStyle',':');
msg = sprintf(' Cutoff = 100 + %.2g \\times 20 / \\surd{n}',zval);
text(cutoff,.15,msg,'Color','b');
text(min(x),conf,sprintf('   %g%% test',100*alpha),'Color','b',...
                         'verticalalignment','top')
palt = normcdf(cutoff,mu1,sig/sqrt(N));
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(h12,'Null hypothesis','Alternative hypothesis');

Этот график показывает вероятность получения значимой статистики (отклонение нулевой гипотезы) для двух различных значений mu, когда N = 16.

Функция мощности определяется как вероятность отклонения нулевой гипотезы, когда альтернатива верна. Это зависит от значения альтернативы и от размера выборки. Мы построим график мощности (которая равна единице минус cdf) как функция N, фиксируя альтернативу на уровне 110. Мы выберем N, чтобы получить мощность 80%. График показывает, что нам нужно около N = 25.

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: \mu = 110')

В этом слишком простом примере существует формула для вычисления требуемого значения непосредственно для получения мощности 80%:

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

N80 =

    25

Чтобы убедиться, что это работает, давайте сделаем симуляцию Монте-Карло и создадим 400 выборок размера 25 как при нулевой гипотезе со средним значением 100, так и при альтернативной гипотезе со средним значением 110. Если мы проверяем каждый образец, чтобы увидеть, является ли его среднее значение 100, мы должны ожидать, что около 5% первой группы и 80% второй группы будут значимыми.

nsamples = 400;
samplenum = 1:nsamples;
N = 25;
h0 = zeros(1,nsamples);
h1 = h0;
for j = 1:nsamples
    Z0 = normrnd(mu0,sig,N,1);
    h0(j) = ztest(Z0,mu0,sig,alpha,'right');
    Z1 = normrnd(mu1,sig,N,1);
    h1(j) = ztest(Z1,mu0,sig,alpha,'right');
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','Location','East')

Тестирование нормального среднего значения с неизвестным стандартным отклонением, двусторонним

Теперь предположим, что мы не знаем стандартное отклонение, и мы хотим выполнить двусторонний тест, то есть тот, который отвергает нулевую гипотезу, является ли среднее значение выборки слишком высоким или слишком низким.

Статистика теста представляет собой статистику t, которая представляет собой разницу между средним и средним испытываемым, деленную на стандартную ошибку среднего. При нулевой гипотезе это имеет распределение Стьюдента с N-1 степенями свободы. Согласно альтернативной гипотезе, распределение является нецентральным t-распределением с параметром нецентральности, равным нормализованной разности между истинным средним и проверяемым средним.

Для этого двустороннего теста мы должны распределить 5% вероятность ошибки при нулевой гипотезе одинаково обоим хвостам и отклонить, если статистика теста слишком экстремальна в любом направлении. Мы также должны рассмотреть оба хвоста в рамках любой альтернативы.

N = 16;
df = N-1;
alpha = 0.05;
conf = 1-alpha;
cutoff1 = tinv(alpha/2,df);
cutoff2 = tinv(1-alpha/2,df);
x = [linspace(-5,cutoff1), linspace(cutoff1,cutoff2),linspace(cutoff2,5)];
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.

N = 16;
x = linspace(90,127);
power = sampsizepwr('t',[100 20],x,[],N);
plot(x,power);
xlabel('True mean')
ylabel('Power')
title('Power function for N=16')

Мы хотим мощность 80%, когда среднее значение 110. Согласно этому графику, наша мощность составляет менее 50% при размере выборки N = 16. Какой размер выборки даст необходимую нам мощность?

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

N =

    34

Нам нужен размер выборки около 34. По сравнению с предыдущим примером нам нужно взять девять дополнительных наблюдений, чтобы разрешить двусторонний тест и компенсировать незнание истинного стандартного отклонения.

Мы можем сделать график функции мощности для различных значений N.

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

И мы можем сделать симуляцию, аналогичную предыдущей, чтобы убедиться, что мы получаем необходимую нам мощность.

nsamples = 400;
samplenum = 1:nsamples;
N = 34;
h0 = zeros(1,nsamples);
h1 = h0;
for j = 1:nsamples
    Z0 = normrnd(mu0,sig,N,1);
    h0(j) = ttest(Z0,mu0,alpha);
    Z1 = normrnd(mu1,sig,N,1);
    h1(j) = ttest(Z1,mu0,alpha);
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','Location','East')

Предположим, мы можем позволить себе взять выборку размером 50. Предположительно, наша мощность для обнаружения альтернативного значения mu = 110 будет больше 80%. Если мы сохраним мощность на уровне 80%, какую альтернативу мы можем обнаружить?

mu1 = sampsizepwr('t',[100 20],[],.8,50)

mu1 =

  108.0837

Тестирование пропорции

Теперь обратимся к проблеме определения размера выборки, необходимого для различения двух пропорций. Предположим, что мы пробуем население, в котором около 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

Усложнение биномиального распределения происходит потому, что оно дискретно. Вероятность превышения значения отсечки не составляет ровно 5%:

1 - binocdf(cutoff, N, p0)

ans =

    0.0340

Еще раз, давайте рассчитаем мощность по альтернативе P = 0,33 для диапазона размеров выборки. Мы будем использовать односторонний (правохвостый) тест, потому что нас интересуют только альтернативные значения, превышающие 30%.

Nvec = 50:50:2000;
power = sampsizepwr('p',p0,p1,[],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,[],'tail','right')

Warning: Values N>200 are approximate.  Plotting the power as a function
of N may reveal lower N values that have the required power. 

approxN =

        1500

Предупреждающее сообщение сообщает нам, что ответ просто приблизительный. Если мы посмотрим на функцию мощности для различных размеров выборки, мы увидим, что функция обычно увеличивается, но нерегулярно, потому что биномиальное распределение дискретно. Рассмотрим вероятность отклонения нулевой гипотезы как для p = 0,30, так и для p = 0,33 в диапазоне размеров выборок от 1470 до 1480.

subplot(3,1,1);
Nvec = 1470:1480;
power = sampsizepwr('p',p0,p1,[],Nvec,'tail','right');
plot(Nvec,power,'ro-',[min(Nvec),max(Nvec)],[DesiredPower DesiredPower],'k:');
ylabel(sprintf('Prob[T>cutoff]\nif p=0.33'))
h_gca = gca;
h_gca.XTickLabel = '';
ylim([.78, .82]);

subplot(3,1,2);
alf = sampsizepwr('p',p0,p0,[],Nvec,'tail','right');
plot(Nvec,alf,'bo-',[min(Nvec),max(Nvec)],[alpha alpha],'k:');
ylabel(sprintf('Prob[T>cutoff]\nif p=0.30'))
h_gca = gca;
h_gca.XTickLabel = '';
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%:

min(Nvec(power>=0.80))

ans =

        1478

Тестирование корреляции

В примерах, которые мы рассматривали до сих пор, мы смогли вычислить предел для тестовой статистики, чтобы достичь определенного уровня значимости, и вычислить вероятность превышения этого предела при альтернативной гипотезе. Для нашего последнего примера давайте рассмотрим проблему, где это не так просто.

Представьте, что мы можем взять выборки из двух переменных X и Y, и мы хотим знать, какой размер выборки нам потребуется, чтобы проверить, являются ли они некоррелированными по сравнению с альтернативой, что их корреляция равна 0,4. Хотя можно проработать распределение корреляции выборки путем преобразования ее в t-распределение, давайте используем метод, который мы можем использовать даже в тех проблемах, где мы не можем выработать распределение тестовой статистики.

Для данного размера выборки мы можем использовать моделирование Монте-Карло, чтобы определить приблизительное значение отсечки для теста корреляции. Давайте сделаем большой прогон моделирования, чтобы мы могли получить это значение точно. Мы начнем с размера выборки 25.

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

cutoff =

    0.3372

Затем мы можем сгенерировать образцы в соответствии с альтернативной гипотезой и оценить мощность теста.

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

power =

    0.6470


powerci =

    0.6165    0.6767

Мы оцениваем мощность в 65%, и у нас есть 95% уверенности, что истинное значение находится между 62% и 68%. Чтобы получить мощность 80%, нам нужен больший размер выборки. Мы можем попытаться увеличить N до 50, оценить значение отсечки для этого размера выборки и повторить моделирование мощности.

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

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

cutoff =

    0.2315


power =

    0.8990


powerci =

    0.8786    0.9170

Этот размер выборки дает мощность лучше, чем наш целевой показатель в 80%. Мы могли бы продолжить эксперименты таким образом, пытаясь найти размер выборки менее 50, который соответствовал бы нашим требованиям.

Заключение

Функции вероятности в инструментарии статистики и машинного обучения могут использоваться для определения размера выборки, необходимого для достижения требуемого уровня мощности в тесте гипотезы. При некоторых проблемах размер выборки можно вычислить непосредственно; в других необходимо выполнять поиск по диапазону размеров выборки до тех пор, пока не будет найдено правильное значение. Генераторы случайных чисел могут помочь проверить соответствие требуемой мощности, а также могут быть использованы для изучения мощности конкретного теста в альтернативных условиях.