Оценка статистической мощности при сравнении долей

В этом сообщении я продолжу тему оценки мощности статистических критериев и покажу, как при помощи 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

Приведенные вычисления размера эффекта можно было значительно сократить, воспользовавшись готовой функцией из пакета pwr -  ES.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() - анализ мощности для одной доли (например, при сравнении ее с каким-либо ожидаемым значением)