Visualização de dados filogenéticos e fósseis no R - BIE5751

Nesse tutorial, vamos relembrar os principais pontos da primeira aula, porém agora observando o que acontece “por trás dos panos” da ferramenta.

O tutorial não contém nenhum muitos comentários sobre os códigos utilizados, pois a intenção aqui não é ser um tutorial de R, mas sim sobre macroevolução. Sendo assim, o tutorial foi escrito de maneira a simplesmente ser necessário copiar e colar os códigos e rodá-los no seu próprio computador. Qualquer dúvida sobre o código, pode perguntar para qualquer professor(a)/monitor(a).

IMPORTANTE: Todas as perguntas a serem respondidas nessa aula estarão em blocos destacados como esse, e as respostas devem vir acompanhadas de figuras que as justifiquem (sejam elas do próprio tutorial, print screens de valores, ou de sua autoria). As respostas devem ser enviadas por e-mail em formato PDF até o início da próxima aula.

Por fim, quase todos os gráficos criados nesse tutorial foram feitos usando o pacote ggplot2. Caso tenha interesse em entender um pouco melhor como ele funciona, a documentação é ótima, e também preparei um tutorial em português com as principais funções.

Sem mais enrolação delongas, vamos à prática!


Carregando os pacotes necessários

library("ape")
library("TreeSim")
library("paleotree")
library("ggplot2")
library("cowplot")
library("phytools")
library("plyr")
library("viridis")

Como vimos nas aulas teóricas, o número de espécies em um cenário de nascimento puro cresce no tempo segundo um modelo exponencial. No entanto, ao simularmos diferentes filogenias com as mesmas taxas e valores de tempo, notamos que há uma variação no número final de espécies obtidos. Isso se deve à estocasticidade do processo, dado que os modelos usados são modelos probabilísticos e não determinísticos.

O lado bom é que por serem modelos estocásticos bem definidos, é possível calcular as médias das características das filogenias resultantes e compará-las com os valores esperados segundo os modelos. É isso que faremos no primeiro item.

Número esperado de espécies

O crescimento exponencial tem a seguinte forma:

\[\begin{equation} N(t) = N_{0}(t) * e^{\alpha t} \end{equation}\]

onde \(\alpha\) é o valor da taxa de acúmulo de espécies. Dependendo da função utilizada, a simulação pode ser iniciada com uma ou duas linhagens, então é importante consultar o help da função sempre que for usá-la.

Podemos então estimar quantas espécies são esperadas após um determinado tempo para uma determinada taxa simplesmente substituindo esses valores diretamente na fórmula:

## Número esperado para uma taxa de 0.2 e um tempo de 15 milhões de anos

expect.N.15 <- 1 * exp(0.2 * 15)
expect.N.15

## Número esperado para uma taxa de 0.2 e um tempo de 30 milhões de anos

expect.N.30 <- 1 * exp(0.2 * 30)
expect.N.30

Vale notar que de maneira similar ao que vimos usando a ferramenta de visualização, ao dobrar o tempo o número de espécies esperado ao fim da simulação aumenta não-linearmente. O mesmo vale comparando simulações que rodam para um mesmo valor de tempo porém com taxas diferentes:

expect.N.02 <- 1 * exp(0.2 * 15)
expect.N.02

expect.N.04 <- 1 * exp(0.4 * 15)
expect.N.04

Mesmo simulando várias vezes o processo é praticamente impossível obter exatamente os valores (arredondados claro) esperados segundo o modelo. Porém, o modelo nos diz a expectativa média do número de espécies, então se pudermos simular várias árvores e calcular a média do número de espécies, provavelmente esse valor se aproximará do número de espécies esperado para o valor taxa utilizado.

Vamos então simular 1000 filogenias com os mesmos valores de taxa e tempo, e depois calcular o número médio de espécies desse conjunto de filogenias. Os comandos abaixo demoram alguns segundos para terminarem de rodar, então não se preocupem.

phylo.list.15 <- sim.bd.age(age = 15, numbsim = 1000, lambda = 0.2, mu = 0, complete = TRUE)
ntip.15 <- unlist(lapply(phylo.list.15, function(x){ifelse(class(x) == "phylo", Ntip(x), 1)}))

