Quando dois é pior do que um, ou: Hurlbert avisou!

Um vento nasceu sobre o Oceano Atlântico. Este vento soprou por cima dos mares e das terras, indo ao norte, em direção à Bahia. Deslocou dunas em Itaúnas, moveu jangadas na RESEX de Canavieiras, derrubou galhos secos na ReBio de Una, e apenas como uma brisa leve chegou a Ilhéus, balançando as folhas das palmeiras. Este vento não era o começo; mas ele era um começo.

E enquanto o vento soprava, eu fazia simulações no meu computador. Em parte inspirado por um post em que Stephen Heard se pergunta se, em trabalhos ecológicos, dois é de fato melhor do que um, resolvi explorar uma situação em que dois pode ser pior do que um. (Código disponível aqui)

O cenário: um ecólogo da Quarta Era interessado em saber se, e como, a abundância de sereias de manguezal está relacionada ao comprimento dos rios em que vivem. Para tal, o ecólogo escolheu quinze rios na Área Protegida Marinha de Tear*, no sul das Terras Ocidentais. Mas, já tendo começado o estudo, percebeu que a amostragem é pouco abrangente, pois todos os rios nesta região eram relativamente **curtos, não ultrapassando os 140 km de comprimento.

Para suprir tal falha e deixar seu estudo mais representativo, o ecólogo tomou coragem, atravessou pântanos impenetráveis e amostrou também quinze rios na Reserva Extrativista de Mayene (assinando um termo em que jurava não revelar a localização dos peixes oleosos, caso por acaso os encontrasse na sua pesquisa), com rios entre 110 e 230 km de comprimento. E finalmente, feito tudo isso, e obtendo as devidas permissões dos Comitês de Ética com Estudo de Seres Humanos, Seres Mais Ou Menos Humanos, e Animais, iniciou a amostragem.

Pois bem. Que problemas haveria em tal pesquisa? O ecólogo tinha duas hipóteses: uma relação linear – quanto maior o comprimento dos rios, mais sereias irão viver neles. Mas poderia também haver um limiar – por exemplo, acima de certo tamanho de rio, o número de sereias se manteria constante, mas decairia rapidamente abaixo deste limiar. E finalmente, poderia não haver relação alguma, caracterizada por um modelo nulo; ou a relação poderia ser mais complexa ainda, podendo ser descrita por um modelo linear ou um modelo aditivo.

Mas há uma outra possibilidade, que não foi prevista pelo ecólogo. Tear é uma nação mais desenvolvida e mais populosa que Mayene, com grande desigualdade social e diversos outros problemas. Enquanto nobres vivem em palácios cercados de muralhas e com ruas pavimentadas entre eles, o povão tem que andar por ruas enlameadas naquele calor insuportável. O que dizer do saneamento básico, então? Os rios de Tear devem ser mais poluídos; e mesmo que não sejam, as sereias de manguezal devem preferir evitar contato excessivo com seres humanos, conhecidos por sua agressividade. Mesmo que os rios de Tear fossem maiores (e eles não o são, não neste estudo), a abundância de sereias de manguezal neles ainda assim seria menor. Como então diferenciar entre o efeito do comprimento dos rios e a variação regional?

Abaixo, um gráfico mostrando essas quatro possibilidades: modelo nulo (sem variação); diferença entre as duas regiões; relação linear; e relacão com limiar. Triângulos representam Tear, círculos mostram Mayene. As linhas representam o modelo (correto) ajustado.

sereias.png

O modelo linear não está muito linear porque, como trabalhamos com dados de abundância, usei um modelo linear generalizado (GLM) com distribuição de Poisson e função de ligação log. Se usasse a função de ligação identidade, talvez o ajuste melhoraria. Um dia escrevo direito sobre isso. Por enquanto, saibam que, ao ajustarmos um modelo linear, ele pode aparecer com uma forma curva por causa de transformações que são feitas nos dados para possibilitar ou facilitar o ajuste (ou algo assim).

Bom, e como podemos determinar qual o melhor modelo estatístico, ou seja, o modelo que melhor se ajusta (sem sobreajuste) a cada conjunto de dados? Podemos usar uma abordagem de seleção de modelos. Não vou entrar em detalhes, mas essa abordagem consiste em ajustarmos diferentes modelos e selecionarmos um deles. Tá, isso não ajudou muito. 🙂 Então, quando falo “modelo”, estou me referindo basicamente à linha que em cada um dos gráficos acima. Só que estas linhas foram colocadas já sabendo qual o modelo correto, ou seja, sabendo como os dados foram gerados. E se não soubermos? Bom, neste caso ajustamos diferentes modelos aos dados e comparamos entre eles. Vou exemplificar com o quarto gráfico, criado a partir do chamado modelo piecewise:

sereias2.png

