Математическая статистика и машинное обучение для наук о жизни

Линейная смешанная модель с нуля

Получите и закодируйте LMM с использованием максимального правдоподобия

Это восемнадцатая статья из колонки Математическая статистика и машинное обучение для наук о жизни, в которой я пытаюсь просто объяснить некоторые загадочные аналитические методы, используемые в биоинформатике и вычислительной биологии. Линейная смешанная модель (также называемая линейной моделью смешанных эффектов) широко используется в науках о жизни, есть много руководств, показывающих, как запускать модель в R, однако иногда неясно, как именно параметры случайных эффектов оптимизирован с помощью процедуры максимизации правдоподобия. В моем предыдущем посте Как работает смешанная линейная модель я представил концепции модели, и в этом руководстве мы выведем и закодируем линейную смешанную модель (LMM) с нуля, применяя максимального правдоподобия (ML), т. е. мы будем использовать простой R для кодирования LMM и сравнения результатов с результатами функций lmer и lme R. Цель этого руководства - объяснить LMM как у моей бабушки, подразумевая, что люди без математического образования должны понимать, что LMM делает под капотом.

Набор данных игрушек

Давайте рассмотрим игрушечный набор данных, который очень прост, но при этом сохраняет все необходимые элементы типичной установки для линейного смешанного моделирования (LMM). Предположим, у нас есть только 4 точки данных / выборки: 2 исходят от лица №1, а другие 2 - от лица №2 . Кроме того, 4 балла распределяются между двумя состояниями: без лечения и с лечением. Предположим, мы измеряем реакцию (Resp) каждого человека на лечение и хотели бы выяснить, привело ли лечение к значительному отклику участников исследования. Другими словами, мы стремимся внедрить что-то подобное парному t-критерию и оценить значимость лечения. Позже мы свяжем результаты LMM и парного t-теста и покажем, что они действительно идентичны. В наборе данных об игрушках 0 в столбце Лечение означает необработанный, а 1 означает обработанный. Во-первых, мы будем использовать наивную методику метода наименьших квадратов (OLS) линейную регрессию, которая не принимает во внимание взаимосвязь между точками данных.

Технически это работает, однако это не подходит, у нас проблема здесь. Обычная линейная регрессия методом наименьших квадратов (OLS) предполагает, что все наблюдения (точки данных на графике) являются независимыми, что должно приводить к некоррелированным и, следовательно, нормально распределенным остаткам. Однако мы знаем, что точки данных на графике принадлежат 2 людям, то есть 2 точкам для каждого человека. В принципе, мы можем подобрать линейную модель для каждого человека отдельно. Однако это тоже не подходит. У нас есть по два очка для каждого человека, поэтому их слишком мало, чтобы подходить для каждого человека. Кроме того, как мы видели ранее, индивидуальные припадки мало что говорят об общем / популяционном профиле, поскольку некоторые из них могут иметь противоположное поведение по сравнению с остальными индивидуальными припадками.

Напротив, если мы хотим сопоставить все четыре точки данных вместе, нам нужно будет как-то учесть тот факт, что они не являются независимыми, т.е. два из них принадлежат индивиду № 1, а два принадлежат индивиду № 2. Это можно сделать в рамках линейной смешанной модели (LMM) или парного теста, например парный t-критерий (параметрический) или знаковый ранговый критерий Вилкоксона (непараметрический).

Линейная смешанная модель с Lmer и Lme

Мы используем LMM, когда между наблюдениями нет независимости. В нашем случае наблюдения группируются внутри отдельных лиц. Давайте применим LMM с фиксированными эффектами для наклонов и пересечений и Случайные эффекты для пересечений, это приведет к добавлению члена (1 | Ind) в формулу Resp ~ Treat:

