Interpretando os resultados de um GLM

Recentemente uma amiga me pediu ajuda para interpretar o resultado de um GLM (eu já escrevi um post explicando o que são GLMs? Não? Poxa. Alguém precisa me lembrar de fazer isso. Haha), e ao dar uma explicação detalhada por email eu pensei, Ei, isso dá um post legal! Tá aqui esse post :-) Agradeço à Bianca Righi por ceder a tabela e os gráficos.

Enquanto eu não escrevo um post sobre o que é um GLM, vai aqui uma explicação rápida: GLM vem de Generalized Linear Model, ou Modelo Linear Generalizado. Um modelo linear é um modelo que modela alguma coisa (nossa variávei Y) como uma função linear de uma outra coisa (nossa variável X). Então um modelo linear pode dizer que quando uma coisa aumenta, a outra também aumenta; ou quando uma coisa aumenta, a outra diminui. E este aumento ou esta diminuição é linear – podemos traçar uma linha reta no gráfico para mostrar isso.

Em um modelo linear podemos também modelar como uma coisa varia em função de duas ou mais outras coisas – como nosso Y varia em função de X1, X2, X3 etc. Isso é feito de forma aditiva: valor de Y = efeito de X1 + efeito de X2 + efeito de X3. Se usarmos uma letra, por exemplo a letra greta β, para representar o efeito de cada coisa, temos que Y = β1X12X23X3. E também podemos incluir interações, que é quando uma coisa afeta o efeito de outra coisa. Por exemplo, pode ser que, sei lá, o metabolismo de dragões vermelhos aumenta com a temperatura em condições de baixa umidade, mas em condições de alta umidade este efeito desapareça (porque a alta umidade atrapalha a geração do fogo interno que estes dragões precisam pra aumentar seu metabolismo) (acho que viajei demais agora…)

Além de assumir relações lineares, modelos lineares assumem que os resíduos do modelo – ou seja, a diferença entre os valores observados e os valores preditos pelo modelo – seguem uma distribuição normal. Aquela distribuição em forma de sino (não confundir com a distribuição paranormal em forma de fantasma). Só que muitas vezes temos dados que não seguem e nem deveriam seguir a distribuição normal – espera-se que dados de contagem, dados de proporção, até mesmo dados de biomassa sigam outras distribuições. Como diretriz geral, dados de contagem (número de indivíduos, número de espécies, número de registros) costumam se ajustar à distribuição de Poisson ou binomial negativa; dados de presença/ausência ou proporção (como frequência – proporção de unidades amostrais em que uma espécie ocorre) seguem distribuição binomial; e biomassa e medidas parecidas podem seguir uma distribuição gamma. E aí que entram os modelos lineares generalizados! O “generalizado” significa que podemos especificar outras distribuições e não nos limitar à distribuição normal.

Dada essa explicação rápida, vamos ao estudo de caso que originou esse post! Foi um trabalho de mestrado sobre ocorrência de baleias-jubarte na costa da Bahia, em dois locais: Abrolhos e Serra Grande. Os dados mostram quantas baleias foram avistadas, por dia, de pontos de observação nestes dois lugares, em diferentes anos (2014, 2015, 2018 e 2019).

SerraGrande

Seria um pecado mortal mencionar Serra Grande e não colocar uma foto.

A análise usada foi um GLM com distribuição de Poisson (porque a variável-resposta é uma variável de contagem, abundância de baleias; neste caso incluindo um parâmetro para modelar a super-dispersão), com três variáveis explanatórias: ano, local, e dia do ano (dia juliano é o número de dias desde o começo do ano). Também foi incluída a interação entre ano e local, significando que a diferença entre os anos pode variar conforme o local e ou que a diferença entre os locais pode variar conforme o ano.

Éis o resultado da análise; preferi apresentar a tabela antes de mostrar figuras porque pode ser um exercício interessante de interpretação:

interpretarGLM_1

Antes de prosseguir com a leitura do texto, sugiro que você gaste um tempo interpretando estes resultados. Mais interessante ainda seria você escrever essa interpretação. Gaste o tempo que quiser, ninguém está com pressa aqui! A não ser que você esteja com pressa. :-) E se quiser colocar sua interpretação nos comentários, eu acharia isso bem interessante!

Bom, agora que você gastou um tempinho interpretando esses resultados, vou dar a minha interpretação/explicação aqui.

