Анализ мощности

Мы рассмотрим

Мощность статистического теста
Способы оценки величины эффекта
A priori анализ мощности
Использование симуляций для анализа мощности
Как влиять на мощность тестов

Мощность статистического теста

Типы ошибок при проверке гипотез

	H0 == TRUE	H0 == FALSE
Отклонить H0	Ошибка I рода α	Верно 1 - α
Сохранить H0	Верно	Ошибка II рода β

Мощность статистического теста

	H0 == TRUE	H0 == FALSE
Отклонить H0	Ошибка I рода α	Верно 1 - α
Сохранить H0	Верно 1 - β	Ошибка II рода β

Мощность статистического теста —

вероятность найти различия там, где они есть \(Power = 1 - \beta\)

Анализ мощности

Какой нужен объем выборки, чтобы найти различия?
Какой величины различия мы можем найти при определенном объеме выборки?

Величина эффекта

\[\frac{\bar \mu_1 - \bar \mu_2}{\sigma}\]

Величина эффекта

\[\frac{\bar \mu_1 - \bar \mu_2}{\sigma}\]

\(\sigma\) — среднеквадратичное стандартное отклонение

\(d = \frac {|\bar x_1 - \bar x_2|} { \sqrt {\frac {s_1^2 + s_2^2 } {2} }}\)

Величина эффекта

\[\frac{\bar \mu_1 - \bar \mu_2}{\sigma}\]

\(\sigma\) — среднеквадратичное стандартное отклонение

\(d = \frac {|\bar x_1 - \bar x_2|} { \sqrt {\frac {s_1^2 + s_2^2 } {2} }}\)

\(\sigma\) — обобщенное стандартное отклонение

\(d = \frac {|\bar x _{1} - \bar x _{2}|} { \sqrt {\frac {(n_{1} - 1) s_1^2 + (n_{2} - 1) s_{2}^2 } {n_{1} + n_{2} - 2} } }\)

Как оценить ожидаемую величину эффекта?

\[\frac{\bar \mu_1 - \bar \mu_2}{\sigma}\]

Варианты:

Пилотные исследования
Литература
Общебиологические знания
Технические требования

Условные уровни величины эффекта

Сильные, умеренные и слабые эффекты (Cohen, 1982)

library(pwr)
cohen.ES(test = "t", size = "large")

## 
##      Conventional effect size from Cohen (1982) 
## 
##            test = t
##            size = large
##     effect.size = 0.8

Задача

Рассчитайте величину умеренных и слабых эффектов для t-критерия при помощи функции cohen.ES()

Решение

# умеренный
cohen.ES(test = "t", size = "medium")

## 
##      Conventional effect size from Cohen (1982) 
## 
##            test = t
##            size = medium
##     effect.size = 0.5

# слабый
cohen.ES(test = "t")

## 
##      Conventional effect size from Cohen (1982) 
## 
##            test = t
##            size = small
##     effect.size = 0.2

Задача

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

Т.е. ожидаемая разница логарифмов экспрессии между опытом и контролем будет примерно \(log2(1) - log2(5) = -2.32\)

Стандартное отклонение этих различий по пилотным данным 1.2

Рассчитайте величину эффекта

Скажите, это будет слабый или сильный эффект?

Решение

abs(-2.32) / 1.2

## [1] 1.93

Расчет величины эффекта по данным пилотного исследования

Величина эффекта по данным пилотного исследования

Представьте себе, что у нас есть результаты пилотного исследования. Мы хотим знать, различается ли уровень клеточной активности у пациентов в зависимости от того, наступила ли ремиссия. (Данные из кн. Freeman 1987).

library(readxl)
rem <- read_excel(path = "data/remission.xlsx", sheet = 1)
head(rem)

## # A tibble: 6 × 3
##      LI     m     r
##   <dbl> <dbl> <dbl>
## 1   0.4     1     0
## 2   0.4     1     0
## 3   0.5     1     0
## 4   0.5     1     0
## 5   0.6     1     0
## 6   0.6     1     0

Это данные о больных раком:

LI — мера клеточной активности
r — индикатор того наступила ли ремиссия (1 — да, 0 — нет)

Убедимся, что данные открылись правильно

str(rem)

## Classes 'tbl_df', 'tbl' and 'data.frame':    27 obs. of  3 variables:
##  $ LI: num  0.4 0.4 0.5 0.5 0.6 0.6 0.6 0.7 0.7 0.7 ...
##  $ m : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ r : num  0 0 0 0 0 0 0 0 0 0 ...

rem$r <- factor(rem$r, levels = c(0, 1), labels = c("no", "yes"))