phylo.list.30 <- sim.bd.age(age = 30, numbsim = 1000, lambda = 0.2, mu = 0, complete = TRUE)
ntip.30 <- unlist(lapply(phylo.list.30, function(x){ifelse(class(x) == "phylo", Ntip(x), 1)}))

Comparando os valores médios obtidos com os esperados:

expect.N.15
mean(ntip.15)

expect.N.30
mean(ntip.30)

Podemos analisar visualmente também:

data.plot <- data.frame(ntip15 = ntip.15, ntip30 = ntip.30)

ggplot(data = data.plot) +
    geom_histogram(aes(x = ntip15)) +
    geom_vline(aes(xintercept = expect.N.15, colour = factor(expect.N.15)), size = 1, linetype = "dashed") +
    geom_vline(aes(xintercept = mean(ntip15), colour = factor(mean(ntip15))), size = 1, linetype = "dashed") +
    labs(x = "Número de Espécies", y = "Frequência", colour = "Estimativas") +
    scale_colour_brewer(palette = "Set1", labels = c("Teórico", "Simulado")) +
    theme(legend.position = "bottom")
    
ggplot(data = data.plot) +
    geom_histogram(aes(x = ntip30)) +
    geom_vline(aes(xintercept = expect.N.30, colour = factor(expect.N.30)), size = 1, linetype = "dashed") +
    geom_vline(aes(xintercept = mean(ntip30), colour = factor(mean(ntip30))), size = 1, linetype = "dashed") +
    labs(x = "Número de Espécies", y = "Frequência", colour = "Estimativas") +
    scale_colour_brewer(palette = "Set1", labels = c("Teórico", "Simulado")) +
    theme(legend.position = "bottom")

Notamos que a média dos valores é muito próxima do número esperado segundo nosso modelo.


Número esperado em um cenário com extinção

Como dito acima, um dos parâmetros do modelo exponencial é a taxa de acúmulo de espécies ao longo do tempo. Esta taxa de acúmulo foi inicialmente representada apenas no argumento lambda (letra grega usada para representar taxas de especiação em filogenias moleculares) nas funções acima, dado que o argumento mu (letra grega usada para representar taxas de extinção em filogenias moleculares) foi ajustado para 0. Portanto, o “saldo” entre ganho e perda de espécies, chamado de taxa líquida de diversificação (net diversification rate em inglês), é determinado apenas pela especiação.

Ajustando ligeiramente nossas simulações acima, vamos ver como o número esperado de espécies se comporta quando a taxa de extinção não é nula. Novamente, as simulações podem demorar alguns poucos minutos para terminar de rodar.

## Número esperado para uma taxa líquida de 0.2 e um tempo de 15 milhões de anos

expect.N.bd.15 <- 1 * exp((0.4 - 0.2) * 15)
expect.N.bd.15

## Número esperado para uma taxa líquida de 0.2 e um tempo de 30 milhões de anos

expect.N.bd.30 <- 1 * exp((0.4 - 0.2) * 30)
expect.N.bd.30

# Simularemos apenas 200 árvores no próximo passo pois é uma etapa lenta computacionalmente
## Simular 200 árvores com especiação = 0.4, extinção = 0.2, tempo = 15
phylo.list.bd.15 <- sim.bd.age(age = 15, numbsim = 200, lambda = 0.4, mu = 0.2, complete = TRUE)

## Obter o número de espécies para as filogenias moleculares que representam o momento final da simulação
ntip.bd.15 <- unlist(lapply(phylo.list.bd.15, function(x){ifelse(class(x) == "phylo", Ntip(drop.fossil(x)), 0)}))

## Simular 200 árvores com especiação = 0.4, extinção = 0.2, tempo = 30
phylo.list.bd.30 <- sim.bd.age(age = 30, numbsim = 200, lambda = 0.4, mu = 0.2, complete = TRUE)

## Obter o número de espécies para as filogenias moleculares que representam o momento final da simulação
ntip.bd.30 <- unlist(lapply(phylo.list.bd.30, function(x){ifelse(class(x) == "phylo", Ntip(drop.fossil(x)), 0)}))
expect.N.bd.15
## [1] 20.08554
mean(ntip.bd.15, na.rm = TRUE)
## [1] 20.945

expect.N.bd.30
## [1] 403.4288
mean(ntip.bd.30, na.rm = TRUE)
## [1] 448.965

E graficamente:

data.plot.bd <- data.frame(ntip15 = ntip.bd.15, ntip30 = ntip.bd.30)

