В этом сообщении я продолжу тему оценки мощности статистических критериев и покажу, как при помощи R можно выполнить анализ мощности при сравнении долей.

Источник: http://apgovernment2010.yolasite.com
Представим, что мы входим в команду кандитата в президенты страны X. Согласно результатам опросов, выполненных командой кандидата-соперника, выяснилось, что популярность нашего кандидата у городских жителей выше, чем у жителей села (28% против 20% среди 100 и 100 опрошенных респондентов соответственно). Безусловно, это важная информация, которая может помочь в планировании агитационных мероприятий (возможно, например, что следует направить больше ресурсов на агитацию среди сельских жителей). Однако, стоит ли доверять информации из лагеря соперника?


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

Начнем с анализа данных, полученных командой соперника. Как уже было показано в одном из предыдущих сообщений, две доли мы можем сравнить при помощи критерия хи-квадрат:

votes <- matrix(c(28, 72, 20, 80), ncol = 2, byrow = T)
votes
     [,1] [,2]
[1,]   28   72
[2,]   20   80
 
chisq.test(votes)
 
 Pearson's Chi-squared test with Yates' continuity correction
 
data:  votes 
X-squared = 1.3432, df = 1, p-value = 0.2465

Выполненный тест хи-квадрат позволяет заключить, что жители города и села статистически не различаются по своим предпочтениям в отношении нашего кандидата (p-value = 0.2465). Казалось бы, это хорошие новости, т.к. нашей команде, по-видимому, нет необходимости предпринимать дополнительные усилия для обеспечения более высокой лояльности среди сельских жителей. Но что если этот вывод неверен, просто в силу недостаточно большого числа опрошенных (т.е. мощность критерия недостаточно велика и мы совершили ошибку второго рода - не отклонили неверную нулевую гипотезу об отсутствии разницы между опрошенными категориями жителей)?

Хотя в базовой комплектации R имеется специальная функция для оценки мощности при сравнении двух долей - power.prop.test(), - мы воспользуемся возможностями специализированного пакета pwr (его можно установить при помощи команды install.packages("pwr")). Реализованные в этом пакете функции следуют формулам из известной книги Cohen (1988). Для анализа мощности критерия хи-квадрат служит функция pwr.chisq.test(), имеющая следующие аргументы:
  • w - размер эффекта
  • N - общее число наблюдений
  • df - число степеней свободы
  • sig.level - уровень значимости
  • power - мощность критерия

Следует пояснить, что подразумевается под размером эффекта w. По Коэну (Cohen 1988), w для \(\chi^2\) рассчитывается как
\[w = \sqrt{\sum_{}\frac{(p0_i - p1_i)^2}{p0_i}},\]

где \(p0_i \) - ожидаемая вероятность для i-й ячейки таблицы сопряженности при верной нулевой гипотезе об отсутствии эффекта, а \(p1_i \) - наблюдаемая вероятность. Иными словами, w измеряет степень отклонения наблюдаемых частот в таблице сопряженности от тех, которые можно было бы ожидать при отсутствии эффекта исследуемого фактора.

Рассчитаем мощность приведенного выше теста \(\chi^2\). Для удобства, сначала сохраним результаты теста в виде самостоятельного объекта (назовем его res):

res <- chisq.test(votes)
 
res
 
 Pearson's Chi-squared test with Yates' continuity correction
 
data:  votes 
X-squared = 1.3432, df = 1, p-value = 0.2465

Объект res представляет собой список из нескольких компонентов, включая таблицы с наблюдаемыми и ожидаемыми частотами. Мы можем легко извлечь эти таблицы и рассчитать соответствующие вероятности:

obs <- res$observed
 
obs
     [,1] [,2]
[1,]   28   72
[2,]   20   80
 
exptd <- res$expected
 
exptd
     [,1] [,2]
[1,]   24   76
[2,]   24   76
 
obs <- obs/200 # здесь и ниже, 200 - общее число опрошенных
 
obs
     [,1] [,2]
[1,] 0.14 0.36
[2,] 0.10 0.40
 
exptd <- exp/200
 
exptd
     [,1] [,2]
[1,] 0.12 0.38
[2,] 0.12 0.38

Используя приведенную выше формулу, рассчитаем размер эффекта:

sqrt(sum((exptd-obs)^2/exptd))
[1] 0.09365858

Приведенные вычисления размера эффекта можно было значительно сократить, воспользовавшись готовой функцией из пакета pwrES.w2(), на которую подается таблица сопряженности с наблюдаемыми вероятностями (в виде матрицы):

library(pwr) 
ES.w2(obs)
[1] 0.09365858

Много это или мало - 0.094? Показатель w изменяется от 0 до 1, и чем ближе его значение к 1, тем сильнее эффект исследуемого фактора. В рассматриваемом случае размер эффекта весьма невысок - место проживания респондентов очень слабо связано с их лояльностью в отношении нашего кандидата. При таком размере эффекта и при общем числе опрошенных респондентов, равном 200, мощность выполненного выше теста \(\chi^2\) составляет

