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 Lineares Generalizados ======


===== GLM: Introdução =====

<WRAP center round box 80%>

{{ youtube>32u0vrXs8AQ |}}


</WRAP>

Os modelos lineares generalizados (**GLMs**) são uma ampliação dos modelos lineares ordinários. Os **GLM's** são usados quando os resíduos (erro) do modelo apresentam distribuição diferente da normal (gaussiana). 
A natureza da variável resposta é uma boa indicação do tipo de distribuição de resíduos que iremos encontrar nos modelos. Por exemplos, variáveis de contagem são inteiras e apresentam os valores limitados no zero. Esse tipo de variável, em geral, tem uma distribuição de erros assimétrica para valores baixos e uma variância que aumenta com a média dos valores preditos, violando duas premissas dos modelos lineares.  Os casos mais comuns de modelos generalizados são de variáveis resposta de contagem, proporção e binária, muito comum nos estudos de ecologia e evolução.

<WRAP center round box 50%>

** Devemos considerar os GLMs principalmente quando a variável resposta é expressa em: **
  * contagens simples
  * contagem expressa em proporções 
  * número de sucesso e tentativa
  * variáveis binárias (ex. morto x vivo)
  * tempo para o evento ocorrer (modelos de sobrevivência)

</WRAP>

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

Uma das formas de entendermos os modelos generalizados é separar o modelo em dois componentes: a relação determinística entre as variáveis (resposta e preditora) e o componente aleatório dos resíduos (distribuição dos erros). Em um modelo linear ordinário a relação entre as variáveis é uma proporção constante, o que define uma relação funcional de uma reta. Quando temos uma contagem, essa relação pode ter uma estrutura funcional de uma exponencial. Para esses casos, os modelos generalizados utilizam uma função de ligação ''log'' para linearizar a relação determinística entre as variáveis.  Portanto, a estrutura determinística dos modelos **GLM's** é definida por um preditor linear, associada à função de ligação. 

O componente aleatório dos resíduos, no caso de uma variável de contagem, segue, em geral, uma distribuição **poisson**. A distribuição **poisson** é uma variável aleatória definida por apenas um parâmetro ($\lambda$), equivalente à média, chamada de ''lambda''. A distribuição **poisson** tem uma característica interessante, seu desvio padrão é igual à média. Portanto, se a média aumenta, o desvio acompanha esse aumento e a distribuição passa a ter um maior espalhamento.



==== Preditor linear e função de ligação ====

O preditor linear está associado à estrutura determinística do modelo e está 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^{-1}(E_{(y)}) $$

Ou seja, nos modelos generalizados não é a variável resposta que tem uma relação linear com a preditora, e sim o preditor linear que tem uma relação linear com as preditoras.

==== Funções de ligações canônicas ====
Para alguns tipos de famílias de variáveis temos funções de ligações padrões. As mais usadas são: 

<WRAP center round tip 60%>

^ Natureza da resposta ^ Estrutura dos resíduos (erro) ^ Função de ligação ^
| contínua | normal | identidade|
| contagem | poisson | log| 
| proporção | binomial| logit| 

</WRAP>





====== GLM Contagem ======
<WRAP center round box 80%>

{{ youtube>RvCOOS7iHDk |}}
</WRAP>

===== Contagem: um exemplo simples =====

Um exemplo, apresentado no livro do Michael Crawley, The R Book, relata a contagem de espécies de árvores em unidades amostrais de florestas com diferentes biomassa e classificadas em três níveis de __ph__ no solo: baixo, médio e alto. O objetivo desse experimento não manipulativo é verificar se há relação entre riqueza de árvores e as preditora biomassa da floresta e __ph__ do solo.

<WRAP center round todo 80%>
__**ATIVIDADE**__

==== Modelo Linear Múltiplo  (LM) ==== 

  -  Importe o arquivo {{ :cursos:planeco:roteiro:species.txt |}} para o Rcmdr. Note que esse arquivo tem como separador de campo a tabulação e decimal como ponto.
  -  Monte o modelo linear clássicos (''lm'') para esse dados, tendo como variável resposta a riqueza de espécies(''Species'') e como preditoras o ''pH'' e ''Biomass'' e as interações possíveis.
  -  Reduza o modelo cheio ao modelo mínimo adequado utilizando como critério de comparação a tabela de anova.
  - Utilizando os coeficientes estimados do modelo, faça a predição do número de espécies para:
      *  um nível alto de ''pH'' com ''Biomass'' de **3.2**
      *  um nível médio de ''pH'' com ''Biomass'' de **15.5**
      *  um nível baixo de ''pH'' com ''Biomass'' de **7.2**