ggplot(data = data.plot.bd) +
    geom_histogram(aes(x = ntip15)) +
    geom_vline(aes(xintercept = expect.N.15, colour = factor(expect.N.15)), size = 1, linetype = "dashed") +
    geom_vline(aes(xintercept = mean(ntip15), colour = factor(mean(ntip15))), size = 1, linetype = "dashed") +
    labs(x = "Número de Espécies", y = "Frequência", colour = "Estimativas") +
    scale_colour_brewer(palette = "Set1", labels = c("Teórico", "Simulado")) +
    theme(legend.position = "bottom")
    
ggplot(data = data.plot.bd) +
    geom_histogram(aes(x = ntip30)) +
    geom_vline(aes(xintercept = expect.N.30, colour = factor(expect.N.30)), size = 1, linetype = "dashed") +
    geom_vline(aes(xintercept = mean(ntip30), colour = factor(mean(ntip30))), size = 1, linetype = "dashed") +
    labs(x = "Número de Espécies", y = "Frequência", colour = "Estimativas") +
    scale_colour_brewer(palette = "Set1", labels = c("Teórico", "Simulado")) +
    theme(legend.position = "bottom")

O número esperado de espécies é muito similar entre os cenários sem e com extinção que simulamos acima, porém os dois cenários resultam em uma coleção de árvores completamente distinta. Qual é essa diferença e o porquê dela? (Não se esqueça de incluir as figuras que justifiquem sua resposta)


Fração de extinção (\(\epsilon\))

Até agora temos olhado para a relação entre especiação e extinção através do \(r\) (diversificação líquida). Porém outra forma de olhar para isso é através da taxa que chamamos de fração de extinção (extinction fraction em inglês). Essa taxa nos indica qual a importância relativa da extinção frente à especiação. Para entendermos melhor o que é a fração de extinção, vamos simular filogenias com o mesmo valor de \(r\), porém com diferentes valores de \(\epsilon\):

lambda.1 <- 0.2
lambda.2 <- 1
mu.1 <- 0.1
mu.2 <- 0.9

## r1
lambda.1 - mu.1

## epsilon 1
mu.1/lambda.1

## r2
lambda.2 - mu.2

## epsilon 2
mu.2/lambda.2


sim.r.eps.1 <- sim.bd.age(age = 15, numbsim = 100, lambda = lambda.1, mu = mu.1)

sim.r.eps.2 <- sim.bd.age(age = 15, numbsim = 100, lambda = lambda.2, mu = mu.2)

## Porcentagem de simulações bem sucedidas

sum(unlist(lapply(sim.r.eps.1, function(x){class(x) == "phylo"})))/100
sum(unlist(lapply(sim.r.eps.2, function(x){class(x) == "phylo"})))/100

## Sorteando 9 filogenias aleatórias de cada cenário para visualização

samp.plot.1 <- sample(which(unlist(lapply(sim.r.eps.1, function(x){class(x) == "phylo"})) == TRUE), 9)
samp.plot.2 <- sample(which(unlist(lapply(sim.r.eps.2, function(x){class(x) == "phylo"})) == TRUE), 9)

## Filogenias completas - Epsilon 1
par(mfrow = c(3, 3))
lapply(sim.r.eps.1[samp.plot.1], function(x){plot(x, show.tip.label = FALSE);axisPhylo()})

## Filogenias completas - Epsilon 2
par(mfrow = c(3, 3))
lapply(sim.r.eps.2[samp.plot.2], function(x){plot(x, show.tip.label = FALSE);axisPhylo()})

## Filogenias moleculares - Epsilon 1
par(mfrow = c(3, 3))
lapply(lapply(sim.r.eps.1[samp.plot.1], drop.fossil), function(x){plot(x, show.tip.label = FALSE);axisPhylo()})

## Filogenias moleculares - Epsilon 2
par(mfrow = c(3, 3))
lapply(lapply(sim.r.eps.2[samp.plot.2], drop.fossil), function(x){plot(x, show.tip.label = FALSE);axisPhylo()})

Por que no segundo cenário (Epsilon 2) apenas cerca de 10%-20% das simulações são bem sucedidas (chegam até o presente), enquanto no primeiro esse número gira em torno de 60%? (Não se esqueça de incluir as figuras que justifiquem sua resposta)