pwr.chisq.test(w = ES.w2(obs), df = 1, N = 200)
 
     Chi squared power calculation 
 
              w = 0.09365858
              N = 200
             df = 1
      sig.level = 0.05
          power = 0.2630843
 
NOTE: N is the number of observations

Рассчитанная мощность (~ 26%) гораздо ниже условно принятого порога в 80%. Таким образом, команда соперника не имеет оснований утверждать, что наш кандидат менее популярен среди сельских жителей - число опрошенных ими респондентов недостаточно для такого вывода при таком небольшом размере эффекта.

Перейдем к планированию нашего собственного, независимого исследования. Поскольку имеющиеся ресурсы ограничены, решено, что результаты этого исследования будут приняты во внимание при планировании агитационных мероприятий только в том случае, если выявленный размер эффекта окажется достаточно большим, и если мощность критерия \(\chi^2\) будет не менее 80% при  уровне значимости 0.05. Однако что значит "достаточно большой" эффект? Согласно классификации, предложенной Коэном (Cohen 1988), w = 0.10 следует считать небольшим эффектом, w = 0.3 - умеренно большим, а w = 0.5 и выше - большим эффектом. Примем, что порог w = 0.15 является достаточно значительным эффектом, чтобы обратить на него внимание в нашей ситуации (должность президента на кону, как никак!). Тогда минимальное число подлежащих опросу респондентов составит (обратите внимание на аргумент N - ему присвоено значение NULL, поскольку рассчитывается именно число наблюдений):

pwr.chisq.test(w = 0.15, N = NULL, df = 1, sig.level = 0.05, power = 0.8)
 
     Chi squared power calculation 
 
              w = 0.15
              N = 348.8382
             df = 1
      sig.level = 0.05
          power = 0.8
 
 NOTE: N is the number of observations

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

В завершение, следует упомянуть еще три функции из пакета pwr, созданные для анализп статистической мощности при работе с долями:
  • pwr.2p.test() - анализ мощности для двух долей, рассчитанных по выборкам одинакового размера
  • pwr.2p2n.test() - анализ мощности для двух долей, рассчитанных по выборкам разного размера
  • pwr.p.test() - анализ мощности для одной доли (например, при сравнении ее с каким-либо ожидаемым значением)


7 Комментарии

Анонимный написал(а)…
Я с большим интересом читаю Ваш блог и очень ценю Ваши усилия, потому что в русскоязычном Интернете совсем мало информации об R. А популяризировать эту программу и всю сферу анализа данных надо. В отношении этого поста у меня есть пара замечаний. Оценка мощности производится, когда количественно сформулирована альтернативная гипотеза. Поскольку конкуренты ее не представили, а читать мысли мы не умеем, подсчет мощности исследования конкурентов невозможен в принципе. Второе, мощность исследования никогда не используется для тестирования гипотез. Гипотезы тестируются с использованием rejection region, который относится к ошибке альфа, но никак не бета. Поэтому сравнение 26% и 80%, чтобы опровергнуть гипотезу конкурентов о разности между городом и селом, неуместно.
Sergey Mastitsky написал(а)…
Уважаемый Анонимный,

Большое спасибо за внимание к блогу и за Ваши полезные замечания по поводу этого поста. Попробую прокомментировать их:

1) "Оценка мощности производится, когда количественно сформулирована альтернативная гипотеза. Поскольку конкуренты ее не представили, а читать мысли мы не умеем, подсчет мощности исследования конкурентов невозможен в принципе".

Действительно, как таковая проверяемая гипотеза в тексте не сформулирована. Однако приведен тест хи-квадрат, который сравнивает две доли, полученные "конкурентами" (28% против 20%). Нулевая гипотеза теста состоит в том, что разницы между этими долями нет, альтернативная - в том, что она есть (без уточнения, в какую сторону и на сколько). Согласно результатам теста, нулевая гипотеза не отклонена.

Вы правы в том, что "конкуренты" действительно не сообщили о размере эффекта, который они имели в виду. Поэтому (в образовательных целях) мы приняли, что имелся в виду максимальный эффект и рассчитали его по приведенной формуле Коэна. Используя полученное значение, была рассчитана так называемая наблюдаемая мощность теста ("observed test power", иногда также "retrospective test power"). Имеет ли смыл ее расчитывать? С точки зрения продемонстрировать принцип - да. В примере как раз и показано, как это можно сделать.

2) Другое дело, что между Р-значем и мощностью теста, как известно, имеется тесная функциональная связь - чем больше Р, тем меньше мощность. Другими словами, получаемые в ходе теста Р-значения > 0.05 (например) по определению указывают на недостаточно высокую мощность теста. В этом смысле Вы совершенно правы, что делать вывод о верности нулевой гипотезы по величине мощности теста неуместно (обсуждению этой проблемы посвящено огромное количество статей; вот, например, одна из наиболее часто цитируемых: http://www.vims.edu/people/hoenig_jm/pubs/hoenig2.pdf). Рассчитав наблюдаемую мощность (26%) и сравнив ее с условным порогом в 80%, я хотел лишь еще раз подчеркнуть взаимосвязь межу величиной эффекта, объемом наблюдений и мощностью теста. Видимо, неудачно сформулировал мысль. Надеюсь, этот комментарий поможет Читалю лучше понять исходную идею.

В дополнение стоит еще раз подчеркнуть, что выполнение анализа мощности имеет смыл только на этапе планирования исследования.
Unknown написал(а)…
Сергей, спасибо за обратную связь. В ответ на Ваши ответы :)