Na figura, a linha cinza contínua é o modelo nulo (sem efeito do comprimento do rio), a linha cinza descontínua é o modelo região (diferença entre as duas regiões), a linha preta contínua é o modelo linear, a linha preta descontínua é o modelo piecewise, e a linha preta pontilhada é o modelo aditivo. Olhando para a figura, já percebemos que os modelos piecewise e aditivo são os que melhor se ajustam aos dados. Isso também pode ser olhado de forma mais objetiva com uma coisa chamada AIC – critério de informação de Akaike. Criado por Hirotugu Akaike, um japa muy phoda estatístico japonês que contribuiu muito para a análise estatística, esta medida pondera quão bem o modelo se ajusta aos dados e a complexidade dele, evitando assim modelos complexos demais – afinal, não queremos incluir parâmetros suficientes para ajustar um elefante!

Usando o AIC – bom, AICc, que é uma correção do AIC para tamanho amostral pequeno, temos o seguinte resultado, fornecido a nós lindamente pelo R:

Modelo dAICc df
aditivo 0.0 4.1
piecewise 0.8 4
linear 67.6 2
região 102.8 2
nulo 147.6 1

O dAICc – delta AICc – é a diferença entre cada modelo e o melhor modelo. O melhor modelo – o que melhor se ajusta aos dados – tem um dAICc = 0. No caso, é o modelo aditivo. O modelo piecewise é só um pouquinho pior, com dAICc de 0.8. Na prática, modelos com dAICc baixo são indistinguíveis um do outro. Existem diferentes recomendações para o que seria um dAICc baixo, variando de 1 a 7; por não interpretarem direito o que Burnham e Anderson escreveram algum motivo, ecólogos costumam usar um dAICc de 2 como limiar para isso. Um dAICc de 2 diz que um modelo é aproximadamente 2.5 vezes mais provável que o outro. Com um dAICc de 4, um modelo é pouco mais de 7 vezes mais provável que outro. Um dAICc de 8 quer dizer um modelo é quase 55 vezes mais provável que outro (de acordo com a página 78 de Burnham e Anderson 2002). Mas, porque estou com preguiça e porque já demorei algumas horas neste post pra seguir o que mais usamos na ecologia, vou usar um dAICc de 2 como critério, o que nos diz que os modelos aditivo e piecewise são igualmente bons.

A outra coluna, df, são os graus de liberdade, ou degrees of freedom. Quanto mais graus de liberdade, mais complexo o modelo. Então, porque normalmente gostamos de simplicidade, vamos escolher o modelo mais simples destes dois – no caso, o piecewise, que tem 0.1 graus de liberdade a menos.*

Até aqui tudo bem. Mas e se tiver diferença entre as regiões, mas nenhum efeito do comprimento do rio?

sereias3.png

Parece que o modelo região é o que melhor se ajusta aos dados, e é isso que o dAIC nos mostra:

Modelo dAICc df
região 0.0 2
aditivo 38.8 5.1
piecewise 48.4 4
linear 52.2 2
nulo 103.9 1

Com um dAICc de 38.8, o modelo região é mais de 22000 vezes mais provável que o modelo aditivo. Até aí, tudo em ordem.
Mas, e se o ecólogo, não estando ciente destas questões, não incluir a região entre as variáveis explanatórias?

sereias4.png

Modelo dAICc df
aditivo 0.0 5.1
piecewise 9.6 4
linear 13.4 2
nulo 65.1 1

Um dAICc de 9.6 corresponde a uma razão de evidência de quase 148. Parece então que agora o modelo aditivo tem um ajuste melhor a estes dados! Seguido do modelo piecewise. Estou um pouco surpreso pelo fato do ponto de quebra no modelo stepwise estar tão distante da diferença entre as duas regiões. Mas, a conclusão se mantém: não incluindo o modelo região, concluiríamos que existe um gradiente, ou um padrão não linear, na relação entre a abundância de sereias e o comprimento do rio, sendo que a explicação real é muito mais simples (e também menos interessante, infelizmente). Também não vi problemas graves ao validar o modelo aditivo graficamente.

Moral da história: ao deixarmos de incluir uma variável importante que pode estar afetando os dados, podemos tirar conclusões bonitinhas mas errôneas. E para não acharem que isso foi fruto do acaso, de uma única simulação, fiz 1000 simulações de cada um dos quatro cenários, verifiquei os modelos selecionados para cada uma delas, e assim calculei a proporção de vezes que o modelo correto foi selecionado. Fiz isso de duas formas: incluindo todos os modelos, e também excluindo o modelo região. Quando mais de dois modelos tinham dAICc < 2, escolhi o com menos graus de liberdade; se mais de um modelo tivessem o mesmo número de graus de liberdade, mantinha os dois. (Tudo isso automatizado num scipt em R, obviamente!)

