Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- historico:2014:ensaios:sugawara [2014/10/25 05:04] – mauro.sugawara
+++ historico:2014:ensaios:sugawara [2022/11/24 14:12] (atual) – edição externa 127.0.0.1
@@ Linha 1: / Linha 1: @@
+====== Modelos em macroevolução (ou porque biólogos devem conduzir estudos de biologia) ======
+=== Mauro Toshiro Caiuby Sugawara ===
+* Mestrando no departamento de Ecologia da USP - Lab. de Macroecologia e Macroevolução
+*maurotcs@gmail.com
+===== "Justificativa" =====
+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
+===== Preâmbulo =====
+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.
+===== Função de Verossimilhança e AICc =====
+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.
+===== Modelagem =====
+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 {{:historico:2014:ensaios:tabela_aicc_mod.csv|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.
+^ Modelo      											^ AICc      ^ dAICc  ^ k
+| Y ~ N(a+b*manchas, sd(Y))					| 448.48	| 0		| 3 |
+| Y ~ N(a*manchasˆb, sd(Y))					| 449.99	| 0.505	| 3 |
+| Y ~ N(a+b*manchas+c*aneis, sd(Y))	| 449.80	| 1.318	| 4 |
+| Y ~ N(a+b*manchas, z*manchasˆw)	| 450.66	| 1.960	| 4 |
+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.
+{{ :historico:2014:ensaios:plot_toshiba.jpg |Figura 1}}
+**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.
+===== Conclusão =====
+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).
+===== Referências bibliográficas =====
+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.
+===== Citação =====
+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.
+===== Código =====
+<code>
+  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")
+</code>