Расчет размера выборки для одномерной логистической регрессии

11

Как рассчитать размер выборки, необходимый для исследования, в котором когорта субъектов будет иметь одну непрерывную переменную, измеренную во время операции, а затем через два года они будут классифицированы как функциональный результат или результат с нарушением.

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

Любые идеи? Любая реализация R.

Фаррел
источник
Ожидаете ли вы каких-либо отсева во время наблюдения? Есть ли другие ковариаты, которые будут включены в вашу модель?
ЧЛ
Позвольте мне высосать уровень отсева из моего большого пальца - 20%. Мы действительно собираем много переменных, например, возраст, оценка травмы, но я хотел, чтобы все было максимально просто для расчета мощности. Я часто находил полезным обсудить первичную модель, а затем вторичные модели, которые загружены с большей утонченностью и нюансами.
Фаррел
Хорошо, но обычно ожидаемый процент выпадения, число ковариат и то, измеряются ли ковариаты с ошибками (см., Например, j.mp/9fJkhb ), введите формулу (во всех случаях это увеличит размер выборки).
ЧЛ

Ответы:

7

Расчет размера выборки для логистической регрессии является сложным. Я не буду пытаться обобщить это здесь. Разумно доступные решения этой проблемы находятся в:

Hsieh FY. Таблицы размеров выборки для логистической регрессии. Статистика в медицине. 1989 июл; 8 (7): 795-802.

Hsieh FY, et al. Простой метод расчета размера выборки для линейной и логистической регрессии. Статистика в медицине. 1998 Jul 30; 17 (14): 1623-34.

Доступное обсуждение проблем с примерами расчетов можно найти в последней главе (раздел 8.5 с. 339-347) « Прикладной логистической регрессии» Hosmer & Lemeshow .

Thylacoleo
источник
7

Обычно мне легче и быстрее запускать симуляции. Бумаги долго читаются, понимают и, наконец, приходят к выводу, что они не применяются в особом случае, в котором они заинтересованы.

Поэтому я бы просто выбрал несколько предметов, смоделировал интересующий вас ковариат (распределенный так, как вы считаете, что он будет), имитировал хорошие / плохие результаты, основываясь на функциональной форме, которую вы поставили (пороговые эффекты ковариации? Нелинейность?) с минимальным (клинически) значительным размером эффекта, который вы хотели бы обнаружить, пропустите результат через свой анализ и посмотрите, найден ли эффект на вашей альфе. Повторите это 10000 раз и посмотрите, нашли ли вы эффект в 80% симуляций (или любую другую необходимую вам мощность). Отрегулируйте количество предметов, повторяйте до тех пор, пока не получите силу, которой вы довольны.

Это имеет то преимущество, что является очень общим, поэтому вы не ограничены конкретной функциональной формой или конкретным числом или распределением ковариат. Вы можете включить отсев, см. Комментарий chl выше, либо наугад, либо под влиянием ковариации или результата. По сути, вы заранее программируете анализ, который вы собираетесь выполнить для окончательной выборки, что иногда помогает мне сосредоточиться на дизайне исследования. И это легко сделать в R (векторизация!).

