Как суммировать достоверные интервалы для медицинской аудитории

21

С помощью пакетов Stan и frontend, rstanarmили brmsя могу легко анализировать данные байесовским способом, как я делал раньше со смешанными моделями, такими как lme. Хотя у меня на столе лежит большая часть книг и статей Крушке-Гельмана-Вагенмейкера и т. Д., Они не говорят мне, как суммировать результаты для медицинской аудитории, разрываясь между Скиллой гнева Байесиана и Харибдой медицинских рецензентов ( «мы хотим значимости, а не этого рассеянного материала»).

Пример: желудочная частота (1 / мин) измеряется в трех группах; здоровые контроли являются эталоном. Есть несколько измерений для каждого участника, поэтому я часто использовал следующую смешанную модель lme:

summary(lme(freq_min~ group, random = ~1|study_id, data = mo))

Немного отредактированные результаты:

Fixed effects: freq_min ~ group 
                   Value Std.Error DF t-value p-value
(Intercept)        2.712    0.0804 70    33.7  0.0000
groupno_symptoms   0.353    0.1180 27     3.0  0.0058
groupwith_symptoms 0.195    0.1174 27     1.7  0.1086

Для простоты я буду использовать 2 * стандартную ошибку как 95% CI.

В частом контексте я бы суммировал это как:

  • В контрольной группе предполагаемая частота составляла 2,7 / мин (возможно, добавьте здесь CI, но иногда я избегаю этого из-за путаницы, создаваемой абсолютной и разностной CI).
  • В группе без симптомов частота была выше на 0,4 / мин, ДИ (от 0,11 до 0,59) / мин, р = 0,006 по сравнению с контролем.
  • В группе с симптомами частота была выше на 0,2 / мин, ДИ (от -0,04 до 0,4) / мин, р = 0,11, чем в контрольной группе.

Речь идет о максимально приемлемой сложности для медицинской публикации, рецензент, вероятно, попросит меня добавить «незначительный» во втором случае.

Здесь то же самое с stan_lmerи по умолчанию приоры.

freq_stan = stan_lmer(freq_min~ group + (1|study_id), data = mo)


           contrast lower_CredI frequency upper_CredI
        (Intercept)     2.58322     2.714       2.846
   groupno_symptoms     0.15579     0.346       0.535
 groupwith_symptoms    -0.00382     0.188       0.384

где CredI - 90% вероятные интервалы (см. виньетку rstanarm, почему 90% используется по умолчанию.)

Вопросов:

  • Как перевести вышеприведенное резюме в байесовский мир?
  • В какой степени требуется предварительное обсуждение? Я совершенно уверен, что статья вернется с обычным «субъективным допущением», когда я упомяну приоры; или, по крайней мере, с "без технической дискуссии, пожалуйста". Но все Байесовские власти требуют, чтобы толкование действовало только в контексте приоры.
  • Как я могу дать некоторую «значимость» суррогата в формулировке, не изменяя байесовским концепциям? Что-то вроде "достоверно другой" (уууу ...) или почти достоверно другой (буоха ..., звучит как "на краю значимости").

Джона Габри и Бен Гудрич (2016). rstanarm: Байесовское моделирование прикладной регрессии через Стэн. Версия пакета R 2.9.0-3. https://CRAN.R-project.org/package=rstanarm

Команда разработчиков Stan (2015). Stan: библиотека C ++ для вероятности и выборки, версия 2.8.0. URL http://mc-stan.org/ .

Пауль-Кристиан Бюркнер (2016). brms: модели байесовской регрессии с использованием Stan. Версия пакета R 0.8.0. https://CRAN.R-project.org/package=brms

Pinheiro J, Bates D, DebRoy S, Sarkar D и R Core Team (2016). nlme: линейные и нелинейные модели смешанных эффектов . Версия пакета R 3.1-124, http://CRAN.R-project.org/package=nlme>.

Дитер Менн
источник
1
У меня нет опыта работы с рецензентами / редакторами медицинских журналов, но, возможно, вы могли бы попытаться сказать, что существует нулевая вероятность того, что перехват отрицательный, нулевая вероятность того, что коэффициент в фиктивной переменной «без симптомов» отрицательный, и вероятность около 5% что коэффициент на фиктивной переменной «с симптомами» отрицателен. Вы можете определить примерно на 5% более точно, выполнив mean(as.matrix(freq_stan)[,"groupwith_symptoms"] < 0).
Бен Гудрич
Мы думали об этом, и 5% звучали хорошо; Исследователи переведут это значение как «значимость», но, поскольку они обычно неправильно понимают значение, они будут правы благодаря двойному отрицанию. «Нулевая вероятность», с другой стороны, убийца: вы бы это приняли? Может быть, <1 / Reff (p <0,001) будет приближением? Но опять же: когда я пишу р <ххх, я нахожусь в мире значимости.
Дитер Менн
Исправьте Reff для n_eff выше.
Дитер Менн
1
Лично я не назвал бы вероятность хвоста «вероятностью менее 1 при n_eff», поскольку n_eff относится к точности, с которой оценивается среднее значение. Возможно, вы могли бы запустить свои цепочки достаточно долго, чтобы получить 1 отрицательное значение для коэффициента, group_nosymptomsа затем сказать, что вероятность того, что он будет отрицательным, равна 1 / draws. Но для перехвата цепочка никогда не собирается уходить в отрицательную территорию для этих данных, поэтому я думаю, вы могли бы сказать, что вероятность меньше, чем 1 / draws.
Бен Гудрич
Я получил несколько полезных советов по включению p-значений для эксперта в области, но не для статистического эксперта-эксперта здесь: stats.stackexchange.com/questions/148649/… . Мы использовали p <минимум (n_eff всех параметров) в качестве консервативной верхней границы, когда p = 0.
stijn

