Введение
В предыдущей статье [1] основное внимание было уделено процедурам определения при помощи дискриминантного анализа параметров электроэнцефалограммы (ЭЭГ), которые вносят наибольший вклад в дифференциацию шизофрении и нормы. Их поиск среди нескольких тысяч параметров осуществлялся путем получения линейных дискриминантных функций (ЛДФ) с наибольшей чувствительностью и специфичностью с тестированием на независимых экзаменующих выборках, которая служила оценкой ошибки классификации для всей генеральной совокупности. Далее рассматриваются более строгие процедуры оценки способности к правильной классификации для новых данных.
Согласно T. Hastie и соавт. [2], существуют две отдельные цели дискриминантного анализа: 1. Выбор модели: оценка эффективности различных моделей для того, чтобы выбрать лучшую. 2. Оценка модели: выбор окончательной модели, оценка ее ошибки прогнозирования (ошибки обобщения, generalization error) по новым данным.
Выбор модели (первая цель) включает в себя нахождение оптимального числа предикторов в ЛДФ и затем поиск тех переменных, которые вносят наибольший вклад в различение двух выборок. Для каждой задачи классификации существует оптимальная сложность модели (т.е. число предикторов), при которой достигается наилучшее качество различения [3, 4], увеличение сложности ведет к эффекту так называемого переобучения (overfitting, overtraining), когда слишком много предикторов дают чрезмерно хорошо подогнанную модель, которая учитывает помимо закономерностей генеральной совокупности и особенности (ошибки измерения, погрешности модели и пр.) в обучающей выборке, отсутствующие в других наборах данных.
Оценка модели (вторая цель) дает правило дискриминации на два класса с приемлемой «обобщающей способностью» (generalization performance) дискриминантной функции, которая определяется как способность к правильному предсказанию на независимых тестовых выборках [2].
При ограниченном объеме данных можно достичь только первой цели, что касается второй цели, то на подтверждение и тестирование обобщающей способности может потребоваться выделение значительной доли материала, по разным оценкам — от одной трети до половины.
Можно предложить различные процедуры отбора параметров ЭЭГ, например отбирать параметры ЭЭГ с наибольшим уровнем значимости различий при сравнении по критерию Стьюдента [5].
Одним из известных способов оценить сложность модели, перспективность дискриминантного анализа и, возможно, получить ЛДФ с приемлемой специфичностью и чувствительностью является метод пошагового включения предикторов [6, 7], в котором пошагово добавляют предиктор из множества исходных параметров ЭЭГ, дающих наилучшее качество различения по различным критериям, при этом поиск прекращается, если качество различения перестало улучшаться. Таким образом, число предикторов в полученных ЛДФ может дать оценку сложности модели, а полученные чувствительность и специфичность наряду с другими параметрами, например лямбдой Уилкса, дадут представление о перспективности дальнейших поисков. Также очень информативен такой показатель, как расстояние Махаланобиса.
Если известны параметры нормальных многомерных распределений обоих классов и для больших выборок [8], для теоретической оценки качества дискриминации в модели Фишера можно использовать квадрат расстояния Махаланобиса D2 (Mahalanobis Squared Distance) между центроидами (центрами распределений) классов в многомерном пространстве. Расстояние Махаланобиса является безразмерной величиной и учитывает коррелированность переменных, а также их единицы измерения и дисперсии [6, 9]. Квадрат расстояния Махаланобиса D2 между центроидами многомерных выборок в матричном виде:
D2=(x1–x2)TC–1(x1–x2),
где C–1 — объединенная выборочная ковариационная матрица, x1–x2 — векторы-столбцы выборочных средних переменных (координаты центроидов), T — знак транспонирования матрицы.
Случайная величина y (дискриминантные баллы) для наблюдений, относящихся к 1-му классу, имеет нормальное распределение с математическим ожиданием D2/2 и стандартным отклонением D (т.е. дисперсией D2), а случайная величина y для наблюдений, относящихся ко 2-му классу, имеет нормальное распределение с математическим ожиданием -D2/2 и также стандартным отклонением D. При этом ошибка классификации Err для обоих классов равна Φ(D/2) [8], где Φ — стандартная функция нормального распределения. Чем больше расстояние Махаланобиса между центроидами классов, тем меньше ошибка классификации. В табл. 1 приведены ошибки различения Err в процентах в зависимости от квадрата расстояния Махаланобиса D2. Таким образом, расстояние Махаланобиса дает представление о возможной ошибке классификации в идеальном случае. Из табл. 1 видно, что если квадрат расстояния Махаланобиса меньше 2, то ошибка разделения на два класса будет недопустимо большой.
Таблица 1. Теоретическая ошибка различения Err (в %) в зависимости от квадрата расстояния Махаланобиса D2 между центроидами выборок в модели Фишера
| D2 | Err | D2 | Err | D2 | Err | D2 | Err |
| 1 | 30,9% | 5 | 13,2% | 9 | 6,7% | 13 | 3,6% |
| 2 | 24,0% | 6 | 11,0% | 10 | 5,7% | 14 | 3,1% |
| 3 | 19,3% | 7 | 9,3% | 11 | 4,9% | 15 | 2,6% |
| 4 | 15,9% | 8 | 7,9% | 12 | 4,1% | 16 | 2,3% |
Также можно при помощи расстояния Махаланобиса оценить ошибку классификации E для нового наблюдения, которая равна вероятности принадлежности к другому классу. Согласно W. Ledermann и E. Lloyd [9], вероятность ошибки классификации нового наблюдения E равна:
E=Φ((|y|+D2/2)/D),
где Φ — стандартная функция нормального распределения, y — дискриминантные баллы, вычисленные при помощи ЛДФ, D — расстояние Махаланобиса между центроидами.
Для небольших выборок также можно воспользоваться этой формулой [9], но с известной осторожностью.
Геометрически вероятность ошибки можно интерпретировать как вероятность нахождения данного наблюдения из первого класса (с дискриминантными баллами Y) от центра распределения случайной величины другого класса. Если величина y близка к нулю (т.е. само наблюдение близко к разделяющей гиперплоскости), то вероятность ошибки достаточно велика.
Можно также выразить расстояние Y от центра распределения другого класса в «сигмах» (т.е. стандартных отклонениях):
S=D2/2+|y|/D,
где |y| — дискриминантные баллы наблюдения, S — расстояние от величины Y до центра распределения другого класса в «сигмах», D — расстояние Махаланобиса.
Как известно, вероятность выхода за пределы двух сигм для нормального распределения составляет 4,55%, трех сигм — 0,27%, четырех сигм — 0,0063%, а так как в нашем случае учитывается только одна сторона, то вдвое меньше. Таким образом, когда удовлетворены условия модели Фишера и известны параметры распределения, если квадрат расстояния Махаланобиса между центроидами D2>4, то вероятность ошибки менее 2,3%, если D2>6, то вероятность ошибки менее 0,14%. Заметим, что при любом распределении Y если D2>6, то, согласно неравенству Чебышева, вероятность ошибки менее 5,6%.
Надо, однако, учитывать, что теоретические оценки качества дискриминации ЛДФ и классификации нового наблюдения достаточно приблизительны в реальных условиях, когда на точность вычислений может влиять множество факторов — от отклонений от модели Фишера до ошибок измерений, артефактов и пр.
Метод пошагового включения предикторов может не дать наилучшую или хотя бы удовлетворительную комбинацию предикторов с точки зрения качества различения классов, для этого необходимо проверить все возможные сочетания, начиная с 2, 3, 4 и т.д. [6], а этот алгоритм требует больших вычислительных и временных ресурсов, поэтому применение его на обычных компьютерах довольно ограниченно. Если сложность модели не превышает 2—3 предикторов, то вполне можно перебрать все возможные комбинации переменных (до 10 000 в рассматриваемом программном комплексе Brainsys), однако для задач дискриминации, требующих получения более сложных моделей, число этих переменных необходимо уменьшить различными способами. Для уменьшения числа переменных в настоящей работе был предложен алгоритм, при котором во множество переменных, из которого затем формируются все возможные комбинации предикторов, включаются переменные, отобранные в процессе процедуры пошагового включения предикторов по нескольким критериям оценки наилучшей дискриминации.
Важно отметить, что лучшим руководством при выборе окончательной модели, которая может хорошо предсказывать новое наблюдение, является процент правильных классификаций, при этом нарушения математических допущений не являются критичными и не наносят большого ущерба [6, 10].
Оценка только по обучающим выборкам будет завышенной («переобучение», «overfitting»), потому что известны только выборочные параметры распределения и дискриминантная функция отражает свойства конкретной выборки более точно, чем свойства всей генеральной совокупности. Одним из способов подтверждения (validation) прогностической способности является метод разделения исходного набора данных на обучающие и подтверждающие выборки [6].
Наряду с методом разделения исходной выборки предложены другие различные методы оценки качества дискриминации. Приведем наиболее известные и простые с точки зрения реализации.
Если данных совсем немного, так что не хватает даже на методику разделения исходной выборки на обучающую и подтверждающую, то можно попробовать использовать скользящий экзамен (leave-one-out cross-validation, LOO CV) [2, 3, 10], который в качестве экзаменующего использует наблюдение, которое временно удаляется из исходной выборки и после пересчета ЛДФ возвращается. Так перебираются все наблюдения из обучающих выборок и вычисляются чувствительность и специфичность ЛДФ. Заметим, что для оценки обобщающей способности полученного правила дискриминации одного скользящего экзамена недостаточно, можно получить только результаты для данной выборки и, соответственно, эффект «переобучения».
Если же исследователь имеет достаточно много исходных данных, то наилучшим подходом, согласно [2], будет являться разделение случайным образом выборки на обучающую, подтверждающую и тестирующую, например в пропорции 50%, 25%, 25%, хотя общепринятых рекомендаций по пропорциям нет. Обучающие выборки (training set) используется для получения ЛДФ, подтверждающие выборки (validation set) — для оценки ошибки предсказания для выбранной модели и выбора наилучшей, тестирующие выборки (test set) — для окончательной оценки обобщающей способности ЛДФ. Тестирующая выборка должна использоваться только в самом конце анализа, иначе она покажет заниженную оценку ошибки обобщающей способности (generalization error) [2]. Разделение исходных данных на три набора (обучающий, подтверждающий и тестирующий, каждый из которых состоит из двух выборок) близко к подходу разбиения исходного материала только на два набора (обучающий и подтверждающий), изложенному W. Klecka [6], и при достаточно большом объеме материала выглядит более предпочтительным. Что касается величины обучающих выборок, то трудно найти общее правило, сколько данных достаточно для дискриминации, это зависит, например, от сложности применяемой модели [2] и других факторов.
При применении методики перебора всех комбинаций предикторов этот подход можно модифицировать, выделив дополнительный набор данных в виде дополнительных экзаменующих выборок, которые отличаются от подтверждающих и тестовых тем, что для них можно установить пороги чувствительности и специфичности. При переборе комбинаций учитываются только те, которые превышают порог чувствительности и специфичности как для обучающих выборок, так и для экзаменующих. Такой способ помогает выделить закономерности, общие для всей генеральной совокупности, уменьшает вычислительные ресурсы, весьма критичные при использовании метода перебора всех комбинаций. Таким образом, исходная выборка разбивается на обучающие, экзаменующие, подтверждающие и, наконец, тестовые. Обучающим и экзаменующим выборкам присваиваются пороги чувствительности и специфичности, а также, при необходимости, устанавливаются пороги квадрата расстояния Махаланобиса и других показателей качества различения, ниже которых ЛДФ с данной комбинацией предикторов не рассматриваются. Подтверждающая выборка служит для выбора лучшей ЛДФ по результатам правильной классификации на всех выборках, и только в самом конце лучшая ЛДФ проверяется на тестовых выборках, которые должны дать окончательную оценку ее обобщающей способности. Разделение на четыре набора данных при всей кажущейся избыточности, однако, может быть полезно в сложных случаях (таких, как в нашем случае получения дискриминантного правила различения нормы и шизофрении), когда процедура пошагового включения не дает удовлетворительных результатов и приходится прибегать к весьма затратному методу перебора всех комбинаций предикторов. Сама методика перебора всех комбинаций требует осторожности, и, если использовать только подтверждающие выборки, то окажется, что они сами участвуют в поиске, что, скорее всего, приведет к эффекту переобучения даже с использованием тестовой выборки.
Также в настоящей работе предложены и другие параметры ЭЭГ помимо использованных ранее нормализованных показателей мощности, относительной мощности и когерентности, которые оказались, в частности, информативнее с точки зрения задачи различения нормы и шизофрении и других сложных случаев. В программный комплекс Brainsys, в частности, были включены действительная и мнимая части когерентности, приведенные к распределению Гаусса при помощи Z-преобразования Фишера (т.е. arth — обратный гиперболический тангенс, арктангенс, ареатангенс), обычно используемому при приведении к нормальному виду линейного коэффициента корреляции Пирсона. Как правило, в научных работах используется так называемый квадрат модуля когерентности (или его квадратный корень — модуль когерентности, коэффициент когерентности), который представляет собой сумму квадратов действительной и мнимой частей комплексной функции когерентности, однако появились работы [11—13], в которых предлагается использовать только мнимую часть когерентности, чтобы избежать так называемого эффекта объемной проводимости. По аналогии можно использовать параметр в виде действительной части функции когерентности. Действительная часть когерентности, по существу, — обычный коэффициент корреляции Пирсона между двумя каналами ЭЭГ, отфильтрованный в заданной полосе частот. Мнимую же часть также можно рассматривать как коэффициент корреляции Пирсона между каналами на заданной частотной полосе, но при этом один из них сдвинут на π/2. В работе S. Valizadeh и соавт. [14] и R. Pascual-Marqui и соавт. [15] квадрат действительной части называется «мгновенная когерентность» (instantaneous coherence), а в первой из этих работ делается вывод о высокой результативности этого параметра при проведении дискриминантного анализа и других алгоритмов машинного обучения. Обозначим действительную и мнимую части функции когерентности, приведенные к нормальному распределению, как arth(R) и arth(I).
В компьютерный программный комплекс топографического картирования электрической активности мозга Brainsys [1, 5, 14] включены программные средства, в которых реализованы вышеуказанные методологические принципы и добавлен необходимый функционал, который позволил получить ЛДФ для дифференциации шизофрении и нормы с удовлетворительной обобщающей способностью.
Материал и методы
Отнесем больных шизофренией к 1-му классу, а обследуемых без психической патологии — ко 2-му классу. Важный момент при поиске ЛДФ состоит в том, что необходимо определить, какую гипотезу исследователю важно проверять с помощью ЛДФ. Так как в нашем случае стоит задача диагностировать пациента с подозрением на шизофрению, то необходимо ориентироваться на получение ЛДФ с наибольшим процентом правильных предсказаний случаев шизофрении у лиц с признаками этого заболевания (т.е. ЛДФ с наибольшей чувствительностью) и не стремиться добиваться высокой специфичности формулы (подтверждения, что здоровый человек здоров).
Были обследованы 220 больных шизофренией и 1400 человек без психической патологии, сопоставимые по основным характеристикам. Возраст больных шизофренией варьировал от 16 до 70 лет (средний возраст 38,6 года, стандартное отклонение 13,3 года), среди них были 163 мужчины и 57 женщин. Возраст обследуемых без психической патологии варьировал от 16 до 79 лет (средний возраст 39,5 года, стандартное отклонение 13,0 года), среди них были 942 мужчины и 458 женщин.
Для указанных групп обследуемых были получены ЛДФ с помощью программы Brainsys, разделяющие ЭЭГ на 2 класса: во всех случаях была записана ЭЭГ с использованием 19 отведений по системе «10—20» и ушных электродов в качестве референтных и проведен спектральный и когерентный анализ c эпохой 4 с при помощи программы Brainsys. Отбор больных шизофренией осуществляли по тем же критериям, что и в предыдущей статье [1].
На первом этапе данные были случайным образом разделены на обучающие (100/100 человек), экзаменующие (40/100), подтверждающие (40/100) и тестовые (40/1100) выборки. Разделение на 4 группы обусловлено сложным алгоритмом поиска ЛДФ и необходимостью получения адекватной оценки ее обобщающей способности (т.е. способности предсказания на генеральной совокупности).
На втором этапе были отобраны независимые переменные (параметры ЭЭГ) — кандидаты в предикторы. Они представляют собой приведенные к нормальному распределению все показатели мощности, относительной мощности, канонограммы (отношение мощностей стандартных диапазонов), когерентности, включая действительную и мнимую части комплексной функции когерентности, их асимметрии и некоторые другие для стандартных диапазонов дельта (0,5—4 Гц), тета (4—8 Гц), альфа (8—13 Гц), бета-1 (13—20 Гц), бета-2 (20—30 Гц). Более узкие полосы (1, 1,5, 2 Гц) не рассматривали, их включение предполагалось в случае неудачи поиска, Были также исключены из рассмотрения показатели, которые могут быть подвержены влиянию артефактов, например параметры дельта-диапазона для лобных отведений. Всего отобрано 1970 переменных, невзирая на уровень значимости различий при сравнении обучающих выборок, чтобы проверить на информативность все возможные варианты методом пошагового включения.
На третьем этапе для отобранных переменных была проведена процедура пошагового включения предикторов по нескольким критериям, которая, однако, не дала удовлетворительных результатов: чувствительность и специфичность — около 70—75%, расстояние Махаланобиса — менее 3,5. Пошаговый поиск остановился на 7 предикторах, что дало грубую оценку оптимального числа предикторов в ЛДФ — не менее 7. В дальнейшем 44 предиктора, которые фигурировали в полученных ЛДФ, были использованы в методе перебора всех возможных комбинаций, так как среди них наиболее вероятны полезные комбинации.
На четвертом этапе, начиная с 7 предикторов, была запущена процедура перебора всех возможных комбинаций из полученного на предыдущем этапе набора переменных, заданы пороги чувствительности для обучающей и экзаменующей групп 1-го класса — 90%, а для обучающей и экзаменующей групп 2-го класса — 75%.
На пятом этапе полученные результаты были проверены на подтверждающих выборках и выбрана наилучшая ЛДФ.
На шестом этапе подтверждена обобщающая способность ЛДФ на тестовой выборке.
В соответствии с задачами исследования были разработаны дополнительные модули дискриминантного анализа к программному комплексу Brainsys (при поддержке научно-медицинской фирмы «Статокин», Россия).
ЭЭГ регистрировали на установке, состоящей из 21-канального усилителя аппаратно-программного комплекса для топографического картирования электрической активности мозга «Нейро-КМ» научно-медицинской фирмы «Статокин» и персонального компьютера. Регистрацию биопотенциалов осуществляли по международной схеме 10—20% от 16 отведений: Fp1, Fp2, F3, F4, F7, F8, СЗ, С4, РЗ, Р4, ТЗ, Т4, Т5, Т6, О1, О2. В качестве референтных использовали ипсилатеральные ушные электроды.
Была получена ЛДФ для восьми предикторов с чувствительностью 90% и специфичностью 80%, уровнем значимости для лямбды Уилкса p<3,9E-28 и расстоянием Махаланобиса между центроидами обучающей выборки 4,6.
Данная формула имеет следующий вид:
Y=1,57·(LRP(T5)4-8)+1,727·(Ln(P)(T5)4-8)–2,244·(arth(R)(T5–O1)8-13)+0,368·(LCoh(C4-FP1)4-8)+5,469·(arth(I)(T3-O1)4-8)+7,674·(arth(I)(F8-F4)4-8)+11,55·(arth(I)(F8-F4)20-30)+2,075·(LRP(C3)20-30)+8,294.
Если Y≥0, то случай следует отнести к классу «шизофрения», если Y<0, то случай относят к классу «норма». Данная формула применима к записям ЭЭГ с ипсилатеральными ушными референтными электродами.
(LRP(T5)4-8) — параметр, вычисляемый по формуле: С=Ln (S/(100–S)), где S — отношение мощности тета-ритма (4—8 Гц) к суммарной мощности биопотенциалов в диапазоне 0,5—30 Гц в процентах на отведении T5;
(Ln(P)(T5)4-8) — натуральный логарифм спектральной мощности тета-ритма (4—8 Гц) в отведении T5;
(arth(R)(T5-O1)8-13) — приведенная к нормальному распределению при помощи Z-преобразования Фишера действительная часть когерентности между отведениями T5-O1 в частотном диапазоне 8—13 Гц (альфа-ритм);
(LCoh(C4-FP1)4-8) — нормализованная когерентность тета-ритма (4—8 Гц) между отведениями C4-FP1;
(arth(I)(T3-O1)4-8) — приведенная к нормальному распределению при помощи Z-преобразования Фишера мнимая часть когерентности между отведениями T3-O1 в частотном диапазоне 4—8 Гц (тета-ритм);
(arth(I)(F8-F4)4-8) — приведенная к нормальному распределению при помощи Z-преобразования Фишера мнимая часть когерентности между отведениями F8-F4 в частотном диапазоне 4—8 Гц (тета-ритм);
(arth(I)(F8-F4)20-30) — приведенная к нормальному распределению при помощи Z-преобразования Фишера мнимая часть когерентности между отведениями F8-F4 в частотном диапазоне 20—30 Гц (бета-2-ритм);
(LRP(C3)20-30) — параметр, вычисляемый по формуле: С=Ln (S/(100–S)), где S — отношение мощности бета-2-ритма (20—30 Гц) к суммарной мощности биопотенциалов в диапазоне 0,5—30 Гц в процентах на отведении C3.
Стандартизованные коэффициенты ЛДФ для различения на классы «Шизофрения» и «Норма» (табл. 2) показывают вклад каждого предиктора. Теоретическая оценка ошибки классификации, полученная по выборочному значению расстояния Махаланобиса (табл. 1), не противоречит полученным результатам.
Таблица 2. Стандартизованные коэффициенты ЛДФ для различения на классы «Шизофрения» и «Норма» и вклад предикторов
| Параметр | Стандартизованные коэффициенты | Вклад предиктора, % |
| (LRP(T5)4-8) | 0,146 | 15 |
| (Ln(P)(T5)4-8) | 0,217 | 22 |
| (arth(R)(T5-O1)8-13) | –0,121 | 12 |
| (LCoh(C4-FP1)4-8) | 0,077 | 8 |
| (arth(I)(T3-O1)4-8) | 0,071 | 7 |
| (arth(I)(F8-F4)4-8) | 0,0618 | 6 |
| (arth(I)(F8-F4)20-30) | 0,112 | 11 |
| (LRP(C3)20-30) | 0,193 | 19 |
Обсуждение
Как видно из табл. 2, большинство предикторов входят в формулу с положительным знаком, т.е. увеличивают вероятность отнесения ЭЭГ к классу «Шизофрения». Среди них наибольший вклад вносят показатели мощности в тета-диапазоне в левой задневисочной области (T5) и бета-2-диапазоне в левой центральной области (C3).
Также увеличивают вероятность отнесения ЭЭГ к классу «Шизофрения» показатели межполушарной когерентности в тета-диапазоне в передних отделах, в проекции C4-FP1 и мнимая часть внутриполушарной когерентности в правом полушарии в передневисочной области F8-F4 в частотных тета- и бета-2-диапазонах и в левом полушарии в затылочно-височной области 3-O1 в тета-диапазоне. Увеличивает вероятность отнесения ЭЭГ к классу «Норма» только один предиктор — arth(R)(T5-O1)8-13, действительная часть когерентности (или коэффициент корреляции между каналами) в альфа-диапазоне в отведениях T5-O1.
Обращает внимание, что при увеличении числа пациентов с 83 (в первой статье) до 220 и усовершенствовании методологии проведения дискриминантного анализа ЭЭГ в два класса с минимизацией эффекта переобучения, выбором оптимальной модели ЛДФ и оценки ее обобщающей способности получены предикторы, сходные с предыдущей формулой, а именно параметры мощности тета- и бета-2-ритмов в левом полушарии и межполушарной когерентности в тета-диапазоне. Такие изменения подтверждают неадекватный уровень внутри- и межполушарной интеграции биоэлектрических процессов преимущественно в частотных тета- и бета-диапазонах при шизофрении.
В результате дискриминантного анализа при помощи программного комплекса Brainsys была получена ЛДФ для различения больных шизофренией и здоровых со способностью предсказывать новые наблюдения (обобщающая способность), с чувствительностью 90% и специфичностью 80%. Для достоверности оценок обобщающей способности и исключения эффекта переобучения была применена методика разделения исходного массива данных на 4 набора: обучающий, экзаменующий, подтверждающий и тестирующий. Эта формула может быть встроена в программный комплекс Brainsys или поставляться в виде отдельной компьютерной программы в качестве помощника при определении диагноза «шизофрения» в спорных случаях и, возможно, с целью контроля лечения. Также реализована оценка вероятности ошибки классификации с использованием выборочного расстояния Махаланобиса обучающих выборок и величины дискриминантных баллов Y нового наблюдения.
Авторы заявляют об отсутствии конфликта интересов.