* Mestrando no departamento de Ecologia da USP - Lab. de Macroecologia e Macroevolução

*maurotcs@gmail.com

Olá!

Como vocês podem ver meu ensaio está bem atrasado. Na verdade, eu estava terminando em cima da hora e quando acabei (~23h20 do dia 24/out) não consegui atualizar a minha página. Fiquei tentando recarregar mas não funcionava, acredito que muitas pessoas estavam fazendo o upload nessa hora. Só agora conseguir atualizar o texto (já tinha colocado uma versão anterior no wiki: Última modificação: 2014/10/24 21:03 por mauro.sugawara). Eu sei que este atraso também é culpa minha, afinal vocês deram um prazo longo para podermos fazer nossa resenha, mas gostaria que pelo menos minha resenha seja avaliada.

Obrigado

A macroevolução é o ramo da biologia evolutiva que estuda os padrões de biodiversidade que ocorrem acima do nível da espécie. Enquanto a microevolução (i.e., ramo da biologia evolutiva que estuda os padrões de biodiversidade dentro da espécie) trata de indivíduos de uma mesma espécie e explica os padrões encontrados em termos de sobrevivência e reprodução diferencial, a macroevolução estuda padrões entre espécies e explica os padrões encontrados em função da variação nas taxas de especiação e extinção. Algumas questões abordadas são: qual a relação entre morfologia e riqueza? Há um limite global para o número de espécies? Porque há tantas espécies de besouros no mundo? Quais condições permitiram o surgimento de novos modos de vida e qual o impacto destes sobre a trajetória da biodiversidade?

Dado a magnitude dos padrões estudados, os trabalhos na área costumam tratar de períodos de milhares de milhões de anos e, em virtude disto, macroevolução é o ramo com abordagem mais histórica dentro de biologia evolutiva. Para poder responder estas questões, os biólogos macroevolutivos dependem fortemente de ajustes e comparações de modelos (outra abordagem com modelos são os estudos com simulações numéricas ou de indivíduo).

Durante muito tempo a macroevolução ficou restrita aos estudos paleontológicos, porém, atualmente as filogenias moleculares datadas (através do relógio molecular) fornecem outra forma de se acessar os fenômenos que moldaram a biodiversidade. Agora os biólogos macroevolutivos podem utilizar dois tipos de dados: fósseis e dados moleculares. Os primeiros estudos baseados apenas em dados moleculares trouxeram resultados completamente discordantes dos já estabelecidos estudos paleontológicos. Uma dessas inconsistências, ainda muito debatida, é o fato da maioria dos estudos com filogenias datadas estimarem uma taxa de extinção igual ou muito próxima a zero (Quental & Marshall, 2009). Para tentar acomodar ambos tipos de dados, os estudos de macroevolução passaram a incorporar modelos cada vez mais complexos (e.g., Silvestro et al., 2014) e há a tendência (pelo menos a meu ver) de se priorizar pesquisadores que saibam mais sobre a parte matemática (método para estudar o padrão) em detrimento daqueles que dominam a teoria biológica (processo que possivelmente gerou tal padrão). Em virtude desta “matematização”, a macroevolução tem atraído cada vez mais pesquisadores da área de exatas (e.g., físicos e matemáticos). Embora, esta interação com pesquisadores de outras não seja em sí algo ruim (muito pelo contrário!), é preciso tomar cuidado para que as análises reflitam fenômenos biológicos.

Neste ensaio pretendo mostrar a importância de desenvolver as hipóteses antes de se utilizar uma abordagem de ajuste e seleção de modelos.

Antes de mostrar os modelos vou apresentar a Função de Verossimilhança e o AICc, abordagem que irei utilizar nas análises.

Considere uma variável X que pode ser descrita por um modelo com função de densidade probabilística $f_x(x | \theta)$, onde $\theta$ representa os parâmetros que controlam o comportamento de X. Para que possamos estimar a densidade probabilística de uma dada observação X = x seria necessário conhecer quais são os parâmetros da função ($\theta$). Contudo, estes parâmetros não são conhecidos. Na verdade, muitas vezes estes parâmetros são o que se pretende estimar com a pesquisa científica! Queremos saber qual a influência de um dado caracter e a riqueza de espécies ou estimar as taxas de especiação, caso já soubéssemos esses valores não seria necessário fazer pesquisa.

