Своими руками: как посчитать квази-биномиальную логистическую регрессию и зачем это вообще нужно

Всем привет!

Я вычислительный биолог, изучаю молекулярные механизмы старения и возрастных заболеваний и работаю с данными секвенирования одиночных клеток (сингл‑селл). Это очень интересный тип данных, поскольку он позволяет смотреть не на «кашу» из органа, а видеть индивидуальные популяции клеток. В частности, можно наблюдать популяции, которых становится больше в старении, или которые наоборот исчезают, или просто что‑то редкое и прикольное. В последние лет 5 сингл‑селл данные дают очень много инсайтов в науке, которые будут безумно полезны для биомедицины. Но как показать, что содержание той или иной популяции меняется при старении или болезни? Давайте разберем.

Итак, допустим, мы отковыряли куски тканей (сердце, легкие, мозг, сосуды, печень — подставьте что нравится) от 10 молодых (или здоровых) и 10 старых (или больных) доноров и хотим сравнить их между собой. Мы провели сингл‑селл, нашли интересующий тип клеток (например, макрофаги или фибробласты) и видим, что они бывают двух «видов»: 

• «Хорошие», которые экспрессируют нормальные гены и явно пребывают в нормальном, гомеостатическом состоянии; 

• «Плохие» — например, хронически воспаленные.

Это небольшое упрощение, но именно такую задачу, например, мне приходилось решать при изучении поведения микроглии (иммунных клеток мозга) при нейродегенеративных заболеваниях (болезни Альцгеймера, Паркинсона и так далее) и старении мозга. Для удобства я придумал данные (табличка на картинке), тут Y1-10 — это молодые доноры, O1-10 — это старые доноры, и указано число хороших/плохих клеток для каждого из них.

Тут можно даже глазами увидеть, что процент плохих как будто бы больше при старении или болезни. Но как это измерить и доказать?

И когда я только начинал работать с сингл‑селл, наивным и самым «логичным» способом сравнить группы мне показался Z‑тест на пропорции (почитать можно тут). Идея простая: суммарно у нас 353 + 461 + … = 2797 хороших клеток у молодых доноров, 2157 хороших клеток у старых доноров, 284 плохих клеток у молодых доноров и 1283 плохих клеток у старых доноров. То есть, если считать пропорции, то p_young = 284 / (2797 + 284) = 0.0922 и p_old = 1283 / (2157 + 1283) = 0.3730. То есть, из всех «молодых» клеток 9.22% — плохие, из всех «старых» — целых 37.30%. Просто проспойлерю, что в конце расчетов Z‑score получается 26.49, а p‑value = 1.24*10^-154. Запредельная достоверность!!! 

Только к сожалению, так делать нельзя. Это проблема носит название псевдорепликация: у нас 6521 клеток, но это псевдореплики — мы наскребли их со всего 20 реальных доноров, биологических повторов. И результаты получаются завышенные, ложноположительные. По противоположной причине тут не стоит использовать тест Манна‑Уитни, чтобы сравнить этих доноров — он слишком слабый. Грубо говоря, информацию «100 из 10000» он будет видеть как «1 из 100», хотя по информативности это небо и земля. Поэтому я посмотрел, как делают в статьях, и нашел, что сейчас такие штуки считают при помощи квази‑биномиальной логистической регрессии. Ну а идея этой статьи у меня родилась из того, что ни на русскоязычных, ни в англоязычных источниках не объяснено нормально, как её считать — обычно объясняют либо простые тесты (t‑критерий, Фишер), либо совсем уже ML модельки, либо это что‑то ультрасложное для математиков. 

Итак, как доказать, что плохих клеток становится больше при старении?

Шаг 1. Сформулируем уравнение модели

Наша модель — это регрессия. То есть функция вида y = kx + b.

В нашем случае y — это логарифм шансов, он же логит (поэтому регрессия называется логистическая), то есть y = ln(p / (1-p)), где p — теоретическая доля (от 0 до 1) плохих клеток для этого вида донора (молодого или старого). Мы строим линейный тренд не для самих долей, а для логарифма шансов, благодаря чему реальные проценты клеток будут красиво укладываться в плавную S‑образную кривую (сигмоиду) и никогда не вылетят за рамки от 0 до 1.

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

x бывает двух видов: 0 для молодых, 1 для старых.

b — это интерцепт, свободный член. Это базовый логарифм шансов найти плохую клетку у молодого донора, когда x = 0.

k — это коэффициент возраста. Он показывает, как изменяется этот логарифм шансов при старении или болезни.

Тогда функция для молодого донора (x = 0) будет иметь вид: ln( p_young / (1 — p_young) ) = b.

А для старого (x = 1) такой: ln( p_old / (1 — p_old) ) = b + k.