Pull of the present

Como discutido anteriormente, em um cenário de nascimento puro o número de espécies cresce exponencialmente com o tempo. Sendo assim, espera-se que o LTT plot de árvores resultantes desse cenário apresentem um crescimento linear quando plotados em escala semi-logarítmica.

par(mfrow = c(1, 1))

pb.trees <- sim.bd.taxa.age(n = 100, numbsim = 100, lambda = 0.2, mu = 0, age = 20)

ltt.plot(pb.trees[[1]], log = "y", xlab = "Tempo", ylab = "Número cumulativo de linhagens")
lapply(pb.trees[2:100], ltt.lines)

Agora repetiremos o processo para árvores simuladas em cenários com extinção. Observe as diferenças em relação aos mesmos plots para um cenário sem extinção (simulações anteriores).

bd.trees <- sim.bd.taxa.age(n = 200, numbsim = 100, lambda = 0.5, mu = 0.3, age = 20)

ltt.plot(bd.trees[[1]], log = "y", xlab = "Tempo", ylab = "Número cumulativo de linhagens")
lapply(bd.trees[2:100], ltt.lines)

Essa diferença é ainda mais acentuada quando a fração de extinção é maior:

bd.trees <- sim.bd.taxa.age(n = 200, numbsim = 100, lambda = 1, mu = 0.8, age = 20)

ltt.plot(bd.trees[[1]], log = "y", xlab = "Tempo", ylab = "Número cumulativo de linhagens")
lapply(bd.trees[2:100], ltt.lines)

Podemos notar que os LTT plots nos cenários com extinção não são iguais aos do cenário de Pure Birth, especialmente próximos do presente. Por que isso acontece? (Não se esqueça de incluir as figuras que justifiquem sua resposta)


Estatística Gamma (\(\gamma\))

Uma das formas que temos para resumir as características de uma filogenia é a estatística Gamma (Pybus & Harvey, 2000). Essa métrica representa, de maneira grosseira, o “centro de massa” da filogenia, ou seja, é um teste para analisar se os nós de uma filogenia estão mais concentrados próximos à raiz ou aos terminais da filogenia quando comparados a um modelo nulo (nascimento puro), onde a expectativa é que os nós estejam espaçados regularmente desde a raiz até o presente.

phylo.gamma.null <- sim.bd.taxa(n = 200, numbsim = 1000, lambda = 0.2, mu = 0)
gamma.null <- unlist(lapply(phylo.gamma.null, gammaStat))

## Simular uma árvore sob Pure Birth

phylo.test.pb <- sim.bd.taxa(n = 200, numbsim = 1, lambda = 0.2, mu = 0)

## Simular uma árvore com desaceleração da diversificação no tempo
set.seed(1234)
b.decl <- function(t){0.5 * exp(-0.2 * t) + 0.1}
phylo.test.decel <- rlineage(birth = b.decl, death = 0, Tmax = 25)

## Comparando a filogenia com desaceleração com filogenias Pure Birth
par(mfrow = c(2, 2))
plot(phylo.test.decel, show.tip.label = FALSE);axisPhylo()
lapply(phylo.gamma.null[sample(1:1000, 3)], function(x){plot.phylo(x, show.tip.label = FALSE);axisPhylo()})

sum(gamma.null <= gammaStat(phylo.test.pb[[1]]))/1000
sum(gamma.null <= gammaStat(phylo.test.decel))/1000

Para variar um pouco e sair do mundo das simulações, vamos agora analisar o valor de \(\gamma\) para uma árvore empírica. Para isso, baixem este arquivo que contém uma filogenia de Viperidae (Alencar et al., 2016, Mol. Phy. Evol.), que vamos assumir que contém 100% da diversidade vivente.

viper.tree <- read.tree("viper_tree.txt")

## Visualizando a árvore
par(mfrow = c(1, 1))
plot(viper.tree, type = "fan", cex = 0.5)

## Calculando o valor de gamma para a árvore empírica
gammaStat(viper.tree)

## Gerando uma distribuição nula de Gamma para testar a significância do valor de gamma da árvore empírica
phylo.gamma.null.viper <- sim.bd.taxa(n = Ntip(viper.tree), numbsim = 1000, lambda = 0.2, mu = 0)

sum(unlist(lapply(phylo.gamma.null.viper, gammaStat)) <= gammaStat(viper.tree))