1. Чтобы понять, что нет свидетельств в пользу различий, нет нужды рассчитывать мощность. Для этого есть данные и нулевая гипотеза. Выбор 28% и 20% для расчета мощности в целом понятен, это все же maximum likelihood estimations, но писать, что они соответствуют максимальному эффекту не корректно (40% и 8% имеют размер эффекта значительно больше).
2. Утверждение "получаемые в ходе теста Р-значения > 0.05 (например) по определению указывают на недостаточно высокую мощность теста" является не точным. Такие значения можно получить в двух случаях: а) если эффекта нет; б) если эффект есть, но нет достаточных свидетельств в его пользу. Утверждение с оговорками относится ко второму случаю.
Sergey Mastitsky написал(а)…
Уважаемый Baurzhan,

1) Пожалуйста, см. мой ответ на Ваш первый комментарий выше - т.н. "наблюдаемая мощность" была рассчитана с целью продемонстрировать принцип работы функции pwr.chisq.test(), а также взаимосвязь между величиной эффекта, объемом наблюдений и мощностью. Вопрос, видимо, надо стативить по-другому: должна ли статистическая программа позволять рассчитывать наблюдаемую мощность теста? Вот здесь [http://research-repository.st-andrews.ac.uk/bitstream/10023/679/5/Thomas-Retrospectivepoweranalysis-postprint.pdf] автор выражает сожаление, что такие "монстры рынка", как SAS, SPSS, и др. (R в их числе) позволяют это делать. С этой точкой зрения, думаю, согласны и Вы, и я.

По поводу "максимального эффекта" - имеется в виду эффект, рассчитанный по имеющимся данным, который составил 0.094. Что если бы "конкуренты" утверждали, что связь между местом проживани респондентов и их лояльностью в отношении "нашего" кандидата составляет, например, 0.04? Этого они не сделали (во всяком случае, у нас такой информации нет и этот эффект мы "протестировать" не может), поэтому и рассчитан действительно наблюдаемый эффект.

2) См. рисунок 1 и соответствующее обсуждение из процитированной ранее статьи ( http://www.vims.edu/people/hoenig_jm/pubs/hoenig2.pdf ). В частности: "Because of the one-to-one relationship between p values and observed power, nonsigni cant p values always correspond to low observed powers"
Unknown написал(а)…
Думаю, что дискуссию на этом можно остановить. Как говорил мой бывший профессор по биостатистике Игорь Зурбенко, ученик и соратник Колмогорова, если хотите рассчитать силу исследования, сделайте симуляцию. Ниже функция, которая рассчитывает мощность исследования по четырем параметрам - количество субъектов в первой группе, количество субъектов во второй группе, ожидаемая частота в первой группе, ожидаемая частота во второй группе. Генерируется 10000 таблиц 2х2 и на основе них рассчитывается сила, то есть вероятность отвергнуть нулевую гипотезу.

prop.power=function(n1,n2,p1,p2) {
twobytwo=matrix(NA,nrow=10000, ncol=4)
twobytwo[,1]=rbinom(n=10000, size=n1, prob=p1)
twobytwo[,2]=n1-twobytwo[,1]
twobytwo[,3]=rbinom(n=10000, size=n2, prob=p2)
twobytwo[,4]=n1-twobytwo[,3]
p=rep(NA, 10000)
chisq.test.v=function(x) as.numeric(chisq.test(matrix(x, ncol=2), correct=FALSE)[3])
p=apply(twobytwo, 1, chisq.test.v)
power=sum(ifelse(p<0.05,1,0))/10000
return(power)
}
Функция дает схожие результаты:
> prop.power(100,100,0.28, 0.20)
[1] 0.2649
Sergey Mastitsky написал(а)…
Спасибо, Baurzhan! Приведенный Вами код - отличное дополнение к посту!
Буду рад Вашим комментариям и в будущем.
Unknown написал(а)…
Спасибо, Сергей!
Я нашел ошибку в коде в строке 'twobytwo[,4]=n1-twobytwo[,3]'
Исправленная функция выглядит так:

prop.power=function(n1,n2,p1,p2) {
twobytwo=matrix(NA,nrow=10000, ncol=4)
twobytwo[,1]=rbinom(n=10000, size=n1, prob=p1)
twobytwo[,2]=n1-twobytwo[,1]
twobytwo[,3]=rbinom(n=10000, size=n2, prob=p2)
twobytwo[,4]=n2-twobytwo[,3]
p=rep(NA, 10000)
chisq.test.v=function(x) as.numeric(chisq.test(matrix(x, ncol=2), correct=FALSE)[3])
p=apply(twobytwo, 1, chisq.test.v)
power=sum(ifelse(p<0.05,1,0))/10000
return(power)
}
Новые Старые