Здесь REML = FALSE просто означает, что мы используем традиционную оптимизацию максимального правдоподобия (ML), а не Ограниченное максимальное правдоподобие (мы поговорим о REML в другой раз). В разделе «Случайные эффекты» выходных данных lmer мы видим оценки для 2 параметров минимизации: остаточная дисперсия, соответствующая стандартному отклонению (Std.Dev.) 4,243, и случайные эффекты ( разделяется между людьми) дисперсия, связанная с точкой пересечения со стандартным отклонением 5,766. Точно так же в разделе «Фиксированные эффекты» выходных данных lmer мы можем видеть две оценки для: 1) точки пересечения, равной 6,5, и 2) наклона / обработки, равной 9. Таким образом, у нас есть 4 параметра оптимизации, соответствующие 4 точкам данных. Значения фиксированных эффектов имеют смысл, если мы посмотрим на самый первый рисунок в разделе «Набор данных игрушек» и поймем, что среднее из двух значений для необработанных образцов составляет (3 + 10) / 2 = 6.5, мы обозначим его как β 1, а среднее значение обработанных образцов составляет (6 + 25) / 2 = 15,5, обозначим его как β 2 . Последнее будет эквивалентно 6,5 + 9, то есть оценке фиксированного эффекта перехвата (= 6,5) плюс оценка фиксированного эффекта наклона (= 9). Здесь мы обращаем внимание на точные значения случайных и фиксированных эффектов, потому что мы собираемся воспроизвести их позже при выводе и кодировании LMM.

По умолчанию пакет lme4 R и функция lmer R не предоставляют меры статистической значимости, такие как p-значения, однако, если вы все же хотите иметь p-value вашего LMM fit, можно использовать функцию lme из пакета nlme R:

Опять же, здесь у нас есть случайные эффекты для перехвата (StdDev = 5,766264) и остатка (StdDev = 4,242649), а также фиксированные эффекты для перехвата (значение = 6,5) и наклона / обработки. (Значение = 9). Довольно интересно, что стандартные ошибки фиксированных эффектов и, следовательно, их t-значения (t-value = 1,5 для Treat) не совпадают между lmer и lme. Однако, если мы требуем REML = TRUE в функции lmer, статистика фиксированных эффектов, включая t-значения, будет идентична для lme и lmer, однако статистика случайных эффектов будет разной.

Это разница между подходами максимального правдоподобия (ML) и Ограниченное максимальное правдоподобие (REML), о которых мы поговорим в следующий раз.

Связь между LMM и парным T-тестом

Ранее мы говорили, что LMM - это более сложная форма простого парного t-критерия. Давайте продемонстрируем, что для нашего набора данных игрушек они действительно дают идентичные результаты. По пути мы также поймем техническую разницу между парными и непарными t-критериями. Давайте сначала запустим парный t-тест между обработанными и необработанными группами образцов с учетом отсутствия независимости между ними:

Мы можем видеть, что значение t = 1,5 и значение p = 0,3743, сообщаемое парным t-критерием, идентично значениям, полученным LMM с использованием функции nlme или lmer с REML = TRUE. Статистика парного t-критерия «среднее значение различий = 9» также согласуется с оценками фиксированного эффекта от lmer и nlme, помните, что у нас была оценка лечения = 9, которая была просто разницей между средними значениями для обработанных и необработанных. образцы.

Теперь, что именно делает парный t-тест? Идея парного t-теста состоит в том, чтобы сделать настройку похожей на одинарный t-тест, где значения в одной группе проверяются на значительное отклонение от нуля. , которое является своего рода средним значением второй группы. Другими словами, мы можем рассматривать парный t-критерий так, как если бы мы сдвигали точки пересечения индивидуальных совпадений (см. Самый первый рисунок) или средние значения необработанной группы до нуля. Проще говоря, это было бы эквивалентно вычитанию необработанных значений Resp из обработанных, т. Е. Преобразованию переменной Resp в Resp_std (стандартизованный ответ), как показано ниже, а затем выполнению непарного t- test на переменной Resp_std вместо Resp:

Мы наблюдаем, что значения Response стали 0 для Treat = 0, то есть для необработанной группы, в то время как значения Response для обработанной группы (Treat = 1) уменьшились на значения для необработанной группы. Затем мы просто использовали новую переменную Resp_std и запустили непарный t-тест, результат эквивалентен запуску парного t-теста для исходной переменной Resp. Таким образом, мы можем резюмировать, что LMM воспроизводит результат парного t-теста, но обеспечивает гораздо большую гибкость, например, не только два (как для t-теста), но и сравнение нескольких групп и т. Д.