Importante ressaltar que, apesar de nos dois casos acima as filogenias geradas para a construção do modelo nulo possuírem o mesmo número de espécies que a filogenia com a qual estamos comparando, isso não é necessário. Para mostrar isso, vamos simular uma distribuição nula de valores de \(\gamma\) com árvores um pouco maiores (demora um pouco também):

phylo.gamma.null.large <- sim.bd.taxa(n = 1000, numbsim = 1000, lambda = 0.2, mu = 0)

## Plotando lado a lado as distribuições nulas (geradas para a árvore empírica e as com um maior número de tips)
gamma.viper <- unlist(lapply(phylo.gamma.null.viper, gammaStat))
gamma.large <- unlist(lapply(phylo.gamma.null.large, gammaStat))

hist(gamma.viper, freq = FALSE, main = NULL,  xlab = "Gamma", ylab = "Densidade")
hist(gamma.large, freq = FALSE, xlab = "Gamma", ylab = "Densidade", add = TRUE, border = "red")
legend("topright", legend = c("Vipers", "Árvores Grandes"), col = 1:2, lty = 1)

Nos casos acima (tanto simulado quanto empírico), há uma premissa importante que assumimos: as filogenias contém 100% das espécies viventes. Na maioria dos casos isso não acontece, e isso pode ter um impacto bastante importante nas conclusões tiradas a partir do valor de \(\gamma\).

## Removendo aleatoriamente 10% das espécies
viper.tree.80 <- drop.random(viper.tree, round(Ntip(viper.tree) * 0.2))

## Removendo aleatoriamente 30% das espécies
viper.tree.60 <- drop.random(viper.tree, round(Ntip(viper.tree) * 0.4))

## Removendo aleatoriamente 50% das espécies
viper.tree.40 <- drop.random(viper.tree, round(Ntip(viper.tree) * 0.6))


## Calculando gamma para os 3 cenários
gamma.viper.80 <- gammaStat(viper.tree.80)
gamma.viper.60 <- gammaStat(viper.tree.60)
gamma.viper.40 <- gammaStat(viper.tree.40)
setNames(c(gammaStat(viper.tree), gamma.viper.80, gamma.viper.60, gamma.viper.40), c("Completa", "80%", "60%", "40%"))

Qual é o efeito do grau de sub-amostragem da filogenia nos valores de \(\gamma\)? Qual seria um possível teste que leve em consideração essa sub-amostragem? (Não se esqueça de incluir as figuras que justifiquem sua resposta)


Taxas variáveis

Até agora, usamos modelos onde as taxas são constantes ao longo do tempo em todas as árvores que simulamos. Em alguns casos essas taxas constantes podem até representar bem a dinâmica do clado (por exemplo quando analisamos a diversificação de um gênero por exemplo, que se dá em “curtas” escalas temporais). No entanto, é possível relaxarmos essa restrição e simular cenários mais realísticos para níveis hierárquicos superiores, onde a premissa de taxas constantes passa a ser pouco plausível. (Para facilitar a análise do resultado, copie e cole o código abaixo inteiro).

set.seed(2018)

## Expansão de Diversidade
b.cresc <- function(t){0.1}

plot(seq(0, 15, length.out = 100), rep(0.4, 100), type = "l", ylim = c(0.05, 0.8))
lines(seq(0, 15, length.out = 100), y = rep(0.1, 100), col = 2)

sim.bd.cresc <- rlineage(birth = 0.4, death = 0.1, Tmax = 15)

## Saturação de Diversidade
b.sat <- function(t){0.8 * exp(-0.2 * t) + 0.1}

plot(seq(0, 20, length.out = 100), b.sat(seq(0, 20, length.out = 100)), type = "l", ylim = c(0, 0.8))
lines(seq(0, 20, length.out = 100), y = rep(b.sat(20), 100), col = 2)

sim.bd.sat <- rlineage(birth = b.sat, death = b.sat(20), Tmax = 20)

## Declínio de Diversidade
b.decl <- function(t){0.8 * exp(-0.2 * t) + 0.1}

plot(seq(0, 25, length.out = 100), b.decl(seq(0, 25, length.out = 100)), type = "l", ylim = c(0.05, 0.8))
lines(seq(0, 25, length.out = 100), y = rep(b.decl(20) + 0.025, 100), col = 2)

