Pure Birth

Vamos iniciar simulando uma árvore sob um modelo simples de Pure Birth. Usaremos um valor alto para o número de espécies para que a estimativa das taxas seja precisa, mas posteriormente veremos como o tamanho da árvore impacta nessas estimativas.

set.seed(42)
sim.time <- 30
tree.pb <- sim.bd.age(age = sim.time, numbsim = 1, lambda = 0.2, mu = 0, complete = FALSE)[[1]]
obj <- ltt(tree.pb)
plotTree(tree.pb, color = makeTransparent("blue", alpha = 50), ftype = "off", add = TRUE, mar = par()$mar)

Se assumimos que o processo de diversificação pode ser descrito como um Processo de Poisson, então temos que a distribuição dos tempos de espera entre sucessivos eventos de especiação deve obedecer a uma distribuição exponencial. Vamos dar uma olhada nesses tempos de espera de nossa filogenia. Além disso, vamos plotar também as curvas esperadas por uma distribuição exponencial com diferentes valores de parâmetros.

wt.pb <- diff(sort(branching.times(tree.pb)))

hist(wt.pb, freq = FALSE, breaks = 100, main = "Histograma de Branching Times", xlab = "Branching Times", ylab = "Densidade")
curve(dexp(x, rate = (0.2 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[1], lwd = 2, lty = 2, add = TRUE)
curve(dexp(x, rate = (0.4 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[2], lwd = 2, lty = 2, add = TRUE)
curve(dexp(x, rate = (1 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[3], lwd = 2, lty = 2, add = TRUE)
curve(dexp(x, rate = (0.05 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[4], lwd = 2, lty = 2, add = TRUE)
curve(dexp(x, rate = (0.01 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[5], lwd = 2, lty = 2, add = TRUE)

Dando um zoom no gráfico, para identificarmos melhor as expectativas das distribuições teóricas com diferentes valores de parâmetro:

hist(wt.pb, freq = FALSE, breaks = 100, xlim = c(0, 2), ylim = c(0, 2), main = "Histograma de Branching Times", xlab = "Branching Times", ylab = "Densidade")
curve(dexp(x, rate = (0.2 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[1], lwd = 3, lty = 2, add = TRUE)
curve(dexp(x, rate = (0.4 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[2], lwd = 3, lty = 2, add = TRUE)
curve(dexp(x, rate = (1 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[3], lwd = 3, lty = 2, add = TRUE)
curve(dexp(x, rate = (0.05 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[4], lwd = 3, lty = 2, add = TRUE)
curve(dexp(x, rate = (0.01 * sum(tree.pb$edge.length))/max(wt.pb)), col = pal[5], lwd = 3, lty = 2, add = TRUE)

Notamos no exemplo acima que é difícil (pra não dizer impossível) estimar qual o valor de parâmetro que melhor se ajustaria à distribuição empírica. Sendo assim, precisamos de uma maneira de quantificar esse ajuste. Para isso, utilizaremos a função de verossimilhança e a usaremos para buscar qual valor de parâmetro provê o melhor ajuste (ou seja, minimiza/maximiza essa função).

## Criando a função de verossimilhança
llik.lambda <- function(lambda){
    -sum(dexp(wt.pb, rate = lambda, log = TRUE))
    }

É importante notar algumas características da função acima. Como mencionado na aula teórica, o produto das probabilidades associadas a cada um dos valores se torna um número muito pequeno muito rapidamente. Esses números costumam ser tão pequenos que fogem à precisão que a maioria dos computadores possui (que geralmente reside na casa de \(10^{-16}\)). Por isso, calculamos o logaritmo natural dessas probabilidades (mas poderia ser em qualquer outra base), e como algebricamente o logaritmo do produto de uma série de números é igual à soma dos logaritmos desses números, calculamos o somatório das log verossimilhanças ao invés do produtório. Por fim, como por definição probabilidades são números entre 0 e 1, os logaritmos dessas probabilidades serão obrigatoriamente negativos, e portanto utilizamos o oposto da soma (isso não é necessário, é só um costume).

Podemos usar a função que criamos acima para calcular a log Verossimilhança negativa de diferentes valores do parâmetro:

## Calculando a verossimilhança de alguns valores de lambda
llik.lambda(0.01)
llik.lambda(0.1)
llik.lambda(0.2)
llik.lambda(0.4)
llik.lambda(1)

Esse procedimento (de testar diferentes valores) é o correto, mas não é nem um pouco prático, especialmente quando usamos distribuições mais complexas (com um maior número de parâmetros). Precisamos então usar funções de otimização, que nada mais são que “robôs” que vão fazer o trabalho por nós de maneira muito mais precisa e eficiente. Um dos pacotes mais utilizados (ao menos em Ecologia) é o pacote bbmle, criado por Ben Bolker. Ele é autor de um ótimo livro sobre o uso de modelos em ecologia.

No nosso caso em específico, usaremos a função mle2 (de Maximum Likelihood Estimation). Essa função necessita de dois elementos: uma função a ser maximizada/minimizada e uma lista com os valores iniciais dos parâmetros a serem estimados. A escolha desses valores iniciais pode ser problemática (especialmente as que possuem muitos parâmetros ou superfícies de verossimilhança muito planas), porém felizmente esse não é nosso caso no momento. O otimizador aceita outros argumentos como o método de otimização, além de outros parâmetros para checagem de convergência por exemplo, mas caso queiram saber mais detalhes consultem o (bom) help da função.

Vamos então colocar o otimizador para funcionar. Passaremos para a função um valor inicial do parâmetro de 0.1, que é um valor de taxa na ordem de grandeza das taxas estimadas para árvores empíricas. Isso já é suficiente para que o nosso otimizador não fique preso em mínimos locais.

## Usando um otimizador para buscar o estimador de máxima verossimilhança de lambda
fit.pb <- mle2(llik.lambda, start = list(lambda = 0.1))

fit.pb

O valor estimado do parâmetro é bem diferente do que usamos para simular a árvore. O que o valor estimado no passo anterior representa? Se multiplicarmos esse valor pela idade da árvore simulada, o que obtemos? Pense em um modo de confirmar sua resposta à essa última pergunta.

Como vimos na aula teórica, ao tentarmos estimar as taxas de diversificação a partir de uma filogenia molecular ajustando os tempos de espera a uma distribuição exponencial, é necessário fazer um truque matemágico uma conversão para que a taxa estimada represente eventos por linhagem por milhão de anos.

fit.pb@coef * (30/sum(tree.pb$edge.length))

Efeito do tamanho da árvore na precisão das estimativas

Vamos agora avaliar brevemente como o tamanho da árvore (= número de espécies) afeta nossas estimativas. Para isso, vamos simular uma árvore com os mesmos parâmetros da árvore anterior, mas com uma menor duração.

set.seed(42)
sim.time.short <- 15
tree.pb.short <- sim.bd.age(age = sim.time.short, numbsim = 1, lambda = 0.2, mu = 0, complete = FALSE)[[1]]
obj <- ltt(tree.pb.short)
plotTree(tree.pb.short, color = makeTransparent("blue", alpha = 50), ftype = "off", add = TRUE, mar = par()$mar)

Vamos repetir os passos anteriores com a nova árvore:

wt.pb.short <- diff(sort(branching.times(tree.pb.short)))

llik.lambda.short <- function(lambda){
    -sum(dexp(wt.pb.short, rate = lambda, log = TRUE))
    }

fit.pb.short <- mle2(llik.lambda.short, start = list(lambda = 0.1))

fit.pb.short@coef * (15/sum(tree.pb.short$edge.length))

Pense em duas possíveis justificativas para as estimativas em árvores pequenas serem piores que para árvores maiores.

Registro Fóssil

Em um cenário de taxas constantes e assumindo um registro fóssil perfeito, estimar tanto as taxas de especiação quanto as de extinção segue o mesmo raciocínio mostrado acima, calculando os tempos de espera entre eventos sucessivos do mesmo tipo (entre especiações sucessivas ou entre extinções sucessivas). É possível repetir o procedimento acima usando os dados da filogenia completa (com as espécies extintas), tratando as espécies que se extinguiram como o registro fóssil do grupo (afinal, elas representam exatamente o registro fóssil perfeito do grupo). Caso você se sinta confortável, tente adaptar o código acima para estimar as taxas de extinção!

Filogenias moleculares

O primeiro passo é obviamente simular uma filogenia em um cenário de birth death. Vamos usar taxas relativamente altas de especiação e extinção afim de facilitar a visualização e explicação do raciocínio.

set.seed(42)
sim.time <- 20
tree.bd <- sim.bd.age(age = sim.time, numbsim = 1, lambda = 0.4, mu = 0.1, complete = FALSE, mrca = TRUE)[[1]]
obj <- ltt(tree.bd)
plotTree(tree.bd, color = makeTransparent("blue", alpha = 50), ftype = "off", add = TRUE, mar=par()$mar)

Como mencionado na aula teórica, os LTT plots de árvores resultantes de um processo birth death costumam apresentar o pull of the present, que é um aumento na inclinação da reta próximo ao presente. Esse pull of the present é mais acentuado quando a taxa de extinção é alta, como na simulação abaixo:

set.seed(42)
#sim.time <- 30
tree.bd.examp <- tree.bd
obj <- ltt(tree.bd.examp, col = 2, lwd = 4)
plotTree(tree.bd.examp, color = makeTransparent("blue", alpha = 50), ftype = "off", add = TRUE, mar=par()$mar)

Essa diferença de inclinação é causada pelas espécies que se originaram próximas ao presente e que ainda não tiveram tempo suficiente para se extinguirem. Desta forma os waiting times próximos ao presente são dominados primordialmente somente pela especiação. Sendo assim, podemos arbitrariamente (porém com palpites embasados) “quebrar” o LTT plot em duas partes: a primeira com uma inclinação menor que é proporcional à diversificação líquida (especiação - extinção), e a segunda com uma inclinação maior que é proporcional somente à especiação.

Vamos então quebrar a árvore nessas duas fases do processo. Visualmente é bastante difícil precisar onde exatamente o pull of the present começa. No entanto, sabemos que a longevidade média de uma espécie é o inverso da taxa de extinção (\(1/\mu\)). Para sermos bem conservadores e eliminar o período de transição entre os dois momentos, vamos pegar apenas os últimos 10% do pull of the present

No nosso caso, com a taxa de 0.1 que usamos na simulação, a longevidade média das linhagens é de 10 milhões de anos. Vamos então obter e ordenar os tempos de espera da filogenia, e depois separá-los em dois grupos: o grupo de eventos que acontece durante o pull of the present (último 1 milhão de anos)e os que acontecem anteriormente a ele.

## Obtendo os branching times
bt.bd <- sort(branching.times(tree.bd))

## Separando nos dois grupos
bt.potp <- bt.bd[which(bt.bd <= 1)]
bt.backg <- bt.bd[which(bt.bd > 1)]

O procedimento agora é similar ao que fizemos na seção anterior: calcularemos os tempos de espera, criaremos as funções de verossimilhança e usaremos um otimizador para estimar os parâmetros que maximizam essa função.

wt.potp <- diff(bt.potp)
wt.backg <- diff(bt.backg)

llik.potp <- function(lambda){
    -sum(dexp(wt.potp, rate = lambda, log = TRUE))
}

llik.backg <- function(netdiv){
    -sum(dexp(wt.backg, rate = netdiv, log = TRUE))
}

fit.potp <- mle2(llik.potp, start = list(lambda = 0.5))
fit.backg <- mle2(llik.backg, start = list(netdiv = 0.5))

Assim como no caso anterior, o valor de parâmetro estimado não representa exatamente a taxa que usamos para simular. A mesma transformação deve ser feita, porém respeitando os intervalos de tempo.

## Obtendo a soma das durações das espécies antes do pull of the present
tree.backg <- treeSlice(tree.bd, slice = sim.time - 1, orientation = "rootwards")

## Obtendo o tempo total das linhagens para a árvore completa
tot.time.full <- sum(tree.bd$edge.length)

## Obtendo o tempo total das linhagens antes do pull of the present
tot.time.backg <- sum(tree.backg$edge.length)

## Obtendo o tempo total das linhagens durante o pull of the present
tot.time.potp <- tot.time.full - tot.time.backg

## Convertendo os parâmetros estimados para as taxas
### Pull of the Present
fit.potp@coef * 1/tot.time.potp

### Background
fit.backg@coef * (sim.time - 1)/tot.time.backg

Compare os resultados acima com as estimativas de lambda obtidas na primeira seção. O que é possível dizer sobre o papel da extinção na estimativa das taxas em filogenias moleculares? Além disso, como que você espera que seja a precisão das estimativas das taxas para o registro fóssil?

Originação

\[\begin{equation} \frac{dD_N (t)}{dt} = \color{orange}{ -( \lambda + \mu ) D_N (t)} + \color{teal}{2 \lambda E(t) D_N (t)} \end{equation}\]

Onde:

• Probabilidade de nada acontecer em \(\delta t\) (parte laranja)

• Probabilidade de uma linhagem especiar em \(\delta t\) condicionada a uma das espécies descendentes se extinguir (parte azul)

Extinção

\[\begin{equation} \frac{dE (t)}{dt} = \color{red}{\mu} \color{orange}{ - (\mu + \lambda) E (t)} \color{teal}{ + 2 \lambda E(t)^2} \end{equation}\]

• Probabilidade de a linhagem se extinguir no \(\delta t\) (parte vermelha)

• Probabilidade de a linhagem sobreviver à esse \(\delta t\) mas se extinguir posteriormente (parte laranja)

• Probabilidade de a linhagem especiar no \(\delta t\) e ambas as linhagens descendentes se extinguirem (parte azul)

Expressão de Verossimilhança para a filogenia

\[\begin{equation} L (\tau) = (n - 1)! \frac{\lambda^{n-2} \bigg[ \displaystyle\prod_{k = 1}^{2n - 2} e^{(\lambda - \mu) (t_{k,b} - t_{k,t})} . \frac{(\lambda - (\lambda - \mu) e^{(\lambda - \mu)t_{k,t}})^2}{(\lambda - (\lambda - \mu) e^{(\lambda - \mu)t_{k,b}})^2} \bigg] }{[1 - E (t_{root})]^2} onde: E(t_{root}) = 1 - \frac{\lambda - \mu}{\lambda - (\lambda - \mu) e^{(\lambda - \mu) t_{root}}} \end{equation}\]