/*
Faça o diagnóstico do resíduos do modelo mínimo adequado.
*/

 
==== Modelo Linear Generalizado (GLM) ====

  - Repita o procedimento de simplificação a partir do modelo cheio, agora com modelo linear generalizado (''glm'') e com ''family = poisson''. 
    * <wrap hi>Caso o Rcmdr não retorne o ''p-valor'' na comparação de modelos por anova, copie a linha de código que foi utilizada com ''anova(...)'' e cole novamente incluindo ''anova(..., test="Chisq")''</wrap>
  - Calcule as mesmas predições acima para o modelo, usando os coeficientes do preditor linear do ''glm''.
  - Transforme os preditos pelo modelo de volta para a escala de observação ((note que é preciso primeiro calcular o predito na escala do preditor linear e depois transformar, o que não é a mesma coisa que transformar os coeficientes e depois calcular o predito)).
  - Faça os gráficos apresentados no tópico [[#gráfico_no_rcmdr| Gráfico no Rcmdr]]


  
<WRAP center round tip 90%>



  * Para a predição no ''glm'' utilize os coeficientes estimados pelo modelo. 
  * Após estimar o predito na escala linear, transforme a predição para a escala de observação. 
  * Como usamos o ''log'' como função de ligação, para retornar a escala da observação devemos utilizar o antilog, no caso, a função exponencial.
 
  
</WRAP>


</WRAP>




===== Gráfico no Rcmdr =====

==== Gráfico dos dados ====

No menu **Graphs**, selecione **XY conditioningh Plot** e selecione as varíáveis, definindo **//ph//** como variável de agrupamento, como no gráfico abaixo.


<WRAP center round box 80%>

{{:cursos:planeco:roteiro:speciesPlot.png?800|}} 

</WRAP>

==== Gráfico dos Modelos ====

No menu **Models>Graphs** selecione **Predict effect plots...** e selecione as variáveis.

<WRAP center round box 80%>
{{:cursos:planeco:roteiro:speciesLMres.png?800|}}
</WRAP>

<WRAP center round tip 80%>
**__Ordenando uma categórica__**


O padrão do R é ordenar as variáveis categóricas por ordem alfabética. No exemplo seria desejável reordenar a variável categórica **//ph//** em uma categórica ordenada  **low>medium>high**. Para reordenar utilize o menu **Data>Manager variable in active data set> Reorder factor levels**. Caso não deseje sobrescrever a variável original, forneça um novo nome para a variável reordenada.

</WRAP>

==== O que preciso entregar ====


<WRAP center round help 100%>
 
 
  * Preencha as perguntas do quadro abaixo ou pelo [[https://forms.gle/Y3Cka1kHQ3to354E9|link do formulário]]




{{url>https://forms.gle/Y3Cka1kHQ3to354E9}}

</WRAP>
===== Contagem: o que faz um aluno faltar às aulas =====

Vamos utilizar um exemplo que está presente no livro de W. Venables e B. Ripley, __Modern Applied Statistics with S-PLUS__((já não tão moderno assim, já que foi publicado pela primeira vez em 1999)), sobre  o número de dias ausentes da escola de crianças na Austrália. 
==== Carregando o pacote MASS ====


No Rcommader (Rcmdr) vá ao menu ** Tools ** > **Load package(s)** e selecione o pacote ''MASS''. 
<wrap hi>

\\

Caso o pacote não apareça listado, significa que ele já está carregado</wrap>, então pule esse passo.

<WRAP center round box 80%>
{{:planeco:roteiro:glm01.png?500|}}

</WRAP>

==== Lendo os dados: quine ====

Em sequida:
  *  abra o menu **Data**  > ** Data in packages** > ** Read data from an attached package...**
  * selecione o pacote **MASS ** e os dados **quine** ((deixe o nome do dado como quine))

<WRAP center round box 80%>
{{:planeco:roteiro:glm02.png?500|}}
</WRAP>


==== Entendendo os dados: quine ====

Os dados estão relacionados ao estudo para entender quais variáveis estão relacionados à ausência (falta) do aluno na escola. A observação está relacionada a alunos amostrados aleatoriamente de escolas na Austrália.


<WRAP center round box 80%>
  * **Days**: variável resposta, número de dias ausente da escola 
  * **Eth**: origem aborígene (A) ou não (N)
  * **Sex**: homem (M) ou mulher (F)
  * **Age**: estágio de educação F0(primário)... quatro níveis.
  * <wrap hi>  **Lrn**: classificação de aprendizado do aluno médio (AL) e fraco (SL)</wrap>((essa variável tem algumas complicações adicionais e por isso vamos deixá-la de lado))
</WRAP>


==== Gráfico dos dados ====


<WRAP center round tip 80%>


O pacote RcmdrPlugin.KMggplot2 é um plugin para Rcmdr que amplia as funções gráficas da interface.
Instale o pacote copiando o comando abaixo no box superior do Rcmdr:

<WRAP center round alert 80%>

  * **guarde os resultados dos modelos fora do Rcmdr pois a instalação e o carregamento do pacote solicita a reinicialização do Rcmdr**
  * **após a instalação e carregamento do pacote, confira se os dados permanecem ativos, confira se precisará carregá-lo novamente **


</WRAP>


<code>install.packages("RcmdrPlugin.KMggplot2") </code>

Em seguida, garanta que o cursor do mouse está na linha de comando e clique no botão **Submit**. Na janela que ira se abrir selecione o repositório **Brasil(SP1)**.

<WRAP center round box 90%>
{{  :cursos:planeco:roteiro:instalpkg.png?700|  }}
</WRAP>
 

Para ativar o plugin selecione o menu **Tools> Load Rcmdr plug-in(s)...** e em seguida selecione o pacote ''RcmdrPlugin.KMggplot2''.

<WRAP center round box 90%>
{{:cursos:planeco:roteiro:loadPluginRcmdr.png?700|}}
</WRAP>

  * clique em sim na janela que solicita a reinicialização do Rcmdr
  * clique na nova opção do menu **KMggplot2 > BoxPlot/...** e selecione as variáveis

<WRAP center round box 90%>
{{:cursos:planeco:roteiro:graficoBoxPlots.png?700|}}
</WRAP>


</WRAP>







==== Ajustando o GLM: dias fora da escola ====

<WRAP center round todo 80%>

**__Atividade__**

Para nosso exercício vamos deixar de lado a variável ''Lrn'' por que há dados faltantes nela com relação a outras variáveis. Vamos construir o modelo cheio com a variável resposta ''Days'' e com  as variáveis preditoras (''Eth'', ''Sex'', ''Age'' ) e todas as possibilidades de interações entre elas. Como estamos tratando de uma variável de contagem podemos partir direto para um modelo **GLM**  indicando a família de distribuição de resíduos  **POISSON** e a função de ligação **//log//**.


  * abra o menu **Statistics** > **Fit model** > ** Generalized Linear Model **
  * construa um modelo cheio com (**Age, Eth e Sex**) e as suas interações possíveis:

<code> Days ~ Eth + Sex + Age + Eth:Sex + Eth:Age + Sex:Age + Eth:Sex:Age </code> 

  * faça a simplificação do modelo para reduzir o modelo ao mínimo adequado
  
 

<WRAP center round box 90%>
{{:cursos:planeco:roteiro:quineGLM.png?700|}}
</WRAP>


</WRAP>

/*
<WRAP center round important 60%>
O nosso modelo cheio não conseguiu estimar alguns dos parâmetros. Isso se deveu ao fato de algumas combinações de níveis de fatores não foram encontradas na amostra. Por exemplo não há nenhum:

  * aluno de etnia Aborígene do sexo feminino no nível máximo de escolaridade com desempenho fraco! 
  * aqui o primeiro passo é jogar fora a interação de quarto nível e prosseguir 
 
</WRAP>
*/



==== Diagnóstico do modelo  ====

Um dos pressupostos do modelo Poisson é que a variância aumenta linearmente com a esperança (média do modelo). Podemos avaliar isso dividindo a ''Residual Deviance '' pelo seu ''degrees of freedom''. Essa razão deve ser próxima a 1. O que não é o caso do nosso modelo. Nesses casos uma das alternativas é:

  * ajustar o modelo usando ** Family**: ''quasipoisson''

 <WRAP center round todo 90%>

==== Ajustando o GLM com sobredispersão ====


 <WRAP center round box 80%>

  * monte o modelo cheio utilizando a família ''quasipoisson'' e 
  * siga em frente simplificando o modelo para o mínimo adequado 
  * interprete o modelo selecionado

</WRAP>


</WRAP>


=== Gráfico do Modelo ===

O gráfico do modelo pode ser obtido no Rcmdr da mesma forma indicada no modelo anterior, no menu:
**Models>Graphs** selecione **Predict effect plots...** e selecione as variáveis.

/*
==== Comparando o GLM-poisson e o LM com dados transformados ====


Agora, vamos comparar os resultados do modelo GLM desenvolvido acima com um modelo linear com distribuição de erro normal mas agora com a variável resposta transformada em log. Ou seja, primeiro vamos transformar a nossa variável resposta em log e depois ajustar um modelo linear com distribuição de erros Normal.

<WRAP center round box 80%>
Faça a comparação dos resultados dos modelos.
Faça também o diagnóstico do modelo com a variável resposta em log.

</WRAP>
*/

==== O que preciso entregar ====


<WRAP center round help 100%>
 
 
  * Preencha as perguntas do quadro abaixo ou pelo [[https://forms.gle/WdAzcHBpF25NeMkd9|link do formulário]]




{{url>https://forms.gle/WdAzcHBpF25NeMkd9}}

</WRAP>
====== GLM Binário ou proporção ======

Em muitos fenomenos naturais a variável de interesse tem apenas dois estados possíveis. Essas variáveis de resposta binária (''presença'' x ''ausência''; ''vivo'' x ''morto'', ''germinou'' x ''não germinou'') tem uma natureza muito distinta das medidas contínuas ou as contagens que vimos anteriormente. Tradicionalmente essa natureza de variável binária foi definida como a probabilidade de sucesso. Neste tópico iremos tratar essas variáveis respostas em **__Modelos Lineares Generalizados Binomiais__** 


<WRAP center round box 60%>

{{ youtube>RabMP1alhe0 |}}

</WRAP>
Os modelos de proporção (ou probabilidade) de sucessos (''sucessos''/''tentativas''), proporção simple ou de resposta binária 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 seu valor difere conforme varia a probabilidade de sucesso. Estas características fazem com que os resíduos apresentem uma estrutura que aumenta e depois diminuí com o aumento da probabilidade de sucessos e o máximo de variância é encontrado nos valores intermediários (probabilidade de sucesso = 0.5).

A função Bernoulli, que é a base para a binomial, é definida pelo parâmetro de probabilidade de sucesso em um evento com apenas duas possibilidades de resultado (binário). O parâmetro da função Bernoulli é a probabilidade de sucesso.  No caso de uma moeda justa seria a probabilidade de ''0,50'' de sair coroa ((ou coroa, dependendo do que chamamos de sucesso)). 

A binomial é uma generalização da Bernoulli, definida pelo número de sucessos em certo número (''n'') de tentativas (número de eventos Bernoulli). Um exemplo de um experimento binomial seria a germinação (sucessos) de sementes em um experimento onde temos 20 sementes (número de tentativas). Neste experimento o que tentamos estimar é a probabilidade de germinação (''sucessos''). Elaborando um pouco mais esse experimento podemos incluir uma variável preditora contínua como a úmidade e/ou uma categórica como o tipo de solo, e perguntar como essas variáveis afetam a ''probabilidade de sucesso'', no caso a germinação.

<WRAP center round tip 80%>

** Conceitos Importantes **

  * n = número de tentativas  
  * s = número de sucessos
  * f = número de falhas

<code> Probabilidade de sucesso </code>

$$p = \frac{s}{n} $$

<code>Probabilidade de falha </code>

$$q = \frac{f}{n} $$

$$q = 1 - p$$


<code> Chance de sucesso (Odds) </code>

$$\text{odds} = \frac{s}{f} $$

$$\text{odds} = \frac{p}{1-p} $$

<WRAP center round box 90%>

Note como a __**chance**__ de ocorrência de um evento é a probabilidade de ocorrência deste evento dividida pela probabilidade da não ocorrência do mesmo evento. 

A __**chance**__  é muito usada em apostas, quando, por exemplo, dizemos que a chance de um time vencer é de ''4:1'' ((ou que está pagando 4:1)), ou seja, a probabilidade de vencer é ''4x'' maior do que a de perder. O conceito de chance é muito importante nos modelos binomiais e devemos evitar confundi-lo com probabilidade. Chance e probabilidade são escalas distintas para medir a ocorrência de sucessos. 

</WRAP>

</WRAP>

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

A estrutura da função de ligação é a mesma para qualquer modelo generalizado, o que muda é o tipo de função aplicada:

O preditor linear $\eta$ está associado à estrutura determinística e relacionado à sua estrutura linear.

$$ \eta = \alpha + \sum\beta_{x_i}$$

Note que: 

$$\alpha + \sum\beta_{x_i}$$

É a estrutura determinística do modelo modelo linear, agora não mais relacionado diretamente à escala da variável resposta.  

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()$ canônica ou padrão para modelos com resposta binária ou proporção é chamada de ''logit'' ou ''logaritmo da chance''((log odds ou log chance)), definida como:


$$ \eta = \log(\frac{p}{1-p})$$

$$ \log(\frac{p}{1-p}) = \alpha + \sum\beta_{x_i} $$

Sendo $\frac{p}{1-p}$  a **chance** ou **odds** em inglês.
 
Para reverter o preditor linear da função logit para a escala de observação usa-se a função inversa:

$$ g^{-1} = \text{logit}^{-1} = \frac{e^{\eta}}{(1+ e^{\eta})} $$


==== Chance e Razão de Chance ====

O predito pelo modelo na escala do preditor linear do modelo binário com função de ligação //**logit**// está na escala de logaritmo da chance ($\log(\frac{p}{1-p})$). 
Dado que, para variáveis categóricas os coeficientes do modelo são relacionados às diferenças entre o nível do tratamento e o controle:


$$\text{exp}(\log(\text{odds}_{\text{trat}}) - \log(\text{odds}_{\text{control}})) = \frac{\text{odds}_{\text{trat}}}{\text{odds}_{\text{control}}}  $$

então, exponenciar os coeficientes do modelo binomial com preditora categórica transforma os coeficientes em razão de chance comparado com o nível basal ((lembre-se que as categóricas são transformadas em variáveis indicadoras ou dummy e um dos níveis é transportado para o intercepto do modelo, sendo esse o nível basal ou controle)).

A //**razão de chance**//  mede o quanto uma chance é proporcionalmente diferente de outra, geralmente comparando com um nível controle. Ou seja, qual a proporção de mudança na chance do tratamento em relação a chance do controle. Pensando em nosso experimento de germinação tendo o ''solo arenoso'' como nível de referência, a razão de chance do ''solo arenoso'' seria o quanto a chance de germinar no ''arenoso'' proporcionalmente maior que a chance de germinar no solo ''arenoso''.  

Parece complicado, mas é apenas por falta de intimidade com essas escalas, a razão de chance é uma medida muito popular em outras áreas da ciência, como medicina. Importante lembrar que a ''razão de chance'' mede o efeito proporcional em relação ao nível de referência. 

No caso de variáveis contínuas a ''razão de chance'' é relacionada à chance de ''x+1'' comparada com ''x'', ou seja, qual a proporção de mudança na chance com o aumento de uma unidade da variável contínua preditora.

Portanto, uma forma de interpretar os coeficientes do modelo binomial é exponenciar e interpretar como razão de chance, sendo o intercepto a chance do nível basal da variável categórica ou a chance quando a variável contínua é zero.  


===== GLM: resposta binária =====

==== Exemplo: pássaro na ilha ====
<WRAP center round box 80%>

{{ youtube>x6eSQ_6HfKo |}}

</WRAP>

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 coletou esses dados foi saber se a ocorrência da ave está relacionada com o isolamento e tamanho da ilha.

<WRAP center round todo 80%>

**__ATIVIDADE__**

    * abra os dados ''isolation.txt'' no Rcmdr (a separação de campo é tabulaçã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 80%>
** Importante: **
  * lembre-se que a ''family'' nesse caso é ''binomial''
  * os modelos com variáveis resposta binárias não tem problema com sobre-dispersão!!!

</WRAP>

</WRAP>

 
==== Interpretação do resultado ====

O modelo prevê a ocorrência da ave na escala de logaritmo da chance (log odds-ratio). Para os coeficientes estimados pelo modelo o melhor é aplicar a função ''exp'' e interpretá-los como ''razão de chance'' entre categorias ou entre ''x+1'' e ''x''. 
Para interpretar 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 ocorrência ($\hat{y}$):

$$\hat{y} = \frac{e^{\hat{\eta}}}{1+e^{\hat{\eta}}} $$ 



 <WRAP center round todo 80%>

**__ATIVIDADE__**

  * calcule o predito pelo modelo na escala de probabilidade de ocorrência para uma ilha de 5.6 Km² e distante 7.2 Km da costa.
  * quanto varia a chance de ocorrência se aumentar 1 Km² no tamanho da ilha?
  * e se aumentar 1 Km no isolamento?
  * faça uma interpretação biológica do modelo selecionado baseado nos seus coeficientes.

</WRAP>

==== O que preciso entregar ====


<WRAP center round help 100%>
 
 
  * Preencha as perguntas do quadro abaixo ou pelo [[https://forms.gle/dUqZad77fU8Hrqrm8|link do formulário]]




{{url>https://forms.gle/dUqZad77fU8Hrqrm8}}

</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>

*/

===== GLM binomial: resposta em proporções =====

<WRAP center round box 80%>

{{ youtube>S3EWkMuV5ng |}}

</WRAP>

==== 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. 

<WRAP center round todo 80%>

    * baixe o arquivo {{ :cursos:planeco:planeco:roteiro:flowering.txt |}}
    * abra os dados no Rcmdr (a separação de campo é tabulaçã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>

</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. 

<WRAP center round todo 80%>

    * 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 80%>
** 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 de contagem, com a diferença que a família aqui é o ''quasibinomial'' 

</WRAP>

</WRAP>

==== Interpretação do resultado ====

Para interpretar os coeficientes use o mesmo procedimento do exercício anterior, que é aplicar a função exponencial (''exp'') nos coeficientes previstos e interpretar como ''chance'' e ''razão de chance''((O Rcmdr apresenta os valores dos coeficientes exponenciados após o resumo do modelo na sua construção )).

Para interpretar 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 todo 80%>

  * calcule o predito pelo modelo e os coeficientes na escala original
  * interprete o efeito da concentração na floração das variedades a partir dos coeficientes do modelo selecionado

</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>


Mas o **Rcmdr** não poderia ficar sem essa funcionalidade para interpretar os valores do predito pelo modelo na escala de observação, para isso utilize o menu **Models> add observation statistic to data...>** e selecione apenas o **Fitted values**. O Rcmdr adiciona uma coluna nos dados chamada ''fitted."nome_do_modelo"'', com os previstos na escala de observação, nesse caso probabilidade.




**__Gráfico para interpretação dos resultados__**

Para um gráfico dos resultados use o menu:

 **Models > Graphs > Predict effect plots...**


</WRAP>






/*




<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>

*/



==== O que preciso entregar ====


<WRAP center round help 100%>
 
 
  * Preencha as perguntas do quadro abaixo ou pelo [[https://forms.gle/81Y14mRvRY6cShfv5|link do formulário]]




{{url>https://forms.gle/81Y14mRvRY6cShfv5}}

</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>


*/
====== Dispersão e acumulo de zeros ======

{{section>cursos:planeco:roteiro:10-glm#Dispersã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 independentes.  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 a 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 dispersões moderadas e não indicam qual a fonte dela.
Há algumas alternativas ao modelo ''quasi'' para a 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 center round box 90%>

<wrap hi>** Variável resposta binária é um caso especial da binomial com apenas uma tentativa, chamado de distribuição de Bernoulli, e não tem problema com sobre-dispersão** 
</wrap>
</WRAP>
 

</WRAP>


====== link para páginas GLMs ======

  * [[cursos:planeco:roteiro:10-glmBinomial|]]
  * [[cursos:planeco:roteiro:10-glmPoisson|]]
Laboratório de Ecologia de Florestas Tropicais

Ferramentas do usuário

Ferramentas do site

Ferramentas da página