sim.bd.decl <- rlineage(birth = b.decl, death = b.decl(20), Tmax = 25)

par(mfrow = c(2, 3))
plot(seq(0, 15, length.out = 100), rep(0.2, 100), type = "l", ylim = c(0.05, 0.8))
lines(seq(0, 15, length.out = 100), y = rep(b.cresc(20), 100), col = 2)
plot(seq(0, 20, length.out = 100), b.sat(seq(0, 20, length.out = 100)), type = "l", ylim = c(0, 0.8))
lines(seq(0, 20, length.out = 100), y = rep(b.sat(20), 100), col = 2)
plot(seq(0, 25, length.out = 100), b.decl(seq(0, 25, length.out = 100)), type = "l", ylim = c(0.05, 0.8))
lines(seq(0, 25, length.out = 100), y = rep(b.decl(20) + 0.025, 100), col = 2)
plot(drop.fossil(sim.bd.cresc), show.tip.label = FALSE);axisPhylo()
plot(drop.fossil(sim.bd.sat), show.tip.label = FALSE);axisPhylo()
plot(drop.fossil(sim.bd.decl), show.tip.label = FALSE);axisPhylo()

Vamos agora comparar os valores de \(\gamma\) para as três filogenias simuladas:

gammaStat(sim.bd.cresc)
gammaStat(sim.bd.sat)
gammaStat(sim.bd.decl)

## Teste de significância

gamma.null.cresc <- sim.bd.taxa(n = Ntip(drop.fossil(sim.bd.cresc)), numbsim = 100, lambda = 0.4, mu = 0)
gamma.null.sat <- sim.bd.taxa(n = Ntip(drop.fossil(sim.bd.sat)), numbsim = 100, lambda = 0.4, mu = 0)
gamma.null.decl <- sim.bd.taxa(n = Ntip(drop.fossil(sim.bd.decl)), numbsim = 100, lambda = 0.4, mu = 0)

sum(gammaStat(sim.bd.cresc) <= unlist(lapply(gamma.null.cresc, gammaStat)))/100
sum(gammaStat(sim.bd.sat) <= unlist(lapply(gamma.null.sat, gammaStat)))/100
sum(gammaStat(sim.bd.decl) <= unlist(lapply(gamma.null.cresc, gammaStat)))/100

O que podemos afirmar com relação à capacidade da estatística \(\gamma\) em distingüir entre os diferentes cenários? (Não se esqueça de incluir as figuras que justifiquem sua resposta)



Simulações do registro fóssil

Como vimos no tutorial da aula passada, os modelos usados para analisar dados do registro fóssil são bastante similares (muitas vezes idênticos) aos modelos usados com filogenias moleculares. Sendo assim, não vamos focar muito nas propriedades desses modelos (estocasticidade, valores esperados, etc.) que já vimos anteriormente, mas vamos avaliar algumas propriedades que são exclusivas dos dados fósseis.

Usaremos as seguintes funções para converter os dados das simulações em dados mais amigáveis para os plots.

fossil.for.plot <- function(x){
    first.last <- plyr::ldply(x, function(y){y$taxa.data})
    sampling.occs <- setNames(plyr::llply(x, function(y){y$sampling.times}), NULL)
    occs <- data.frame(samp.times = 0, taxon.id = 0) #matrix(NA, ncol = length(sampling.occs[which.max(unlist(lapply(sampling.occs, length)))][[1]]))
    for(i in 1:dim(first.last)[1]){
        #print(i)
        if(length(sampling.occs[[i]]) == 0){
            occs <- rbind(occs, data.frame(samp.times = c(first.last$orig.time[i], first.last$ext.time[i]), taxon.id = rep(i, 2)))
        } else {
            temp <- sampling.occs[[i]]
            occs <- rbind(occs, data.frame(samp.times = temp, taxon.id = rep(i, length(temp))))
        }
    }
    return(list(first.last, occs[-1, ]))
}

taxa.for.plot <- function(x){
    res.df <- data.frame(occs = NA, taxon.id = NA)
    for(i in 1:length(x)){
        res.df <- rbind(res.df, data.frame(occs = x[[i]], taxon.id = rep(i, length(x[[i]]))))
    }
    res.df <- res.df[-1,]
    return(res.df)
}
        

Efeitos da qualidade do registro

Registro “perfeito”

