====== Modelos Lineares Múltiplos III ======

===== Interação entre preditoras e Colinearidade =====

<WRAP center round box 80%>
{{ youtube>62cRvHQzIFc |}}

</WRAP>
Neste terceiro e último roteiro sobre Modelos Lineares Múltiplos vamos trabalhar com conjuntos de dados com **variáveis preditoras contínuas**. Inicialmente, iremos avaliar uma importante premissa dos modelos associados a esse tipo de dado, aprendendo a identificar colinearidade entre variáveis, entendendo os efeitos sobre a seleção de modelos e interpretando os coeficientes do modelo selecionado. Depois, vamos utilizar um conjunto de dados mais complexo com **variáveis contínuas e categóricas** para reforçar o procedimento de seleção de modelos e posteriormente obter e interpretar, por meio dos resumos dos modelos, os coeficientes dos parâmetros, incluindo as interações. Nessa segunda parte vamos exercitar todos os conceitos importantes de Modelos Lineares Múltiplos e fechar esse módulo da disciplina.

\\


===== Colinearidade entre variáveis =====


{{:cursos:planeco:roteiro:Figura matriz correlacao.jpg?nolink&200 |}}

Uma importante premissa de modelos lineares múltiplos é que as variáveis preditoras sejam independentes entre si. Entretanto, em estudos observacionais ou exploratórios é relativamente comum que as variáveis preditoras não sejam independentes. Quando duas variáveis preditoras estão correlacionadas e estão explicando a mesma porção da variância da variável resposta estamos diante de um problema de colinearidade. Nos casos mais extremos, a colinearidade pode afetar a significância de algumas variáveis e até mesmo o sinal do efeito.

Existem várias formas de lidar com a colinearidade, mas vamos focar nossa atividade em identificar e remover variáveis que estejam inflando as estimativas de variação. Para isso, vamos usar um índice chamado de Variance Inflation Factor (VIF), que é calculado a partir dessa equação:

$$ VIF_{i} = \frac{1}{1 - R_{i}^{2}} $$


/*

{{:cursos:planeco:roteiro:VIF_equation.png?nolink&200|}}

VIF<sub>X</sub><sub>1</sub> = 1 / ( 1 - R<sup>2</sup><sub>X</sub><sub>1</sub>)

Essa equação estava dando problemas em alguns browsers, por isso troquei por uma figura

*/

Esse R<sup>2</sup> é obtido ajustando um modelo linear múltiplo para analisar a relação entre cada variável preditora, por exemplo //i// = X<sub>1</sub>, e todas as outras preditoras no modelo de interesse (X<sub>2</sub>, X<sub>3</sub>,..., X<sub>n</sub>). Fazendo isso para todas as preditoras teremos um ''VIF'' para cada uma delas. Um alto R<sup>2</sup> significa que grande parte da variação na preditora em questão é compartilhada pelas outras variáveis. Veja no exemplo abaixo um passo-a-passo para calcular o ''VIF''.

Quanto maior for o valor do ''VIF'', mais os valores de erro padrão dos parâmetros do modelo serão inflados e mais dificilmente um efeito será detectado. Além da imprecisão nas estimativas dos parâmetros colineares, um outro problema que pode emergir é o modelo mínimo adequado ser diferente, dependendo da ordem da simplificação do modelo cheio. Um valor frequentemente usado para definir um limite aceitável de ''VIF'' é 4,0, acima desse valor as estimativas do modelo podem ficar comprometidas.

Nessa abordagem, após identificar quais são as variáveis com maiores valores de VIF, elas serão removidas sequencialmente. A cada variável retirada verifica-se novamente se os valores de VIF diminuíram ou se ainda precisam ser retiradas outras variáveis colineares. 


<WRAP center round important 90%>

É importante entender que a escolha de qual variável retirar vai depender também do sentido biológico/ecológico de cada variável. Em alguns casos, pode valer a pena manter uma variável cujo VIF é levemente mais alto, pois o mecanismo de explicação pode ser mais explícito para essa variável.
</WRAP>

Vamos ver como isso funciona na prática abaixo.
\\


====Biomassa de manguezais e variáveis ambientais====