Para poder estimar os parâmetros do modelo que melhor se ajustam aos dados coletados utilizamos a função de verossimilhança ($\mathcal{L}$):

$$\mathcal{L} (\theta | X=x) = f_x (\theta | X = x)$$

Esta função estima o quão verossímil é um modelo e seus respectivos parâmetros, dado o conjunto de observações. Dessa maneira, o que varia na função de verossimilhança são os parâmetros do modelo, enquanto os dados são mantidos constantes. É importante ressaltar que, diferentemente de valores de probabilidade, o valor de verossimilhança por si só não é interpretável. Os valores de verossimilhança só podem ser utilizados para fazer inferência quando comparados com outros valores de verossimilhança.

Para tornar mais palpável, podemos pensar em um baralho com cartas pretas e vermelhas. Podemos considerar duas hipóteses a priori: Ha = cartas pretas e vermelhas ocorrem na mesma proporção ($\theta_a=0.5$); e Hb = há três vezes mais cartas pretas do que vermelhas ($\theta_b=0.75$). Vamos supor que eu sorteei 10 cartas e saíram 6 cartas pretas. Neste caso temos:

$$Ha: \mathcal{L} (0.5 | X=6) = 0.828$$

e

$$Hb: \mathcal{L} (0.75 | X=6) = 0.224$$

Para poder decidir qual dos parâmetros melhor se ajusta aos dados utilizamos a razão de verossimilhança, na qual se divide o maior valor de verossimilhança pelo menor. Neste caso: 0.828/0.224 = 3.70. Portanto podemos dizer que Ha ($\theta=0.5$) é 3.70 vezes mais plausível do que Hb ($\theta=0.75$). Dessa maneira, ao testar diversos valores é possível estimar os valores que maximizam a verossimilhança (MLE, do inglês “Maximum Likelihood Estimate”). Para o exemplo descrito acima o MLE seria $\theta=0.6$.

Na prática nunca lidamos com um único valor e a multiplicação ($\mathcal{L} (\theta | X=x_1) \times \mathcal{L} (\theta | X=x_2)\ldots$) de diversos números pequenos rapidamente tende a zero. Para contornar esse problema utilizamos o logaritmo (na base e) da verossimilhança (logLik, do inglês “Log Likelihood”). Como o logaritmo de valores menores que zero é negativo, é comum utilizar o negativo do logLik.

Uma vez que já se tenha estimado os parâmetros do modelo (os MLE), pode-se comparar o ajuste de diferentes modelos. O método mais comumente utilizado para este fim é o AIC (do inglês “Akaike Information Criterion”):

$AIC = -2 \times logLik + 2 \times k$,

onde k representa o número de parâmetros do modelo. Quanto menor o AIC melhor o modelo. Uma vez calculado o AIC de todos os modelos concorrentes, estima-se a diferença entre os AIC e o AIC do melhor modelo:

$$\delta_i = AIC_i - min(AIC)$$

Por convenção, todos os modelos com $\delta$AIC menor ou igual a 2 são considerados igualmente plausíveis.

O AIC é sensível a amostras muito pequenas (i.e., número de observações para cada parâmetro do modelo menor do que 40). Portanto, em geral em biologia utilizamos o AIC corrigido para amostras pequenas (AICc):

$AICc = -2 \times logLik + 2 \times k \times \fraq{n}{n-k-1}$,

onde n representa o número de observações.

Neste ensaio utilizei três séries temporais disponíveis no pacote “datasets” (versão 2.15.3) do R, a saber: “discoveries”, número de “grandes” invenções e descobertas científicas por ano; “sunspot.year”, número de manchas solares por ano; e “treering” comprimento de anéis de crescimento. Os dados analisados compreendem 100 anos, desde 1860 a 1959. Escolhi estas séries temporais por não apresentarem nenhuma relação causal entre si (alguém poderia argumentar que o número de manchas pode influenciar os anéis de crescimento, mas me parece muito improvável dado o efeito que outras variáveis como clima possuem sobre a formação dos anéis).

