Терапевтическая эквивалентность: размер выборки

Одним из типов подтверждающих клинических исследований являются исследования с целью доказательства терапевтической эквивалентности. Методологически они похожи на исследования неуступающей/превышающей эффективности, но с одной особенностью: вместо статистический значимости различий интерес представляет ширина доверительного интервала для величины различий по главной(-ым) переменной(-ым). Доверительный интервал не должен выходить за границы зоны терапевтической эквивалентности, которые зачастую принимаются равными величине различий +/- 20%:

Источник: Sample Size Tables for Clinical Studies (Machin D., 2009)

Использование дихотомической переменной с категориями "терапия эффективна" и "терапия не эффективна" является очевидным и предпочтительным подходом и при изучении эффективности как таковой, и при доказательстве терапевтической эквивалентности. Для дихотомической переменной оценка размера выборки может быть выполнена с использованием простой формулы (из книги Sample size calculations in clinical research):

$$n_1=n_2=\frac{\left ( z_\alpha + z_\beta \right )^{2}\left ( p_1\left ( 1-p_1 \right ) + p_2\left ( 1-p_2 \right ) \right )}{\left ( \left | \delta \right |-\varepsilon\right )^{2}},$$

где $z_\alpha$ и $z_\beta$ - процентные точки стандартного нормального распределения (берется двусторонняя вероятность для ошибок первого и второго рода), $p_1$ и $p_2$ - доля положительных результатов в первой и во второй группах, $\delta=p_1-p_2$, $\varepsilon$ - величина предельно допустимых различий.

Но тут нужно вспомнить, что при такой оценке размера выборки предполагается дальнейшее сравнение  с помощью критерия хи-квадрат Пирсона или эквивалентное такому сравнению построение доверительного интервала. А критерий хи-квадрат Пирсона дает корректные результаты только тогда, когда минимальная ожидаемая частота в таблице сопряженности больше или равна 5 (иными словами - только для случаев, когда ожидаемые доли в пропорциях не слишком близки к 0 или 1). Соответственно, для случая, когда предполагаемые доли положительных исходов в группах составляют, например, 97% и 98%, нужно использовать иной подход и при оценке размера выборки, и при обработке результатов.

Для более корректного расчета воспользуемся формулой с использованием метода максимального правдоподобия, приведенной в книге Sample Size Tables for Clinical Studies (Machin D., 2009):
$$n_1=n_2=\frac{z_\alpha \sqrt{\bar{p}_{1D}\left ( 1-\bar{p}_{1D} \right ) + \bar{p}_{2D}\left ( 1-\bar{p}_{2D} \right )} + z_\beta \sqrt{p_1\left ( 1-p_1 \right ) + p_2\left ( 1-p_2 \right )}}{\left ( \left | \delta \right | - \varepsilon \right )^{2}},$$

где все обозначения идентичны предыдущей формуле, а $\bar{p}_{1D}$ и $\bar{p}_{2D}$ - оценки максимального правдоподобия для истинных значений долей $p_1$ и $p_2$ в предположении, что различия между ними равны $\varepsilon$.

$\bar{p}_{1D}$ является корнем уравнения $ax^3+bx^2+cx+d$, в котором:
$$a=2,$$  
$$b=-\left ( 2+p_1+p_2+3\varepsilon \right ),$$
$$c=\varepsilon^2+\varepsilon\left (2p_1+2\right )+p_1+p_2,$$
$$d=-p_1\varepsilon\left (1+\varepsilon \right )$$
Также может быть показано, что $\bar{p}_{1D}=2u\cos \left ( w \right ) - b/3a$ и $\bar{p}_{2D}=\bar{p}_{1D}-\varepsilon$, где $$w=\frac{ \pi+\arccos \left ( v/u^3 \right )} {3},$$ $$v=\frac{b^3} {(3a)^3}-\frac{bc} {6a^2}+\frac{d} {2a},$$ $$u=sgn(v)\sqrt{\frac{b^2}{(3a)^2}-\frac{c}{3a}}$$
Реализация на R выглядит следующим образом:

propSample <- function(p1, p2, epsilon = 0.2, alpha = 0.05, beta = 0.2) {
 
    # p1, p2 - предполагаемая доля положительных результатов (событий и т.д.)
    # epsilon -  величина предельно допустимых различий, по умолчанию 0.2
    # alpha - уровень значимости (двусторонняя вероятность), по умолчанию 0.05
    # beta - вероятность ошибки второго рода (двусторонняя), по умолчанию 0.2
 
    z_alpha <- qnorm(1 - alpha / 2)
    z_beta <- qnorm(1 - beta / 2)
 
    delta <- p1 - p2
 
 
    if (p1 >= 0 && p2 >= 0 && p1 <= 1 && p2 <= 1 && 
        epsilon >= 0 && epsilon <= 1 && delta < epsilon ) {
 
        delta <- p1 - p2
 
        a <- 2
        b <- -(2 + p1 + p2 + 3 * epsilon)
        C <- epsilon^2 + epsilon * (2 * p1 + 2) + p1 + p2
        d <- -p1 * epsilon * (1 + epsilon)
 
        v <- b^3 / (3 * a)^3 - b * C / (6 * a^2) + d / (2 * a)
        u <- sqrt(b^2 / (3 * a)^2 - C / (3 * a))
        u <- ifelse(v >= 0, u, -u)
        w <- (pi + acos(v / u^3)) / 3
 
        p1D <- 2 * u * cos(w) -b /(3 * a)
        p2D <- p1D - epsilon
 
        n <- (z_alpha * sqrt(p1D * (1 - p1D) + p2D * (1 - p2D)) + 
              z_beta * sqrt(p1 * (1 - p1) + p2 * (1 - p2)))^2 / 
              (abs(delta) - epsilon)^2
        return(ceiling(n))
 
    } else if (delta >= epsilon) {
 
        stop("p1 - p2 больше величины предельно допустимых различий")
 
    } else {
 
        stop("Доля меньше 0 или больше 1; проверьте правильность ввода данных")
 
    }
 
}

Рассчитаем размер каждой группы для p₁=0,97 и p₁=0,98:

> propSample(0.97, 0.98)
[1] 35

Для сравнения, "обычная" формула даст n=15, т.е. требуемый размер каждой группы будет недооценен более чем вдвое. Продолжение следует.
P.S. Обозначения в формулах приведены к единообразному виду, а сами формулы упрощены для случая, когда n₁=n₂. Варианты формул для произвольного соотношения размеров первой и второй группы см. в первоисточниках.