Aqui você vê as diferenças entre duas revisões dessa página.
Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior Próxima revisão Ambos lados da revisão seguinte | ||
cursos:planeco:roteiro:10-glmbinomial [2020/04/14 11:19] adalardo [Função de ligação] |
cursos:planeco:roteiro:10-glmbinomial [2020/04/19 14:18] adalardo [GLM: componentes] |
||
---|---|---|---|
Linha 1: | Linha 1: | ||
- | ====== Modelos Generalizados: binomial ====== | + | ====== Modelos Lineares Generalizados: binomial ====== |
===== GLM: introdução ===== | ===== GLM: introdução ===== | ||
Linha 13: | Linha 13: | ||
{{section>cursos:planeco:roteiro:10-glm#glm:_componentes}} | {{section>cursos:planeco:roteiro:10-glm#glm:_componentes}} | ||
+ | |||
====== GLM: binomial ====== | ====== GLM: binomial ====== | ||
- | Os modelos de proporção de sucessos (sucessos/tentativas), proporção simple (%) ou de resposta binária (presença/ausência, vivo/morto) são modelados, normalmente, com estrutura do erro binomial. Nesses casos os limites dos valores da variável resposta é bem definido: entre 0 e 1. Além disso, a variância não é constante e varia conforme a média. Essas características fazem com que os resíduos apresentem uma estrutura que aumenta e depois diminuí, e normalmente o máximo de desvios é encontrado nos valores intermediários. | + | {{section>cursos:planeco:roteiro:10-glm#glm_binário_ou_proporção}} |
- | ===== Função de ligação ===== | + | ===== GLM: dispersão e acumulo de zeros ====== |
- | A estrutura da função de ligação é a mesma para qualquer modelo: | + | {{section>cursos:planeco:roteiro:10-glm#dispersão_e_acumulo_de_zeros}} |
- | + | ||
- | O preditor linear está associado à estrutura determinística do modelo e relacionado à linearização da relação, aqui definido como $\eta$: | + | |
- | + | ||
- | $$ \eta = \alpha + \beta x$$ | + | |
- | + | ||
- | A função de ligação é o que relaciona o preditor linear com a esperança do modelo: | + | |
- | + | ||
- | $$ \eta = g(E_{(y)}) $$ | + | |
- | + | ||
- | A função de ligação $g()$ para modelos com resposta binária ou proporção é chamada de ''logit'' ou ''log odds''((chance = $\frac{p}{1-p}$)), definida como: | + | |
- | + | ||
- | $$ \eta = \log(\frac{p}{1-p})$$ | + | |
- | + | ||
- | Para reverter o preditor linear da função logit para a escala de observação usa-se a função inversa: | + | |
- | + | ||
- | $$ \text{logit}^{-1} = \frac{e^{\eta}}{(1+ e^{\eta})} $$ | + | |
- | + | ||
- | ===== Resposta: proporções ===== | + | |
- | + | ||
- | ==== Exemplo: floração ==== | + | |
- | + | ||
- | {{:cursos:planeco:roteiro:FloradaIpesPantanal.png?200 | }} | + | |
- | Mais um exemplo apresentado no livro do Michael Crawley, //The R Book//. | + | |
- | Neste experimento o objetivo foi avaliar a floração de 5 variedades de plantas tratadas com hormônios de crescimento (6 concentrações). Depois de seis semanas as plantas foram classificadas em floridas ou vegetativas. | + | |
- | + | ||
- | + | ||
- | <WRAP center round box 60%> | + | |
- | + | ||
- | **//__Conjunto de Dados__//**: ''flowering.txt'' | + | |
- | * **//flowered//**: número de plantas que floresceram | + | |
- | * **//number//**: número de plantas acompanhadas | + | |
- | * **//dose//**: concentração da dose de hormônio | + | |
- | * **//variety//**: variedade da planta (categórica 5 níveis) | + | |
- | </WRAP> | + | |
- | + | ||
- | ==== Hipótese ==== | + | |
- | + | ||
- | O objetivo do estudo que gerou esses dados é saber se o evento de floração é influenciado pelo dose de hormônio e a variedade da planta. | + | |
- | * baixe o arquivo {{ :cursos:planeco:planeco:roteiro:flowering.txt |}} | + | |
- | * abra os dados no Rcmdr (a separação de campo é espaço) com o nome ''flower'' | + | |
- | * crie a variável ''prop'' pelo menu **Data> Manage variables in active data set> Compute new variable...**, colocando no campo **Expression to compute**: | + | |
- | + | ||
- | <code> | + | |
- | cbind(sucess = flowered, fail = number - flowered) | + | |
- | + | ||
- | </code> | + | |
- | + | ||
- | <WRAP center round box 80%> | + | |
- | {{:cursos:planeco:roteiro:newVarFlower.png?700|}} | + | |
- | </WRAP> | + | |
- | + | ||
- | Esse comando acima cria uma nova variável nos dados **flower** chamada **prop**. Essa nova variável tem duas colunas (**sucess e fail**) contendo o número de plantas floridas e o número de plantas que não floresceram, respectivamente. | + | |
- | + | ||
- | + | ||
- | * use a variável ''prop'' como resposta (sucessos, falhas) | + | |
- | * monte o modelo cheio com todas a variáveis preditoras e interações | + | |
- | * simplifique o modelo para o mínimo adequado | + | |
- | + | ||
- | <WRAP center round important 60%> | + | |
- | ** Use os mesmos passos do modelo anterior no Rcmdr ** | + | |
- | * lembre-se que a ''family'' nesse caso é ''binomial'' | + | |
- | * o procedimento para a sobre-dispersão é o mesmo que no exemplo anterior | + | |
- | + | ||
- | </WRAP> | + | |
- | + | ||
- | ==== Interpretação do resultado ==== | + | |
- | + | ||
- | Para interpretar tanto os coeficientes quanto os valores previsto é necessário aplicar a função inversa do ''logit'', ou seja, nosso modelo faz previsões na escala de log(odds-ratio), nosso preditor linear $\eta$, e precisamos retornar para a escala de observação que é a probabilidade de florescer ($\hat{y}$): | + | |
- | + | ||
- | $$\hat{y} = \frac{e^{\hat{\eta}}}{1+e^{\hat{\eta}}} $$ | + | |
- | + | ||
- | <WRAP center round box 80%> | + | |
- | + | ||
- | * calcule o predito pelo modelo e os coeficientes na escala original | + | |
- | * interprete o efeito da concentração na floração das variedades | + | |
- | </WRAP> | + | |
- | + | ||
- | + | ||
- | <WRAP center round tip 80%> | + | |
- | + | ||
- | **__Transformar os coeficientes e valores preditos pelo GLM:__** | + | |
- | + | ||
- | Para transformar o valor predito pelo modelo (log(odds-ratio)) na escala de medida (proporção) é preciso transformar os preditos pelo modelo. Para predizer na escala de medida usamos a função ''predict'', como no código abaixo. O predito pelo modelo, está na escala do preditor linear, portanto devemos transformar essa medida com a função inversa da logit, como no código abaixo. <wrap hi>Lembre-se de mudar, no código, o "nomedomodelo"</wrap> pelo nome que usou quando construiu o glm. | + | |
- | + | ||
- | <code> | + | |
- | + | ||
- | (preditoLinear <- predict("nomedomodelo")) | + | |
- | (preditoProp <- exp(preditoLinear)/(1+ exp(preditoLinear))) | + | |
- | + | ||
- | </code> | + | |
- | + | ||
- | A própria função ''predict'', também faz o serviço completo se colocarmos o argumento ''type="response"'', como abaixo: | + | |
- | + | ||
- | <code> | + | |
- | + | ||
- | predito <- predict("nomedomodelo", type = "response") | + | |
- | predito | + | |
- | + | ||
- | </code> | + | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | </WRAP> | + | |
- | + | ||
- | === Gráfico e interpretação dos resultados === | + | |
- | + | ||
- | Para um gráfico dos resultados use o menu: | + | |
- | + | ||
- | **Models > Graphs > Predict effect plots...** | + | |
- | + | ||
- | <WRAP center round todo 60%> | + | |
- | <wrap hi> | + | |
- | A partir dos gráficos e do modelo selecionado faça um relato (5 linhas) das interpretações biológicos. Esse relato, junto ao resultado e gráficos, deve ser enviado aos professores ao final da atividade.</wrap> | + | |
- | + | ||
- | + | ||
- | </WRAP> | + | |
- | + | ||
- | /* | + | |
- | + | ||
- | <WRAP center round tip 80%> | + | |
- | + | ||
- | **__Um gráfico | + | |
- | + | ||
- | ** Verifique se o pacote ''ggplot2'' está disponível para carregar em:** | + | |
- | No menu **Tools**> **Load package(s)...** | + | |
- | + | ||
- | * caso não esteja copie e cole o código abaixo para instalá-lo: | + | |
- | + | ||
- | install.packages(ggplot2) | + | |
- | + | ||
- | * construa a variável de proporção: | + | |
- | + | ||
- | flower$pf <- flower$flowered/flower$number #Proporção de flores que floresceram | + | |
- | + | ||
- | * construa o gráfico: | + | |
- | + | ||
- | <code> | + | |
- | ggplot(data = flower, aes(x = dose, y = pf, fill = variety)) + | + | |
- | geom_point(shape = 21, size = 2) + | + | |
- | scale_x_continuous(name = "Dose de hormônio de crescimento", | + | |
- | breaks = c(0,10,20,30), labels = c("0","10","20","30")) + | + | |
- | scale_y_continuous(name = "Proporção de indivíduos florescendo", | + | |
- | breaks = c(0.00,0.25,0.50,0.75,1.00), labels = c("0.00","0.25","0.50","0.75","1.00")) | + | |
- | + geom_smooth(data = flower, aes(colour = variety), | + | |
- | method = "glm", method.args = list(family = quasibinomial), se = FALSE) | + | |
- | + | ||
- | </code> | + | |
- | + | ||
- | + | ||
- | </WRAP> | + | |
- | + | ||
- | */ | + | |
- | + | ||
- | ===== Resposta: binária ===== | + | |
- | + | ||
- | ==== Exemplo: pássaro na ilha ==== | + | |
- | + | ||
- | O conjunto de dados que vamos usar, {{ :planeco:roteiro:isolation.txt |}} tem como variável: | + | |
- | <WRAP center round box 60%> | + | |
- | **__//Conjunto de dados//__**: ''isolation.txt'' | + | |
- | + | ||
- | * **//incidence//**: presença/ausência da espécie de ave (reprodução) | + | |
- | * **//area//**: área total da ilha ($km^2$) | + | |
- | * **//isolation//**: distância do continente (km) | + | |
- | </WRAP> | + | |
- | + | ||
- | + | ||
- | ==== Hipótese ==== | + | |
- | + | ||
- | O objetivo do estudo que gerou esses dados é saber se a ocorrência da ave (reprodução) está relacionada com o isolamento e tamanho da ilha. | + | |
- | + | ||
- | * abra os dados ''isolation.txt'' no Rcmdr (a separação de campo é espaço) | + | |
- | * monte o modelo cheio com todas a variáveis preditoras e interações | + | |
- | * simplifique o modelo para o mínimo adequado | + | |
- | + | ||
- | <WRAP center round important 60%> | + | |
- | ** Use os mesmos passos do modelo anterior no Rcmdr ** | + | |
- | * lembre-se que a ''family'' nesse caso é ''binomial'' | + | |
- | * o procedimento para a sobre-dispersão é o mesmo que no exemplo anterior | + | |
- | + | ||
- | </WRAP> | + | |
- | + | ||
- | + | ||
- | ==== Interpretação do resultado ==== | + | |
- | + | ||
- | O modelo prevê a ocorrência da ave na escala de logaritmo da chance (log odds-ratio). Para interpretar tanto os coeficientes quanto os valores previsto é necessário aplicar a função inversa do ''logit'', como no exercício anterior: | + | |
- | + | ||
- | <WRAP center round box 80%> | + | |
- | + | ||
- | * calcule o predito pelo modelo e os coeficientes na escala original | + | |
- | * interprete o efeito do tamanho e distância na ocorrência da espécie | + | |
- | + | ||
- | </WRAP> | + | |
- | + | ||
- | /* | + | |
- | <WRAP center round tip 100%> | + | |
- | + | ||
- | Caso já tenha o pacote ''ggplot2'' instalado, carregue-o usando o menu **Tools**> **Load package(s)...**. Caso não tenha, instale com os seguinte comando: | + | |
- | + | ||
- | install.packages("ggplot2") | + | |
- | + | ||
- | Em seguida retorne ao menu para carregá-lo. | + | |
- | + | ||
- | </WRAP> | + | |
- | + | ||
- | <WRAP center round tip 100%> | + | |
- | + | ||
- | === Fazendo o gráfico === | + | |
- | + | ||
- | == Ocorrência dependendo da Área == | + | |
- | + | ||
- | <code> | + | |
- | ggplot(isolation, aes(x = area, y= incidence))+ | + | |
- | geom_point(shape=21, fill="white") + | + | |
- | geom_smooth(method = "glm", method.args = list(family = binomial), se = FALSE) | + | |
- | + | ||
- | </code> | + | |
- | + | ||
- | == Ocorrência dependendo do Isolamento == | + | |
- | + | ||
- | <code> | + | |
- | ggplot(isolation, aes(x = isolation, y= incidence))+ | + | |
- | geom_point(shape=21, fill="white") + | + | |
- | geom_smooth(method = "glm", method.args = list(family = binomial), se = FALSE) | + | |
- | + | ||
- | </code> | + | |
- | + | ||
- | </WRAP> | + | |
- | + | ||
- | */ | + | |
- | + | ||
- | <WRAP center round todo 80%> | + | |
- | + | ||
- | <wrap hi> | + | |
- | + | ||
- | **__ O que deve entregar?__** | + | |
- | + | ||
- | + | ||
- | </wrap> | + | |
- | + | ||
- | Para cada exercício feito, deve ser entregue, em um único arquivo: | + | |
- | + | ||
- | * o resultado do modelo mínimo adequado | + | |
- | * os coeficientes estimados, na escala de observação | + | |
- | * gráficos que apresentem os resultados principais | + | |
- | * um relato de no máximo 5 linhas, ou em tópicos, da interpretação biológica dos resultados | + | |
- | + | ||
- | + | ||
- | </WRAP> | + | |
- | + | ||
- | ====== Sobredispersão e acumulo de zeros ====== | + | |
- | + | ||
- | Os modelo GLM poisson e binomial apresentam a variância acoplada à média dos valores, diferentemente dos modelos com distribuição normal onde a média e a variância são idependentes. Caso haja uma variação maior ou menor nos dados do que o previsto por essas distribuições, o modelo não consegue dar conta. Essa sobre-dispersão ou sub-dispersão dos dados indica que temos mais ou menos variação do que é predito pelos modelos. Isso pode ser decorrência de vários fontes de erro na definição do modelo, alguns exemplos são: | + | |
- | * o resíduo dos dados pode não ter sido gerado por um processo aleatório poisson ou binomial | + | |
- | * há mais variação do que predito pela ausência de preditoras importantes | + | |
- | * muitos zeros, além do predito pelas distribuições, em decorrência de diferentes processos: um que gera a ausência e outro que gera a variação nas ocorrências de sucesso | + | |
- | + | ||
- | <WRAP center round tip 100%> | + | |
- | **__ Soluções para a sobre-dispersão e acumulo de zeros__** | + | |
- | + | ||
- | + | ||
- | A solução mais simples para lidar com sobre-dispersão são os modelo quasipoisson e quasibinomial, que estimam um parâmetro a mais, relacionando a média à variância, o parâmetro de dispersão. Entretanto, os modelos ''quasi'' dão conta apenas de sobre-dispersões moderadas e não indicam qual a fonte dela. | + | |
- | Há algumas alternativas ao modelo ''quasi'' para a sobre-dispersão dos dados, alguns deles estão listados abaixo: | + | |
- | * modelo binomial negativo | + | |
- | * modelo de mistura, considerando dois processos distintos | + | |
- | * modelos mistos, considerando a ausência de independência das observações | + | |
- | * modelos com acúmulos de zeros (Zero Inflated Models). | + | |
- | + | ||
- | Não é objetivo deste curso mostrar todas essas alternativas, mas caso se deparem com esse problema, muito frequente na área da biológica, saibam que existem alternativas robustas para solucioná-lo. | + | |
- | + | ||
- | + | ||
- | </WRAP> | + | |