## Пропущенные значения
colSums(is.na(rem))

## LI  m  r 
##  0  0  0

## Объем выборки
table(rem$r)

## 
##  no yes 
##  18   9

Задача

Визуализируйте при помощи бокс-плотов различия клеточной активности в двух группах пациентов.

Используйте геом geom_boxplot()

Обозначьте группу при помощи заливки.

Подпишите оси при помощи элемента labs()

Решение

library(ggplot2)
ggplot(rem, aes(x = r, y = LI)) + 
  geom_boxplot(aes(fill = r)) + 
  labs(x = "Remission", y = "Cell activity", fill = "Remission")

Давайте оценим величину эффекта

library(effsize)
effect <- cohen.d(rem$LI, rem$r)
effect

## 
## Cohen's d
## 
## d estimate: -1.32 (large)
## 95 percent confidence interval:
##    inf    sup 
## -2.235 -0.398

Задача

Посмотрите на структуру объекта effect при помощи функции str() и добудьте из него средствами R значение величины эффекта

Решение

str(effect)

## List of 6
##  $ method    : chr "Cohen's d"
##  $ name      : chr "d"
##  $ estimate  : Named num -1.32
##   ..- attr(*, "names")= chr "no"
##  $ conf.int  : Named num [1:2] -2.235 -0.398
##   ..- attr(*, "names")= chr [1:2] "inf" "sup"
##  $ conf.level: num 0.95
##  $ magnitude : Ord.factor w/ 4 levels "negligible"<"small"<..: 4
##  - attr(*, "class")= chr "effsize"

effect$estimate

##    no 
## -1.32

Анализ мощности в R

Анализ мощности t-критерия в R

Давайте рассчитаем, какой будет нужен объем выборки, чтобы показать различия между группами с вероятностью 0.8 (т.е. чтобы мощность теста была 80%)

library(pwr)
pwr.t.test(n = NULL, d = effect$estimate, power = 0.8, 
           sig.level = 0.05, type = "two.sample", 
           alternative = "two.sided")

## 
##      Two-sample t test power calculation 
## 
##               n = 10.1
##               d = 1.32
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

Задача

Экстрагируйте из результатов pwr.t.test() величину объема выборки

Решение

R1 <- pwr.t.test(n = NULL, d = effect$estimate, power = 0.8,
                 sig.level = 0.05, type = "two.sample", 
                 alternative = "two.sided")
str(R1)

## List of 7
##  $ n          : num 10.1
##  $ d          : Named num 1.32
##   ..- attr(*, "names")= chr "no"
##  $ sig.level  : num 0.05
##  $ power      : num 0.8
##  $ alternative: chr "two.sided"
##  $ note       : chr "n is number in *each* group"
##  $ method     : chr "Two-sample t test power calculation"
##  - attr(*, "class")= chr "power.htest"

R1$n

## [1] 10.1

Округляем в большую сторону

ceiling(R1$n)

## [1] 11

Задача

А если группы будут разного размера?

Пациентов с ремиссией всего 5.

Сколько нужно обследовать пациентов без ремиссии, чтобы обнаружить различия клеточной активности между группами?

Используйте функцию pwr.t2n.test

Решение

R2 <- pwr.t2n.test(n1 = 5, n2 = NULL, d = effect$estimate, 
                   power = 0.8, sig.level = 0.05, 
                   alternative = "two.sided")
ceiling(R2$n2)

## [1] 68

Другой пример

Данные анализа мочи. В некоторых пробах обнаружены кристаллы оксалата кальция. (Данные из сборника датасетов Andrews Herzberg 1985)

Представим, что это данные пилотного исследования.

ur <- read_excel(path = "data/urine.xlsx", sheet = 1)
head(ur)

## # A tibble: 6 × 7
##       r gravity    ph  osmo  cond  urea  calc
##   <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     0    1.02  4.91   725    NA   443  2.45
## 2     0    1.02  5.74   577  20.0   296  4.49
## 3     0    1.01  7.20   321  14.9   101  2.36
## 4     0    1.01  5.51   408  12.6   224  2.15
## 5     0    1.00  6.52   187   7.5    91  1.16
## 6     0    1.02  5.27   668  25.3   252  3.34

r — индикатор присутствия кристаллов оксалата кальция
calc — концентрация кальция

Убедимся, что данные открылись правильно

str(ur)