Primeiro: o que são as colunas? No caso, a primeira coluna é a variável explanatória: o intercepto (que é como o intercepto de uma regressão linear – o valor de Y quando X for igual a zero), ano, dia, local, e interação entre local e ano. A segunda coluna (Estimate) é a estimativa do efeito: valores positivos significam que o efeito é positivo, valores negativos significam que o efeito é negativo. Quanto mais positivo ou mais negativo o valor, mais intenso é o efeito. A terceira coluna (p-value) é o p-valor; quanto mais baixo for o p-valor, mais evidência temos contra a hipótese nula (a hipótese de que não há um efeito). 2e-16 significa 2*10-16, ou seja, 0.0000000000000006. Onde está <2e-16 significa portanto que é um valor muito, muito, muito baixo, mesmo. E a última coluna são estrelinhas resumindo os níveis de significância(menor que 0.001, entre 0.001 e 0.01, entre 0.01 e 0.05, e entre 0.05 e 0.1. Se for maior que 0.1 não ganha estrelinha).

Segundo: as variáveis explanatórias local e ano são variáveis categóricas. Teoricamente o ano poderia ser tratado como uma variável contínua, mas como são apenas quatro anos isso não seria muito interessante. Temos então quatro anos: 2014, 2015, 2018 e 2019; e temos dois locais: Abrolhos e Serra Grande.

Então por que abóboras que 2014 e Abrolhos não aparecem na tabela? Por que só aparece Serra Grande e os anos 2015, 2018 e 2019? (Essa é uma das perguntas que mais recebo sobre interpretação de GLM). O motivo é que para variáveis categóricas não faz muito sentido estimar o intercepto e um valor para cada nível – isso implicaria em estimar um valor a mais do que o necessário e do que é possível. Por exemplo, digamos que você compara a abundância de dragões vermelhos em cavernas naturais e em antigos reinos anânicos abandonados. Digamos que você encontra em média dois dragões por caverna natural e cinco por reino abandonado (talvez em Krynn algo assim possa ocorrer). Você pode representar isso como Ndragões = (2 se caverna natural ou 5 se reino abandonado). Ou você pode representar isso como Ndragões = 2 + (3 se reino abandonado). Neste último caso, nós temos o intercepto (2) e temos o efeito do reino abandonado (3); se o ambiente for caverna natural, não teremos o efeito do reino abandonado, e então Ndragões = 2.

Outra forma de pensar isso é usando um dummy variable – uma variável, que vou chamar de X, que assume o valor de 0 se o local for caverna natural e 1 se for reino abandonado. Então, Ndragões = 2 + 3*X. Se o ambiente for caverna natural, X=0 e Ndragões=2 + 3*0 = 2. Se o ambiente for reino abandonado, X=1 e Ndragões=2 + 3*X = 5.

Por que não incluir o intercepto e os dois coeficientes? Bom, isso deixa o modelo mais complicado: teríamos por exemplo as dummy variables X1 = 1 se caverna natural e 0 se reino abandonado e X2 = 0 se caverna natural e 1 se reino abandonado. É redundante. Além disso, incluindo o intercepto, podemos ter infinitas soluções: Ndragões = 1 + 1*X1 + 4*X2; ou Ndragões = -1 + 3*X1 + 6*X2; Ndragões = 3 + (-1)*X1 + 2*X2; enfim. E como é mais simples trabalhar com um intercepto e um coeficiente relacionado à variável (embora talvez menos intuitivo), o coeficiente das cavernas naturais é omitido e apenas o do reino abandonado é utilizado: Então, Ndragões = 2 + 3*X; X=1 se reino abandonado e X=0 caso contrário.

E é por isso que não temos coeficiente pra 2014 e pra Abrolhos na nossa tabela de resultados para as baleias! Estes valores já estão implícitos no intercepto. Se nada for especificado, o primeiro nível de cada variável categórica (em ordem alfabética) é usado como nível-base, e os efeitos (Estimate) mostrados mostram o quanto cada outro nível difere deste nível base. Eu acho sensacional que Bianca já colocou isso embaixo da tabela: “Diferença de 2014” e “Diferença do Arquipélago dos Abrolhos”.

Sendo assim. o intercepto se refere ao valor no dia 0 (dia é variável contínua), para o ano de 2014 (nível-base da variável Ano) e para o Arquipélago de Abrolho (nível-base da variável Local). Aí na segunda linha vemos que o Estimate para ano 2015 é -0.24. É um valor negativo – ou seja, em 2015 teve menos baleias do que em 2014. Em 2019 teve menos ainda. O Estimate negativo para Serra Grande indica que em Serra Grande tem menos baleias que em Abrolhos.