Шаг 2. Найдем коэффициенты k и b

Это сделать не очень сложно. Шансы (odds) считаются так: 

Odds_young = p_young / (1 — p_young) = 284 / 2797 = 0.1015. 

Odds_old = p_old / (1 — p_old) = 1283 / 2157 = 0.5948.

То есть тогда как у молодых на одну плохую клетку приходится почти 10 хороших, у старых — меньше 2 хороших. 

Далее b = ln(Odds_young) = -2.2873. Ну а k = ln(Odds_old) — b = -0.5185 + 2.2873 = 1.7678.

То есть, мы вывели тот самый закон природы, как % плохих клеток изменяется при старении. А именно: 

ln( p / (1-p) ) = -2.2873 + 1.7678 * x, где x = 0 для молодых и x = 1 для старых. Простыми словами, для молодого донора шансы найти старую клетку равны e^-2.2873 = 0.1015 (мы это уже считали), а при старении эти шансы растут в e^1.7678 = 5.8580 раз! 

Ну, осталось самое сложное – доказать статистическую значимость наших находок. 

Шаг 3. Расчет ожидаемых значений и хи-квадратов для всех образцов

Партия сказала надо, в следующем шаге поймете, зачем.

Для каждого донора мы считаем ожидаемое количество плохих и хороших клеток. Процент плохих клеток у молодых = p_young = 0.0922, процент хороших = 1 — p_young = 0.9078. У стариков процент плохих клеток p_old = 0.3730, процент хороших = 1 — p_old = 0.6270.

Для расчета ожидаемого количества плохих и хороших клеток мы общее число клеток каждого донора умножаем на его проценты, при этом до целого не округляем. То есть для первого молодого донора (Y1) мы ждем 370 * 0.0922 = 34.114 плохих и 370 * 0.9078 = 335.886 хороших клеток, для Y2 489 * 0.0922 = 45.086 плохих и 489 * 0.9078 = 443.914 хороших клеток; для первого старого (O1) 334 * 0.3730 = 124.582 плохих и 334 * 0.6270 = 209.418 хороших клеток. И так далее, масштабировать можно в экселе или похожих программах.

Далее считаем хи-квадраты. В данном случае для каждого образца χ² = ((реальные плохие — ожидаемые плохие)²) / (ожидаемые плохие * доля хороших в группе). То есть для Y1 будет χ² = ((17 — 34.114)²) / (34.114 * 0.9078) = 9.458, а для O1 верно χ² = ((115 — 124.582)²) / (124.582 * 0.6270) = 1.175. Опять-таки, масштабируется экселем.

Шаг 4. Оценка параметра сверхдисперсии

Теперь просуммируем получившиеся хи‑квадраты по всем 20 образцам. У меня получилось 682.093 по значениям экселя, если считать на калькуляторе, будет чуть по‑другому из‑за округления.

И прикинем число степеней свободы. Для этого надо из числа доноров (20) вычесть число параметров (их два, k и b), которые мы считаем. То есть, df = 20 — 2 = 18.

Из этого можно посчитать параметр (или коэффициент) сверхдисперсии φ = сумма хи‑квадратов / df = 682.093 / 18 = 37.894.

Что это вообще такое? Тут раскрывается смысл ещё одного магического термина «квази‑биномиальная». В обычной биномиальной регрессии дисперсия = (p * (1-p)) / n. Ну то есть если бы мы подбрасывали абсолютно сферических коней в вакууме, а все наши 20 доноров были бы идеальными биологическими роботами. В реальности же у кого‑то скрытый воспалительный процесс, у кого‑то генетика кривая, какой‑то образец просто на столе протух — в биологии всё сложно.

И поэтому из 682.093 единиц разброса только 18 (по центральной предельной теореме) может объяснить технический шум. Остальные — это биология, и наша модель, которая теперь называется квази‑биномиальная регрессия, определяет дисперсию как φ * (p * (1-p)) / n. И собственно, этот φ мы нашли.

Шаг 5. Раздувание стандартной ошибки и финальная оценка p-value

Доказать статистическую значимость в регрессии — это значит, доказать, что коэффициент b или k не равен нулю. И собственно, этим мы сейчас займемся. Сначала для b. 

Посчитаем информационный вес для молодых доноров: W_young = общее число клеток молодых доноров * процент плохих * процент хороших = 3081 * 0.0922 * 0.9078 = 257.81. Старые нам тут не нужны, так как b работает для x = 0. 

Теперь стандартная ошибка для b. SE_binomial(b) = sqrt( 1 / W_young ) = 0.0623. И вот теперь мы вспоминаем, что наши данные упоротые, имеют сверхдисперсию, и стандартную ошибку нам придется раздуть. Это плохо для стат. значимости, но зато очень хорошо для честности и адекватной оценки того, что происходит в образце. А поможет нам тут как раз коэффициент φ. SE_quasi(b) = SE_binomial(b) * sqrt(φ) = 0.3834.

