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 | ||
cursos:planeco:roteiro:10-glmbinomial [2020/04/14 10:59] adalardo |
cursos:planeco:roteiro:10-glmbinomial [2020/04/19 14:21] adalardo |
||
---|---|---|---|
Linha 1: | Linha 1: | ||
- | |||
- | |||
- | |||
- | {{section>cursos:planeco:roteiro:10-glm#glm:_introdução}} | ||
- | |||
====== Modelos Lineares Generalizados: binomial ====== | ====== Modelos Lineares Generalizados: 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. | + | ===== GLM: introdução ===== |
- | + | <WRAP center round tip 60%> | |
- | ===== Função de ligação ===== | + | Essa introdução aos GLM é a mesma do tutorial [[cursos:planeco:roteiro:10-glm|]], caso já tenha feito, pode passar diretamente para o tópico [[cursos:planeco:roteiro:10-glmbinomial#glm:_binomial|GLM: binomial]] |
- | + | ||
- | A estrutura da função de ligação é a mesma para qualquer modelo: | + | |
- | + | ||
- | 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-ratio'', 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: | + | |
- | + | ||
- | $$ 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> | </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. | + | {{section>cursos:planeco:roteiro:10-glm#glm:_introdução}} |
- | * 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> | + | ===== GLM: componentes ===== |
- | <WRAP center round box 80%> | + | {{section>cursos:planeco:roteiro:10-glm#glm:_componentes}} |
- | {{: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. | ||
+ | ===== GLM: dispersão e acumulo de zeros ====== | ||
- | * use a variável ''prop'' como resposta (sucessos, falhas) | + | {{section>cursos:planeco:roteiro:10-glm#dispersão_e_acumulo_de_zeros}} |
- | * 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 ==== | + | ====== GLM: binomial ====== |
- | 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}$): | + | {{section>cursos:planeco:roteiro:10-glm#glm_binário_ou_proporção}} |
- | $$\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> | ||