Quer dizer, seria só isso não tivesse interações no modelo. Como temos interações, essa diferença entre os anos diz respeito a Abrolhos, e essa diferença entre locais diz respeito a 2014. Em outros anos a diferença entre locais pode ter sido diferente e em Serra Grande a diferença entre os anos pode ter sido diferente.

Mas antes de falar das interações, tem um outro detalhe que ainda não explorei: o I(dia juliano^2). Que abóboras é isso?, você talvez tenha se perguntado (talvez não com essas palavras). Bom, pra quem conhece as baleias na costa da Bahia, sabe que elas não ficam o tempo todo aqui – chegam por volta de julho, ficam alguns meses e voltam pra se alimentar nas Ilhas Geórgia e Sanduiche do Sul (que ficam tipo a meio caminho entre Brasil e Antârtida). Então não faria sentido modelar sua abundância como uma relação linear do tempo, faria? O que precisamos é um modelo que aumenta e depois diminui – por exemplo uma parábola! E como modelamos uma parábola? Usando uma equação quadrática!

E é para isso que servem as equações de segundo grau que aprendemos na escola :-) (Entre outras coisas) Onde está escrito dia juliano^2 quer dizer dia juliano elevado ao quadrado (o ^ significa “elevado a”). E o I que vem antes serve basicamente pro R pegar o dia juliano e elevar ao quadrado, ao invés de fazer alguma outra coisa. E assim teremos uma equação do segundo grau, que modela algo que aumenta e depois diminui (ou algo que diminui e depois aumenta).

Mas Pav…
– Sim?
Você não disse que são modelos lineares?
– Disse!
Mas este modelo aí não é um modelo parabólico?
– É sim! E também é um modelo linear :-)
…Say whaaaat?

Explicando… Lembram que falei que em um modelo linear o efeito de diferentes coisas é modelado de forma aditiva? Y = β1X12X2, etc. Pois então, neste caso, X1 é o dia juliano, e X2 é o dia Juliano elevado ao quadrado. De modo que Y = β1Dia+β2Dia2. O resultado é uma parábola; mas funciona como um modelo linear.

Como então interpretar esse resultado? Bom, primeiro vemos que o efeito de I(dia juliano^2) foi significativo (p<2*10-16) – de modo que a relação é de fato quadrática, e não linear. Segundo que o coeficiente é negativo. Em uma equação quadrática, quando o elemento quadrático tem sinal negativo, isso significa que a curva sobe e depois desce (visualize mentalmente os valores de -x2 perto de zero para entender melhor). Se o sinal for positivo, a curva desce e depois sobe. (E confesso que precisei perguntar isso para meu amigo Marcos, que mantém este blog bemmm dahora sobre coisas matemáticas legais). Como o sinal é negativo, concluímos que a abundância de baleias é menor no começo do ano, depois sobre, atinge um pico, e desce. O que está de acordo com o que já sabemos sobre esses serzinhos lindos.

Então para avaliar como a abundância de baleias varia ao longo do tempo, podemos fazer uma curva mostrando essa relação quadrática. Vamos usar o intercepto e os coeficientes do dia juliano e do dia juliano^2. Assim estaremos modelando a variação na abundância de baleias ao longo do ano de 2014 em Abrolhos; se quiséssemos modelar o ano 2015, por exemplo, teríamos que incluir o valor correspondente.

A equação (usando dois dígitos significativos) fica assim: Y = -25.82 + 0.25*X – 0.00054*X2, sendo Y o número de baleias registradas e X o dia do ano. Para fazer este gráfico em R, usamos o comando curve: curve(-25.82 + 0.25*x - 0.00054*x^2, from=1, to=365, xlab="Dia juliano", ylab="Algo relacionado à abundância").

interpretarGLM_2

Sim, é uma parábola… Mas percebem que tem algo muito, muito, muito estranho nela? Especificamente no eixo Y… Ele mostra que a maior parte do ano temos um número negativo de baleias! Oxente, é possível uma coisa dessas? Será que temos anti-baleias habitando na costa da Bahia, e cada anti-baleia conta como abundância negativa? Bom, talvez na verdade é mais simples do que isso. Um GLM usa uma coisa chamada link function, ou função de ligação, para transformar os valores em algo que se aproxima de uma distribuição normal e facilitar o ajuste. Porque sim. :-) No caso do GLM com distribuição de Poisson, a função de transformação é logarítmica – portando precisamos de uma transformação exponencial para voltar à escala original. A fórmula então fica como Y=e-25.82 + 0.25*X – 0.00054*X^2; em R fica como curve(exp(-25.82 + 0.25*x - 0.00054*x^2), from=1, to=365, xlab="Dia juliano", ylab="Abundância").