Bom… Os resultados****:

Quando o modelo nulo é o verdadeiro:

Incluindo região na análise:

Modelo Proporção de seleção
aditivo 0.8%
linear 0.8%
linear ou região 2.5%
nulo 92.9%
região 1.9%
piecewise 1.1%

Sem incluir região na análise:

Modelo Proporção de seleção
aditivo 0.8%
linear 0.8%
nulo 94.9%
piecewise 1.3%

O modelo nulo foi selecionado em 93 a 95% das vezes; a inclusão da região aumentou o erro tipo I (achar um efeito pelo acaso) em apenas 2%, ou seja, nenhum problema aí.

Quando o modelo região é verdadeiro:

Incluindo região na análise:

Modelo Proporção de seleção
aditivo 0.2%
região 99.6%
piecewise 0.2%

Sem incluir região na análise:

Modelo Proporção de seleção
aditivo 59.4%
linear 26.2%
piecewise 14.4%

E aqui que está o problema! Incluindo a região, temos 99% de acerto – um resultado excelente. Mas esquecendo desta variável, podemos inferir que é uma relação linear, uma relação piecewise, ou uma relação não-linear mais complexa. Além das instabilidades na resposta, iremos corroborar a nossa hipótese, de que o comprimento dos rios afeta a abundância das sereias de menguezal, que é totalmente falsa neste caso.

Bom, mas e as desvantagens? Será que incluir o modelo região dificultaria a gente encontrar a resposta certa nos outros dois casos? Vamos ver…

Quando o modelo linear é verdadeiro:

Os resultados incluindo ou deixando de incluir região foram idênticos:

Modelo Proporção de seleção
aditivo 59.7%
linear 9.2%
piecewise 31.1%

Assustador, né? O modelo correto, que é o linear, só foi selecionado 9% das vezes! Acho que isso indica que precisamos tomar mais cuidado ao usar modelos com limiares ou modelos aditivos (por mais que me doa falar assim – adoro modelos aditivos!), pois eles podem ser selecionados mesmo quando o padrão subjacente é linear. Ainda preciso explorar em detalhes isso. Talvez usando a função de ligação “identidade” o problema se resolveria – mas, como decidir isso com dados reais?

E quando o modelo piecewise é verdadeiro:
Independentemente de incluir ou não o modelo região, o resultado é:

Modelo Proporção de seleção
aditivo 23.1%
piecewise 76.9%

Com 77% de acerto (e mais 23% de acerto parcial, pois um modelo aditivo mostra também uma relação não-linear), estamos bem. Ou seja, quando a relação é do tipo piecewise, parece que ela é relativamente fácil de ser detectada. Mas quando a relação subjacente é linear, uma relação piecewise também é relativamente fácil de ser detectada.

Bom, tais problemas talvez possam ser resolvidos graficamente. Pegando os dados da primeira simulação, pro modelo linear, temos:

Modelo dAICc df
piecewise 0.0 4
aditivo 4.3 5.1
linear 26.6 2
região 97.0 2
nulo 252.3 1

Colocando o piecewise e o aditivo no gráfico:

sereias5.png

Ou seja, interpretando a figura, podemos corrigir a nossa conclusão e falar que a relação é linear apesar do modelo linear não ter sido selecionado.

Resumindo, morais da história:

  • Olhar os gráficos, não confiar cegamente nos resultados do AIC;
  • Tome cuidado com variáveis que podem estar afetando seus resultados sem que você perceba! Não queremos detectar relações espúrias, queremos?
  • Caso você ainda não leu o artigo clássico do Hurlbert, leia – eu apenas fiz simulações para mostrar uma pequena parte do que ele discute;
  • Se for estudar sereias de manguezal, não esqueça dos comitês de ética!

* Sim, estou ciente que Robert Jordan não falou nada de sereias em Tear, nem de APMs. Mas ele também não falou que elas não existem! E se falou, pode ter se enganado também. Quanto às APMs – do jeito que Rand era progressivo, bem capaz que tenha deixado instruções para isso também. 🙂 Quanto a manguezais – de fato, talvez Tear e Mayene sejam lugares frios demais para tal vegetação, mas não me parece impossível que manguezais tenham expandido sua zona no tempo entre a Terceira e a Quarta Eras.

** Não trabalho com rios, não sei dizer o que é curto ou longo. 🙂

*** Graus de liberdade costumam ser números inteiros, exceto para modelos aditivos, nos quais os chamados graus de liberdade efetivos podem ser valores decimais também.

**** Houve problemas de convergência em 20 modelos, de todos os ajustados (excetuando os _piecewise_); e às vezes os modelos _piecewise_ não convergiram – nestes casos não os incluí nos cálculos do AIC.

Anúncios

Um pensamento sobre “Quando dois é pior do que um, ou: Hurlbert avisou!

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s