Laboratório de Ecologia de Florestas Tropicais

Essa página está em modo somente de leitura. Você pode visualizar a fonte, mas não alterá-la. Informe-se com o administrador do Wiki, caso você ache que isso está incorreto.
====== Modelos Generalizados: binomial ======

===== GLM: introdução =====
<WRAP center round tip 60%>
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]]
</WRAP>


{{section>cursos:planeco:roteiro:10-glm#glm:_introdução}}


===== GLM: componentes =====

{{section>cursos:planeco:roteiro:10-glm#glm:_componentes}}
====== GLM: binomial ======

{{section>cursos:planeco:roteiro:10-glm#glm:_binomial}}


====== GLM: binomial OLD ======

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.

===== Função de ligação =====

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''((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>
Laboratório de Ecologia de Florestas Tropicais

Ferramentas do usuário

Ferramentas do site

Ferramentas da página