interpretarGLM_3

Agora sim, não temos abundâncias negativas :-) A forma da curva mudou por causa da função exponencial; no fim do post tem os gráficos com os dados originais mostrando que o ajuste nem tá tão perfeito assim… Porque não se pode fazer coisas super sofisticadas com modelos quadráticas. Mas o momento do pico foi modelado até que com certa precisão.

E se quisermos ver a variação nas baleias neste mesmo ano em Serra Grande? Bom, aí adicionamos o efeito de Serra Grande a esta última curva: Y=e-25.82 + 0.25*X – 0.00054*X^2 – 1.75. Mais pro final do post mostro o gráfico. O que quero explorar agora é a parte das interações… Que é a parte muitas vezes mais interessante e mais complicadinha de explorar!

Então, acima eu falei dos efeitos dos anos: cada ano mostra o quanto aquele ano difere do ano-base, que é 2014. Vemos que em todos anos tivemos menos baleias, e este efeito foi mais intenso em 2019… Para Abrolhos! E também vimos que em Serra Grande tivemos menos baleias do que em Abrolhos… Em 2014! Agora para vermos as diferenças entre os locais para outros anos, ou as diferenças para os outros anos para Serra Grande, precisamos olhar também para as interações.

Por exemplo: o efeito de 2015 foi -0.24 (em um modelo linear significaria que tivemos 0.24 baleias a menos; em um GLM o efeito é negativo mas interpretar sua intensidade é um pouco mais complicado…), e o efeito de Serra Grande foi de -1.75. Mas o efeito combinado de 2015 e Serra Grande foi de 0.65! Então juntando tudo isso, para comparar as baleias de Serra Grande em 2015 com as de Abrolhos em 2014, fazemos -0.24 (efeito de 2015) – 1.75 (efeito de Serra Grande) + 0.65 (interação Serra Grande com 2015). A interação é positiva, o que significa que teve mais baleias do que seria esperado sem a interação. Isso é mais visível para 2019: o 1.76 da interação basicamente cancela o -1.75 do efeito do local, indicando que neste ano a quantidade de baleias nos dois locais foi similar.

Então olhando para os coeficientes, eu tiraria as seguintes conclusões gerais: A quantidade de baleias em Abrolhos foi diminuindo de 2014 pra 2019, com maior diminuição em 2019 (mas nenhuma diferença significativa entre 2014 e 2015); por outro lado, ela foi aumentando de 2014 pra 2019 em Serra Grande – somando o efeito do ano com a interação, temos efeitos positivos; e em Serra Grande no geral temos menos baleias que em Abrolhos, mas a diferença foi reduzindo com o tempo e em 2019 a quantidade foi similar (a interação e o efeito geral se cancelam).

Então, seguindo o modelo do código acima, vou fazer curvas para os anos de 2014 e 2019 para os dois locais:

par(mfrow=c(2,2), mar=c(3,3,2,2), oma=c(2,2,2,2)) # Cria espaço para 4 gráficos


curve(exp(-25.82 + 0.25*x - 0.00054*x^2), from=1, to=365, xlab="", ylab="", main="Arquipélago dos Abrolhos (2014)", ylim=c(0,25))


curve(exp(-25.82 + 0.25*x - 0.00054*x^2 -1.06), from=1, to=365, xlab="", ylab="", main="Arquipélago dos Abrolhos (2019)",
ylim=c(0,25))

curve(exp(-25.82 + 0.25*x - 0.00054*x^2 -1.75), from=1, to=365, xlab="", ylab="", main="Serra Grande (2014)", ylim=c(0,25))

curve(exp(-25.82 + 0.25*x - 0.00054*x^2 -1.75 -1.06 + 1.76), from=1, to=365, xlab="", ylab="", main="Serra Grande (2019)", ylim=c(0,25))

mtext(side=1,text="Dia juliano", outer=T, cex=1.5) # Legenda do eixo X para todos os gráficos
mtext(side=2, text="Número de baleias", outer=T, cex=1.5) # Legenda do eixo Y para todos os gráficos