Inicialmente vamos simular um cenário ideal: todos os fósseis que aparecem na história do grupo são preservados, e sabemos os momentos de surgimento e desaparecimento verdadeiros das espécies. O comando a seguir simula o registro, bem como plot as curvas de diversidade real e amostrada.

record.full <- simFossilRecord(p = 0.3, q = 0.2, r = sProb2sRate(R = 0.99999), nTotalTaxa = c(30, 50), nExtant = c(0, 10))

record.full.plot <- fossil.for.plot(record.full)

## Plotando os ranges verdadeiros das espécies

ggplot(record.full.plot[[1]]) +
    geom_segment(aes(x = -orig.time, xend = -ext.time, y = factor(taxon.id), yend = factor(taxon.id), colour = factor(taxon.id)), size = 1) +
    labs(x = "Tempo", y = "Espécie") +
    theme(legend.position = "none")

## Plotando as ocorrências  das espécies

recfull.plot <-
    ggplot(record.full.plot[[2]]) +
    geom_point(aes(x = -samp.times, y = factor(taxon.id), colour = factor(taxon.id)), alpha = 0.5, size = 1) +
    labs(x = "Tempo", y = "Espécie") +
    theme(legend.position = "none") +
    xlim(-max(record.full.plot[[1]]$orig.time), 0)
recfull.plot

É possível notar que os dois gráficos acima mostram praticamente a mesma informação, já que quando a amostragem é quase perfeita, os momentos de aparecimento e desaparecimento das espécies são muito próximos (se não idênticos) aos “verdadeiros”. Isso também é ressaltado ao plotarmos as curvas de diversidade real e amostrada:

## Plotando as curvas de diversidade
divCurveFossilRecordSim(record.full)

Veremos o que acontece quando o registro não é perfeito no próximo passo.


Degradando o registro perfeito

Agora vamos usar os dados simulados acima com preservação perfeita e simular cenários onde a preservação é mediana e baixa.

record.full.taxa <- fossilRecord2fossilTaxa(record.full)

record.50 <- sampleRanges(record.full.taxa, r = sProb2sRate(0.5), ranges.only = FALSE)
record.10 <- sampleRanges(record.full.taxa, r = sProb2sRate(0.1), ranges.only = FALSE)

record.50.rang <- setNames(ldply(record.50, range)[, c(1, 3, 2)], c("taxon.id", "orig.time", "ext.time"))
record.10.rang <- setNames(ldply(record.10, range)[, c(1, 3, 2)], c("taxon.id", "orig.time", "ext.time"))

record.50.plot <- taxa.for.plot(record.50)
record.10.plot <- taxa.for.plot(record.10)

---

## Ocorrências para 50% de preservação

rec50.plot <-
    ggplot(data = record.50.plot) +
    geom_point(aes(x = -occs, y = factor(taxon.id), colour = factor(taxon.id)), alpha = 0.5, size = 1) +
    labs(x = "Tempo", y = "Espécie") +
    theme(legend.position = "none") +
    xlim(-max(record.full.plot[[1]]$orig.time), 0)

## Ocorrências para 10% de preservação

rec10.plot <-
    ggplot(data = record.10.plot) +
    geom_point(aes(x = -occs, y = factor(taxon.id), colour = factor(taxon.id)), alpha = 0.5, size = 1) +
    labs(x = "Tempo", y = "Espécie") +
    theme(legend.position = "none") +
    xlim(-max(record.full.plot[[1]]$orig.time), 0)

plot_grid(recfull.plot,
          rec50.plot,
          rec10.plot,
          ncol = 3,
          align = 'h'
          )

Ao plotarmos o range das espécies, essas diferenças ficam mais evidentes:


range.full.plot <-
    ggplot(record.full.plot[[1]]) +
    geom_vline(xintercept = 0, colour = "darkgrey", linetype = "dashed") +
    geom_segment(aes(x = -orig.time, xend = -ext.time, y = factor(taxon.id), yend = factor(taxon.id), colour = factor(taxon.id)), size = 1) +
    labs(x = "Tempo", y = "Espécie") +
    theme(legend.position = "none")