## Classes 'tbl_df', 'tbl' and 'data.frame':    79 obs. of  7 variables:
##  $ r      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ gravity: num  1.02 1.02 1.01 1.01 1 ...
##  $ ph     : num  4.91 5.74 7.2 5.51 6.52 5.27 5.62 5.67 5.41 6.13 ...
##  $ osmo   : num  725 577 321 408 187 ...
##  $ cond   : num  NA 20 14.9 12.6 7.5 25.3 17.4 35.9 21.9 25.7 ...
##  $ urea   : num  443 296 101 224 91 252 195 550 170 382 ...
##  $ calc   : num  2.45 4.49 2.36 2.15 1.16 3.34 1.4 8.48 1.16 2.21 ...

ur$r <- factor(ur$r, levels = c(0, 1), 
               labels = c("absent", "present"))

## Пропущенные значения
colSums(is.na(ur))

##       r gravity      ph    osmo    cond    urea    calc 
##       0       0       0       1       1       0       0

## Объем выборки
table(ur$r)

## 
##  absent present 
##      45      34

Задача

Постройте гистограмму распределения концентрации кальция. Заливкой обозначьте присутствие кристаллов оксалата кальция в пробе.

Решение

ggplot(ur, aes(x = calc, fill = r)) + 
  geom_histogram(binwidth = 1, position = "dodge") + 
  labs(x = "Calcium, mmol/l", y = "N", fill = "Crystals of calcium oxalate")

Задача

Сколько нужно анализов, чтобы показать, что содержание кальция в этих пробах достоверно отличается при уровне значимости \(\alpha = 0.001\)?

При равном объеме выборок
Если будет только 16 проб с кристаллами оксалата кальция

Решение

effect_ur <- cohen.d(ur$calc, ur$r)
R3 <- pwr.t.test(n = NULL, d = abs(effect_ur$estimate),
                 power = 0.8, sig.level = 0.001, 
                 type = "two.sample", alternative = "two.sided")
ceiling(R3$n)

## [1] 24

R4 <- pwr.t2n.test(n1 = 16, n2 = NULL, d = effect_ur$estimate, 
                   power = 0.8, sig.level = 0.001, 
                   alternative = "two.sided")
ceiling(R4$n2)

## [1] 43

Анализ мощности t-критерия при помощи симуляции

При помощи простой симуляции можно проверить, будет ли достаточно 24 наблюдений в группе

Сначала готовим данные

# средние значения из пилотного исследования
mu <- tapply(X = ur$calc, INDEX = ur$r, FUN = mean)

# обобщенное стандартное отклонение sigma в пилотном исследовании
ns <- table(ur$r)
var <- tapply(X = ur$calc, INDEX = ur$r, FUN = var)
sigma <- sqrt(((ns[1] - 1) * var[1] + (ns[2] - 1) * var[2])/(length(ur$r) - 2))

n <- 24 # проверяемое значение объема выборки
reps <- 1000 # число повторов в симуляции

ВАРИАНТ 1 - Используем цикл

pvals <- rep(NA, reps) # Место для результата
set.seed(42) # зерно генератора случайных чисел

# в цикле  многократно симулируем данные 
for (i in 1:reps){
  # генерируем две случайные выборки 
  # из нормального распределения
  x1 <- rnorm(n, mu[1], sigma)
  x2 <- rnorm(n, mu[2], sigma)
  # сравниваем их средние значения t-критерием
  t_res <- t.test(x = x1, y = x2)
  # добываем уровень значимости
  pvals[i] <- t_res$p.value
}
# Доля уровней значимости меньше критического уровня (здесь p < 0.001)
mean(pvals < 0.001)

## [1] 0.818

ВАРИАНТ 2 - Используем `replicate()`

# функция, генерирующая одну симуляцию
# n - объем выборки
my_sim <- function(n) {
  x1 <- rnorm(n, mu[1], sigma)
  x2 <- rnorm(n, mu[2], sigma)
  t_res <- t.test(x = x1, y = x2)
  return(t_res$p.value)
}
# my_sim(24) # пример употребления - одна симуляция

# Повторяем симуляцию 1000 раз
set.seed(42)
pvals <- replicate(1000, my_sim(24))
mean(pvals < 0.001)

## [1] 0.818

Как проверить несколько возможных объемов выборки?

Вариант 2 легче масштабируется.

Достаточно дописать еще одну функцию, и мы сможем оценить мощность теста при разных объемах выборки.

# функция, возвращающая долю симуляций с p < 0.001
# x - объем выборки
my_sim_range <- function(x){
  pvals <- replicate(1000, my_sim(x))
  pw <- mean(pvals < 0.001)
  return(pw)
}
set.seed(42)
sapply(X = 20:25, FUN = my_sim_range)