interpretarGLM_4

No último gráfico, temos o efeito do local (Serra Grande), do ano (2019), e a interação entre local e ano (ano2019:localSerraGrande). Reparem como o pico de abundância sempre é mostrado no mesmo momento. Isso acontece porque não foi incluída no modelo interação do dia juliano com local ou ano. Incluindo essa interação, a forma da curva (não apenas seu tamanho) também poderia mudar.

E aqui os gráficos (feitos pela Bianca!) mostrano os números reais, para avaliarmos se as nossas curvas fazem sentido.

interpretarGLM_5

Bom, é isso por hoje e espero que tenham etendido e que tenha feito sentido e se quiserem perguntar algo (ou corrigir algo), comentem! :-)

8 pensamentos sobre “Interpretando os resultados de um GLM

  1. Olá!! Muito bom seu post, parabéns!! Me ajudou muito. Estou trabalhando com GLM agora e sofrendo um pouco, em especial com o uso de uma variável categórica com muitos níveis (22 para ser exato rsrs). Infelizmente não consigo agrupar muito para diminuir essa quantidade de níveis… e está muito complicado a interpretação dos resultados. Alguma dica? É possível trabalhar com tantos níveis? Obrigado!!

    Curtir

    • Oi Ina! Legal que você gostou do post :-) Olha, tem um pacote chamado emmeans que parece que faz comparações par-a-par. Tipo o teste de Tukey pra ANOVA, sabe? Mas pra GLM e similares. Eu nunca usei, mas talvez seja uma boa pra você! Fora isso, faça gráficos, gráficos lindos sempre ajudam :-)

      Curtir

  2. Boa sir,

    os exemplos com dragões, creio que possam interessar muito a este blog: https://www.blogs.unicamp.br/nasasasdodragao/

    Acho que é cabível investigar mais sobre as anti-baleias antes de presumirmos que o gráfico parabólico tenha pontos do domínio que sejam incoerentes…

    Apenas para esquentar mais um pouco a discussão, você sabe o que é um hiperplano? É um nome bonito para chamar em um espaço de N dimensões, um espaço de N-1 dimensões. Por exemplo, o hiperplano de um espaço de 4 dimensões, seria o espaço de 3 dimensões. O hiperplano de um espaço de 3 dimensões seria um espaço de 2 dimensões (ou o que entendemos por plano). O hiperplano de um espaço de 2 dimensões seria um espaço de 1 dimensão (ou o que entendemos por reta). Então se quiser dizer de forma chique, podemos dizer que este modelo das baleias pode ser expresso como um hiperplano no espaço bidimensional :P

    Curtir

    • Obrigado! Então, o que sei sobre a distribuição Gamma é que ela é usada pra variáveis contínuas com valores positivos (ou seja, maiores do que zero). Pode ser usada, por exemplo, pra descrever a distribuição de biomassa de indivíduos (nós podemos ter indivíduos com biomassa pequena, mas não com biomassa nula, porque aí não haveria um indivíduo né…) E tem a distribuição Tweedie, que é tipo uma mistura de Gamma com Poisson – pode ser usada pra modelar a distribuição total de biomassa em uma área, onde pode haver indivíduos ou não. Um detalhe sobre a distribuição Gamma é que ela usa a transformação inversa, então um coeficiente negativo representa uma relação positiva e vice-versa.

      Curtir

  3. Olá, muito boa sua explanação. Me ajudou muito. Mas gostaria de saber se vc conhece os software statistica ou jamovi. Como organizo a planilha para fazer a análise? Tenho Número de ovos de msoquitos para tres armadilhas em tres períodos verão, outono e inverno. Tenho também o número de ovos das tres armadilhas e dados do clima ( chuva, UR, temperatura máx, méd e mínima.

    Curtir

    • Obrigado!
      Eu conheço superificialmente o Jamovi. Fiz um vídeo introdutório: https://www.youtube.com/watch?v=LsbrcQaaI4o
      Eu imagino que para um GLM a organização da planilha seja a mesma que para uma regressão. Ou seja, cada coluna seria uma variável: identidade da armadilha, período, número de ovos, chuva, UR etc. Acho que tem post sobre organização de dados aqui, mas tem um post muito mais legal sobre organização de dados no blog do Marco Mello (blog Sobrevivendo na Ciência), procura lá ;-)

      Curtir

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 )

Foto do Google

Você está comentando utilizando sua conta Google. 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 )

Conectando a %s