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

Одним из типов подтверждающих клинических исследований являются исследования с целью доказательства терапевтической эквивалентности. Методологически они похожи на исследования неуступающей/превышающей эффективности, но с одной особенностью: вместо статистический значимости различий интерес представляет ширина доверительного интервала для величины различий по главной(-ым) переменной(-ым). Доверительный интервал не должен выходить за границы зоны терапевтической эквивалентности, которые зачастую принимаются равными величине различий +/- 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₂. Варианты формул для произвольного соотношения размеров первой и второй группы см. в первоисточниках.