Nesta análise pretendo ver a influência dos anéis de crescimento e manchas solares sobre o número de descobertas. Para tanto, vou utilizar modelos de regressão. Estes modelos são muito comuns na estatística frequentista, contudo, utilizando a abordagem de verossimilhança você pode construir seus próprios modelos, indo além dos modelos de regressão clássicos. Rapidamente construí 10 modelos diferentes para tentar explicar o número de descobertas por ano (ver arquivo Tabela AICc). Os melhores modelos podem ser vistos na Tabela 1.

Tabela 1. Os quatro modelos com melhor ajuste (i.e., diferença de AICc menor ou igual a 2). dAICc, diferença de AICc; k, número de parâmetros.

Apenas um dos quatro melhores modelos (i.e., igualmente plausíveis) inclui os anéis de crescimento como variável preditora, enquanto todos os modelos (incluindo o modelo que leva em consideração os anéis) incluem o número de manchas solares. Dessa forma podemos concluir que o efeito do número de manchas solares sobre o número de descobertas é muito mais forte do que o efeito dos anéis de crescimento, correto? O único problema é que estas variáveis não tem nenhum motivo para estarem correlacionadas. Na verdade, não consigo pensar nem sequer em uma variável que esteja influenciando ambas concomitantemente.

Figura 1. Relação entre número de descobertas e número de manchas solares entre os anos de 1860 e 1959. A linha vermelha representa o modelo de regressão linear com melhor ajuste.

Além disso, se olharmos para o gráfico do melhor modelo, veremos que mesmo sendo o melhor ele não explica grande parte da variação (Figura 1). Isto é um indicativo de que há um modelo, possivelmente com outras variáveis preditoras (e.g., verba alocada para pesquisas nos três anos anteriores), que seja muito melhor para explicar o número de descobertas por ano. Em outras palavras, a abordagem de ajuste e seleção de modelos é muito poderosa pois permite ao pesquisador testar diversas hipóteses ao mesmo tempo (ao contrário da comparação entre hipótese nula e alternativa, da abordagem frequentista), contudo é possível ajustar qualquer modelo aos seus dados e é extremamente fácil fugir da realidade, construindo modelos sem nenhum realismo. Portanto, os modelos devem sempre estar embasados em teorias biológicas, a conexão entre o modelo e a hipótese deve estar muito clara. De modo que, ao aceitar um modelo como o mais plausível para explicar um padrão, seja possível identificar qual o processo associado.

Esta abordagem de ajuste e seleção de modelos (explicação detalhada abaixo) possui diversas vantagens, especialmente nos estudos macroevolutivos em que se tem diversas hipóteses concorrentes (Johnson & Omland, 2004), e a interação entre biólogos e pesquisadores de outras áreas com certeza é muito vantajosa. Entretanto, a capacidade de construir modelos rebuscados não deve ficar em primeiro plano. O conhecimento sobre os processos biológicos é imprescindível para a análise adequada dos dados, bem como para a interpretação dos resultados.

"E lembre-se, com grandes poderes, vem grandes responsabilidades..."

Sheldon Cooper. The Bing Bang Theory - The Panty Piñata Polarization (Temp.2, ep.7).

Johnson, J. B. & K. S. Omland, 2004. Model selection in ecology and evolution. Trends in Ecology & Evolution 19:101–108.

Quental, T. B. & C. R. Marshall, 2009. Extinction during evolutionary radiations: reconciling the fossil record with molecular phylogenies. Evolution 63:3158–3167.

Silvestro, D., J. Schnitzler, L. H. Liow, A. Antonelli & N. Salamin, 2014. Bayesian estimation of speciation and extinction from incomplete fossil occurrence data. Systematic biology 63(3): 349-367.

