Анализ мощности для порядковой логистической регрессии
12
Я ищу программу (в R или SAS или автономную, если бесплатно или по низкой цене), которая будет выполнять анализ мощности для порядковой логистической регрессии.
Я предпочитаю проводить силовой анализ за рамками моделирования. С предварительно упакованными пакетами я никогда не совсем уверен, какие предположения делаются.
Имитация мощности довольно проста (и доступна) с использованием R.
решить, как вы думаете, как должны выглядеть ваши данные и как вы будете их анализировать
написать функцию или набор выражений, которые будут имитировать данные для данного отношения и размера выборки и выполнять анализ (функция предпочтительна в том смысле, что вы можете преобразовать размер выборки и параметры в аргументы, чтобы упростить выбор различных значений). Функция или код должны возвращать значение p или другую статистику теста.
используйте replicateфункцию для запуска кода сверху несколько раз (обычно я начинаю примерно со 100 раз, чтобы понять, сколько времени это занимает, и получить правильную общую область, затем до 1000, а иногда и 10000 или 100000 для окончательные значения, которые я буду использовать). Доля раз, когда вы отвергли нулевую гипотезу, является силой.
Повторите вышеизложенное для другого набора условий.
+1, это очень надежный, универсальный подход. Я часто использовал это. Я хотел бы предложить еще одну особенность: вы можете сгенерировать данные для максимального значения вы бы рассмотрели, а затем подогнать модель под пропорции этих данных, последовательно подбирая 1-е n из них через равные промежутки времени до (например, n = 100 120, 140, 160, 180 и 200). Вместо сохранения значения p из каждого сгенерированного набора данных вы можете сохранить строку значений p. Усреднение по каждому столбцу дает вам быстрое и грубое представление о том, как меняется мощность без , и помогает быстро найти подходящее значение. N NNNN
gung - Восстановить Монику
2
@gung: ваш комментарий имеет смысл, не могли бы вы добавить свои коды, чтобы люди с опытом в R могли извлечь из этого пользу? спасибо
1
Я смотрю на это снова, и у меня есть пара вопросов: 1) Почему форма х на 1:10? 2) Как бы вы обобщили это для более чем 1 независимой переменной?
Питер Флом - Восстановить Монику
1
@PeterFlom, x должен был быть чем-то, поэтому я выбрал (произвольно), чтобы он был равномерным в диапазоне от 0 до 10, это также могло бы быть нормально, гамма и т. Д. Лучше всего было бы выбрать что-то, похожее на то, что мы ожидаем от реального х переменных, чтобы выглядеть. Чтобы использовать более 1 предикторную переменную, генерируйте их независимо (или из многовариантной нормали, связки и т. Д.), А затем просто включите их все в часть eta1, например eta1 <- beta0 + beta1*x1 + beta2*x2 + beta3*x3.
Грег Сноу,
1
replicatemeanα = 0,05
3
Я бы добавил еще одну вещь к ответу Сноу (и это относится к любому анализу мощности с помощью симуляции) - обратите внимание на то, ищете ли вы 1 или 2 хвостатых теста. В популярных программах, таких как G * Power, по умолчанию используется односторонний тест, и если вы пытаетесь проверить, соответствуют ли им симуляции (это всегда хорошая идея, когда вы изучаете, как это сделать), вам следует сначала проверить это.
Чтобы Сноу выполнял односторонний тест, я бы добавил параметр «tail» к входам функции и поместил что-то вроде этого в саму функцию:
#two-tail test
if (tail==2) fit$stats[5]
#one-tail test
if (tail==1){
if (fit$coefficients[5]>0) {
fit$stats[5]/2
} else 1
Версия с 1 хвостом в основном проверяет, является ли коэффициент положительным, а затем сокращает значение p пополам.
Помимо отличного примера Сноу, я полагаю, что вы также можете выполнить имитацию мощности путем повторной выборки из существующего набора данных, который имеет свой эффект. Не совсем самозагрузки, так как вы не выборки-с-заменой той же п , но та же идея.
Итак, вот пример: я провел небольшой самоэксперимент, который привел к положительной точечной оценке, но, поскольку он был небольшим, не был почти статистически значимым в порядковой логистической регрессии. С этой точной оценкой, какой большой n мне понадобится? Для различных возможных значений n я много раз генерировал набор данных, запускал порядковую логистическую регрессию и видел, как мало было p- значение:
library(boot)
library(rms)
npt <- read.csv("http://www.gwern.net/docs/nootropics/2013-gwern-noopept.csv")
newNoopeptPower <- function(dt, indices) {
d <- dt[sample(nrow(dt), n, replace=TRUE), ] # new dataset, possibly larger than the original
lmodel <- lrm(MP ~ Noopept + Magtein, data = d)
return(anova(lmodel)[7])
}
alpha <- 0.05
for (n in seq(from = 300, to = 600, by = 30)) {
bs <- boot(data=npt, statistic=newNoopeptPower, R=10000, parallel="multicore", ncpus=4)
print(c(n, sum(bs$t<=alpha)/length(bs$t)))
}
В этом случае при n = 600 мощность составила 32%. Не очень обнадеживает.
(Если мой подход к симуляции неверен, пожалуйста, кто-нибудь, скажите мне. Я собираюсь написать несколько медицинских статей, в которых обсуждается силовая симуляция для планирования клинических испытаний, но я не совсем уверен в своей точной реализации.)
Я до сих пор не уверен, как должно выглядеть моделирование с более (в частности, тремя) независимыми переменными. Я понимаю, что должен «включить их всех в часть eta1, например, eta1 <- beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3» (как упоминалось выше). Но я не знаю, как настроить остальные параметры в функции. Может ли кто-нибудь помочь мне с этим?
eta1 <- beta0 + beta1*x1 + beta2*x2 + beta3*x3
.replicate
mean
Я бы добавил еще одну вещь к ответу Сноу (и это относится к любому анализу мощности с помощью симуляции) - обратите внимание на то, ищете ли вы 1 или 2 хвостатых теста. В популярных программах, таких как G * Power, по умолчанию используется односторонний тест, и если вы пытаетесь проверить, соответствуют ли им симуляции (это всегда хорошая идея, когда вы изучаете, как это сделать), вам следует сначала проверить это.
Чтобы Сноу выполнял односторонний тест, я бы добавил параметр «tail» к входам функции и поместил что-то вроде этого в саму функцию:
Версия с 1 хвостом в основном проверяет, является ли коэффициент положительным, а затем сокращает значение p пополам.
источник
Помимо отличного примера Сноу, я полагаю, что вы также можете выполнить имитацию мощности путем повторной выборки из существующего набора данных, который имеет свой эффект. Не совсем самозагрузки, так как вы не выборки-с-заменой той же п , но та же идея.
Итак, вот пример: я провел небольшой самоэксперимент, который привел к положительной точечной оценке, но, поскольку он был небольшим, не был почти статистически значимым в порядковой логистической регрессии. С этой точной оценкой, какой большой n мне понадобится? Для различных возможных значений n я много раз генерировал набор данных, запускал порядковую логистическую регрессию и видел, как мало было p- значение:
С выходом (для меня):
В этом случае при n = 600 мощность составила 32%. Не очень обнадеживает.
(Если мой подход к симуляции неверен, пожалуйста, кто-нибудь, скажите мне. Я собираюсь написать несколько медицинских статей, в которых обсуждается силовая симуляция для планирования клинических испытаний, но я не совсем уверен в своей точной реализации.)
источник
Ссылаясь на первое моделирование (предложено Snow; /stats//a/22410/231675 ):
Я до сих пор не уверен, как должно выглядеть моделирование с более (в частности, тремя) независимыми переменными. Я понимаю, что должен «включить их всех в часть eta1, например, eta1 <- beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3» (как упоминалось выше). Но я не знаю, как настроить остальные параметры в функции. Может ли кто-нибудь помочь мне с этим?
источник