## [1] 0.666 0.735 0.730 0.780 0.787 0.844

Задание

По итогам предыдущей симуляции постройте кривую зависимости мощности от объема выборки (от 5 до 35)

Решение

set.seed(42)
sim <- sapply(X = 5:35, FUN = my_sim_range)

pwr_dat <- data.frame(N = 5:35, Power = sim)
ggplot(pwr_dat, aes(x = N, y = Power)) + 
  geom_point(colour = "green3") + 
  geom_line(aes(group = 1), colour = "green3") +
  geom_hline(yintercept = 0.8, linetype = "dashed")

Задание

При помощи симуляций проведите анализ мощности для данных о клеточной активности у больных раком

Постройте кривую зависимости мощности теста от объема выборки

Решение: подготовка

# средние значения из пилотного исследования
mu <- tapply(X = rem$LI, INDEX = rem$r, FUN = mean)

# обобщенное стандартное отклонение sigma в пилотном исследовании
ns <- table(rem$r)
var <- tapply(X = rem$LI, INDEX = rem$r, FUN = var)
sigma <- sqrt(((ns[1] - 1) * var[1] + (ns[2] - 1) * var[2])/(length(ur$r) - 2))

reps <- 1000 # число повторов в симуляции

Решение: симуляция

# функция, генерирующая одну симуляцию
my_sim <- function(n) {
  x1 <- rnorm(n, mu[1], sigma)
  x2 <- rnorm(n, mu[2], sigma)
  t_res <- t.test(x = x1, y = x2)
  return(t_res$p.value)
}

# функция, возвращающая долю симуляций со значением p меньше критического
my_sim_range <- function(x){
  pvals <- replicate(1000, my_sim(x))
  fraction <- mean(pvals < 0.001)
  return(fraction)
}

set.seed(42)
sim <- sapply(X = 5:15, FUN = my_sim_range)

Решение: график зависимости мощности от объема выборки

pwr_dat <- data.frame(N = 5:15, Power = sim)
ggplot(pwr_dat, aes(x = N, y = Power)) + 
  geom_point(colour = "steelblue") + 
  geom_line(aes(group = 1), colour = "steelblue") +
  geom_hline(yintercept = 0.8, linetype = "dashed")

Как влиять на мощность теста?

Мощность зависит

от объема выборки
от величины эффекта
от уровня значимости

Чем больше объем выборки—тем больше мощность

Чем больше уровень значимости—тем больше мощность

Чем больше величина различий—тем больше мощность

Скажите

Какие из факторов, влияющих на мощность теста, мы не можем контролировать?

Мы не можем контролировать внешние факторы
величину эффекта (\(ES\))
фоновую изменчивость (\(\sigma^2\))

Каким образом можно повлиять на мощность теста?

Мощность теста можно регулировать, если
изменить число повторностей
выбрать другой уровень значимости (\(\alpha\))
определиться, какие эффекты действительно важны (\(ES\))

Take home messages

Способность выявлять различия зависит
от объема выборки,
от уровня значимости
от величины эффекта

Дополнительные ресурсы

Quinn, Keough, 2002, pp. 164-170
OpenIntro Statistics

Мы рассмотрим

Мощность статистического теста

Типы ошибок при проверке гипотез

Мощность статистического теста

Мощность статистического теста —

Анализ мощности

Величина эффекта

Величина эффекта

Величина эффекта

Величина эффекта

Как оценить ожидаемую величину эффекта?

Условные уровни величины эффекта

Задача

Решение

Задача

Решение

Расчет величины эффекта по данным пилотного исследования

Величина эффекта по данным пилотного исследования

Убедимся, что данные открылись правильно

Задача

Решение

Давайте оценим величину эффекта

Задача

Решение

Анализ мощности в R

Анализ мощности t-критерия в R

Задача

Решение

Задача

Решение

Другой пример

Убедимся, что данные открылись правильно

Задача

Решение

Задача

Решение

Анализ мощности t-критерия при помощи симуляции

Анализ мощности t-критерия при помощи симуляции

ВАРИАНТ 1 - Используем цикл

ВАРИАНТ 2 - Используем replicate()

Как проверить несколько возможных объемов выборки?

Задание

Решение

Задание

Решение: подготовка

Решение: симуляция

Решение: график зависимости мощности от объема выборки

Как влиять на мощность теста?

Мощность зависит

Чем больше объем выборки—тем больше мощность

Чем больше уровень значимости—тем больше мощность

Чем больше величина различий—тем больше мощность

Скажите

Take home messages

Дополнительные ресурсы

ВАРИАНТ 2 - Используем `replicate()`