Теперь считаем t‑критерий по формуле коэффициент делить на стандартную ошибку: t = b / SE_quasi(b) = -2.2873 / 0.3834 = -5.9658. Это много по модулю, скажу я Вам, и по таблице/экселю/питону t = -5.9658 при 18 степенях свободы дает p‑value = 1.2*10^-5.

То есть, b прошел порог статистической значимости. Если бы не прошел, то это значило бы, что логарифм шансов для молодого донора найти плохую и хорошую клетку достоверно не отличается от нуля, и p_young достоверно не отличается от 0.5. Ну а у нас отличается: всё‑таки у молодых доноров, как мы помним, только 10% клеток — плохие.

Теперь то же самое для k. Нам ещё понадобится W_old = общее число клеток старых доноров * процент плохих * процент хороших = 3440 * 0.3730 * 0.6270 = 804.49. W_young берем из предыдущих расчетов.

SE_binomial(k) = sqrt( [1 / W_young] + [1 / W_old] ) = 0.0716. Снова раздуваем ошибку: SE_quasi(k) = SEbinomial(k) * sqrt(φ) = 0.4406. 

Аналогично считаем t = k / SE_quasi(k) = 1.7678 / 0.4406 = 4.0127, что при 18 степенях свободы дает нам p-value = 8.177*10^-4. Снова достоверно! 

Это уже доказывает, что при старении действительно статистически значимо меняется процент плохих клеток. Но обратите внимание – p-value = 0.0008 это уже почти “на грани”, особенно если будут ещё какие-то поправки на множественное сравнение (ведь обычно мы тестируем не одну, а десятки популяций!).

Вот только если бы мы не раздували стандартные ошибки, а использовали SE_binomial(b) и SE_binomial(k), то для b мы бы получили p-value 2.68*10^-295, а для k — 1.02*10^-134. То есть как будто бы супер круто и биологически достоверно, но на деле — тотальный научный терроризм и так считать нельзя. А раздувание ошибок позволяет нам не получить ложноположительный результат и выдать его за реальную биологию.

Решение в гугл таблицах смотреть тут. В питоне — тут.

Если отрисовать в Seaborn, выглядеть это будет вот так:

Итоги и применение

Ну а теперь о том, где это применяется. Как и обещал в начале, приведу два хороших примера из мира сингл‑селл биологии.

В первой статье Раймондо и соавторы изучали механизм развития псориатического артрита, правда, у мышей. Это когда оно начинается на коже как псориаз, а дальше болезнь поражает суставы и развивается уже как артрит. В общем, выделили «плохую» популяцию CD2⁺MHC‑II⁺CCR2⁺ миелоидных предшественников (плюсики означают экспрессию соответствующих генов), которая увеличивается в коже, в крови и в синовиальной ткани сустава при болезни. В суставах эти ребята превращаются в макрофагов и начинают стимулировать Т лимфоциты производить всякие провоспалительные молекулы (интерлейкин 17, например), что и приводит к печально известным последствиям [1].

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

Во второй в лаборатории Манолиса Келлиса отсеквенировали больше миллиона ядер клеток мозга в здоровом состоянии и при болезни Альцгеймера. При помощи квази‑биномиальной регрессии обнаружили, что селективно умирают именно те нейроны, в которых активен сигнальный путь Reelin (а именно гены RELN и DAB1). То есть тут ситуация обратная: нашли «хорошую» популяцию, которая исчезает при заболевании [2].

Зачем нам эта информация? Путь Reelin — как минимум молекулярная сигнатура, а может быть, и причина уязвимости. Зная это, мы можем пытаться воздействовать на него, чтобы такие нейроны умирали медленнее, и пытаться тормозить развитие болезни Альцгеймера. Конечно, в составе комплексной терапии, но данные очень ценные.

P. S. Картинка в начале статьи собрана из этих двух публикаций.

Ссылки

[1] Raimondo, Maria G., Hashem Mohammadian, Mario R. Angeli, Stefano Alivernini, Vladyslav Fedorchenko, Kaiyue Huang, Richard Demmler et al. “Skin‑derived myeloid precursors and joint‑resident fibroblasts spread psoriatic disease from skin to joints.” Nature Immunology (2026): 1–13.

[2] Mathys, Hansruedi, Carles A. Boix, Leyla Anne Akay, Ziting Xia, Jose Davila‑Velderrain, Ayesha P. Ng, Xueqiao Jiang et al. “Single‑cell multiregion dissection of Alzheimer”s disease.’ Nature 632, no. 8026 (2024): 858–868.