Ответы:

16

Быстрые мысли:

1) Ключевой вопрос заключается в том, какой прикладной вопрос вы пытаетесь ответить для своей аудитории, потому что это определяет, какую информацию вы хотите получить из своего статистического анализа. В этом случае мне кажется, что вы хотите оценить величину различий между группами (или, возможно, величину соотношений групп, если эта мера более знакома вашей аудитории). Величина различий напрямую не указывается в анализе, который вы представили в вопросе. Но получить байесовский анализ просто так: вы хотите последующее распределение различий (или отношений). Затем из апостериорного распределения различий (или отношений) вы можете сделать прямое утверждение вероятности, такое как это:

«Наиболее вероятные различия в 95% находятся между [низким пределом ИЧР в 95%] и [высоким пределом ИЧР в 95%]» (здесь я использую интервал с высокой плотностью 95% [ИЧР] в качестве вероятного интервала, и потому что определение значений параметров наивысшей плотности, которые они обозначают как «наиболее достоверные»)

Аудитория медицинского журнала интуитивно и правильно поняла бы это утверждение, потому что это то, о чем обычно думает аудитория, это значение доверительного интервала (даже если это не означает доверительный интервал).

Как вы получаете различия (или отношения) от Стэна или JAGS? Только путем постобработки завершенной цепочки MCMC. На каждом этапе цепочки вычислите соответствующие различия (или соотношения), затем изучите последующее распределение различий (или соотношений). Примеры приведены в DBDA2E https://sites.google.com/site/doingbayesiandataanalysis/ для MCMC, как правило, на рисунке 7.9 (с. 177), для JAGS на рисунке 8.6 (с. 211) и для Stan в разделе 16.3 (с. 468) и т. Д.!

2) Если вы по традиции вынуждены делать заявления о том, отклоняется или нет нулевая разница, у вас есть два байесовских варианта.

2А) Один из вариантов - сделать вероятностные заявления относительно интервалов, близких к нулю, и их связи с ИЧР. Для этого вы устанавливаете область практической эквивалентности (ROPE) около нуля, которая является просто порогом принятия решения, подходящим для вашей прикладной области - насколько большая разница тривиально мала? Установление таких границ обычно выполняется, например, в клинических тестах на неполноценность. Если у вас есть мера «размера эффекта» в вашей области, могут быть соглашения для «маленького» размера эффекта, и пределы ROPE могут быть, скажем, половиной небольшого эффекта. Затем вы можете сделать прямые вероятностные заявления, такие как эти:

«Только 1,2% заднего распределения различий практически эквивалентно нулю»

и

«95% наиболее вероятных различий практически не эквивалентны нулю (т. Е. 95% ИЧР и ВЕРЕВКА не перекрываются), и поэтому мы отвергаем ноль». (обратите внимание на различие между утверждением вероятности и последующим распределением и последующим решением, основанным на этом утверждении)

Вы также можете принять нулевую разницу для практических целей, если все 95% наиболее вероятных значений практически эквивалентны нулю.

2B) Второй байесовский вариант - это проверка байесовской нулевой гипотезы. (Обратите внимание, что метод выше не былназывается «тестирование гипотез»!) Байесовское тестирование нулевых гипотез делает сравнение байесовской модели предыдущего распределения, которое предполагает, что разница может быть только нулевой, по сравнению с альтернативным предварительным распределением, которое предполагает, что разница может быть некоторой разбросанной областью возможностей. Результат такого сравнения моделей (обычно) очень сильно зависит от конкретного выбора альтернативного распределения, и поэтому необходимо сделать тщательное обоснование для выбора альтернативы. Лучше использовать как минимум умеренно информированные априоры как для нулевого, так и для альтернативного, чтобы сравнение моделей было действительно значимым. Обратите внимание, что сравнение моделей дает иную информацию, нежели оценка различий между группами, потому что сравнение моделей касается другого вопроса. Таким образом, даже при сравнении моделей,