Математика за линейной смешанной моделью

Давайте теперь попробуем вывести несколько уравнений LMM, используя наш игрушечный набор данных. Мы снова рассмотрим 4 точки данных и сделаем некоторые математические обозначения, учитывающие эффекты лечения, β, что не что иное, как фиксированные эффекты , и блочная структура u из-за кластеризации точек в пределах двух индивидов, что на самом деле вклад Случайные эффекты. Мы собираемся выразить координаты отклика (отклика) y через β и u параметры.

Здесь β 1 - это реакция людей в состоянии без лечения, а β 2 - это Ответ на Лечение. Можно также сказать, что β 1 - это среднее значение необработанных образцов, а β 2 - среднее значение обработанные образцы. Переменные u 1 и u 2 - это блочные переменные, учитывающие эффекты, характерные для индивидуума №1 и индивидуума №2. , соответственно. Наконец, ϵij ∼ N (0, σ ²) - это Остаточная ошибка. , т. е. ошибку, которую мы не можем смоделировать, и можем только попытаться минимизировать ее как цель задачи оптимизации Максимального правдоподобия. Следовательно, мы можем записать переменную ответа y как комбинацию параметров β , u, т.е. фиксированные и случайные эффекты и ϵ как уравнение. (1). В общем виде эту систему алгебраических уравнений можно переписать как Ур. (2), где индекс i = 1, 2 соответствует лечебным эффектам, а j = 1, 2 описывает индивидуальные эффекты. Мы также можем выразить эту систему уравнений в матричной форме Eq. (3). Таким образом, мы приходим к следующей известной матричной форме LMM, которая показана во всех учебниках, но не всегда должным образом объясняется, уравнение. (4).

Здесь X называется матрицей дизайна, а K называется блочной матрицей, он кодирует взаимосвязь между точки данных, то есть исходят ли они от связанных лиц или даже от одного и того же человека, как в нашем случае. Важно отметить, что лечение моделируется как фиксированный эффект, потому что уровни обработанный-необработанный исчерпывают все возможные результаты лечения. Напротив, блочная структура данных моделируется как случайный эффект, поскольку отдельные лица были выбраны из популяции и могут неправильно представлять всю популяцию людей. Другими словами, есть ошибка, связанная со случайными эффектами, например uj ∼ N (0, σs ²) , в то время как фиксированные аффекты считаются безошибочными. Например, секс обычно моделируется как фиксированный эффект, потому что обычно предполагается, что он имеет только два уровня (мужчины, женщины), в то время как групповые эффекты в науках о жизни должны моделироваться как случайные эффекты, поскольку потенциально дополнительные экспериментальные протоколы или лаборатории будут производить гораздо больше, то есть много уровней, систематических различий между образцами, которые затрудняют анализ данных. На практике можно подумать, что фиксированные эффекты не должны иметь много уровней, в то время как случайные эффекты обычно представляют собой многоуровневые категориальные переменные, где уровни представляют собой лишь выборку всех возможностей, но не все из них.

Перейдем к вычислению среднего и дисперсии точек данных Y. Поскольку и ошибка случайного эффекта, и остаточная ошибка происходят из нормального распределения с нулевым средним, а ненулевой компонент в E [Y] происходит из фиксированного эффекта, мы можем выразить ожидаемое значение Y как уравнение. (5). Затем дисперсия члена фиксированного эффекта равна нулю, поскольку предполагается, что фиксированные эффекты не содержат ошибок, поэтому для дисперсии Y мы получаем уравнение. (6).

Уравнение (6) получено с учетом того, что var (K u) = K * var (u) * K ^ T и var ( ϵ) = σ ² * I и var (u) = σ s² * I , где I - это единичная матрица размером 4 x 4. Здесь σ ² - это остаточная дисперсия (немоделированная / невосстановленная ошибка), а σ s ² - это дисперсия случайных эффектов перехвата (общая для всех точек данных). Матрица перед σ s ² называется матрицей родства и задается уравнением. (7). Матрица родства кодирует все родства между точками данных. Например, некоторые точки данных могут быть получены от генетически связанных людей, географических местоположений в непосредственной близости, некоторые точки данных могут быть получены из технических копий. Эти отношения закодированы в матрице родства. Таким образом, матрица дисперсии-ковариации точек данных из уравнения. (6) принимает окончательную форму уравнения. (8). Как только мы получили матрицу ковариации дисперсии, мы можем продолжить процедуру оптимизации функции максимального правдоподобия, которая требует ковариации дисперсии.