range.50.plot <-
    ggplot(data = record.50.rang) +
    geom_vline(xintercept = 0, colour = "darkgrey", linetype = "dashed") +
    geom_segment(aes(x = -orig.time, xend = -ext.time, y = factor(taxon.id), yend = factor(taxon.id), colour = factor(taxon.id)), size = 1) +
    geom_point(data = record.50.rang[(record.50.rang$orig.time - record.50.rang$ext.time) < 0.000001, ], aes(x = -orig.time, y = factor(taxon.id), colour = factor(taxon.id)), size = 1) +
    labs(x = "Tempo", y = "Espécie") +
    xlim(-max(record.full.plot[[1]]$orig.time), 0) +
    theme(legend.position = "none")

range.10.plot <-
    ggplot(data = record.10.rang) +
    geom_vline(xintercept = 0, colour = "darkgrey", linetype = "dashed") +
    geom_segment(aes(x = -orig.time, xend = -ext.time, y = factor(taxon.id), yend = factor(taxon.id), colour = factor(taxon.id)), size = 1) +
    geom_point(data = record.10.rang[(record.10.rang$orig.time - record.10.rang$ext.time) < 0.000001, ], aes(x = -orig.time, y = factor(taxon.id), colour = factor(taxon.id)), size = 1) +
    labs(x = "Tempo", y = "Espécie") +
    xlim(-max(record.full.plot[[1]]$orig.time), 0) +
    theme(legend.position = "none")

plot_grid(range.full.plot,
          range.50.plot,
          range.10.plot,
          ncol = 3,
          align = 'h'
          )

De posse dos registros com diferentes graus de preservação, podemos também traçar as curvas de diversidade para cada um dos casos:

## Preparando os dados para o gráfico
rang.data <- list(full = record.full.taxa, pres50 = record.50.rang[,2:3], pres10 = record.10.rang[,2:3])
rang.list <- list()
for(i in 1:3){
    rang.list[[i]] <- as.data.frame(taxicDivCont(rang.data[[i]], int.length = 1))
    rang.list[[i]]$foo <- c(NA, rang.list[[i]]$int.div[-length(rang.list[[i]]$int.div)])
    rang.list[[i]]$name <- rep(names(rang.data)[i], dim(rang.list[[i]])[1])
}

div.data <- dplyr::bind_rows(rang.list)

ggplot(data = div.data) +
    geom_segment(aes(x = -int.start, xend = -abs(int.end), y = factor(int.div), yend = factor(int.div), colour = name), size = 1, alpha = 0.75) +
    geom_segment(aes(x = -int.start, xend = -int.start, y = factor(int.div), yend = factor(foo), colour = name), size = 1, alpha = 0.75) +
    scale_colour_brewer(palette = "Set1", labels = c("Perfeito", "50%", "10%")) +
    labs(x = "Tempo", y = "Número de Espécies", colour = "Preservação") +
    theme(legend.position = "bottom")

Por fim, vamos avaliar como a resolução temporal do registro fóssil pode alterar a reconstrução da diversidade ao longo do tempo:


rang.list.res <- list()
for(i in 1:3){
    rang.list.res[[i]] <- as.data.frame(taxicDivCont(record.full.taxa, int.length = c(5, 1, 0.1)[i]))
    rang.list.res[[i]]$foo <- c(NA, rang.list.res[[i]]$int.div[-length(rang.list.res[[i]]$int.div)])
    rang.list.res[[i]]$name <- rep(paste0("Resolução = ", c(5, 1, 0.1)[i], " MYr"), dim(rang.list.res[[i]])[1])
}

div.data.res <- dplyr::bind_rows(rang.list.res)
div.data.res$name <- factor(div.data.res$name, levels = unique(div.data.res$name)[c(3, 2, 1)])

ggplot(data = subset(div.data.res, int.start < 37)) +
    geom_segment(aes(x = -int.start, xend = -abs(int.end), y = int.div, yend = int.div, colour = name), size = 1, alpha = 0.75) +
    geom_segment(aes(x = -int.start, xend = -int.start, y = int.div, yend = foo, colour = name), size = 1, alpha = 0.75) +
    scale_colour_brewer(palette = "Set1") +
    labs(x = "Tempo", y = "Número de Espécies", colour = "Preservação") +
    theme(legend.position = "bottom") +
    #xlim(-37, 0) +
    facet_grid(name  ~ .)

Avalie o efeito da preservação e da resolução temporal na caracterização da trajetória da diversidade no tempo. (Não se esqueça de incluir as figuras que justifiquem sua resposta)

see ya'