{{:cursos:planeco:roteiro:MaguezalMarajo2.jpg?nonlink&300|http://www.imagensdobrasil.art.br/produtos/3148/3/15/Mangue_de_Maraj%C3%B3#.YlmY-1zMKV4}}
\\

O objetivo dessa pesquisa foi avaliar quais variáveis ambientais predizem melhor a biomassa acima do solo (Aboveground Biomass-AGB) de manguezais em diferentes locais do mundo((dados fictícios, mas baseados em valores reais)). Foram utilizadas 3 variáveis ambientais que são facilmente obtidas em bases de dados mundiais.
\\

1. Baixe o conjunto de dados {{ :cursos:planeco:roteiro:mangrove.csv |mangrove.csv}}, importe para o Rcmdr, usando __vírgula__ como separador de campo, e visualize os dados para entender o arquivo.

2. Entenda as variáveis do arquivo:  
  * variável resposta: 
             * **AGB_carbon** = Estoque de Carbono estimado a partir da biomassa acima do solo (AGB) de árvores de manguezais (MgC/ha) no hemisfério sul
  * variáveis preditoras:
             * **lat** = latitude (em graus)((Os dados são predominantemente do hemisfério sul e por isso seriam esperados valores negativos de latitude, porém foram transformados em valores positivos para facilitar a interpretação. Alguns valores do hemisfério norte aparecem como negativos na planilha)) 
             * **temp** = Temperatura média anual (em graus Celsius)
             * **ppt** = Precipitação anual (mm)


3. Inspecione a correlação entre todas as variáveis preditoras contínuas:

Para fazer isso no //Rcmdr//, você tem duas opções:

  * Uma opção é avaliar numericamente as correlações. Para isso, entre em **Statistics-> Summary -> Correlation Matrix**, selecione todas as variáveis preditoras e clique em OK. Você verá os valores de correlação de todos os pares de variáveis. Observando os valores mais altos de correlação, você já pode ter uma ideia se existem variáveis com potencial para apresentar colinearidade. 

  * Outra opção é avaliar graficamente as correlações entre as variáveis. Para isso, entre em **Graphics -> Dispersion Matrix** e selecione todas as variáveis preditoras contínuas. Na aba **Options** selecione "Minimum Square Line" e clique em OK. Na figura que foi gerada, você poderá avaliar quais pares de variáveis parecem ter uma maior correlação entre elas.

<WRAP center round important 90%>
Esse procedimento de analisar a correlação entre todas as nossas variáveis preditoras contínuas deveria ser sempre realizado antes de fazermos nossas análises. 
</WRAP>


4. Ajuste um modelo, relacionando AGB-carbon com todas as variáveis preditoras, mas ainda sem incluir as interações. Nomeie esse modelo como "carbon1". No //summary// do modelo, repare nos efeitos e na significância de cada um dos parâmetros.


5. Calcule os VIFs((A função utilizada no //Rcmdr// para calcular o ''VIF'' utiliza uma variante chamada ''GVIF'', uma generalização que pode ser aplicada também para variáveis categóricas com a mesma interpretação colocada acima o ''GVIF'' foi desenvolvida pelo John Fox, mesmo autor do Rcmdr. Veja o artigo no link [[https://www.tandfonline.com/doi/abs/10.1080/01621459.1992.10475190]])) para as variáveis incluídas no modelo

Para isso, entre em **Models -> Numerical diagnostics -> Variance-inflation factors**. O primeiro resultado apresentado é uma linha com os valores de VIF para cada parâmetro do modelo. O segundo resultado apresentado é uma matriz de correlação das __estimativas dos parâmetros__. Note que os valores são diferentes das correlações feitas diretamente para as variáveis (item 3, acima).

<WRAP center round important 90%>
Importante: Como o valor de VIF de cada parâmetro depende de quais outros parâmetros estão sendo incluídos no modelo, só é possível calcular os VIFs depois de ter ajustado um modelo. Ao usar o //Rcmdr//, fique sempre atento(a) se o modelo ativo é realmente o modelo para o qual você quer calcular os VIFs.
</WRAP>


<WRAP center round box 70%>

**PAUSA OPCIONAL** caso você queira aprender a calcular manualmente os valores de VIF:

Em primeiro lugar, reveja a equação de cálculo de VIF apresentada acima.

Agora, vamos calcular manualmente o valor de VIF para a variável preditora **lat** e comparar com o valor obtido acima no //Rcmdr//. Para isso, precisamos calcular o R<sup>2</sup> da relação entre essa variável preditora e todas as outras preditoras que estavam no modelo completo, sem as interações (carbon1). Para isso, vamos criar um novo modelo no qual a variável para a qual estamos interessados em calcular o VIF (**lat**) passará agora a ser a __variável resposta__ desse novo modelo que criaremos. 

Entre em **Statistics -> Fit models -> Linear model**. Coloque **lat** como __variável resposta__ na caixa da esquerda da equação e coloque as outras 2 variáveis preditoras (**temp** + **ppt**) na caixa da direita da equação. Defina o nome desse modelo como "viflat". No //summary// do modelo será apresentado o valor de R<sup>2</sup> Múltiplo (//Multiple R-square//). Utilize esse valor na equação de cálculo de VIF e veja se o resultado é igual ao valor de VIF calculado pelo R Commander para a variável **lat** a partir do modelo "carbon1" feito acima. Deveria ser. Se não foi, peça ajuda a alguém da equipe.

Repita o mesmo procedimento para outra variável de sua escolha. Você pode fazer isso para todas as variáveis do modelo, se quiser.

</WRAP>

Continuando nossa análise sobre o estoque de Carbono em manguezais:

6. Após analisar os valores dos VIFs do modelo "carbon1", se houver alguma variável com valor maior que 4, remova a variável com o maior VIF e ajuste um novo modelo. Coloque "carbon2" como nome desse modelo. Olhe para o //summary// desse modelo e para as variáveis que permaneceram nele. Cheque os valores dos coeficientes e a significância de cada variável em relação ao modelo "carbon1". **Houve alguma alteração? Alguma variável deixou de ser significativa? Alguma variável passou a ser significativa? O sinal do efeito mudou?** 

7. Calcule os VIFs das variáveis do modelo "carbon2" usando o caminho **Models -> Numerical diagnostics -> Variance-inflation factors** e veja se ainda tem alguma variável com VIF maior que 4.

8. Repita os procedimentos anteriores até não haver nenhuma variável com VIF maior que 4.

9. E possível que algumas das variáveis remanescentes, mesmo que não sejam colineares entre si, não sejam relevantes para definir o estoque de carbono em manguezais. Então, para iniciar o procedimento de seleção de modelos, crie um modelo completo  que inclua as variáveis remanescentes **e suas interações** (nomeie como "carbon_int"). Analise o //summary//do modelo.

10. Realize o procedimento de seleção do modelo mínimo plausível pelo método de simplificação para o mínimo adequado, conforme explicado no item [[cursos:planeco:roteiro:09-lm02#simplificando_modelos|Simplificando modelos]] do roteiro I de Modelos Lineares Múltiplos]] 

11. Analise os resultados do modelo final.
\\
\\
\\



===== Modelos Lineares Múltiplos: preditoras contínuas e categóricas =====

{{:cursos:planeco:roteiro:plotBabies01.png?300  |}}

\\
\\
\\
\\
Nesse último tópico do bloco vamos resgatar os principais conceitos que emergiram com a generalização do modelo linear, agora com múltiplas preditoras, a partir de um exemplo que tem duas variáveis preditoras contínuas e duas categórica. Acreditamos que esse exemplo incorpora as complexidades tratadas e ajuda a agrupar os tópicos que devem ficar atentos nos modelos com múltiplas preditoras.     
\\
\\
\\ 
\\
\\
\\
\\

<WRAP center round box 60%>

==== Desafios dos modelos com múltiplas preditoras ====


Ao final desta seção é desejável que tenha compreendido nos modelos lineares múltiplos:

  * compreender a partição da variância do modelo;
  * interpretar a ''tabela de anova'' na comparação de dois modelos;
  * entender o procedimento da ''anova'' para simplificação do modelo;
  * saber interpretar os gráficos diagnósticos do modelo;
  * avaliar a colinearidade entre variáveis no modelo;
  * interpretar os coeficientes estimados;
  * entender quais níveis estão representados no intercepto do modelo;
  * compreender os termos de interação;
  * compor o predito pelo modelo a partir dos coeficientes;
  * interpretar biologicamente o resultado do modelo.


</WRAP>
\\
\\

==== VIF e as interações ====

No //Rcmdr// o  ''VIF'' é aplicado ao modelo ativo pelo menu ''Models > Numerical diagnostics > Variance-inflation factors''), calculando o valor para todos os termos do modelo, inclusive as interações. Como interações e as variáveis isoladas compartilham parte da variação explicada, a correlação entre eles é esperada.  Ou seja, não é possível fazer a avaliação do VIF das variáveis em modelos com interação diretamente. Uma soluçao é fazer modelos sem as interações como fizemos anteriormente. 
Uma outra forma de contornar esse problema é fazer uma transformação simples nas variáveis contínuas, centralizando a média em ''zero'', subtraindo o valor observado da média( $x_i - \bar{x}$ ). 

Com essa transformação o valor ''0'' passa a representar a média e os valores positivos o aumento em relação a média e negativos a diminuição, na mesma unidade de escala da variável original. 
A **centralização** das variáveis contínuas é uma transformação corriqueira pois não dificulta a interpretação e ao contrário, evita muitos problemas analíticos e de interpretação. Entre as vantagens da centralização está a possibilidade de interpretar o ''VIF'' diretamente no modelo selecionado e incorporar uma interpretação biológica para o valor do intercepto, onde muitas vezes não existia.  

\\
\\
\\

==== Peso de bebês ao nascer ====

{{:cursos:planeco:roteiro:cegonha final_133910457_dreamstime.jpg?nolink&250 |}}
O objetivo dessa pesquisa foi saber quais fatores afetam o tamanho de bebês ao nascer, de modo que fosse possível orientar campanhas de conscientização para evitar o nascimento de bebês com baixo peso, uma vez que isso pode implicar em maiores custos e muitos riscos ao bebê devido à permanência no hospital. Três variáveis preditoras (explicadas abaixo) foram consideradas relevantes para essa pesquisa, mas também havia um interesse genuíno em saber se alguma das variáveis poderia interferir no efeito das outras. Como a variável resposta, peso do bebê ao nascer, foi medida em ''onças'' vamos primeiro transformar em uma escala de medida que temos mais facilidade para  interpretar, multiplando essa variável por ''0.02835'' para transformar em ''kg''.



<WRAP center round box 80%>

  * Abra o arquivo {{ :planeco:roteiro:babies.csv |}} no Rcmdr, usando __tabulação(Tabs)__ como separador de campo
  * Garanta que os dados foram lidos corretamente>
  * Abra a janela para criar uma nova variável no menu ''Data > Manage variables in active data set > Compute a new variable'';
  * Na caixa ''New variable name'' nomei a nova variável como ''pesoKg'';
  * Na caixa ''Expression to compute'' coloque a expressão: ''bwt * 0.02835'';
  * Ajuste um modelo contendo __apenas as variáveis indicadas abaixo__ e todas as interações entre elas:
  * variável resposta:  ''pesoKg'' = peso do bebê (medido em ''kg'') 
  * preditoras:
            * ''gestation'' = tempo de gestação (dias)
            * ''age'' = idade da mãe
            * ''smoke'': ''FALSE'' mãe não fumante; ''TRUE'' mãe fumante

  *  Selecione o modelo mínimo plausível pelo método de simplificação para mínimo adequado (ver roteiro I de MLM)

  * Calcule o VIF do modelo selecionado pelo menu ''Models > Numerical diagnostics > Variantion Inflation Factor''
  * Guarde o resultado dos VIF destes modelos;
  * Crie uma nova variável pelo menu : ''Data > Manage variable in active data set > Computer new variable'';
  * Na janela que se abre coloque em ''New variabel name'' o nome ''ageCenter'' e em ''Expression to compute''   inclua a expressão ''age - mean(age)'';
<WRAP center round box 60%>
{{:cursos:planeco:roteiro:centerAgeBabies.png?400|}}
</WRAP>
  *  Faça o mesmo para uma nova variável com o nome ''gestationCenter'' usando a expressão '' gestation - mean(gestatation)'';
  * Construa o modelo selecionado utilizando estas novas variáveis contínuas centralizadas em substituição às originais;
  * Refaça o calculos dos VIFs para esse novo modelo com as variáveis selecionadas. Guarde o resultado.
  * Para o modelo final selecionado, com as variáveis preditoras contínuas centralizadas:
    * avalie os gráficos diagnósticos; 
    * faça a avaliação da colinearidade entre os termos do modelo;
    * identifique qual(is) nível(is) está(ão) representado(s) no intercepto;
    * interprete cada um dos parâmetros do modelo, incluindo interações, se houver;
  * A partir dos resultados do modelo proponha uma campanha para evitar que bebês nasçam com baixo peso.


</WRAP>
\\

<WRAP center round tip 60%>
Retorne à [[cursos:planeco:roteiro:09-lm02b#desafios_dos_modelos_com_multiplas_preditoras|lista de desafios dos modelos com múltiplas preditoras]] do início desta seção e avalie se todos os pontos foram compreendidos. 
</WRAP>




===== Exercício=====

<WRAP center round help 90%>

Responda o  [[https://forms.gle/qZg13LsMEG2xHRAe6|o formulário MLM III]] incluindo arquivos de resultados e figuras quando solicitado. 


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

</WRAP>



/*
- subir um arquivo com os resumos de alguns modelos, incluindo o modelo final selecionado e os gráficos diagnósticos referentes aos dados babies.csv, interpretar o modelo final e responder as perguntas propostas.

- subir um arquivo com os resultados referentes às análises de colinearidade (por meio dos VIFs) para os dados birds_clim.csv e responder as perguntas propostas

- subir um arquivo com a seleção de modelos dos dados de birds.csv, a partir do modelo completo com as variáveis que permaneceram após a remoção daquelas com altos VIFs. Porém, para esse exercício, caso tenham permanecido três ou mais variáveis, faça o modelo completo contendo as variáveis que permaneceram, mas **apenas as __interações duplas__**((ou seja, não precisa incluir as  interações triplas, quádruplas ou quíntuplas)) entre elas. Esse arquivo deve ter os resumos dos modelos e as comparações entre modelos até chegar ao modelo final selecionado. Para o modelo final, apresente também os gráficos diagnósticos. 

- interpretar o modelo final selecionado.




ATIVIDADES ANTIGAS QUE FORAM EXCLUÍDAS OU MODIFICADAS EM 2020:


<WRAP center round box 80%>
  * abra o arquivo {{ :planeco:roteiro:babies.csv |}} no R ou Rcmdr
  * faça um modelo plausível com as variáveis:
    * resposta **bwt** : peso do bebê ao nascer
    * preditoras:
        * gestation: tempo de gestação (dias)
        * age: idade
        * weight: peso
        * smoke: 0 não fumante; 1 fumante
        * interações: até a 3 ordem para aquelas que são plausíveis e passíveis de interpretação biológica
    * selecione o modelo mínimo plausível pelo método de simplificação para mínimo adequado 
    * interprete o resultado
</WRAP>


<WRAP center round box 80%>
  - abra o arquivo {{ :planeco:roteiro:davis.csv |}} no R ou no Rcmdr
  - monte o modelo para peso com relação às variáveis preditoras: altura, sexo e suas interações.
  - a partir  do modelo mais cheio, simplifique o modelo até o mínimo adequado
  - apresente o resultado em um gráfico

</WRAP>



===== Exercício Extra =====

Um artigo recente (Fisher, R. & Ai C. 2018)((note que temos dois homônimos aqui, o mais famoso estatístico e um critério de seleção de modelos juntos, como autores!! Os outros nomes foram omitidos propositadamente.)) sobre métodos de regressões múltiplas, apresenta dados sobre peixes de recifes de corais. Entre as questões apresentadas pelos autores originais do trabalho estava se a biomassa de diferentes guildas de peixes em zonas protegidas ou não de recifes. Aqui vamos usar apenas a guilda de peixes que se alimentam de plânctons e apenas parte das preditoras. O método desenvolvido no artigo é uma forma de automatizar a seleção de preditoras em modelos com muitas variáveis potenciais de influenciar a resposta.

  * faça uma análise exploratória das variáveis do dados {{ :cursos:planeco:roteiro:plankivore.csv |}} e suas relações;  
  * monte o modelo mínimo adequado, partindo de todas as variáveis preditoras e suas interações;
  * verifique se não há nenhuma variável com relação não linear com a biomassa, e se houver use mais uma variável representada pelo quadrado dessa variável;
  * ao final faça o diagnóstico do modelo e veja se o resíduo cumprem com as premissas do modelo linear;
  * caso diagnostique problema no modelo, transforme a variável resposta usando o logaritmo natural;
  * refaça a seleção e o diagnóstico com a resposta na escala log;
  * interprete o resultado.
 

{{ :cursos:planeco:roteiro:Fisher_et_al-2018-Ecology_and_Evolution.pdf | Fisher, R.; Ai C. L. et al. 2018. A simple function for full-­subsets multiple regression in ecology with R. Ecology and Evolution 8: 6104-6113}}



=== Gerando dados com colinearidade ===

Gerar dados com colinearidade.

*/