Modelos Lineares Múltiplos I [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 Múltiplos I ======

{{:cursos:planeco:roteiro:datamodel.jpg?300|}}


Uma extensão do modelo linear simples ((modelo linear com apenas uma preditora)) são os modelos lineares com mais de uma preditora, aqui definido como modelos múltiplos. Quando temos mais de uma preditora o modelo aumenta em complexidade com mais parâmetros para estimar. Além disso, a estrutura mais complexa do modelo gera desafios para a interpretação e dificulta a avaliação da adequação do modelo aos dados. Uma primeira complexidade está relacionada a como simplificar a estrutura do modelo com a finalidade de facilitar a interpretação e melhorar a estimação dos parâmetros. A tomada de decisão sobre quais variáveis devemos reter em nosso modelo e quais podem ser retiradas, por não terem efeito na variável resposta, pode ser feita utilizando diferentes critérios e técnicas. A seguir apresentamos uma das técnicas utilizadas para essa tomada de decisão e que iremos utilizar ao longo desse curso. Outros critérios ou técnicas podem ser utilizadas com vantagens ou desvantagens em relação ao que utilizaremos. Não é objetivo desse curso se debruçar sobre essas diferentes técnicas. 


===== Duas preditoras categóricas =====

O primeiro exemplo que iremos trabalhar é baseado nos dados utilizados para exemplificar o [[cursos:planeco:roteiro:07b-anovarcmdr#particao_da_variancia|teste de Anova]]. Vamos criar um experimento plausível a partir dele.


======= Simulando um experimento plausível =======

Vimos que existe um efeito do tipo de solo na produção de um cultivar. Uma expectativa plausível é que a adição de adubo também tenha efeito na produtividade. Ou seja, os tipos de solo tem produtividade diferente, assim como o adubo aumenta a produtividade.

Nos dados originais do exercício de ANOVA a produtividade média nos solos foi de:

  * arenoso:   9.9
  * argiloso: 11.5
  * humico:   14.3

  
Vamos, a partir dessa informação, criar um experimento onde, além da diferença do solo, metade dos cultivos foram tratados com adubo orgânico.

  * 1. Abra o arquivo {{ :cursos:planeco:roteiro:cropMult.xlsx |cropMulti}} em uma planilha eletrônica:

<WRAP center round box 80%>

{{ :cursos:planeco:roteiro:cropMult01.jpg?500 |}}
</WRAP>
 
  *2. Preencha a coluna ''efeitoAdubo'' com o valor de  ''1.2'' para todas as parcelas  adubadas ((coluna ''adubo'' igual a ''sim'')) e ''0'' para aquelas que não foram ((coluna ''adubo'' igual a  ''não'' )).
  *3. Preencha a célula ** E2** da coluna ''desvios normal'' com a fórmula //** = INV.NORM.N(ALEATÓRIO(); 0 ; 1.5)**//((Essa expressão retorna valores associados a uma distribuição normal com média 0 e desvio padrão 1.5. Para libreoffice use '' = NORM.INV(RAND(), 0, 1.5)'')).
  *4. Some os valores em uma mesma linha 

Ao final sua planilha deve estar preenchida como a que segue, apenas com os valores da coluna resíduo diferentes:

<WRAP center round box 80%>

{{ :cursos:planeco:roteiro:adubo02.png?500 |}}

</WRAP>
 
<WRAP center round todo 80%>
__**Procedimentos**__

  *  Salve a planilha com o nome ''soloAduboAditivo.csv'' em formato texto com campos separados por vírgula;
  *  Abra os dados no Rcmdr;
  *  Produza um modelo chamado ''mlSolo_Adubo_Aditivo'' da seguinte forma:
    ''prodCampo ~ solo + adubo'';
  * Avalie o modelo pelo seu sumário e pela tabela de Anova;
  * Faça uma interpretação biológica do resultado do modelo.

</WRAP>


=====  Modelos Plausíves =====



O nosso modelo tem duas preditoras e pode ser simplificado. Nesse caso, como temos poucas possibilidades de comparação, podemos comparar os modelos plausíveis, desde que sejam aninhados. O que produzimos acima tem o efeito de solo e de adubo, podemos pensar em mais algumas possibilidades de modelo:

  * **__mlSolo__** só com o efeito do solo:

''prodCampo ~ solo''

  * **__mlAdubo__** só com efeito do adubo:
 
''prodCampo ~ adubo''

  * **__mlNull__** sem efeito de solo ou adubo:

''prodCampo ~ 1''

O valor ''1'' na última formula indica que o modelo não tem nenhuma variável preditora ((o valor 1  indica que a resposta é predita apenas pela sua própria média))

/*
Antes de construir modelos alternativos e compará-los, vamos complicar um pouco os nosso dados simulados incluindo algo que é a interação entre as variáveis preditoras. 

*/ 


==== Interação entre preditoras ====

Nos modelos acima, desconsideramos um elemento importante que emerge quanto temos mais de uma preditora, a possibilidade de uma variável preditora interferir no efeito de outra, efeito esse chamado de interação.
A interação é um elemento muito importante quando temos mais de uma preditora, pois desconsiderá-la pode limitar o entendimento dos processos envolvidos. Um exemplo cotidiano da interação é visto no uso de medicamentos e o alerta da bula sobre interação medicamentosa ou efeitos colaterais para pessoas portadoras de doenças crônicas. Dizemos que um medicamento tem interação com outra substância quando o seu efeito é modificado pela presença de outra substância, como por exemplo a ingestão de álcool junto com muitos medicamentos. Nos modelos a interação tem uma interpretação similar, a resposta pelo efeito de uma variável preditora se altera com a presença de outra preditora. Muitas vezes a interação pode ser o efeito de interesse do estudo, como na pergunta: //O efeito de solo na produtividade agrícola  depende da quantidade de adubo orgânico adicionado?// Ou em outras palavras: //O efeito da adubação orgânica depende do tipo de solo?//  Note que nestas perguntas o foco não é se há ou não efeito do adubo ou solo, mas se a presença de uma variável afeta o efeito de outra.  


<WRAP center round todo 80%>

  * No conjunto de modelos acima, não incluímos o termo da interação. Produza o modelo abaixo incluindo o termo da interação e avalie esse modelo e seus coeficientes.
 
    ''prodCampo ~ solo + adubo + solo:adubo''
 

</WRAP>


Não é esperado encontrar interação entre as preditoras nos dados simulados da maneira como fizemos, ele pode emergir por acaso, apenas porque temos uma variável aleatória ((se o termo da interação foi significativo, confira os cálculos e mantenha o resultado como está, esse resultado emerge com baixa frequência, simplesmente por acaso. )). Da maneira como simulamos os dados temos duas preditoras que tem efeitos  aditivos onde não há interação. Uma outra forma de dizer isso é que o efeito do ''adubo'' não interfere no efeito do ''solo'', ou que esses efeitos são independentes. A interpretação biológica nesse caso também pode ser feita independentemente. 



======= Simulando dados com interação =======

Seguindo a mesma abordagem anterior, vamos produzir dados simulando a interação entre as variáveis ''solo'' e ''adubo''. Para isso precisamos produzir dados em que o efeito do adubo depende do tipo de solo.


  - Abra o arquivo {{ :cursos:planeco:roteiro:cropMult.xlsx |cropMulti}} em uma planilha eletrônica:

<WRAP center round box 80%>

{{ :cursos:planeco:roteiro:cropMult01.jpg?500 |  }}
</WRAP>
 
  - Preencha a coluna ''efeitoAdubo'' com os valores:
    * 2.7 para ''arenoso'' com ''adubo'' igual a ''sim'' 
    * 0.7 para ''argiloso'' com ''adubo'' igual a ''sim'' 
    * 0.2 para ''humico'' com ''adubo'' igual a ''sim'' 
  - O campos da coluna ''efeitoAdubo'' onde ''adubo'' é igual a ''não'' devem ser preenchidos com ''0'' 
  - Preencha a célula ** E2** da coluna ''desvios normal'' com a fórmula //** = INV.NORM.N(ALEATÓRIO(); 0 ; 1.5)**//((Essa expressão retorna valores associados a uma distribuição normal com média 0 e desvio padrão 1.5. Para versões antigas do libreoffice a função pode ser '' = NORM.INV(RAND(), 0, 1.5)'')), as atuais utilizam a mesma que o excel.
  - Some na coluna ''prodCampo'' os valores ''prodSolo + efeitoAdubo + desviosNormal'' 

Ao final sua planilha deve estar preenchida como a que segue, apenas com os valores da coluna resíduo diferentes:

<WRAP center round box 80%>

{{ :cursos:planeco:roteiro:cropMult02.jpg?500 |}}
</WRAP>
 
<WRAP center round tip 80%>
__**Procedimentos**__

  - Salve a planilha com o nome ''soloAduboInteracao.csv'';
  - Importe os dados para o Rcmdr. **Atenção nomeie os dados na aba de importação com o nome ''soloAduboInt'', em alguns casos o Rcmdr não importa se a planilha e os dados importados tiverem o mesmo nome de uma importação anterior**
  - Confira se os dados foram lidos corretamente, inclusive se a decimal é ''.'';
  - Produza o modelo cheio ''mlSolo_AduboAll'' com a seguinte formula:
     * ''prodCampo ~ solo + adubo + solo:adubo''
     * interprete o resumo, comparando com o resumo do modelo similar proveniente da planilha de dados anterior

</WRAP>


===== Simplificando Modelos =====

Durante o curso usaremos o procedimento de simplificar o modelo a partir do modelo cheio.
O procedimento consiste em comparar modelos aninhados((o modelo mais simples está contido no mais complexo)), dois a dois, retendo o que está mais acoplado aos dados. Para comparar os modelos utilizaremos o procedimento da partição da variância baseado na tabela de anova. Quando os modelo comparados são diferentes retemos o mais complexo, pois explica mais variação dos dados ((Este é um atributo associado aos modelos aninhados: aquele que tem mais variáveis ou parâmetros só pode explicar mais ou a mesma quantidade de variação do mais simples, já que todos os parâmetros do modelo mais simples estão contidos no mais complexo)). Por outro lado, quando os modelos não são diferentes no seu poder explicativo, retemos o modelo mais simples, apoiados no princípio da parcimônia. Para tomar a decisão se os modelos são iguais ou diferentes utilizamos a estatística ''F'' da tabela de anova. 


<WRAP center round box 80%>

==== Princípio da parcimônia (Navalha de Occam) ====

  * número de parâmetros menor possível
  * linear é melhor que não-linear
  * reter menos pressupostos
  * simplificar ao mínimo adequado
  * explicações mais simples são preferíveis

</WRAP>

<WRAP center round box 80%>

==== Método do modelo cheio ao mínimo adequado ====

  - ajuste o modelo máximo (cheio)
  - simplifique o modelo:
    * inspecione os coeficientes (summary)
    *  remova termos não significativos ((um de cada vez))
  - ordem de remoção de termos: 
    * interações não significativas (primeiro as de maior ordem)
    * termos quadráticos ou não lineares
    * variáveis explicativas não significativas
  - caso faça sentido, agrupe níveis de fatores sem diferença
  - verifique se a ordem de remoção não interfere na seleção do modelo
    * retorne ao modelo cheio
    * retire as variáveis que não foram retidas no outro procedimento em outra ordem
    * confirme que o modelo mínimo adequado é o mesmo 
  - Faça o diagnóstico do modelo mínimo adequado
  - Interprete o modelo selecionado 
    

</WRAP>


<WRAP center round box 80%>

==== Tomada de decisão  ====
<WRAP center round tip 80%>
** A diferença não é significativa: **
    * retenha o modelo mais simples
    * continue simplificando

</WRAP>

<WRAP center round tip 80%>

**A difereça é significativa: **
    * retenha o modelo complexo
    * verifique se existe termo que pode e ainda não foi retirado 
    * caso não haja nenhum termo que possa ser retirado, este é o modelo __MINÍMO ADEQUADO__

</WRAP>

</WRAP>








<WRAP center round box 80%>

==== Interpretando Variáveis Indicadoras (Dummy) ====


As variáveis indicadoras devem ser interpretadas com cuidado. No exemplo do modelo cheio acima ((aquele que inclui a interação entre ''solo'' e ''adubo'')), o modelo pode ser descrito da seguinte forma:


$$ y_{tr} =  \alpha + \beta_1 * arg + \beta_2 * hum + \beta_3 * adubo + \beta_4 * arg * adubo + \beta_5 * hum * adubo $$


As variáveis __arg__, __hum__ e __adubo__ são dummy ou indicadoras, representadas por 1 quando presente e 0 quando ausentes. $\alpha, \beta_i$ representam as estimativas do modelo e estão relacionados, nesse caso, ao efeito de cada tratamento.


Para calcular o valor predito para o tratamento no solo arenoso com adubo, temos:

$$ y_{arenAdubo} =  \alpha + \beta_3 * adubo $$

Isso em decorrência do tratamento **arenoso sem adubo** estar representado pelo intercepto ($\alpha$) do modelo.


Para o tratamento de solo **argiloso com adubo** o predito é:

$$ y_{argAdubo} =  \alpha + \beta_1 * arg + \beta_3 * adubo + \beta_4 * arg * adubo $$

E assim por diante, usando as variáveis indicadoras e os coeficientes estimados para o cálculo do predito pelo modelo.

 
</WRAP>




<WRAP center round todo 80%>
**__Procedimento__**

  - Faça a seleção do modelo mínimo adequado para o conjunto de dados da última planilha, partindo do modelo com a interação, simplificando até o modelo mínimo adequado. Utilize o procedimentos de comparação de modelo pela partição de variância;
  - Avalie o modelo selecionado pelo sumário e pela tabela de Anova. Reconheça os valores utilizados para gerar os dados a partir das estimativas do modelo.
  - Preencha a aba  ''cropIntera'' da planilha [[https://docs.google.com/spreadsheets/d/1ISjeCppXA4kP-l_6uHPLst7GecvXRak6vryS9UqY0Rw/edit?usp=sharing|lmCrop2pred]] com os resultados do modelo selecionado
  - Na planilha onde os dados foram gerados, calcule, a partir dos coeficientes estimados, os valores preditos pelo modelo para cada uma das observações, coloque esses valores em uma coluna nomeada de ''predito''. Veja como calcular os valores preditos no quadro [[cursos:planeco:roteiro:09-lm02#interpretando_variáveis_indicadoras_(dummy)|]]
  - Calcule os resíduos do modelo ((diferença entre observado e o predito pelo modelo)) em uma coluna denominada ''residuos''
  - Eleve o valor dos resíduos ao quadrado em uma coluna denominada ''resQuad''. A  soma destes valores representa a variabilidade não explicada pelo modelo
  - Calcule a média da variável resposta e calcule a diferença deste valor para todas as observações e eleve ao quadrado e armazene em uma coluna ''desvQuadTotal'', a soma destes valores representa a variabilidade total dos dados
  - Calcule o R² do modelo, baseado no ''resQuad'' e no ''desvQuadTotal'' ((O R² é a razão entre  (''desQuadTotal'' - ''resQuad'')  sobre a ''desvQuadTotal''. Ou seja, quanto da variação dos dados é explicada pelo modelo em relação ao total de variação dos dados))

</WRAP>

===== Modelos Lineares Multiplos: ANCOVA  =====

Entres os testes frequentistas clássicos existe um que se presta para avaliar o efeito de uma preditora categórica e uma contínua influenciando uma variável respostas contínua, denominado de ''Análise de Covariância''. Na unificação existente nos modelos lineares, não precisamos nos preocupar com essa nomenclatura antiga, podemos simplesmente incorporar essa variáveis indistintamente no modelo e ele automaticamente irá transformar os níveis da variável categóricas em variáveis indicadoras. Entretanto, precisamos entender essa transformação para poder interpretar o modelo resultante e ficar atentos, especialmente, quando misturamos variáveis preditoras contínuas com categóricas.

Vamos agora retomar os procedimentos que vimos no inicio desse roteiro e tratar de um exemplo onde as variáveis preditoras apresentam essa complexidade.  

==== Exemplo: peso, altura e sexo ====

Vamos utilizar nesse exemplo dados parciais de uma publicação (Davis 1990((Weight and diet preoccupation and addictiveness: The role of exercise.Personality and Individual Differences 11(8): 823-827))) que contém informação sobre peso, altura e sexo de indivíduos.  

<WRAP center round todo 80%>
  * baixe o arquivo de dados {{ :cursos:planeco:roteiro:altPeso.txt |}}; 
  * importe os dados para o {{:cursos:planeco:planeco:logorcmdr01.png?20|}} com as opções:
      * nome das variáveis no arquivo
      * separador de campo tabulação (''Tabs'')
      * separador de decimais como ponto
{{:cursos:planeco:roteiro:importAltPeso.png?400|}}

  * explore os dados fazendo um gráfico de dispersão pelo menu ''Gráficos'' > ''Diagrama de Dispersão'', selecionando as variáveis ''weight'' , ''height'' e na opção ''Gráfico por grupos...'' selecione a variável ''sex'' 
  * salve o gráfico para postar no formulário

</WRAP>


====  Modelos Plausíves ====

<WRAP center round box 80%>
{{ youtube>89hcuYp9oWo |}}
</WRAP>

Nesse exemplo estamos interessados em entender a influência de altura e sexo no peso dos indivíduos. Neste caso temos alguns modelos plausíveis concorrentes que iremos nomear como:


  * **__lmNull__** sem efeitos de altura ou sexo:

''weight ~ 1''

  * **__lmSex__** só com efeito do sexo:
 
''weight ~ sex''


  * **__lmAlt__** só com o efeito do altura :

''weight ~ height''


  * **__lmAltSex__** efeito aditivo de altura e sexo:

''weight ~ height + sex''

  * **__lmAltXSex__** efeito de interação entre altura e sexo:

''weight ~ height * sex''


O valor ''1'' na última fórmula indica que o modelo não tem nenhuma variável preditora ((o valor 1  indica que a resposta é predita apenas pela sua própria média)).  




<WRAP center round todo 80%>
{{:cursos:planeco:roteiro:ancova01.png?400|}}

  * construa cada um desses modelos no {{:cursos:planeco:planeco:logorcmdr01.png?20|}} pelo menu ''Estatísticas'' > ''Ajuste de Modelo'' > ''Modelo Linear...'';
  * não esqueça de nomear os modelos de forma diferente para que não sejam sobrescritos;
  * faça o ''summary'' de cada um deles e salve o resultado em um arquivo;
  * descreva em uma frase a conclusão biológica associada a cada um desses modelos;
  * indique qual o valor que indica a precisão da conclusão para cada um dos modelos.
</WRAP>


==== Modelo lmAltXSex ====

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

</WRAP>

Vamos agora explorar o nosso modelo mais complexo, aquele que além da ''height'' e ''sex'' tem também a interação entre essas duas variáveis ''height:sex''. 


<WRAP center round todo 60%>
Para o modelo com a interação ''lmAltXsex'':
  * utilize os coeficientes do modelo para:
      * estimar o peso de homens com alturas de: 1.5; 1.7 e 1.9 metros;
      * estimar o peso de mulheres com alturas de 1.5; 1.7 e 1.9 metros;
      * descreva a equação, com os valores dos coeficientes que usou para o cálculo;
  * anote os valores para postar no formulário final.
</WRAP>


==== Seleção do Modelo ====

Agora, vamos usar o procedimento de [[cursos:planeco:roteiro:09-lm02#simplificando_modelos| simplificação do modelo]] que usamos no exemplo anterior para chegar ao modelo mínimo adequado. 

<WRAP center round box 80%>
{{ youtube>5av4ffv89A0 |}}

</WRAP>

<WRAP center round todo 60%>

  * inicie com o modelo cheio e compare, dois a dois, modelos aninhados;
  * a cada etapa de comparação, salve a ''tabela de Anova'';
  * a cada etapa, descreva a decisão tomada;
  * descreva o resultado do modelo selecionado.

</WRAP>


/*
===== 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 com mais variáveis preditoras contínuas e  categóricas.


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


Além da interação entre as preditoras, tratada no exemplo anterior, a inclusão de múltiplas preditoras eleva a complexidade do modelo e dificulta sua interpretação. Por exemplo, um modelo com três variáveis preditoras pode potencialmente ter quatro termos de interação((três interações duplas e uma tripla)), além dos termos dos efeitos isolados. Com quatro preditoras já são 11 potenciais termos de interações possíveis. Além do número elevado de termos possíveis, as interações de níveis mais elevados são difíceis de interpretar. Essa é uma das razões para buscarmos a simplificação dos modelos antes de interpretá-los. Um outro ponto importante é nunca incluir termos no modelo que não temos capacidade para interpretar, a menos que seu modelo seja apenas preditivo. Para modelos que se propõem apenas a fazer predições, um termo complexo não impede sua aplicação. Entretanto, se a proposta do modelo é entender os efeitos causais e os processos subjacentes, um termo complexo, que não pode ser interpretado, inviabiliza o modelo. Nesses casos, é importante apenas incluir termos no modelo que são passíveis de serem interpretados. 

Uma outra questão que emerge dos modelos com múltiplas preditoras contínuas é a multicolinearidade. Quando duas variável preditoras apresentam muita sobreposição de variação explicada, ou seja, ambas explicam a mesma porção considerável da variação na resposta, o modelo pode ter problemas na estimativa dos coeficientes. Além disso, modelos diferentes podem ser selecionados na busca do mínimo adequado, dependendo da ordem em que os termos são retirados ou incluídos. Variáveis com muito correlação entre elas tem grande potencial de apresentar multicolinearidade no modelo. Uma boa abordagem para quem está iniciando na modelagem estatística é investigar a correlação entre as preditoras e, caso a correlação seja muito forte, ficar atento com as estimativas do coeficientes e com a ordem de retirada dos termos do modelo ao buscar o modelo mínimo adequado. Se possível, descartar uma das preditoras colineares é uma das estratégias, já que ambas estão explicando a mesma porção da variação e são em grande parte redundantes.  

Caso tenha interesse sobre outra técnica de diagnóstico de colinearidade, consulte o roteiro [[cursos:planeco:roteiro:09-lm02b|]].  




Quando usamos variáveis preditoras contínuas nos modelos múltiplos devemos avaliar se as variáveis contínuas apresentam muita sobreposição de variação explicada. O compartilhamento de variação explicada entre preditoras torna as estimativas dos coeficientes dos modelos problemática. 



<WRAP center round box 60%>



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

  * 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;
  * 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;
  * fazer o diagnóstico 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 muitos riscos ao bebê  maiores custos devido à permanência no hospital. As variáveis preditoras consideradas para essa pesquisa estão listadas abaixo, mas também havia um interesse genuíno em saber se o efeito de uma variável poderia interferir no efeito das outras. 



A descrição destes dados pode ser consultada em [[https://www.stat.berkeley.edu/users/statlabs/labs.html#babies]]. 

<WRAP center round box 80%>
=== Descrição dos dados ===
    * variável resposta 
        * ''bwt'' : peso do bebê ao nascer em onças(oz)
    * preditoras:
        * ''gestation'': tempo de gestação (dias)
        * ''age'': idade da mãe (anos)
        * ''height'': altura da mãe (polegadas)
        * ''weight'': peso da mãe (libras)
        * ''smoke'': 0 não fumante; 1 fumante
   
</WRAP>

Notem que as preditoras estão relacionadas a características da mãe: dias de gestação, idade, peso, altura e se ela é fumante ou não. Como a variável resposta, peso do bebê ao nascer, foi medida em ''onças'', vamos transformar em uma escala de medida que temos mais facilidade para  interpretar, multiplicando essa variável por ''0.02835'' para transformar em ''kg''. Vamos fazer o mesmo para o peso da mãe que foi medido em libras e deve ser multiplicada por ''0.453592'' para ser transformada em ''kg''.

<WRAP center round box 80%>

  * Abra o arquivo {{ :planeco:roteiro:babies.csv |}} no Rcmdr((os campos neste arquivo são separados por tabulação));
  * 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 ''bebeKg'';
  * Na caixa ''Expression to compute'' coloque a expressão: ''bwt * 0.02835'';
  * Em seguida crie a variável ''maeKg'' multiplicando ''weight'' por ''0.453592''.
</WRAP>



Agora que temos uma variável que caracteriza o tamanho do bebê ao nascer em uma escala que faz mais sentido (kg), vamos começar a nossa avaliação de quais características da mãe tornam os bebês mais susceptíveis a nascerem com pesos baixo.

<WRAP center round important 90%>
Para simplificar nosso exemplo, vamos deixar de lado duas variáveis preditoras que foram coletadas nesse estudo: ''parity'' e ''height''. Ao final ficamos com quatro variáveis preditoras:
 ''gestation'', ''age'', ''maeKg'' e ''smoke''      

</WRAP>



===== Um mal começo =====

Um procedimento de modelagem ineficiente é colocar todas as variáveis preditoras e suas interações e torcer para ter alguma resultado interessante. No caso do peso dos bebês, que é razoavelmente simples com quatro variáveis preditoras sendo uma delas categórica, o modelo resultante é bastantes complexo. Vamos construir esse modelo e verificar o resultado:


<WRAP center round box 80%>
No menu ''Estatísticas > Ajustes de modelos > Modelo linear...'', construa o modelo ''lmFull'' com todas as quatro preditoras ((''gestation'', ''age'', ''maeKg'' e ''smoke'')) e suas interações. O modelo resultante tem a seguinte expressão:

<code> bebeKg ~ age * gestation * maeKg * smoke</code>

Na linguagem R os símbolos de ''*'' em expressões de notação de modelos representa tanto a variável isoladamente quanto as interações possíveis. 
Verifique o resumo do modelo gerado. 

</WRAP>


Um modelo como o construído acima tem muitos problemas. Nesse caso, foram estimados 16 coeficientes um para cada termo do modelo. Note também que nenhum dos coeficientes estimados é significativo nos testes marginais, ou seja, as incertezas nas estimativas dos coeficientes não nos permite afirmar que nenhum deles seja diferente de zero. Portanto, temos um modelo estatístico complexo, com 16 termos, sendo que todos os termos são multiplicados por coeficientes não distintos de zero! 


<WRAP center round box 80%>

==== Complexidade das interações ====

Vamos calcular o número de interações possíveis em  modelo com 4 preditoras:

<WRAP center round tip 80%>


=== Combinatória ===

Combinatória é a operação matemática para calcular de quantas maneiras conseguimos organizar ou combinar um conjuntos de elementos. No nosso caso, temos 4 variáveis e podemos nos perguntar de quantas formas podemos combinar esses elementos em diferentes conjuntos de dois a dois (interação dupla) ou três (interação tripla) e quatro (interação quadrupla)((Essa operação é chamada de combinatória simples, pois a ordem dos elementos nas combinações não importa)). A expressão matemática para essa operação é:

$$C^{r}_{n} =  \frac{n!}{r! (n-r)!} $$

onde:
 n: número total de elementos
 r: número de elementos combinados

A expressão ''!'' é a operação fatorial na matemática.

Para calcular o número de combinações possíveis de grupos de 2 e 3 elementos para um total de 4, temos:

$$C^{2}_{4} =  \frac{4!}{2! (4-2)!} = 6$$
$$C^{3}_{4} =  \frac{4!}{3! (4-3)!} = 4$$
$$C^{4}_{4} =  \frac{4!}{4! (4-4)!} = 1$$ 

</WRAP>

A lista de interações triplas e quadruplas tem essas expressões:


  * ''gestation:age:maeKg:smokeTRUE''
  * ''age:gestation:maeKg''
  * ''age:gestation:smoke''
  * ''age:maeKg:smoke''
  * ''gestation:maeKg:smoke''

As interações duplas:

  * ''age:gestation'' 
  * ''age:maeKg''
  * ''age:smoke''
  * ''gestation:maeKg''
  * ''gestation:smoke'' 
  * ''maeKg:smoke'' 

Além dessas, temos os termos isolados para cada variável e o intercepto do modelo, totalizando os 16 termos do modelo ''lmFull''. 

Se tivéssemos mais uma preditora contínua em nosso modelo o número de interações possíveis subiria para ''26'' e o total de termos do modelo para ''31''. 


Note que mesmo com poucas variáveis as possibilidade de interações e grande para decidirmos o que será incluído no modelo. Esse processo não é trivial e deve estar embasado no conhecimento prévio do sistema e na teoria para definir aquilo que faz sentido e pode ser interpretado, caso seja mantido no modelo. 

</WRAP>




Vamos agora comparar esse modelo cheio com o modelo sem nenhuma preditora. 

<WRAP center round box 80%>

===== Comparação com o modelo nulo =====

  * Produza o modelo nulo chamado ''lmNull'' com a formula:
      * ''bebeKg ~ 1''
  * compare com o modelo cheio produzido acima ''lmFull'' utilizando a função:
      * ''anova(lmNull, lmFull)''

</WRAP>


Nosso modelo ''lmFull'', apesar de não ter nenhum termo com coeficientes significativamente diferentes de zero, explica uma porção razoável e "significativa" da variação dos dados. Por volta de 25% da variação dos dados é explicada pelo modelo, por outro lado, não temos nenhum termo do modelo que valha a pena interpretar. Um poder de predição razoável mas com nenhum poder de interpretação dos processos subjacentes. 


Faça uma tentativa de interpretar o que o termo mais complexo deste modelo, a interação que aparece na última linha dos coeficientes estimados, está informando:


''age:gestation:maeKg:smoke[T.TRUE]''


O que significa uma interação entre as características da mãe: idade, dias de gestação, altura e que é fumante? Normalmente, conseguimos interpretar interações de segundo nível (levando em conta duas variáveis), algumas poucas vezes há sentido em interpretar a interação de três preditoras. Então o primeiro passo é retirar qualquer termo de interação de ordem mais alta. 

===== Definindo as interações =====

Para modelos com várias variáveis preditoras é sempre complicado definir quais interações devem ser contempladas. E como vimos acima, a estimativa e a interpretação dos coeficientes do modelo ficam comprometidos. Se a intenção do modelo é entender os processos causais da variação na variável resposta, é importante iniciarmos a seleção do modelo mínimo adequado com um modelo cheio que contemple apenas os termos que são importantes para o processo em questão.  







===== Modelo Cheio =====

Depois de definir quais são os termos que queremos incluir no nosso modelo (variáveis simples e interações), podemos iniciar o procedimento de modelagem seguindo algum protocolo para chegar ao modelo mínimo adequado. No nosso caso, iremos partir do modelo cheio, simplificar até o mínimo adequado. Como não temos experiência prévia do sistema e não temos "muita" experiência sobre tamanhos de bebês ao nascer, "consultamos vários especialistas na área" e eles chegaram à conclusão que os termos que deveriam ser contemplados são: 

 ''bebeKg ~ gestation + age + maeKg + smoke + gestation:age + gestation:smoke + 
            age:maeKg + age:smoke + maeKg:smoke + gestation:age:smoke ''


<WRAP center round box 80%>

==== Seleção do mínimo adequado ====

  * No menu ''Estatísticas > Ajustes de modelos > Modelo linear...'', construa o modelo com a seguinte expressão:


<code> bebeKg ~ gestation + age + maeKg + smoke + 
            gestation:age + gestation:smoke +  
            age:maeKg + age:smoke + maeKg:smoke + 
            gestation:age:smoke
</code>

  * simplifique esse modelo até o mínimo adequado;

<WRAP center round tip 80%>

Durante o processo de simplificação, quando nos defrontamos com vários termos de mesma ordem não significativos, um bom procedimento é retirar um deles de cada vez e, mesmo que o termo não seja retido no modelo, retorna-lo ao modelo antes de retirar o outro. Caso a ordem de retirada não torne nenhum dos termos significativos, ambos podem ser retirados. No caso de serem mais do que dois termos de mesma ordem, é importante também testar a retirada de dois a dois termos depois do procedimento de retirada de um a um não ter tornado nenhum termo significativo. Isso garante que a ordem de retirada não define o termo que será retido no modelo. Lembre-se que um termo "não significativo" em um modelo mais complexo pode se tornar "significativo" em um modelo mais simples. 

</WRAP>

  * qual a predição do modelo selecionado do peso do bebê das seguintes mães:
      * 26 anos, 280 dias de gestação, peso 55 kg e não fumante
      * 26 anos, 280 dias de gestação, peso 55 kg e fumante
      * 26 anos, 210 dias de gestação, peso 55 kg e não fumante
      * 26 anos, 210 dias de gestação, peso 55 kg e fumante
      * 26 anos, 280 dias de gestação, peso 70 kg e não fumante
      * 26 anos, 280 dias de gestação, peso 70 kg e fumante
      * 26 anos, 210 dias de gestação, peso 70 kg e não fumante
      * 26 anos, 210 dias de gestação, peso 70 kg e fumante
  * faça as mesmas predições acima, agora para uma mãe com 35 anos;
  * coloque as características da mãe e as predições do modelo para os casos acima em um planilha ((essa planilha pode ser feita no R como um ''data.frame'', exportando depois como ''csv''));
  * interprete o resultado do modelo biologicamente.

</WRAP>



==== Tamanho e sobrepeso da mãe ====

Um "outro especialista", ao analisar o modelo mínimo adequado selecionado acima, afirmou que a variável ''height'' deveria ser incluida para controlar o peso do bebê ser maior apenas porque a mãe era grande.   Além disso, sugeriu incluir uma nova variável que indicasse a quanto a mãe tem de sobrepeso ou subpeso antes de engravidar. Indicou ainda que uma forma de criar essa variável seria utilizar os resíduos de uma regressão simples da variável peso como resposta e a variável altura como preditora. Todos os "outros especialistas" concordaram que  eram uma boa sugestão, já que essa nova variável representa o quanto a mãe tem mais peso ou menos peso do que esperado para uma mãe com a mesma altura.  
 
<WRAP center round box 80%>
  * inclua a variável ''height'' no modelo mínimo selecionado anteriormente e verifique se deve se essa variável deve ser retida no modelo, com a comparação dos modelos pela partição da variação; 
  * crie uma nova variável ''resPeso'' com a expressão ''residuals(lm(weight ~height))'';
{{  :cursos:planeco:roteiro:resPeso.png?700  |}}
  * construa um modelo que inclua a variável ''resPeso'' como o primeiro termo no modelo mínimo adequado selecionado no procedimento anterior;
  * cronstrua um segundo modelo que inclua ''resPeso'' agora como o último termo no modelo mínimo adequado selecionado nos passos anteriores;
  * compare os resumos de ambos os modelos e anote as diferenças encontradas;
  * no menu ''Estatísticas'' > ''Resumos'' > ''Teste de correlação'' selecione as variáveis ''resPeso'' e ''maeKg'', anote o valor da correlação entre essas variáveis;
{{  :cursos:planeco:roteiro:corPeso.png?700  |}}
  * a partir da avaliação da correlação acima tome a decisão de reter ambas, uma ou nenhuma das duas variáveis do passo acima, justifique sua decisão e construa o modelo resultante;
  * simplifique o modelo do item anterior para o mínimo adequado, caso necessário;
  * faça uma interpretação biológica do modelo final.
</WRAP>












===== Um procedimento =====


Antes de montar o modelo precisamos nos perguntar e avaliar na literatura quais efeitos e interações fazem sentido e devem ser incorporadas no modelo. Vamos olhar, por exemplo, a primeira interação dupla que aparece no modelo que construímos acima:

''age:gestation'' 

Faz sentido e há alguma evidência que o peso do bebê depende da idade da mãe e que esse efeito da idade seja diferente dependendo do tempo de gestação? Caso faça sentido e/ou tenhamos evidências de fontes externas que essa interação é plausível, mantemos ela no nosso modelo. Um bom exercício é tentar explicar a interação em termos biológicos e isso fazer sentido. Por exemplo:

<WRAP center round tip 60%>
==== Interação   ====

Uma interpretação plausível para a interação ''age:gestation'', poderia ser: 


"A idade da mãe afeta negativamente o peso do bebê ao nascer, entretanto, se o tempo de gestação for longo, o efeito negativo se perde" 

</WRAP>


Vamos avaliar outra interação dupla:

''height:weight''

Faz sentido esperar que pessoas maiores tenham bebês maiores, além disso, faz sentido esse efeito ser modulado pelo peso da mãe, já que se ela tem sobrepeso ou subpeso, pode afetar o desenvolvimento do feto. 



Sem conhecer bem o tema, imaginamos porque a altura da mãe teria um efeito no peso do bebê, e ainda por cima modulada pelo peso da mãe? Faria mais sentido dizer que o peso (altura) isoladamente não afetam o peso do bebê ao nascer, mas sim o sobrepeso ou o subpeso da mãe. Mais ainda, diria que o peso e altura da mãe são variáveis correlacionadas e devem ter alguma colinearidade na explicação do peso do bebê.






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


*/

===== Formulário de resposta =====

<WRAP center round help 90%>

Responda o  [[https://docs.google.com/forms/d/e/1FAIpQLSewq_huLqaDEA-Hu6hxGC-J1tOXE7k4wcHojIbNjiF4yC5LLg/viewform?usp=sharing|o formulário MLM I]] incluindo arquivos de resultados quando solicitado. 


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

*/



/*

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


<WRAP center round help 100%>
 
  * 1. As estimativas dos modelos devem ter sido incluídas nas planilhas quando foram solicitados ao longo do roteiro 
  * 2. Preencha as perguntas do quadro abaixo ou pelo [[https://forms.gle/LvN2j8iE7JWpeEBr9|link do formulário]]


</WRAP>

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

*/