4.3: Оцінка ставок з використанням максимальної ймовірності

Last updated
Save as PDF

Page ID: 4223

$ \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } $ $ \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} $$\newcommand{\id}{\mathrm{id}}$ $ \newcommand{\Span}{\mathrm{span}}$ $ \newcommand{\kernel}{\mathrm{null}\,}$ $ \newcommand{\range}{\mathrm{range}\,}$ $ \newcommand{\RealPart}{\mathrm{Re}}$ $ \newcommand{\ImaginaryPart}{\mathrm{Im}}$ $ \newcommand{\Argument}{\mathrm{Arg}}$ $ \newcommand{\norm}[1]{\| #1 \|}$ $ \newcommand{\inner}[2]{\langle #1, #2 \rangle}$ $ \newcommand{\Span}{\mathrm{span}}$ $\newcommand{\id}{\mathrm{id}}$ $ \newcommand{\Span}{\mathrm{span}}$ $ \newcommand{\kernel}{\mathrm{null}\,}$ $ \newcommand{\range}{\mathrm{range}\,}$ $ \newcommand{\RealPart}{\mathrm{Re}}$ $ \newcommand{\ImaginaryPart}{\mathrm{Im}}$ $ \newcommand{\Argument}{\mathrm{Arg}}$ $ \newcommand{\norm}[1]{\| #1 \|}$ $ \newcommand{\inner}[2]{\langle #1, #2 \rangle}$ $ \newcommand{\Span}{\mathrm{span}}$

Ми також можемо оцінити швидкість еволюції, знайшовши значення параметра максимальної правдоподібності для броунівської моделі руху, придатних до наших даних. Нагадаємо, що значення параметрів ML - це ті, які максимізують ймовірність даних, заданих нашою моделлю (див. Главу 2).

Ми вже знаємо, що за броунівською моделлю руху стани символів верхівки черпаються з багатовимірного нормального розподілу з матрицею дисперсійно-коваріації C, яка обчислюється на основі довжин гілок і топології філогенетичного дерева (див. Розділ 3). Ми можемо обчислити ймовірність отримання даних за нашою броунівською моделлю руху, використовуючи стандартну формулу ймовірності чергування з багатовимірного нормального розподілу:

\[ L(\mathbf{x} | \bar{z}(0), \sigma^2, \mathbf{C}) = \frac {e^{-1/2 (\mathbf{x}-\bar{z}(0) \mathbf{1})^\intercal (\sigma^2 \mathbf{C})^{-1} (\mathbf{x}-\bar{z}(0) \mathbf{1})}} {\sqrt{(2 \pi)^n det(\sigma^2 \mathbf{C})}} \label{4.5} \]

Тут наші параметри моделі - σ ² і$\bar{z}(0)$, значення кореневої ознаки. x - вектор n × 1 значень ознак для n видів верхівки в дереві, з видами в тому ж порядку, що і C, а 1 - вектор стовпців n × 1. Зверніть увагу, що (σ ² C) ⁻¹ - матриця, обернена матрицею σ ² C

Як приклад, за даними ссавців можна обчислити ймовірність для моделі зі значеннями параметрів σ ² = 1 і$\bar{z}(0) = 0$. Нам потрібно працювати з ln-ймовірностями (LnL), як тому, що значення тут настільки маленьке, так і для полегшення майбутніх розрахунків, так:

\[lnL(\mathbf{x} | \bar{z}(0), \sigma^2, \mathbf{C}) = -116.2.\]

Щоб знайти оцінки ML параметрів нашої моделі, нам потрібно знайти значення параметрів, які максимізують цю функцію. Один (не дуже ефективний) спосіб зробити це - обчислити ймовірність в широкому діапазоні значень параметрів. Потім можна візуалізувати отриману поверхню правдоподібності та визначити максимум функції правдоподібності. Наприклад, поверхня ймовірності для даних розміру тіла ссавців за даними броунівської моделі руху показана на малюнку 4.3. Зверніть увагу, що ця поверхня має пік близько σ ² = 0,09 і$\bar{z}(0) = 4$. Оглядаючи матрицю значень ML, знаходимо найвищу ln-правдоподібність (-78.05) при σ ² = 0,089 і$\bar{z}(0) = 4.65$.

Малюнок 4.3. Поверхня ймовірності еволюції маси тіла ссавців за даними Garland (1992). Зображення автора, може бути використано повторно за ліцензією CC-BY-4.0.

Описаний вище розрахунок неефективний, тому що доводиться розраховувати ймовірності при широкому діапазоні значень параметрів, далеких від оптимальних. Краща стратегія передбачає використання алгоритмів оптимізації, добре розвиненої області математичного аналізу (Nocedal and Wright 2006). Ці алгоритми відрізняються своїми деталями, але ми можемо проілюструвати, як вони працюють на загальному прикладі. Уявіть, що ви знаходитесь поблизу гори. Сент-Хеленс, і вам покладено завдання знайти вершину цієї гори. Туманно, але ви можете побачити область навколо ніг і мати точний висотомір. Одна зі стратегій - просто подивитися на схил гори, де ви стоїте, і піднятися в гору. Якщо схил крутий, ви, ймовірно, все ще далеко від вершини, і повинні швидко підніматися; якщо схил неглибокий, ви можете бути біля вершини гори. Може здатися очевидним, що це доставить вас до місцевої вершини, але, можливо, не найвищої вершини гори. Сент-Хеленс. Математичні схеми оптимізації також мають цю потенційну складність, але використовуйте деякі хитрощі, щоб стрибати в просторі параметрів і спробувати знайти загальний найвищий пік, коли вони піднімаються. Деталі реальних алгоритмів оптимізації виходять за рамки цієї книги; докладнішу інформацію див. Nocedal and Wright (2006).

Один простий приклад базується на методі оптимізації Ньютона [як реалізований, наприклад, функцією r nlm ()]. Ми можемо використовувати цей алгоритм, щоб швидко знайти точні оцінки ML ².

За допомогою алгоритмів оптимізації ми знаходимо$\hat{\sigma}_{ML}^2 = 0.08804487$ розв'язок ML при і$\hat{\bar{z}}(0) = 4.640571$, з l n L = −78.04942. Важливо, що рішення можна знайти лише за допомогою 10 обчислень ймовірності; це значення хороших алгоритмів оптимізації. Я побудував шлях через простір параметрів, взятий методом Ньютона при пошуку оптимального на малюнку 4.4. Зверніть увагу на дві речі: по-перше, що функція починається в якийсь момент і прямує вгору на поверхні ймовірності, поки не буде знайдений оптимальний; і по-друге, що цей розрахунок вимагає набагато менше кроків (і набагато менше часу), ніж обчислення ймовірності для широкого діапазону значень параметрів.

Малюнок 4.4. Поверхня ймовірності еволюції маси тіла ссавців за даними Garland (1992). Тут показаний шлях, пройдений алгоритмом оптимізації, щоб знайти пік поверхні ймовірності. Останні п'ять кроків цього десятикрокового алгоритму занадто близькі один до одного, щоб їх можна було побачити на цьому малюнку. Зображення автора, може бути використано повторно за ліцензією CC-BY-4.0.

Використання алгоритму оптимізації також має додаткову перевагу надання (приблизних) довірчих інтервалів для значень параметрів на основі Гессіана поверхні правдоподібності. Такий підхід передбачає, що форма поверхні ймовірності в безпосередній близькості від піку може бути наближена квадратичною функцією, і використовує кривизну цієї функції, визначену Гессіаном, для наближення стандартних похибок значень параметрів (Burnham and Anderson 2003). Якщо поверхня сильно пікована, ТО СЕС буде мало, в той час як якщо поверхня дуже широка, то СЕС будуть великими. Наприклад, поверхня ймовірності навколо значень ML для еволюції розміру тіла ссавців має Гессіан:

\[ H = \begin{bmatrix} -314.6 & -0.0026\\ -0.0026 & -0.99 \\ \end{bmatrix} \label{4.6}\]

Це дає стандартні помилки 0,13 (для $\ hat {\ sigma} _ {ML} ^2$) і 0.72 [для $\ hat {\ bar {z}} (0) $]. Якщо припустити, що похибка навколо цих оцінок є приблизно нормальною, ми можемо створити достовірні оцінки шляхом додавання та віднімання вдвічі більшої стандартної помилки. Потім отримано 95% CI 0,06 − 0,11 (для$\hat{\sigma}_{ML}^2$) та 3,22 − 6,06 [for$\hat{\bar{z}}(0)$].

Небезпека в алгоритмах оптимізації полягає в тому, що іноді можна застрягти на локальних піках. Більш складні алгоритми, що повторюються для декількох відправних точок, можуть допомогти вирішити цю проблему, але не потрібні для простого броунівського руху по дереву, як розглянуто тут. Чисельна оптимізація є складною проблемою в філогенетичних порівняльних методах, особливо для розробників програмного забезпечення.

У конкретному випадку підгонки броунівського руху до дерев виявляється, що навіть наш швидкий алгоритм оптимізації виявився непотрібним. При цьому оцінка максимальної ймовірності для кожного з цих двох параметрів може бути розрахована аналітично (O'Meara et al. 2006).

\[ \hat{\bar{z}}(0) = (\mathbf{1}^\intercal \mathbf{C}^{-1} \mathbf{1})^{-1} (\mathbf{1}^\intercal \mathbf{C}^{-1} \mathbf{x}) \label{4.7}\]

і:

\[ \hat{\sigma}_{ML}^2 = \frac {(\mathbf{x} - \hat{\bar{z}}(0) \mathbf{1})^\intercal \mathbf{C}^{-1} (\mathbf{x} - \hat{\bar{z}}(0) \mathbf{1})} {n} \label{4.8}\]

де n - кількість таксонів у дереві, C - матриця дисперсії-коваріації n × n при броунівському русі для верхівкових символів, заданих філогенетичним деревом, x - вектор значень ознак n × 1 для видів верхівок у дереві, 1 - вектор-стовпчик n × 1 одиниць, $\ hat {\ bar {z}} (0) $ - розрахунковий кореневий стан символу, а $\ hat {\ sigma} _ {ML} ^2$ - розрахункова чиста швидкість еволюції.

Застосовуючи такий підхід до розміру тіла ссавців, отримано оцінки, точно такі ж, як наші результати числової оптимізації: $\ hat {\ sigma} _ {ML} ^2 = 0.08$ і$\hat{\bar{z}}(0) = 4.64$.

Рівняння (4.8) є упередженим і буде послідовно оцінювати темпи еволюції, які є занадто малими; неупереджена версія, заснована на обмеженій максимальній ймовірності (REML) і використовується Garland (1992) та іншими:

\[ \hat{\sigma}_{REML}^2 = \frac {(\mathbf{x} - \hat{\bar{z}}(0) \mathbf{1})^\intercal \mathbf{C}^{-1} (\mathbf{x} - \hat{\bar{z}}(0) \mathbf{1})}{n-1} \label{4.9}\]

Ця корекція змінює нашу оцінку швидкості розмірів тіла ссавців від$\hat{\sigma}_{ML}^2 = 0.088$ до$\hat{\sigma}_{REML}^2 = 0.090$. Рівняння\ ref {4.8} точно ідентичне розрахунковій швидкості еволюції, розрахованої за допомогою середнього квадрата незалежного контрасту, описаного вище; тобто$\hat{\sigma}_{PIC}^2 = \hat{\sigma}_{REML}^2$. Насправді PIC є формулюванням моделі REML. «Обмежена» частина REML посилається на той факт, що ці методи обчислюють ймовірності на основі перетвореного набору даних, де ефект параметрів неприємності був видалений. В даному випадку параметром nuisance є розрахунковий стан кореня$\hat{\bar{z}}(0)$ ³.

Для прикладу розміру тіла ссавців ми можемо додатково вивчити різницю між REML та ML з точки зору статистичних інтервалів довіри, використовуючи ймовірності, засновані на контрастах. Знову ж таки, ми припускаємо, що всі контрасти взяті з нормального розподілу із середнім 0 та невідомою дисперсією. Якщо ми знову використаємо метод Ньютона для оптимізації, ми знайдемо максимальну ймовірність логу REML -10.3 at$\hat{\sigma}_{REML}^2 = 0.090$. Це повертає матрицю 1 × 1 для Гессіана зі значенням 2957,8, що відповідає SE 0,018. Це трохи більший SE відповідає 95% CI для$\hat{\sigma}_{REML}^2$ 0,05 − 0,13.

У контексті порівняльних методів REML має дві основні переваги. По-перше, PIC розглядають стан кореня дерева як параметр неприємності. Зазвичай ми маємо дуже мало інформації про цей кореневий стан, так що це може бути перевагою підходу REML. По-друге, PIC легко обчислити для дуже великих філогенетичних дерев, оскільки вони не вимагають будівництва (або інверсії!) будь-яких великих дисперсійно-коваріаційних матриць. Це важливо для великих філогенетичних дерев. Уявіть, що у нас було філогенетичне дерево всіх хребетних (~60 000 видів) і хотіли обчислити швидкість еволюції розмірів тіла. Щоб використовувати стандартну максимальну правдоподібність, ми повинні обчислити C, матрицю з 60 000 × 60 000 = 3,6 мільярда записів, і інвертувати її, щоб обчислити C ⁻¹. Для обчислення ПІК, навпаки, нам потрібно провести лише близько 120 000 операцій. На щастя, зараз існують алгоритми обрізки, щоб швидко розрахувати ймовірність великих дерев під різними моделями (див., наприклад, FitzJohn 2012; Freckleton 2012; і Ho and Ané 2014).