Стефан Коласса
источник
У вас есть работающий случай в R?
Фаррел
1
@Farrel - вот очень короткий сценарий, который предполагает [0,1] - равномерно распределенные ковариаты, ИЛИ = 2 между первым и третьим квартилем ковариатного и стандартного нормального шума, что приводит к мощности .34 для n = 100. Я бы поиграл с этим, чтобы увидеть, насколько все чувствительно к моим предположениям: пробеги <- 1000; nn <- 100; set.seed (2010); обнаружения <- репликация (n = прогоны, expr = {ковариация <- runif (nn); результат <- runif (nn) <1 / (1 + exp (-2 * log (2) * ковариация + rnorm (nn)) ); summary (glm (исход ~ ковариат, семейство = "биномиальное")) $ coefficients ["covariate", "Pr (> | z |)"] <.05}) cat ("Power:", sum (обнаружения) /
Run
1
Вы можете прикрепить свой код как pastie ( pastebin.com ) или Gist ( gist.github.com ), если считаете, что это более удобно, и сослаться на него в своем комментарии.
ЧЛ
@chl: +1, спасибо большое! Вот суть: gist.github.com/607968
Стефан Коласса
Отличный код, но есть проблема. Я не такой умный, как ты. Мне нужно, чтобы это было разбито по ступеням. Я так понимаю, работает ли количество симуляций? Что такое nn? Это количество предметов в исследовании? Затем я вижу, что вы создали распределение ковариат и заставили их определять да или нет в зависимости от порога.
Фаррел
4

После публикации Стефана Коласса (я не могу добавить это в качестве комментария), у меня есть альтернативный код для симуляции. Здесь используется та же базовая структура, но она разбита немного больше, так что, возможно, ее немного легче читать. Он также основан на коде Кляйнмана и Хортона для моделирования логистической регрессии.

nn - число в образце. Ковариата должна быть непрерывно нормально распределена и стандартизирована, чтобы означать 0 и sd 1. Мы используем rnorm (nn), чтобы сгенерировать это. Мы выбираем соотношение шансов и сохраняем его в odds.ratio. Мы также выбираем номер для перехвата. Выбор этого числа определяет, какая доля образца испытывает «событие» (например, 0,1, 0,4, 0,5). Вы должны играть с этим числом, пока не получите правильную пропорцию. Следующий код дает пропорцию 0,1 с размером выборки 950 и ИЛИ 1,5:

nn <- 950
runs <- 10000
intercept <- log(9)
odds.ratio <- 1.5
beta <- log(odds.ratio)
proportion  <-  replicate(
              n = runs,
              expr = {
                  xtest <- rnorm(nn)
                  linpred <- intercept + (xtest * beta)
                  prob <- exp(linpred)/(1 + exp(linpred))
                  runis <- runif(length(xtest),0,1)
                  ytest <- ifelse(runis < prob,1,0)
                  prop <- length(which(ytest <= 0.5))/length(ytest)
                  }
            )
summary(proportion)

Резюме (пропорция) подтверждает, что пропорция составляет ~ 0,1

Затем, используя те же переменные, мощность рассчитывается за 10000 прогонов:

result <-  replicate(
              n = runs,
              expr = {
                  xtest <- rnorm(nn)
                  linpred <- intercept + (xtest * beta)
                  prob <- exp(linpred)/(1 + exp(linpred))
                  runis <- runif(length(xtest),0,1)
                  ytest <- ifelse(runis < prob,1,0)
                  summary(model <- glm(ytest ~ xtest,  family = "binomial"))$coefficients[2,4] < .05
                  }
            )
print(sum(result)/runs)

Я думаю, что этот код верен - я проверил его по примерам, приведенным в Hsieh, 1998 (таблица 2), и, похоже, он согласен с тремя приведенными там примерами. Я также проверил это на примере 342 - 343 Хосмера и Лемешоу, где он нашел мощность 0,75 (по сравнению с 0,8 в Хосмере и Лемешоу). Так что может случиться так, что в некоторых обстоятельствах этот подход недооценивает власть. Однако, когда я запустил тот же пример в этом онлайн-калькуляторе , я обнаружил, что он согласен со мной, а не с результатами в Hosmer и Lemeshow.

Если кто-нибудь может сказать нам, почему это так, мне было бы интересно узнать.

Андрей
источник
У меня есть 2 вопроса, если вы не возражаете. 1) Функция пропорциональности просто для того, чтобы получить правильный перехват? 2) какова логика использования ytest (сравнение prob со случайным uni draw)?
B_Miner
@B_Miner 1) С другой стороны - чтобы получить правильную пропорцию, вам нужно правильно установить перехват - так что отрегулируйте перехват, пока не получите пропорцию, которую вы ожидаете. 2) Логика ytest заключается в том, что нам нужно получить дихотомический 0 или 1 результат. Таким образом, мы сравниваем каждую выборку из равномерного распределения с вероятностью (prob), чтобы получить наш дихотомический результат. 'Runis' не должен быть взят из случайного равномерного распределения - биномиальное или другое распределение может иметь больше смысла для ваших данных. Надеюсь, это поможет (извините за задержку с ответом).
Андрей
3

простой вопрос о размере выборки: каков размер выборки, чтобы получить 95% доверительный интервал не более 2d для [неизвестного] среднего значения распределения данных. Другой вариант: насколько большой образец должен иметь мощность 0,9 при при тестировании H . Похоже, вы не указали критерий выбора размера выборки.0 : θ = 0θ=10:θ=0

на самом деле, звучит так, что ваше исследование будет проводиться последовательно. в этом случае он может заплатить, чтобы сделать это явной частью эксперимента. последовательная выборка часто может быть более эффективной, чем эксперимент с фиксированным размером выборки [в среднем требуется меньше наблюдений].

Фаррел: я добавляю это в ответ на ваш комментарий.

чтобы получить размер выборки, обычно указывается какой-то критерий точности для оценки [такой как длина CI] ИЛИ мощности при определенной альтернативе теста, который должен быть выполнен на данных. Вы, кажется, упомянули оба этих критерия. в принципе, в этом нет ничего плохого: вам просто нужно сделать два вычисления размера выборки - одну для достижения желаемой точности оценки, а другую - для получения желаемой мощности при указанной альтернативе. тогда требуется больший из двух размеров выборки. [кстати - кроме того, что вы говорите о 80% мощности - вы, кажется, не упомянули, какой тест вы планируете провести - или альтернативу, при которой вы хотите 80% мощности.]

что касается использования последовательного анализа: если предметы включаются в исследование одновременно, то имеет смысл фиксированный размер выборки. но если предметов мало и далеко друг от друга, может потребоваться год или два [или больше], чтобы получить необходимое количество зачисленных. таким образом, испытание может продолжаться в течение трех или четырех лет [или более]. в этом случае последовательная схема дает возможность остановить раньше, чем это - если ожидаемый эффект становится статистически значимым в начале испытания.

ronaf
источник
Критериями будет разница в 10% в вероятности хорошего против плохого результата. Или, скажем, так как это будет логистическая регрессия, отношение шансов = 2. альфа = 0,05, мощность = 80%, я пока не знаю, что такое объединенная дисперсия для непрерывной переменной, но давайте предположим, что стандартное отклонение составляет 7 мм рт. Последовательный анализ был бы хорош, но окончательный результат - через два года после измерения.
Фаррел