LMM из принципа максимального правдоподобия (ML)

Почему мы потратили так много времени на вывод ковариационной матрицы и какое отношение она имеет к линейной регрессии? Оказывается, что вся концепция подбора линейной модели, как и многие другие, если не все концепции традиционной статистики Frequentist, исходит из максимального правдоподобия (ML) принцип. Для этого нам нужно максимизировать многомерную гауссову функцию распределения по параметрам β 1, β 2, σs ² и σ ², уравнение . (9).

Здесь | Σ y | обозначает определитель ковариационно-дисперсионной матрицы. Мы видим, что обратная матрица и определитель ковариационной матрицы явно включены в функцию правдоподобия, поэтому нам пришлось вычислить ее выражение через дисперсию случайных эффектов σs ² и остаточная дисперсия σ ². Максимизация функции правдоподобия эквивалентна минимизации логарифмической функции правдоподобия, уравнение. (10).

Нам нужно будет выполнить утомительный символический вывод определителя матрицы ковариации дисперсии, обратной матрицы ковариации дисперсии и произведения обратной дисперсии- ковариационная матрица с Y - X β условиями. По моему опыту, это сложно сделать в R / Python, однако мы можем использовать Maple (или аналогично Mathematica или Matlab) для символьных вычислений и вывести выражения для определителя и обратная матрица ковариации дисперсии:

Используя Maple, мы можем получить определитель ковариационно-дисперсионной матрицы как уравнение. (11). Далее, последний член в формуле. (10) для логарифма правдоподобия принимает форму уравнения. (12).

Теперь мы готовы выполнить числовую минимизацию функции логарифма правдоподобия относительно β 1, β 2, σs ² и σ ² с использованием optim R:

Мы видим, что алгоритм минимизации успешно сходимся, так как мы получили сообщение «конвергенция = 0». В выходных данных σ = 4,242640687 - это стандартное остаточное отклонение, которое в точности воспроизводит результат из lme и lmer. (с REML = FALSE). По аналогии, σs = 5.766281297 - это общее стандартное отклонение, которое снова точно воспроизводит соответствующие выходные данные перехвата случайных эффектов от lme и lmer (с REML = ЛОЖЬ) функций. Как и ожидалось, фиксированный эффект β 1 = 6,5 представляет собой среднее значение для необработанных образцов в соответствии с оценкой фиксированного эффекта перехвата от lmer и lme. . Затем β 2 = 15,5 - это среднее значение обработанных образцов, которое представляет собой оценку фиксированного эффекта пересечения (= 6,5) плюс оценка фиксированного эффекта наклона / обработки (= 9). из функций lmer и lme R.

Отличная работа! Мы успешно воспроизвели выходные данные «Фиксированные эффекты» и «Случайные эффекты» из функций lmer / lme путем получения и кодирования смешанной линейной модели (LMM) с нуля!

Резюме

В этой статье мы узнали, как вывести и закодировать смешанную линейную модель (LMM) с фиксированными и случайными эффектами на набор данных игрушек. Мы рассмотрели связь между LMM и парным t-тестом и воспроизвели параметры фиксированных и случайных эффектов из функций lmer и lme R.

В комментариях ниже дайте мне знать, какие аналитические методы из наук о жизни кажутся вам особенно загадочными, и я постараюсь осветить их в будущих публикациях. Проверяйте коды из поста на моем Github. Подписывайтесь на меня на Medium Николай Осколков, в Twitter @NikolayOskolkov и подключайтесь в Linkedin. В следующем посте мы рассмотрим разницу между подходами максимального правдоподобия (ML) и Ограниченное максимальное правдоподобие (REML), следите за обновлениями.