В связи с недавними волнениями на тему «Баги в программном обеспечении для МРТ-сканеров», мы попросили прокомментировать ситуацию Екатерину Печенкову, главного редактора «Российского журнала когнитивной науки» (http://www.cogjournal.ru) и руководителя группы фМРТ головного мозга Центра лучевой диагностики (Лечебно-реабилитационный центр, Москва).
Как многие, наверное, уже знают, недавно в PNAS вышла статья Андерса Эклунда, Томаса Николса и Ханса Кнутсона «Cluster failure: Why fMRI inferences for spatial extent have inflated false-positive rates» [1]. Эта статья наделала неожиданно много шума в прессе. Я всех призываю прочесть статью самостоятельно (имеет смысл читать также приложение к статье, где подробно изложена методика и собственно результаты [2]), но многие просили меня ее прокомментировать, поскольку тем, кто работает в другой области, в некоторых вещах, изложенных в статье, достаточно трудно разобраться.
Начну с некоторых предварительных сведений. Если разбираться с самого начала, то в ходе фМРТ мы регистрируем серию объемных изображений головного мозга. Регистрация одного измерения — объема головного мозга — занимает обычно 2−3 секунды. В течение нескольких минут мы регистрируем несколько десятков или сотен измерений. При обычной фМРТ, связанной с заданием, в это время происходят различные события — предъявления стимулов и ответы испытуемых.
Каждый объем головного мозга — это трехмерное изображение, состоящее из объемных элементов, которые называются вокселы. Это аналог с picture elements — пикселами — для двумерного компьютерного изображения. Каждый воксел в каждом функциональном объеме мозга имеет определенное значение яркости, которое можно выразить числом, а можно визуализировать на картинке. Яркость зависит от того, какая ткань находится в соответствующей части мозга, а также от той активности, которая происходит в нервной ткани. При увеличении активности происходят изменения в крови, которые на функциональном МРТ-изображения отображаются как небольшое изменение яркости (изменение МР-сигнала, полученное за счет изменений насыщенности крови кислородом; это называется BOLD-сигнал). При увеличении активности нервной ткани яркость увеличивается, при снижении метаболизма — уменьшается. Зная, в какие моменты времени испытуемым предъявлялись стимулы и когда они давали ответы, мы можем найти изменения яркости отдельных вокселов, коррелирующие с этими событиями.
Таким образом, яркость в каждом вокселе можно считать отдельной зависимой переменной в эксперименте. В упрощенном виде можно считать, что эффект, оказываемый независимыми переменными на каждую из таких зависимых переменных, оценивается с помощью статистического критерия. В результате на выходе мы получаем карту значения статистических критериев (например, значений t-критерия) по отдельным вокселам. Это и есть карта активации.
Поскольку зависимых переменных более одной, то результаты нужно подвергать поправке на множественные сравнения. Вообще-то более одной — это очень мягко сказано. Вокселов может быть несколько десятков тысяч и даже более сотни тысяч (например, 36 срезов по 64×64 воксела — это 147 456 вокселов), поэтому возникает проблема со слабо выраженными эффектами, каковых в современных когнитивных исследованиях большинство. Чтобы набрать необходимое количество измерений, нужно столько держать человека в томографе, что он перестанет выполнять задачу и просто-напросто заснет. Поэтому изыскиваются всякие способы сделать поправки на множественные сравнения менее консервативными и более рациональными.
Для функциональных изображений мозга обычно размер воксела составляет от 2x2x2 до 4x4x4 мм. При этом большие зоны занимают много вокселов, но часто нейроны внутри такой функциональной зоны ведут себя примерно одинаково на всем ее протяжении. Поэтому получается, что в соседних вокселах изменения яркости изображения, происходящие со временем, будут коррелировать. Предварительная обработка функциональных изображений перед созданием карт активации (так называемый препроцессинг) приводит к тому, что эта корреляция возрастает. По скольким именно вокселам будет наблюдаться корреляция, будет зависеть в числе прочего и от того, какой размер воксела задал экспериментатор. Поэтому поправка на множественные сравнения производится обычно не из расчета количества вокселов, а из расчета количества реальных независимых измерений, которые были произведены.
Поскольку это количество неизвестно, то оно оценивается с помощью специального математического аппарата, основанного на теории гауссовых случайных полей (random field theory, RFT). Эмпирически рассчитывается количество так называемых резелов (resolution elements) — таких элементарных единиц объема нервной ткани, внутри каждой из которых яркость вокселов изменяется согласованно. Для того, чтобы можно было корректно рассчитать количество резелов, на этапе препроцессинга необходимо подвергнуть данные пространственному сглаживанию, которое может производиться с разными параметрами.
Дальше можно ввести поправку Бонферрони на множественные сравнения, просто разделив требуемый уровень значимости на количество резелов. В реальности в большинстве статистических пакетов все происходит несколько сложнее и с более удовлетворительным для экспериментатора результатом. Однако исходно теория случайных полей применялась к данным нейровизуализации (сначала ПЭТ, а затем и фМРТ) именно для внесения поправки на множественные сравнения, основанной на FWER (групповой вероятности ошибки первого рода). Подробнее про все это в достаточно доступной форме можно почитать в главе, написанной Мэтью Бреттом и соавторами [3].
Кроме того, поправка на множественные сравнения часто вводится не на уровне вокселов, а на уровне кластеров — предполагается, что если мы получили несколько активированных вокселов, находящихся по соседству, то это скорее подлинная активация, а если несколько разрозненных вокселов в разных частях мозга, то это скорее шум (хотя шум также будет коррелировать в пространстве). Рассчитывается вероятность ошибки при нахождении активированного кластера заданного размера. Чем крупнее кластер, тем меньше вероятность ошибиться.
Все расчеты и поправки, о которых я только что сказала, производятся с помощью программного обеспечения для постобработки. Есть несколько популярных пакетов, основные среди которых — SPM, FSL и AFNI. Все эти пакеты имеют открытый исходный код. Работе этих пакетов и посвящена статья Eklund et al (2016). Начиная с получения самих функциональных изображений, программное обеспечение томографа никак в обработке не участвует. Поэтому если вы слышали что-то вроде «найден баг в программном обеспечении томографов» — это полная ерунда.
На этом мы заканчиваем с введением и переходим к разбору статьи.
Начну с того, что как минимум первый автор статьи имеет непосредственное отношение к разработке программного обеспечения для непараметрического анализа фМРТ-данных, которое называется BROCCOLI (см. статью 2014-го года [4]). И по результатам, которые докладывают Eklund et al. в статье 2016-го года, это программное обеспечение участвует во всех сравнениях и «выигрывает» по точности. Причем, как оказалось, Андерс Эклунд уже в течение пары последних лет в рассылке пользователей SPM анонсировал, что они с коллегами пишут статью, где BROCCOLI будет сравниваться с другими пакетами. Ни о каких революционных открытиях речь не шла. Держим это в голове и двигаемся дальше.
Обычно алгоритмы, разрабатываемые для статистической обработки функциональных изображений, тестируются на так называемых синтетических данных. Берутся изображения мозга, а изменения яркости, связанные с активацией мозга, моделируются путем искусственного случайного повышения или понижения значений. Плюс делаются «добавки», имитирующие шум. Схема может быть и более сложной, но в любом случае данные будут зависеть от представлений исследователя о характере реальных данных. Дальше смотрят, насколько хорошо тот или иной алгоритм отыскивает «загаданную» разработчиком активацию в таких синтетических данных.
Eklund et al. справедливо такой подход критикуют, так как в реальных данных могут быть спрятаны более сложные закономерности, которые не будут учтены при моделировании и создании синтетических данных. Они предлагают вместо этого использовать реальные данные для оценки уровня ошибки, которую дают различные методы обработки. Однако здесь тоже возникает проблема, т.к. непонятно, для какого набора реальных данных мы могли бы точно делать утверждения о наличии или отсутствии активации, вызванной тем или иным событием в тот или иной момент времени.
Авторы статьи предположили, что реальная активация будет отсутствовать в данных фМРТ покоя. Это функциональные объемы головного мозга, которые регистрируются в течение нескольких минут, когда испытуемый бодрствует, но при этом не занят никакой систематической мыслительной деятельностью (часто с закрытыми глазами). Все основные выводы статьи основаны на предположении, что или (а) систематическая активация тех или иных зон головного мозга в условиях фМРТ покоя отсутствует; или (б) в условиях фМРТ, связанной с задачей, происходит простая суммация спонтанной активности в условиях фМРТ покоя, и активации, связанной с задачей.
Разумеется, любой исследователь, работающий с фМРТ покоя, скажет, что первое предположение не соответствует действительности. В данных фМРТ покоя имеет место ритмическая спонтанная активность на низких частотах (порядка 0.1 Гц), благодаря которой в основном мы и выделяем сети покоя, в том числе знаменитую default mode network.
Второе предположение имеет под собой некоторые основания, в частности, были получены некоторые подтверждения этого предположения в работе Fair et al. (2007) [5], однако предполагать идентичность фоновых колебаний сигнала в активном состоянии и в покое было бы опрометчиво. Есть стойкое ощущение, что это утверждение приведет к чему-то сродни поискам выраженного альфа-ритма на ЭЭГ в процессе активного выполнения задачи на том основании, что альфа-ритм присутствует при закрытых глазах в покое.
Итак, Eklund et al. взяли данные фМРТ покоя около 500 человек из проекта Human Connectome Project (данные получены в нескольких странах на нескольких томографах). И они обработали эти данные в четырех пакетах (SPM, FSL, AFNI, BROCCOLI), как если бы это были данные фМРТ, связанной с задачей (например, предположили, что активное условие и покой чередовались каждые 10 секунд или каждые 30 секунд, а также как если бы присутствовал план фМРТ-эксперимента, связанный с событиями).
Поскольку предполагалось, что в заданные моменты времени никаких реальных стимулов не предъявлялось, вся полученная активация была сочтена ошибочной (ложно-положительной).
Как нетрудно догадаться, если мы ищем в данных фМРТ покоя события, происходящие раз в 60 секунд, мы их находим, т.к. это будет просто каждое пятое событие, происходящее на частоте 0.08 Гц или каждое шестое, происходящее на частоте 0.1Гц, а таких периодических колебаний в данных фМРТ покоя немало. Собственно, это прекрасно видно по предыдущей статьей Eklund et al. (2012) [6], где они ровно такие результаты и получили. С планом, связанным с событиями, все несколько сложнее, но там тоже иногда выделяются низкочастотные периодические колебания.
При подсчете статистики по отдельным вокселам SPM, FSL, AFNI и BROCCOLI детектируют эти низкочастотные колебания примерно одинаково плохо (т.е. никакой «активации» не выявляют, и авторы заявляют о низком уровне ошибки).
Теперь переходим к кластерам. Авторы берут три способа выделения кластеров и применения к ним поправок на множественные сравнения. При этом заявляют, что это три распространенных способа. Первый способ — отбор кластеров из вокселов, активированных на уровне p<0.01 без поправки на множественные сравнения и последующая поправка FWE на уровне кластеров. Второй способ — то же самое, но отбор производится из вокселов, активированных на уровне p<0.001. Третий способ — отбор из вокселов, активированных на уровне p<0.001, но без поправки на множественные сравнения на уровне кластеров, а просто отбираются все кластеры, в которые входят 10 и более вокселов (или объемом более 80 мм3).
Я не знаю, почему авторы считают первый способ распространенным. Я никогда в жизни, наверное, в реальном исследовании его не видела. Даже при рецензировании плохих статей. Третий способ тоже не особо распространен для групповых данных. На индивидуальных часто берут уровень значимости p<0.001 uncorr и k>5 (k — это размер кластера в вокселах), иногда в качестве порогового выбирают тот же объем кластера 80 мм³ (см., например, знаменитую статью про лосося [7]), но откуда 80 мм³ для групповых данных? Второй способ, примененный Eklund et al., действительно очень распространен, но гораздо чаще в этом случае применяется поправка FDR (доля ложных отклонений нуль-гипотез), чем FWE (групповая вероятность ошибки первого рода).
Теперь — внимание! — «ошибки» под 70% и более (те самые, растиражированные в прессе) получены в основном третьим способом. Для первого способа ошибки тоже достаточно высоки, но тут ничего удивительного. Им потому никто и не пользуется, что получаются неинформативные кластеры, занимающие половину объема мозга. Просто уже даже на уровне здравого смысла становится очевидным, что модель может в этом случае не работать, а в ряде пакетов и вовсе учитываются ограничения теории случайных полей для подобных случаев. Для второго же способа подсчета, например, для SPM «ошибки» в основном получились порядка 8%-10% вместо 5% расчетных с отдельными выбросами по вариантам экспериментального плана или подвыборкам.
То есть мы просто получаем еще один аргумент в пользу того, что нужно применять поправки на множественные сравнения.
Разумеется, все «ошибки» получены только для SPM, FSL и AFNI. BROCCOLI остается примерно в пределах заданных 5%.
Теперь про 40 тысяч статей. Откуда взялась эта цифра, я не поняла. Это прямо как сорок тысяч братьев из «Гамлета». Если имеются в виду все исследования с использованием фМРТ, проведенные за 25 лет, то утверждение про необходимость пересчета результатов 40 тысяч исследований некорректно: часть фМРТ-исследований вообще использует только анализ по зонам интереса, и их все сказанное выше вовсе не касается; из оставшейся части очень многие вводили поправки FWE или FDR на уровне вокселов и не делали поправку на уровне кластеров; наконец, из оставшихся часть применяла поправку на уровне кластеров, но FDR, а не FWE (хотя это не слишком принципиально, т.к. речь все равно идет о random field theory), или вовсе даже основанную на random permutations, и уже из оставшихся работ подавляющее большинство пользовалось достаточно безопасным вторым способом применения поправок по Eklund et al.
То есть если взять способы обработки данных, которыми мало кто пользуется, и данные, природу которых еще толком никто не понимает, и показать, к каким страшным последствиям это может привести, а затем написать броские выводы по статье, то можно проснуться автором единственного в мире непогрешимого пакета для обработки фМРТ. Вся эта история очень напоминает историю с книжкой Минского и Пейперта «Перцептроны» 1969-го года, благодаря которой почти на 20 лет остановилось развитие моделирования познания с помощью нейронных сетей.
Но, разумеется, вся статья не сводится к этому. Авторы проделали большую работу, и в статье есть целая россыпь рациональных зерен, с которыми дальше необходимо разбираться. Есть некая непонятная вещь с функцией пространственной автокорреляции (и вот она-то, похоже, и связана напрямую с томографами). Есть зависимость «ошибок» от применяемого пространственного сглаживания на этапе препроцессинга. Есть совершенно здравое замечание, что резелы не могут быть однородными и иметь одинаковый объем по всему мозгу (но это замечание уже неоднократно высказывалось ранее, и предложено даже как минимум одно софтверное решение [8]). Есть не менее здравое замечание, что параметрические пакеты не тестируют свои исходные предположения о нормальности распределения, и что неплохо было бы тестировать и шире использовать непараметрические методы. Но все это достаточно локальные рабочие вещи, которые действительно интересны специалистам внутри области и не должны привлекать чрезмерного внимания прессы.
[1] Eklund et al. (2016) Cluster failure: Why fMRI inferences for spatial extent have inflated false-positive rates http://www.pnas.org/content/113/28/7900.full.pdf
[2] Приложение к статье Eklund et al. http://www.pnas.org/content/suppl/2016/06/27/1 602 413 113.DCSupplemental/pnas.1 602 413 113.sapp.pdf
[3] Brett et al. (2003) An Introduction to Random Field Theory http://www.fil.ion.ucl.ac.uk/spm/doc/books/hbf2/pdfs/Ch14.pdf
[4] Eklund et al. (2014) BROCCOLI: Software for fast fMRI analysis on many-core CPUs and GPUs http://journal.frontiersin.org/article/10.3389/fninf.2014.24/full
[5] Fair et al. (2007). A method for using blocked and event-related fMRI data to study «resting state» functional connectivity. Neuroimage 35, 396−40 510.1016/j.neuroimage.2006.11.051 http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2563954/
[6] Eklund et al. (2012) Does parametric fMRI analysis with SPM yield valid results?—An empirical study of 1484 rest datasets http://www.sciencedirect.com/science/article/pii/S1053811912003825
[7] Bennett at al. (2011) Neural Correlates of Interspecies Perspective Taking in the Post-Mortem Atlantic Salmon: An Argument For Proper Multiple Comparisons Correction http://jpeelle.net/reprints/Bennett-2011-Neural_correlates_of_interspecies_perspective_taking.pdf
[8] Hayasaka et al. (2004) Nonstationary cluster-size inference with random field and permutation methods http://www.fil.ion.ucl.ac.uk/spm/doc/papers/HayasakaNichols_NonstCls.pdf
Примечание от TCTS, комментарии Эклунда и Николса (в частности, откуда взялось 40к публикаций, и насколько это корректно):
[9] Is the bulk of fMRI data questionable? http://retractionwatch.com/2016/07/12/is-the-bulk-of-fmri-data-questionable/
[10] Bibliometrics of Cluster Inference http://blogs.warwick.ac.uk/nichols/entry/bibliometrics_of_cluster/