Могут быть способы сделать тест байесовской нулевой гипотезы из вывода Stan / JAGS / MCMC, но я не знаю в этом случае. Например, можно попробовать приближение Сэвиджа-Дики к коэффициенту Байеса, но это будет зависеть от знания предыдущей плотности различий, что потребует некоторого математического анализа или некоторого дополнительного приближения MCMC из предыдущего.

Два метода принятия решения о нулевых значениях обсуждаются в гл. 12 из DBDA2E https://sites.google.com/site/doingbayesiandataanalysis/ . Но я действительно не хочу, чтобы это обсуждение было отвлечено дискуссией о «правильном» способе оценки нулевых значений; они просто разные и предоставляют разную информацию. Основным пунктом моего ответа является пункт 1 выше: посмотрите на последующее распределение различий между группами.

Джон К. Крушке
источник
3
Добро пожаловать на наш сайт! Здорово, что ты стал частью нашего сообщества!
Тим
Если вы хотите объединить свой аккаунт с этим stats.stackexchange.com/users/16592 (который, похоже, тоже ваш), вы можете сделать это автоматически через stats.stackexchange.com/contact .
говорит амеба: восстанови Монику
Вы можете выполнить проверку гипотезы, описанную здесь, используя brms. Смотрите: github.com/paul-buerkner/brms
bjw
3

Следуя SO этикету, это должно было быть написано как комментарий к @John K. Kruschke, но более длинные комментарии трудно структурировать. Сожалею.

  • @ Джон К. Крушке пишет: «Просто постобработкой законченной цепочки MCMC ...»

lower_CredIи upper_CredIв исходном посте были вычислены, как вы упомянули из полных цепочек MCMC, и только слегка переформатированы для лучшего сравнения с lmeвыходными данными. Хотя вы предпочитаете ИЧР, это простые квантили; с симметричным апостериором в этом примере это не имеет большого значения.

  • ВЕРЕВКА и размер эффекта

Я видел заявления в комитеты по этике, в которых статистическая мощность была рассчитана без указания предположения о величине эффекта. Даже для случая, когда нет никакого способа определить «клинически значимый эффект», трудно объяснить эту концепцию медицинским исследователям. Для испытаний, не связанных с неполноценностью, это немного проще, но они не так часто становятся предметом исследования.

Поэтому я совершенно уверен, что введение ROPES не будет приемлемым - другие предположения, люди не могут иметь в виду более одного числа. Байесовские факторы могут сработать, потому что есть только одно число, которое можно вернуть домой, как p-значения ранее

  • Приоры

Я удивлен, что ни @Джон К. Крушке, ни @Бен Гудрич из команды Стэна не упоминают приоры; большинство работ по этому вопросу требуют подробного обсуждения предварительной чувствительности при представлении результатов.

Было бы неплохо, если бы в следующем издании вашей книги - надеюсь, со Стэном - вы могли бы добавить блоки «Как опубликовать это (в нестатистическом документе) со 100 словами» для выбранных примеров. Когда я взял бы вашу главу 23.1 словом, типичный медицинский исследовательский документ был бы длиной в 100 страниц и цифрами ...

Дитер Менн
источник
* Главное было посмотреть на последующее распределение различий (между группами, между комбинациями групп). Это то, что нуждается в последующей обработке цепочки MCMC.
Джон К. Крушке
* ВЕРЕВКА: Вы «совершенно уверены, что ВЕРЕВКИ не будут приемлемы» и «это трудно объяснить концепцию медицинским исследователям». Тогда я не понимаю, как байесовские факторы будет легче объяснить или принять, поскольку байесовский фактор требует еще более детального объяснения и обоснования какого-то определенного порогового значения BF для принятия решения !! Мне кажется, вы предположили, что ваша аудитория постоянно окостеневала в частых рамках; если это так, просто используйте статистику по частоте или отправьте свою работу в более просвещенный журнал.
Джон К. Крушке
* Вы сильно преувеличиваете в отношении рекомендаций главы 23.1, которые на самом деле могут быть кратко изложены в небольшом тексте, особенно для простых моделей, которые вы здесь используете. Продолжение в следующем комментарии ...
Джон К. Крушке
1
(I) мотивировать использование байесовского - он дает вам богато информативные апостериорные распределения. (ii) Объясните модель и ее параметры, что легко в этом случае. (iii) Обоснуйте предыдущее - в этом случае опять тривиально, просто чтобы сказать, что вы использовали диффузные априоры, которые практически не влияют на заднюю часть. (Но НЕ, если вы используете байесовские факторы, для которых важен предыдущий.) (Iv) Сообщить о плавности цепочки MCMC - просто сказать, что ESS составляла около 10 000 для всех параметров и различий. Продолжение в следующем комментарии ...
Джон К. Крушке
1
(v) Интерпретировать апостериор: просто укажите центральную тенденцию (например, модус) апостериора и его 95% ИЧР для каждой разницы интересов. Это не так коротко, как твит, но это всего лишь пара абзацев.
Джон К. Крушке