Este ensaio é um produto de disciplina da pós-graduação da Universidade de São Paulo. Para citá-lo:

Sugawara, M.T.C. 2014. Modelos em macroevolução (ou porque biólogos devem conduzir estudos de biologia). In: Prado , P.I & Batista, J.L.F. Modelagem Estatística para Ecologia e Recursos Naturais. Universidade de São Paulo. url: http://cmq.esalq.usp.br/BIE5781/doku.php?id=historico:2014:ensaios:sugawara.

  data(discoveries)
  data(sunspot.year)
  data(treering)
  dados <- data.frame(matrix(data=NA, nrow=100, ncol=4))
  colnames(dados) <- c("Ano", "Descobertas","ManchasSolares","AnelCresc")
  dados$Ano <- 1860:1959
  dados$Descobertas <- discoveries
  dados$ManchasSolares <- sunspot.year[160:259]
  dados$AnelCresc <- treering[7861:7960]
  dados
  # MODELOS
  require(bbmle)
  nLLlinear1 <- function(a=5, b=-0.05, sigma=sd(dados$Descobertas)){
    mu <- a + b*dados$ManchasSolares
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  (linear1 <- mle2(nLLlinear1))
  nLLlinear2 <- function(a=5, b=-0.05, z=1, w=-0.05){
    mu <- a + b*dados$ManchasSolares
    sigma <- z * dados$ManchasSolares^w
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  (linear2 <- mle2(nLLlinear2))
  nLLlinear3 <- function(a=1, b=-0.05, sigma=sd(dados$Descobertas)){
    mu <- a * dados$ManchasSolares^b
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  (linear3 <- mle2(nLLlinear3))
  nLLlinear4 <- function(a=5, b=-0.05, z=1, w=-0.05){
    mu <- a * dados$ManchasSolares^b
    sigma <- z * dados$ManchasSolares^w
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  (linear4 <- mle2(nLLlinear4))
  nLLlinear5 <- function(a=5, b=-0.05, sigma=sd(dados$Descobertas)){
    mu <- a + b*dados$AnelCresc
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  (linear5 <- mle2(nLLlinear5))
  nLLlinear6 <- function(a=5, b=-0.05, z=1, w=-0.05){
    mu <- a + b*dados$AnelCresc
    sigma <- z * dados$ManchasSolares^w
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  (linear6 <- mle2(nLLlinear6))
  nLLlinear7 <- function(a=5, b=-0.05, sigma=sd(dados$Descobertas)){
    mu <- a * dados$AnelCresc^b
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  (linear7 <- mle2(nLLlinear7))
  nLLlinear8 <- function(a=5, b=-0.05, c=-0.05, sigma=sd(dados$Descobertas)){
    mu <- a + b*dados$ManchasSolares + c*dados$AnelCresc
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  (linear8 <- mle2(nLLlinear8))
  nLLlinear9 <- function(a=5, b=-0.05, c=-0.05, z=1, w=-0.05){
    mu <- a + b*dados$ManchasSolares + c*dados$AnelCresc
    sigma <- z * dados$ManchasSolares^w
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  linear9 <- mle2(nLLlinear9)
  linear9
  nLLlinear10 <- function(a=5, b=-0.05, c=-0.05, z=1, w=-0.05){
    mu <- a + b*dados$ManchasSolares + c*dados$AnelCresc
    sigma <- z * dados$AnelCresc^w
    -sum(dnorm(x=dados$Descobertas, mean=mu, sd=sigma, log=TRUE))
  }
  (linear10 <- mle2(nLLlinear10))
  # tabela com o aic
  aic <- AICctab(linear1, linear2, linear3, linear4, linear5, linear6, linear7, linear8, linear9, linear10, base=TRUE, weights=TRUE, nobs=nrow(dados))
  plot(x=dados$ManchasSolares, y=dados$Descobertas, xlab="Numero de Manchas Solares", ylab="Numero de Descobertas")
  abline(a=coef(linear1)[1], b=coef(linear1)[2], lwd=2, col="red")