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 Mistos (LMM) ====== 


<WRAP center round info 80%>
__**Objetivo desse roteiro**__ 

Ao final que você seja capaz de: 

  * 1. Entender o que são efeitos fixos e aleatórios.  
  * 2. Compreender a estrutura básica de um modelo linear misto.  
  * 3. Fazer uma análise de modelo misto no **Rcmdr** usando o pacote ''lme4'' (Bates et al. 2014)   
  * 4. Entender o //output// da função ''lmer''.  
  * 5. Decidir quais efeitos aleatórios e fixos manter no seu modelo final.
  * 6. Interpretar o resultado de um modelo misto selecionado.  

</WRAP>


 

/*
<WRAP center round tip 100%>

===== Antes de começar =====

O Rcmdr não tem uma interface para a instalação de pacotes e, normalmente, precisa ser reiniciado para carregar o  pacote instalado. Como isso pode levar a que o trabalho feito anteriormente seja perdido, sugerimos fortemente que instale os seguintes pacotes antes de abrir o Rcmdr, ou seja, logo no inicio da sessão do R. Para isso, copie e cole as seguintes linhas de comando:

<code>

pkgs <- c("lme4", "RVAideMemoire")
pkgsNot <- pkgs[!(pkgs %in% rownames(installed.packages()))]
install.packages(pkgsNot)

</code>


</WRAP>

*/



===== Modelos Lineares Mistos (LMM) =====

Os modelos Lineares Mistos são uma generalização de modelos lineares que afrouxam uma premissa básica de todos os modelos que vimos até agora: <wrap em>a independência entre as observações</wrap>. Essa falta de independência emerge por vários motivos em experimentos e muitas vezes são desejáveis pois controlam fatores de confusão. Entre as principais fontes de onde emergem essas dependências estão:

  * espacial
  * temporal
  * biológicas

No primeiro caso, é esperado que objetos que estão mais próximos sejam mais parecidos entre si do que objetos mais longe, pois compartilham muitas condições associadas ao mesmo espaço. Dependências temporais emergem principalmente quando medimos um mesmo objetos de estudo em diferentes momentos. As biológicas estão relacionadas ao compartilhamento de caracteres ou funcionalidades entre indivíduos mais relacionados geneticamente, por exemplo. 

Essa dependência portanto, cria agrupamentos de observações que compartilham similaridades que não são compartilhadas com observações de outros grupos. 

Em geral, não estamos interessados nas diferenças que existem entre esses grupos, mas precisamos contemplar a falta de independência de alguma forma. Os **Modelos Lineares Mistos** justamente lidam com essas falta de independência desse agrupamentos de observações.  Para isso, os **LMM** incorporam uma estrutura aletoria adicional nos modelos lineares para acomodar esses agrupamentos definidos agora como **variáveis categóricas preditoras aleatórias**.  Diferente das **variáveis categóricas preditoras fixas** onde o interesse é interpretar as diferenças entre os níveis, nas aleatórias queremos apenas estimar a variabilidade associada aos agrupamentos.

<WRAP center round box 80%>

==== Variáveis categóricas fixas ou aleatórias ====

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

Uma forma de diferenciar as variáveis categóricas fixas e aleatórias é imaginar a realização de múltiplos experimentos, que está por trás da teoria da estatística frequentista. Em um experimento usual com a variável categórica fixa estamos interessados na estimativa das médias para cada nível da variável e  buscamos saber se esses níveis são diferentes em média, como exemplificado na primeira realização do experimento (**1**)  no painel da esquerda (**a**) da figura acima. Ao refazer o mesmo experimento, iriamos recoletar os mesmos níveis (experimentos **2** e **3** ) e estimar valores de médias muito próximos, como exemplificado no painel **a** da figura acima no seu canto direito. Por outro lado, quando temos uma variável categórica aleatória não estamos interessados em nenhum nível específico e outro experimento iria coletar outros níveis aleatórios ((por isso o nome!)) dentro todos os níveis possíveis dessa variável categórica, como mostra o painel **b** da figura acima. Neste caso, estamos interessados, não na diferença das médias entre os níveis e sim, na variabilidade associada à população de níveis dessa variável categórica, como mostrado no canto direito do painel **b**.    

</WRAP>



O vídeo abaixo apresenta uma breve introdução sobre os modelos lineares mistos .

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

</WRAP>

===== Crescimento de Árvores =====


Vamos contextualizar os **LMM** a partir de um exemplo e suas diferença com os modelos **LM**. 
O interesse do estudo do exemplo a seguir foi entender se a **deciduidade** (perda de folha síncrona) de árvores está relacionada ao crescimento do indivíduo. Algumas árvores perdem as folhas de forma sincronizada (decíduas), enquanto outras trocam de folhas ao longo do tempo sem haver essa sincronicidade (perenes).  

A pergunta: árvores decíduas crescem menos que árvores não decíduas? Tem seus desafios e fontes de confusão. Um complicador é que o crescimento está diretamente relacionado com o tamanho da árvore. Por isso, os pesquisadores decidiram incluir o tamanho da árvore como uma preditora fixa e com isso  condicionaram o crescimento à variável tamanho. Além disso, a perda de folha síncrona é característica da espécie e eles precisaram incluir essa variável no desenho experimental.
Ao final, os pesquisadores selecionaram ''5'' espécies decíduas e outras ''5'' perenes e, para cada espécie anotaram o crescimento de ''10'' indivíduos, totalizando 100 indivíduos. 
<WRAP center round box 80%>
**__Deciduidade em Árvores__**


{{youtube>EN2_5l8FXQk  |}}

</WRAP>


Vamos iniciar nosso exemplo com uma análise de modelo linear simples, sempre testando a hipótese que o crescimento de árvores decíduas e perenes é diferente.

<WRAP center round box 80%>
Cresce árvores: modelo linear

  - baixe os dados do experimento {{ :cursos:planeco:roteiro:alturaCresc.txt |cresce arvores}};
  - faça o modelo ''cresc ~ alt + dec + alt:dec '';
  - simplifique até o mínimo adequado;
  - interprete o modelo selecionado.

</WRAP>
 

==== Procedimento em duas etapas: cresce árvores ====

O problema com o modelo que fizemos acima é que desrespeitamos uma pressuposto muito importante para os modelos lineares: <wrap em> a independência entre as observações </wrap>. Os indivíduos de uma mesma espécie podem compartilhar características que poderiam fazer o crescimento de todos os indivíduos desta espécies crescer mais ou menos, independente se a espécie é decídua ou não. Por essa razão precisamos tratar essa dependência entre as observações associada aos indivíduos serem da mesma espécie. Antes de termos a instrumentação dos modelos mistos ou outras técnicas para lidar com essas dependências, uma forma de tratar as dependências entre observações era a redução da informação, por exemplo a média, para cada nível de agrupamento desta dependência, no caso aqui ilustrado, para cada uma das espécies. Como temos o crescimento dependendo do tamanho do indivíduo, a média aqui não seria uma boa redução de informação pois poderíamos ter amostrado indivíduos maiores em um espécie e menores em outra, simplesmente por acaso. Para esse caso, o que podemos fazer é a redução da relação ''cresc ~ alt'' para cada espécie e em seguida testar se as inclinações das espécies decícuoas são diferentes das espécies perenes. Se as inclinações não forem diferentes, ainda resta testar se os interceptos são diferentes. Vamos fazer esse teste. 

<WRAP center round box 80%>
  - para cada uma das 10 espécies faça o modelo ''cresc ~ alt'';
  - guarde os coeficientes de inclinação e intercepto para cada uma das espécies junto a informação se é decídua ou perene;
  - faça o modelo que testa se a inclinação é  diferente entre decíduas e perenes;
  - caso o a inclinação não seja significativamente diferente entre os dois grupos, faça o mesmo para a intercepto;
  - interprete o resultado.
</WRAP>

 

==== LMM: cresce árvores ====

Vamos agora modelar esses dados com o modelo linear mistos. O **misto** vem da ideia de incorporar variáveis fixas e aleatórias como preditoras. Nestes modelos a estrutura aleatória preditora pode influenciar o intercepto ou a inclinação ((a modelagem da inclinação pressupõem o intercepto)). 
A influência no intercepto significa que a estimativa de interceptos para as espécies provêm de uma distribuição normal de interceptos definidos por um intercepto médio e um dado desvio padrão. A expressão geral para um modelo com uma preditora e a variável aleatória afetando apenas o intercepto é:  

$$ y_{ij} = (\bar{\alpha} + \epsilon_j) + \beta x_{ij}  + \epsilon_{ij} $$




$$ \epsilon_j = N(0, \sigma_{entre}) $$

$$ \epsilon_{ij} = N(0, \sigma_{intra}) $$

O $\bar{\alpha}$ é o intercepto médio e o $\sigma_{entre}$ é a estimativa do desvio padrão associado à distribuição de interceptos para a variável aleatória. Ou seja, no nosso caso, a estimativa da variabilidade de interceptos associado às espécies.  

O modelo em que a variável aleatória influencia a inclinação tem a seguinte expressão: 

$$y_{ij} = (\bar{\alpha} + \epsilon_{\alpha_j}) + (\bar{\beta} + \epsilon_{\beta_j})  x_{ij} + \epsilon_{\text{res}_{ij}}$$

Note que não ao influenciar a inclinação, por consequência o intercepto também é influenciado pela variável aleatória.  Neste caso o modelo estima o intercepto médio, a inclinação média e os seguintes desvios padrão:

$$ \epsilon_{\alpha} = N(0, \sigma_{\alpha}) $$
$$ \epsilon_{\beta} = N(0, \sigma_{\beta}) $$
$$ \epsilon_{\text{res}} = N(0, \sigma_{\text{res}}) $$


<WRAP center round tip 80%>


=== LMM no Rcmdr ===

Existem vários pacotes disponíveis no R para realizar análises de modelos mistos. Neste roteiro usaremos o ''lme4'' (Bates et al. 2014), que possui funções para analisar modelos lineares mistos, modelos lineares mistos generalizados e modelos mistos não lineares. 

Desde o ano de 2020, os modelos mistos passaram a ser contemplados no Rcmdr. Por isso, podemos utilizar essa interface para fazer os modelos. Nas versões recentes do Rcmdr o pacote ''lme4'' é carregado por padrão quando se constroi um modelo misto. Caso ele não seja carregado automaticamente, é possível carregá-lo. Para isso, vá em ''Tools -> load packages'' e escolha ''lme4''((Caso queira usar as funções que utilizaremos a seguir no console do ''R'' diretamente, você precisará instalar e carregar o mesmo pacote)). 


</WRAP>



Para criar o nosso modelo cheio no R com a variável aleatória influenciando a inclinação, a notação é:

<code> lmer(cresc ~ alt + dec + alt:dec + (alt|sp)) </code>

A expressão ''(alt|sp)'' está indicando que a variável aleatória ''sp'' irá modelar a inclinação da variável contínua ''alt''.

O modelo em que a variável ''sp'' influência apenas o intercepto é construído com a seguinte expressão: 

<code> lmer(cresc ~ alt + dec + alt:dec + (1|sp)) </code>

O código ''(1|sp)'' indica que a variável ''sp'' vai influenciar apenas o intercepto. 

<WRAP center round box 90%>
==== Selecionando a estrutura aleatória do modelo ====

Assim como na estrutura fixa do modelo podemos reduzir a estrutura aleatória ao mínimo adequado para interpretar. Para isso, seguimos o procedimento descrito em Zuur et al. 2009 ((veja referência ao final do roteiro)).
  - defina estrutura cheia das variáveis fixas;
  - inicie com a forma mais complexa da estrutura aleatória;
  - compare modelos aninhados por simplificando a estrutura aleatória;
  - utilize a função ''anova'' com o argumento ''refit = FALSE'' quando comparar modelos com diferentes estruturas aleatórias;
  - finalize a simplificação da estrutura aleatória ((a equipe da disciplina acredita que ao menos a influência no intercepto deve ser mantida para contemplar a dependência do desenho experimental)),
  - faça a simplificação da estrutura fixa do modelo como em um modelo ''lm'', mantendo a estrutura aleatória da forma que foi selecionada nos passos acima. 

<WRAP center round todo 80%>
==== Atividade ====

  * baixe o arquivo de {{ :cursos:planeco:roteiro:crescArv.txt |dados de crescimento de árvores}} ((apresenta campos separados por tabulação));
  * leia o arquivo de dados no Rcmdr;
  * construa o modelo cheio;
  * faça a simplificação do modelo seguindo os passos indicados acima;
  * interprete o modelo selecionado.

</WRAP>

</WRAP>


<WRAP center round help 90%>

**PARA ENTREGAR ANTES DA PRÓXIMA AULA**
Acesse o [[https://forms.gle/Uh7QZvrKmcFFNkHX6|formulário]] e responda às questões propostas referentes aos dados sobre crescimento de árvores. 


</WRAP>





===== Riqueza em praias =====

<WRAP center round box 60%>

{{https://d3i71xaburhd42.cloudfront.net/e00b5c0fea6fb05b1370cc3dc537ed2e821033b7/5-Figure1-1.png?200  |}}


Esse exemplo no roteiro é baseado em material desenvolvido por alunos de nosso programa de PG ((A versão original está disponível neste [[https://melinaleiteblog.netlify.com/2017/09/19/introducao-aos-modelos-mistos/|site]]))[[https://melinaleite.weebly.com/|Melina Leite]], [[https://mariliagaiarsa.weebly.com/|Marília Gaiarsa]] e [[http://www.guimaraes.bio.br/people.html|Lucas Medeiros]] e vem sendo modificado, ao longo dos anos, para adaptá-lo ao curso. 

</WRAP>


Para exemplificar o que são modelos mistos e sua importância em ecologia, usaremos como exemplo um conjunto de dados presente no [[cursos:planeco:roteiro:11-lmm#referências_e_recomendações|capítulo 5 de Zuur et al. (2009)]]. Esses dados são referentes a um estudo sobre a riqueza de espécies da macrofauna em 9 praias na costa da Holanda. Em cada uma das praias os autores coletaram dados em cinco localidades diferentes. Para cada localidade existe informação sobre a altura da estação de amostragem em relação à altura média da maré (''NAP'', variável contínua) e também um índice de exposição da praia (''Exposure'', variável categórica). 

Vamos supor que estamos interessados em verificar se a altura em relação à altura média da maré (variável ''NAP'') influencia a riqueza de espécies nessas praias, deixando de lado a variável de exposição nesse momento. Podemos, por exemplo, construir um **modelo linear simples** da seguinte forma:

$$\text{riqueza} = \alpha + \beta * \text{NAP} + \varepsilon$$

sendo que $\varepsilon$ é o resíduo com distribuição normal (ou gaussiana) de média zero e um dado desvio padrão:

$$ \varepsilon = N(0, \sigma) $$

Aqui estamos modelando a riqueza de cada ponto de amostragem em função da variável ''NAP''. O coeficiente $\alpha$ é o intercepto do modelo, $\beta$ é a inclinação (coeficiente angular ) de ''NAP'', que definem a esperança (média) do modelo.  O resíduo do nosso modelo ($\varepsilon$) é a variação da riqueza que não conseguimos explicar com nossa variável preditora ''NAP''. Este é o  escopo dos modelos lineares que estudamos até o momento.

No entanto, neste exemplo o modelo viola uma premissa fundamental de modelos lineares, a de que **as observações que geram os dados são independentes umas das outras** (Winter, 2013). Os dados obtidos em uma mesma praia não são independentes entre si (dependência espacial). É possível imaginar diversas características particulares de cada praia, como o tipo de grão de areia ou a inclinação da praia, que podem influenciar a riqueza de espécies e que não estão contempladas nesse modelo. Deste modo, os valores de riqueza obtido pelas unidades amostrais de uma mesma praia tendem a ser mais parecidos entre si, independente da relação com as variáveis preditoras que estamos analisando. Veja o gráfico abaixo, em que cada cor representa uma praia, para avaliar visualmente se a relação entre riqueza e ''NAP'' é diferente entre praias:

<WRAP center round box 90%>

**__Modelos Lineares para cada praia__**

{{:cursos:planeco:planeco:roteiro:lmm1.png?650|}}

</WRAP>

Nesta figura temos as retas do modelo linear simples entre ''Riqueza'' e ''NAP'' aplicado para cada praia . Entretanto, a pergunta do estudo não diz respeito a cada praia separadamente, o objetivo é modelar a relação para todas as praias amostradas, sem especificar o que acontece para cada praia isoladamente. Por outro lado, se considerarmos cada ponto de amostragem como se fossem independente entre si, estaremos incorrendo em um erro denominado de pseudo-replicação ou falta de independência das observações((podemos dizer também que as observações da mesma praia são correlacionadas)). Ou seja, queremos considerar essa dependência sem falar especificamente de nenhuma praia e, ao mesmo tempo, avaliar a variação associada às praias. Além disso, queremos que as praias amostradas representem uma amostra aleatória de todas as praias que poderíamos  contemplar no nosso estudo.

==== Efeito Aleatório: riqueza em praias ====

Para isso, precisamos de um modelo linear que incorpore o fato de que nossos dados estão agrupados em praias. Chamamos de **efeito aleatório** uma variável que agrupa nossos dados e que, geralmente, seu efeito sobre a variável resposta não nos interessa diretamente ((nesse exemplo não interessa, mas dependendo da análise, pode interessar - veja discussão em McGill 2015)). Nesse exemplo, os autores amostraram ''9'' praias, mas poderiam ter feito o mesmo estudo amostrando outras ''9'' praias em uma amostra aleatória. Nesse sentido, as ''9'' praias da amostra são apenas ''9'' possíveis praias amostradas aleatoriamente de todas as possíveis ''9'' que poderiam ter sido amostradas ao acaso. Por esta razão, ''praias'' é um variável aleatória.

O efeito aleatório ''praia'' __organiza__ parte da variação nos nossos dados que o modelo não  explicaria com as variáveis preditoras e que estaria presente no resíduo($\varepsilon$) do modelo linear simples. Por outro lado, as __variáveis preditoras__ usuais em modelos lineares são chamadas de **efeitos fixos**. No exemplo das praias, a variável de diferença do nível da maré (''NAP'') é uma variável de efeito fixo e estamos diretamente interessados em seu efeito sobre a variável resposta (''Richness'').  O nome **misto** refere-se ao fato de  existirem tanto efeitos fixos, quanto aleatórios nas variáveis preditoras do modelo. Esses modelos são chamados também de **__hierárquicos__** ou **__multinível__**. 

Em ecologia é comum encontrarmos delineamentos amostrais ou experimentais que geram dados com algum tipo de agrupamento (delineamento hierárquico ou aninhado). Por exemplo, quando uma amostragem é feita por parcelas ou quando um experimento é separado em diferentes blocos. Nesses casos, não podemos tratar nossos dados como independentes e uma abordagem estatística adequada é a de **modelos mistos**. 

Inicialmente, a forma mais simples de incorporar o efeito aleatório em nosso exemplo é definir que cada praia apresenta um intercepto diferente no modelo. Ou seja, queremos ajustar uma reta para nossa relação entre ''riqueza'' e ''NAP'', mas permitir que cada praia tenha seu próprio intercepto, que é relacionado ao intercepto da reta **principal** (efeito fixo). 

Neste caso, o modelo terá a seguinte estrutura:

$$\text{riqueza}_{i} = \alpha + \epsilon_{\text{praia}_i} + \beta * NAP + \epsilon_{\text{res}}$$

sendo:

$$\epsilon_{\text{praia}} = N(0, \sigma_{\text{praia}})$$

$$ \epsilon_{\text{res}} = N(0, \sigma_{\text{res}})$$

A riqueza da praia ''i'' é explicada pelo efeito fixo $\beta*NAP$ e pelo efeito aleatório $\epsilon_{\text{praia}_i}$, que se soma ao valor do intercepto fixo do modelo, $\alpha$, para formar o intercepto da praia ''i''. Veja que o índice ''i'' está atrelado a $\epsilon_{\text{praia}_i}$ e, portanto, o modelo permite que cada uma das praias apresente um intercepto diferente com mesmas inclinações $\beta$. O $\beta$ faz parte do efeito fixo e não possui nenhum índice ''i'' atrelado a ele, portanto, é o mesmo para todas as praias. Finalmente, $\epsilon_{\text{res}}$ representa o resíduo associado aos desvios dos valores observados em relação aos preditos pela reta da praia relativa à observação.  

Ambos, $\epsilon_{\text{praia}}$ e $\epsilon_{\text{res}}$, são variáveis aleatórias com média ''zero'' e desvios padrão $\sigma_{\text{praia}}$ e  $\sigma_{\text{res}}$, respectivamente. 
 
====  LMM: riqueza em praias =====


/*

Antes de tudo, precisamos baixar e instalar o pacote ''lme4'':

<code>
##caso não tenha instalado o pacote tire o comentário da linha seguinte
#install.packages("lme4")
library(lme4)
</code>

*/


/*

Nos dados disponibilizados aqui no nosso roteiro criamos uma variável categórica ''fExp'', com dois níveis categóricos "low" e "high")), mas para esse roteiro recomendamos que os dados sejam obtidos pelo link abaixo, pois já fizemos alguns ajustes necessários no arquivo do link. 

*/




==== Lendo os dados praias no R ====

<WRAP center round box 100%>


  * baixe o arquivo {{ :cursos:planeco:roteiro:praia01.txt | vamos à praia}};
  * no Rcmdr, importe o arquivo, que está formatado com campos separados por tabulação e decimal como ponto;
  * atribua o nome de ''praia'' ao conjunto de dados;
  * visualize o conjunto de dados para se familiarizar com as variáveis.

</WRAP>

/*
<code>

dados <- read.table("praia.txt", header = TRUE, sep = "\t" , as.is = TRUE)
head(dados) #observando as primeiras linhas de dados
str(dados) # vendo se a estrutura dos dados está ok!

</code>
##### RETIREI ESSA PARTE TODA EM 2020, POIS O ARQUIVO "praia.txt" JÁ TEM A VARIÁVEL fExp COM "low" e "high".

*/
/*
<code>

dados$fExposure <- factor(dados$fExp, levels = c("low", "high"))

</code>
*/

==== Um modelo plausível: riqueza em praias ====

Vamos construir um modelo simmples que tem apenas ''NAP'' como variável fixa preditora e a variável ''Beach'' como aleatória apenas no intercepto.  Esse modelo pode ser descrito no ''lme4'' da seguinte forma((pacotes podem ter diferentes sintaxes para descrever modelos mistos)):


<code>
Richness ~ NAP +  (1 | Beach)
</code>


<WRAP center round box 100%>


No Rcmdr siga os passos:

  * abra o menu ''Statistics'' > ''Fit models'' > ''Linear mixed models...''
  * nomeie o modelo como ''lmm.riq'' e coloque a fórmula acima nos campos adequados

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

</WRAP>



Em seguida, vamos dar uma olhada no resultado deste modelo que aparece na janela inferior do Rcmdr.

==== Interpretando o modelo: riqueza em praias ====

Na figura abaixo temos o resultado resumido do modelo ''lmm.riq'':

<WRAP center round box 100%>
{{  :cursos:planeco:roteiro:lmmriqSummary.png?500  |}}
</WRAP>

Vamos começar olhando a parte dos efeitos aleatórios (''Random effects''). O desvio padrão (''Std.Dev.'') é uma medida do quanto a variabilidade da nossa variável dependente  ''Richness'' é particionada aos efeitos aleatórios que estamos analisando (''Beach'') e o quanto não foi explicado (''Residual''). Podemos ver em ''Beach'' o desvio padrão associado às diferenças de intercepto entre praias ( $\sigma_{\text{praia}}$ do modelo associado ao $\epsilon_{\text{praia}}$). A última linha apresenta o resíduo, que indica o quanto da variabilidade não é prevista pela variável aleatória  ''praia'' nem pela variável fixa ''NAP'', ou seja, o que não é explicado por nenhum termo do nosso modelo. 

Na seção ''Fixed effects'' são apresentados os coeficientes estimados para cada um dos fatores que estamos considerando como fixos. No caso, temos o intercepto médio((todas as praias)) associado à riqueza quando o ''NAP'' é zero (''6.5819'') e o coeficiente angular de ''NAP'' (''-2.5684''), que indica que há uma diminuição de riqueza da ordem de mais de ''2.5'' espécies a cada unidade de ''NAP'' acrescida.

Note que nesse modelo a inclinação é a mesma para todas as praias e que o intercepto da estrutura fixa é o intercepto médio dos coeficientes da estrutura aleatória do modelo.

A correlação dos efeitos fixos é o quanto os coeficientes fixos estimados são correlacionados. Mais formalmente é uma estimativa do quanto teríamos de correlação entre o intercepto e a inclinação se refizessemos o experimento muitas vezes((No caso, a correlação esperada entre os valores de intercepto e inclinação seria por volta de ''15.7%''. Essa correlação é calculada à partir da matriz de variância/covariância do modelo e pode ser calculada para qualquer modelo com estrutura fixa. Uma confusão comum é achar que esse valor está associada necessariamente à multicolinearidade das variáveis preditoras, o que não é o caso.)). Apesar de ser apresentada por padrão no resumo do ''lme4'' essas correlações não tem uma interpretação direta e muito menos intuitiva ((Segundo alguns autores é útil apenas em casos especiais e não deveria ser apresentada por padrão. Para retirar essa correlação do resumo do modelo é necessário usar o argumento ''correlation = FALSE'' na função ''summary'')).  


/*

<code>

mean(randMM[,1])

</code>

*/

==== Predito pelo Modelo: riqueza em praias ====

Uma boa maneira de interpretar os resultados é a partir do gráfico do predito pelo modelo, se possível junto aos valores das observações dos dados. Temos apresentado gráficos com a experança do modelo ao longo de todo o curso na intenção de auxiliar na interpretação do modelo. Além disso, calcular o predito através dos coeficientes do modelo ajuda a entender a estrutura atrás dos coeficientes.

Para fazer o gráfico do modelo ajustado, com a reta média predita pela estrutura fixa do modelo e as retas da estrutura aleatória, preditas para cada praia, vamos criar um objeto no Rcmdr com todos os coeficientes necessários. Copie as linhas de código a seguir no painel superior do Rcmdr (o painel chamado //Rscript//), em seguida, **selecione todas as linhas copiadas** e clique no botão ''Submit'' ((Obs.: As linhas que começam com ''#'' não são lidas pelo Rcmdr e podem ser usadas para comentar as linhas de código, isso ajuda a lembrar o que a linha de código faz)).

<code>
#criando objeto com os coeficientes do modelo (efeitos fixos)

fixLMM <- fixef(lmm.riq)
fixLMM

#criando objeto com os coeficientes do modelo (efeitos aleatorios)

randLMM <- ranef(lmm.riq)$Beach
randLMM

</code>

<WRAP center round box 80%>

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


</WRAP>


O predito pela estrutura das variáveis fixas do modelo estimam o **modelo médio**, que seria a //esperança// do modelo para uma praia média, ponderada pelas observações de cada uma das praias. Além disso, a estrutura aleatória do modelo prediz um intercepto para cada praia das nossas observações a partir da variabiabilidade existente entre praias, definido pelo efeito aleatório $\epsilon_{praia}$.

 
<WRAP center round important 80%>
Note que nesse modelo temos apenas um coeficiente de inclinação, pois não modelamos a inclinação na estrutura aleatória do modelo, apenas o intercepto. Mais adiante faremos um modelo incluindo inclinação para a variável aleatória e poderemos comparar os resultados.
</WRAP>

==== Gráfico do Modelo: riqueza em praias ====


Vamos primeiro construir um gráfico base, com os valores observados. Para tornar o gráfico mais informativo vamos colocar cores para identificar em que praia cada observação foi coletada. 
No código abaixo temos as seguintes sequência de eventos: 

  * ''cores...'': cria um objeto com nove cores tiradas de uma paleta das cores do arco-íris
  * ''par(...) '': modifica parâmetros da janela gráfica, no caso abaixo modificamos as margens do gráfico
  * ''plot(...)'': cria o gráfico 
      * 1° linha: indica quais variáveis devem estar no eixo Y e X, e em qual objeto estão os dados
      * 2° e 3° linhas: modificações de parâmetros do gráfico, como símbolos, letras, etc

Então, copie os códigos abaixo para a janela superior do Rcmdr (RScript), selecione essas linhas e aperte o botão ''Submit'':

<code>
cores <- rainbow(9)
par( mar=c(5,5,2,2))
plot(Richness ~ NAP ,data = praia,
     pch = 19, col = cores[as.factor(Beach)] , las = 1,
     cex=1.5, cex.lab= 1.7, cex.axis = 1.5
     )

</code>

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



=== Incluindo o modelo no gráfico ===

Vamos usar os valores dos coeficientes (aqueles que calculamos anteriormente com as funções ''fixef'' e ''ranef'') para representar as esperanças do modelo e criar as retas para cada praia e a reta média. 

Primeiro o gráfico com as retas de cada praia:

Copie os códigos abaixo para a janela superior do Rcmdr (RScript), selecione essas linhas e aperte o botão ''Submit'':

<code>

abline(randLMM[1,'(Intercept)']+ fixLMM[1], fixLMM[2], col=cores[1])
abline(randLMM[2,'(Intercept)']+ fixLMM[1], fixLMM[2], col=cores[2])
abline(randLMM[3,'(Intercept)']+ fixLMM[1], fixLMM[2], col=cores[3])
abline(randLMM[4,'(Intercept)']+ fixLMM[1], fixLMM[2], col=cores[4])
abline(randLMM[5,'(Intercept)']+ fixLMM[1], fixLMM[2], col=cores[5])
abline(randLMM[6,'(Intercept)']+ fixLMM[1], fixLMM[2], col=cores[6])
abline(randLMM[7,'(Intercept)']+ fixLMM[1], fixLMM[2], col=cores[7])
abline(randLMM[8,'(Intercept)']+ fixLMM[1], fixLMM[2], col=cores[8])
abline(randLMM[9,'(Intercept)']+ fixLMM[1], fixLMM[2], col=cores[9])

</code>


Para finalizar o gráfico, vamos colocar o predito pelo modelo médio e uma legenda:

Repita o mesmo procedimento de copiar, selecionar e "submeter":
<code>

abline(fixLMM[1], fixLMM[2], col = "black", lwd=3, lty= 2)
legend("topright", c( "praia média",paste("praia", 1:9)) , lty=c(2, rep(1,9)), col= c(1,cores), bty="n", cex=1.2)


</code>

Ao final o nosso gráfico deve ser parecido com o que está a seguir:

<WRAP center round box 90%>

{{  :cursos:planeco:roteiro:praiasLMM.png?600  |}}

</WRAP>

/* 

<code>
# Primeiro criamos um novo conjunto de dados com as características no antigo
novo <- expand.grid(Beach = unique(dados$Beach),
                    NAP = seq(-1.3,2.2,0.5))

#agora calculamos os valores preditos para esse novo conjunto de dados
preditos <- predict(modelo.riqueza, newdata = novo)

#guardando tudo em um novo data.frame
dados.preditos <- data.frame(pred = preditos, novo)
</code>


Agora podemos plotar os dados com o ajuste do nosso modelo ((estamos usando o pacote `ggplot2` para fazer o gráfico, que tem uma sintaxe um pouco diferente de um gráfico do pacote base)).

{{:cursos:planeco:planeco:roteiro:lmm2.png?650|}}



{{:cursos:planeco:roteiro:LMMModel.png?600|}}


==== Código do Gráfico ====

Vamos usar um pacote do R chamado **ggplot2**. Para usá-lo é preciso instalar (''install.packages'') e carregar o pacote (''library'')

<code>
#não esqueça de instalar o pacote,caso ainda não esteja instalado, tire o sustenido da linha logo abaixo. 
#install.packages("ggplot2")

#library(ggplot2)

funi <- function(x){cof[1] + cof[2]*x} # para plotar a predição dos efeitos fixos
dados$Praia = as.factor(dados$Beach)# transformando praia em fator

ggplot(data = dados, aes(x = NAP, y = Richness, color = Praia)) + # dados e eixos
  geom_point(size = 3, shape = 19) +  # plotando os pontos das praias 
  geom_line(data = dados.preditos, aes(y = pred, x = NAP, 
                            col = as.factor(Beach))) + # retas de cada praia
  stat_function(fun = funi, col = "black", size = 2) +  # reta do modelo fixo 
  scale_color_brewer(palette = "Set1") + # a partir daqui estética do gráfico
  theme_bw() +
  theme(axis.text = element_text(size = 13), 
        axis.title = element_text(size = 15),
        axis.text.x = element_text(size = 11),
        axis.text.y = element_text(size = 11),
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 13)) +
  xlab("NAP") +
  ylab("Riqueza")

</code>

*/

Podemos ver nessa figura a predição do nosso modelo em relação aos parâmetros fixos (reta tracejada em preto) e as predições para cada praia separadamente. Como o efeito aleatório do nosso modelo estava apenas variando os valores de intercepto das praias, as retas para cada praia são paralelas.



/*


==== Inferência e diagnóstico do modelo: riqueza em praias ====

Nesta parte, depois de já escolhermos a estrutura aleatória do nosso modelo, podemos averiguar qual a real influência dos efeitos fixos sobre a riqueza de espécies, utilizando o procedimento que temos adotado nessa disciplina, ou seja o **teste de hipóteses por meio da comparação de modelos pela tabela de ANOVA**. Depois de selecionarmos o modelo mínimo adequado, que melhor se ajusta aos dados, vamos fazer o diagnóstico dos resíduos deste modelo para ver se ele atende às premissas de um modelo linear misto.


==== Seleção de modelo: riqueza em praias ====

O procedimento para verificar se termos do modelo são ou não significativos é a construção de modelos aninhados (ou seja, retirando um parâmetro do modelo com mais parâmetros) e a comparação por uma tabela de ''anova''. 

**Muita atenção aqui:**

No caso da **comparação da __estrutura fixa__ de modelos mistos com mesmo efeito aleatório** precisamos reajustar nossos modelos para que **as estimativas sejam feitas por ML (Máxima Verossimilhança)**. Mas, relembrando o que já foi dito anteriormente, a função ''anova'', por padrão, já reajusta os modelos para ML quando compara dois modelos. Portanto, não precisamos nos preocupar em adicionar nenhum termo, mas é importante saber que esses procedimentos estão sendo feitos. 

Então vamos começar usando nosso modelo cheio ''lmmint'' na comparação.

=== llmmint: ===


<code>

Richness ~ fExposure + NAP + fExposure:NAP + (1|Beach)

</code>


Vamos agora construir o modelo __sem interação__ entre ''fExposure'' e ''NAP'', denominando-o de **lmm01**.
Depois, comparamos esses modelos com o teste de hipótese por meio da função ''anova'' pelo menu do Rcmdr **Models>Hipothesis tests>Compare two models…**. E assim por diante, simplificando os modelos como já feito nos roteiros anteriores até chegar ao modelo mínimo adequado. Importante lembrar que a **estrutura aleatória** deverá ser sempre a mesma em todos os modelos que serão comparados!


/*

### Retirei toda essa parte para deixar que os alunos quebrem um pouco a cabeça para fazer essa comparação de modelos e chegar ao modelo mínimo. Eles já fizeram isso muitas vezes, então, será tranquilo)

Como o resultado da comparação entre os modelos foi significativo, nós retemos o modelo mais complexo, ou seja, o modelo com interação. O que podemos fazer é assegurar que o modelo com interação é também diferente do modelo nulo (sem efeitos fixos):

<code>
# modelo nulo

m6 <- lmer(Richness ~ 1 + (1|Beach), data = dados, REML = F)

anova(m1, m6)
</code>

Então podemos concluir que tanto ''Exposure'' quanto ''NAP'' influenciam na riqueza de espécies. Caso não tivesse havido diferença entre o modelo com interação e sem, nós deveríamos prosseguir com a seleção de modelos, fazendo a comparação entre o modelo com as duas variáveis e modelos com cada variável separadamente.

Vamos interpretar o resumo do nosso modelo selecionado, ele aparece no janela inferior de **Output** do Rcmdr quando construimos o modelo. Para solicitar novamente o resumo do modelo insira a linha de comando abaixo na janela superior, selecione essa linha e clique em **Submit**. 

<code>
summary(lmmint)
</code>

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



<WRAP center round box 80%>
Após definir o modelo mínimo adequado, analise o "summary" do modelo e avalie:
  - na estrutura aleatória o valor do desvio padrão da variável aleatória é maior, menor ou equivalente ao desvio padrão dos resíduos?
  - o que essa comparação representa?
  - o que está representado no intercepto da estrutura fixa do modelo?
  - qual a interpretação biológica dos coeficientes dos parâmetros da estrutura fixa do modelo?
</WRAP>

*/


==== Seleção de LMM: riqueza em praias ====

Uma questão central em planejamento e análise de dados é decidir como compor um conjunto de modelos plausíveis a partir das hipóteses e um conjunto de variáveis. Além disso, a partir deste conjunto de modelos plausíveis, é preciso tomar a decisão sobre qual((ou quais)) é o melhor para representar nossos dados. Nesse curso utilizamos a partição da variação para comparar modelos aninhados. Vamos continuar a usar o mesmo procedimento, ou seja, fazer o teste de hipótese para comparar modelos dois a dois a partir do modelo cheio para chegar ao mínimo adequado. Esse procedimento se torna mais complexo nos modelos mistos ao envolver a simplificação da estrutura aleatória, além da estrutura fixa, como fizemos nos modelos anteriores.  


Existem diversas possibilidades de modelos com esses mesmos dados. Por exemplo, podemos adicionar o nível de exposição de cada praia (''Exposure'') como uma variável preditora fixa categórica e/ou podemos inserir o efeito aleatório da variável ''Beach'' na inclinação da variável ''NAP'' (efeito fixo). Essa inclusão da interação entre a variável aleatória (''Beach'') e a preditora fixa contínua (''NAP'') faz com que cada praia possa também ter uma relação específica com ''NAP'' (inclinações de reta diferentes para cada praia), mas sempre ponderadas pela relação média.



<WRAP center round todo 80%>

=== Ajustando os dados:  riqueza em praias ===


Os dados que estamos usando estão disponíveis no [[https://github.com/jmadin/zuur/blob/master/data/RIKZ.txt|github do livro do Zuur et al. (2009)]]. A variável ''Exposure'' originalmente tem ''3'' níveis: ''8'', ''10'' e ''11''. Como o nível ''8'' foi observado apenas em uma praia, reclassificamos esta praia para o nível seguinte, no caso o ''10'' (Zuur et al. 2009).  


A partir da variável ''Exposure'', criamos a ''fExp'' que contempla apenas os valores **10** e **11** para designar os dois níveis de exposição das praias. Para tornar essa variável mais explicita com relação ao seu significado vamos transformar os valores **10** e **11** de ''fExp'' em fatores **low** e **high**, nesta ordem, criando uma nova variável ''fExposure'':

  * abra o menu ''Data'' > ''Manage variables in active data set'' > ''Convert numeric variables to factors...'';
  * coloque o nome da nova variável como ''fExposure'' no campo ''<same as variable>'';
  * deixe selecionada a opção ''Supply level names'';
  * no quadro que irá se abrir digite ''low'' para o valor ''10'' e ''high'' para o valor ''11'';
  * após criar verifique se a variável foi corretamente criada clicando no botão ''View data set''.




</WRAP>


==== Seleção dos efeitos aleatórios: riqueza em praias ====

Existem modelos e, portanto, perguntas e delineamentos amostrais que requerem apenas um efeito aleatório para indicar o agrupamento dos dados. Entretanto, como colocado na seção anterior, há também modelos que podem incluir mais de um efeito aleatório. Esse é o caso da interação entre ''praia'' e ''NAP'' mencionada acima. 

/*
Se olharmos para o primeiro gráfico que fizemos para esse conjunto de dados ((no início desse tutorial)), veremos que cada praia parece ter seu próprio intercepto e inclinação, o que torna plausível pensarmos que o modelo com a interação se ajuste melhor a esses dados.
*/

Para conjunto de modelos plausíveis com diferentes efeito aleatório, Zuur et al. (2009) sugerem um protocolo para a escolha da melhor estrutura de efeitos aleatórios, que vamos seguir neste curso. Os passos principais deste procedimento são listados a seguir. 

<WRAP center round box 80%>

I. Incluir no modelo cheio todos os efeitos fixos de interesse, incluindo termos de interação e qualquer ou estrutura não aleatória, no caso teremos os efeitos fixos:
  * ''Richness ~ fExposure + NAP + fExposure:NAP''

II. associar a esta estrutura fixa cheia os efeitos aleatórios plausíveis em sua estrutura mais complexa. No nosso exemplo, temos  duas possibilidades de efeito aleatório: 
  * apenas variações no intercepto entre praias (''(1|Beach)''); 
  * contendo a interação entre praia e NAP, resultando __também__ na variação de inclinação entre praias (''(NAP|Beach)'');
  
<WRAP center round important 80%>
Os modelos mistos serão ajustados utilizando uma forma diferente de estimar os parâmetros. Ao invés da máxima verossimilhança (//ML - Maximum Likelihood//) que utilizamos nos modelos lineares os modelos mistos usam a **máxima verossimilhança restrita** (//REML - Restricted Maximum Likelihood//). Isso acontece porque a ''ML'' é enviesada para as estimativas de variância e a ''REML'' corrige esse viés. Na prática, nós não precisamos fazer nada, pois a função ''lmer'' do pacote ''lme4''utiliza a ''REML'' (argumento ''REML = TRUE'') como padrão.

</WRAP>


III. utilizar ''anova'' para comparar modelos com estruturas aleatórias aninhadas e mesma estrutura fixa cheia, buscando a simplificação da estrutura aleatória. Para que a comparação seja feita com o ''REML'' é preciso usar o argumento da função ''anova'' como ''refit = FALSE'';

IV. Mantenha a estrutura aleatória selecionada e siga simplificando a estrutura fixa, comparando modelos aninhados por ''anova''. Utilize o padrão da ''anova'' com o argumento ''refit = TRUE''.

</WRAP>

<WRAP center round info 80%>
Poderíamos também ajustar um modelo sem efeito aleatório (usando a função ''lm'') e ver se a inclusão do efeito aleatório resulta em um melhor ajuste do modelo, entretanto, a dependência entre as amostras de uma mesma praia são parte do desenho experimental dos dados e acreditamos que essa característica deve ser mantida no modelo. Além disso, é importante lembrar que a retirada de todos os efeitos aleatórios superestima os graus de liberdade do modelo, ou seja o modelo linear estaria inflado no seu poder de explicação da variabilidade dos dados. 
</WRAP>



A sequência do procedimento de simplificação da estrutura aleatória para o nosso exemplo, seria:


  * Criar  o modelo ''lmmfull'' com todos os efeitos aleatórios, pelo menu do Rcmdr ''Statistic'' > ''Linear mixed models...''  deixando o argumento com o padrão ''REML=TRUE'' e inserindo os termos do modelo como indicado abaixo:

<code>
 Richness ~ fExposure + NAP + fExposure:NAP + (NAP|Beach)
</code>

  * Em seguida faça o mesmo procedimento para criar o modelo ''lmmint'' com a variável aleatória somente no intercepto:

<code>
 Richness ~ fExposure + NAP + fExposure:NAP + (1|Beach)
</code>

  * Compare os modelos anteriores com o comando ''anova'' buscando a simplificação do modelo ao mínimo adequado((note que estes modelos são aninhados, a estrutura aleatória no inclinação ''NAP'' incorpora o intercepto também, apesar de não estar explicito na nossa formulação)), usando o menu do Rcmdr ''Models > Hipothesis tests > Compare two models...'' e selecionando os modelos ''lmmint'' e ''lmmfull''. 

<WRAP center round important 90%>
Um problema com a interface do Rcmdr é que precisamos incluir o argumento ''refit = FALSE'' na função ''anova'' para garantir que o procedimento irá manter a estimativa por ''REML'' e não irá reajustar os modelos por ''ML'' antes de comparar os modelos. O padrão da função ''anova'' quando aplicada à comparação de dois modelos __mistos__ é reajustar para ''ML''. Como a interface não oferece essa opção precisamos editar a linha de comando do painel superior. 

Quando você clicou ''ok'' na função ''anova'' de comparação de modelos, apareceu esse código na janela superior do Rcmdr:

<code>
anova(lmmint, lmmfull)
</code>

copie esse código novamente no painel superior, mas agora inclua a expressão ''refit = FALSE'', como indicado abaixo:

<code>
anova(lmmint, lmmfull, refit= FALSE)

</code>

Selecione essa linha e clique no botão ''Submit''. O resultado deve ser como o apresentado abaixo:

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

</WRAP>


Note que as estimativas do primeiro resultado (ou seja, da anova com reajuste para ML) e do segundo resultado (ou seja, sem reajuste para ML, usando ''refit = FALSE'') são diferentes. É importante saber que **no processo de simplificação da __estrutura aleatória__ do modelo** devemos ignorar os primeiros resultados, pois __não queremos que seja feito um reajuste para ML__, ou seja, **queremos que a comparação das estruturas aleatórias seja feita usando as estimativas por REML**. 

O "p-valor" observado nessa comparação é alto, indicando que os modelos não apresentam diferenças marcantes, sendo assim, devemos reter o modelo mais simples e seguimos simplificando. Como só há mais um termo na estrutura aleatória modelando o intercepto, podemos finalizar a seleção da estrutura aleatória aqui. 

<WRAP center round info 80%>

Alguns autores ((inclusive o próprio Zuur do capítulo que indicamos nesse tutorial como leitura)) advogam que se nenhum termo da estrutura aleatória for significativo, deve-se abandonar o modelo misto e partir para o modelo linear sem estrutura aleatória. De fato, isso facilitaria muito a interpretação e a apresentação dos resultados do modelo. Porém, outros autores que seguiremos aqui, indicam que devemos manter a coerência do delineamento experimental/amostral e portanto devemos contemplar a dependência das observações de uma mesma praia, utilizando o modelo com a variável aleatória  ''Beach'' no intercepto.   
</WRAP>

/*

<code>
anova(m1, m0, refit= FALSE)

</code>


Aqui, por outro lado, o valor do ''p'' é baixo e retemos o modelo mais complexo, ou seja aquele com o variável ''praia'' no efeito aleatório do intercepto.

*/

Após simplificar a estrutura de efeito aleatório o resultado é que as praias diferem na riqueza média, mas essa variação das praias não influencia a relação de ''Richness'' com o ''NAP''.  

Note que quando observamos o primeiro gráfico desse tutorial, parecia que a inclusão de diferentes inclinações para cada praia no modelo seria importante. Porém, o resultado desse teste acima nos indica que, apesar de algumas praias terem inclinações aparentemente diferentes, não ganhamos muita informação adicional relevante ao incluir inclinações diferentes para todas as nove praias.

==== Estrutura fixa do modelo: riqueza de praias ====


O próximo passo é prosseguir com a verificação e simplificação dos efeitos fixos através do mesmo método que utilizamos até então. Relembrando, estamos interessados em verificar se existe um efeito de ''NAP'',  ''fExposure'' e da interação ''NAP:fExposure'' na riqueza de espécies. Começamos com o nosso modelo com todos os efeitos fixos e a estrutura aleatória já simplificada: 

<code>
 Richness ~ fExposure + NAP + fExposure:NAP + (1|Beach)
</code>

Primeiro iremos simplificar retirando o termo fixo mais complexo, a interação ''fExposure:NAP'':

<code>
Richness ~ fExposure + NAP +  (1 | Beach)
</code>

Caso o modelo mais simples tenha um poder de explicação da variação dos dados similar ao mais complexo, devemos continuar a simplificação com os modelos só com ''NAP'' ou só com ''fExposure''. Nos dois casos iremos comparar com o modelo que contém ambos termos ''fExposure + NAP'',  para que os modelos comparados sejam aninhados ((não é possível comparar o modelo só com ''NAP'' com o modelo só com ''fExposure''!)). 

Siga a simplificação até chegar ao modelo mínimo adequado para esse conjunto de dados. 







<WRAP center round box 80%>

==== Diagnósticos dos modelos ====

Depois de selecionado o modelo que melhor se ajusta aos nossos dados é sempre necessário avaliar o ajuste deste modelo a suas premissas. O ''Rcmdr'' faz alguns gráficos diagnósticos e de efeito do modelo que podem ser explorados no menu: ''Models'' > ''Graph''.


Para a construção de outros gráficos diagnósticos podemos carregar o pacote ''lattice'', selecionando-o pelo menu ''Tools'' > ''Load package(s)...''.


Depois, copie os códigos abaixo no painel superior (RScript) do Rcmdr, mas, antes de apertar o botão ''Submit'' **substitua o argumento NOME_DO_MODELO pelo nome do seu modelo mínimo adequado((caso não tenha nomeado os modelos, precisa encontrar o nome que o ''Rcmdr'' designou ao modelo))**. Depois disso pode selecionar e submeter as linhas de comando.

**Atenção: Faça um gráfico de cada vez e salve o ''pdf'' de cada gráfico antes de começar a fazer o próximo**

<code>
## 1) gráfico quantil-quantil (normalidade)
qqmath(NOME_DO_MODELO,id=0.05) 
</code>

<code>
## 2) gráfico de valores ajustados x resíduos (homocedasticidade)
plot(NOME_DO_MODELO,type=c("p","smooth")) 
</code>

<code>
## 3) gráfico de valores ajustados x resíduos padronizados (homocedasticidade)
plot(NOME_DO_MODELO,sqrt(abs(resid(.)))~fitted(.), type=c("p","smooth")) 
</code>



<WRAP center round alert 90%>
**OPS!** Olhando o gráfico diagnóstico dos resíduos, parece que os dados não são tão homocedásticos como deveriam, pois vemos algo parecido com um funil se abrindo da esquerda para a direita, o que indica que o modelo viola esta premissa. Isso não deveria ser um
=== Ajustando os dados ===


Os dados que estamos usando estão disponíveis no [[https://github.com/jmadin/zuur/blob/master/data/RIKZ.txt|github do livro do Zuur et al. (2009)]]. A variável ''Exposure'' originalmente tem ''3'' níveis: ''8'', ''10'' e ''11''. Como o nível ''8'' foi observado apenas em uma praia, reclassificamos esta praia para o nível seguinte, no caso o ''10'' (Zuur et al. 2009).  


A partir da variável ''Exposure'', criamos a ''fExp'' que contempla apenas os valores **10** e **11** para designar os dois níveis de exposição das praias. Para tornar essa variável mais explicita com relação ao seu significado vamos transformar os valores **10** e **11** de ''fExp'' em fatores **low** e **high**, nesta ordem, criando uma nova variável ''fExposure'':

  * abra o menu ''Data'' > ''Manage variables in active data set'' > ''Convert numeric variables to factors...'';
  * coloque o nome da nova variável como ''fExposure'' no campo ''<same as variable>'';
  * deixe selecionada a opção ''Supply level names'';
  * no quadro que irá se abrir digite ''low'' para o valor ''10'' e ''high'' para o valor ''11'';
  * após criar verifique se a variável foi corretamente criada clicando no botão ''View data set''.



a surpresa já que a variável ''Richness'' é uma contagem e como tal tem a variância acoplada à média.

**E agora??!** Bem, quase sempre existe um caminho! O problema aqui é que assumimos que a riqueza de espécies, uma variável de contagem, poderia ser modelada como uma distribuição normal. Entretanto, dados de contagem geralmente são melhor modelados usando a distribuição de ''Poisson''. Já aprendemos isso no tutorial [[cursos:planeco:roteiro:10-glmpoisson]].  

Então, para fazermos a modelagem correta dos nossos dados teremos que usar um modelo com distribuição de ''Poisson'' dentro do contexto de modelo misto. Nesse caso teríamos que usar uma classe de modelos que junta o GLM com o LMM que se chama **Modelo Generalizado Misto**.   

Mas isso é tema para outro roteiro... [[cursos:planeco:roteiro:12-glmm|]]
</WRAP>

</WRAP>
\\
\\







/*
Para isso usamos a função `plotresid` do pacote `RVAidememoire`:

{{:cursos:planeco:planeco:roteiro:lmm3.png?600|}}

<code>
library(RVAideMemoire)

plotresid(m1, shapiro = T)
</code>

Olhando o gráfico da direita (o qqplot), temos que os valores extremos dos dados não se comportam muito bem como uma distribuição normal (as linhas em vermelho no gráfico indicam a área em que os pontos deveriam estar para que pudéssemos considerar os resíduos como normalmente distribuídos). Isso fica evidente quando fazemos o teste de Shapiro e encontramos que os dados são significativamente diferentes de uma distribuição normal (nesse teste se dá significativo é que não é normal).
*/

<WRAP center round help 90%>

**PARA ENTREGAR ANTES DA PRÓXIMA AULA**
Acesse o [[https://forms.gle/AEk1zdTyYzXu2MnV7|formulário]] e responda às questões propostas referentes aos dados sobre a riqueza em praias. 


</WRAP>



===== Cultivando aveia =====

Nesse experimento clássico descrito por Yates em 1935, foram utilizadas três variedades de aveia tratadas com 4 níveis de adubação com nitrogênio. O desenho experimental deste estudo((''split-plot'')) é bastante complexo e contempla diferentes níveis de dependência das observações que precisam ser incorporadas no modelo como variáveis aleatórias.  Vamos entender esse estudo.

O primeiro passo é fazer a leitura dos dados do objeto ''oats'' que se encontra no pacote ''MASS'':

  
<WRAP center round box 80%>

- Carregue o pacote ''MASS'' através do menu ''Tools'' > ''Load packages...''

{{  :cursos:planeco:roteiro:loadMass.png?700  |}}


- Entre no menu ''Data'' > ''Data in packages'' > ''Read data set on attached package...'' 

- Selecione na janela ''Package'': ''MASS'' e na  ''Data set'':  ''oats''

{{  :cursos:planeco:roteiro:oatsData.png?700  |}}



</WRAP>
 

==== Entendendo Aveia ====

 
Para entender os dados utilize o menu ''Help'' > ''Help on active data set (if available)'' e leia a documentação que irá se abrir onde  há uma descrição das variáveis:

<WRAP center round box 80%>

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

O desenho experimental contempla ''6'' localidades (''B'' de //Blocks//) distintas, sendo que cada uma  apresenta três campos de cultivo (''plot'') para cada variedade de aveia (''V'' de //Varieties//).  Cada campo, por sua vez, é dividido em quatro partes (''subplot'') com os diferentes tratamentos com adubo nitrogenado (''N'' de //Nitrogen treatment//). No desenho experimental temos um aninhamento de ''bloco'' com ''plot'' e ''subplot''. Para auxiliar no entendimento deste desenho experimental, preparamos a figura abaixo:


{{ :cursos:planeco:roteiro:Figura Split Plot oats.png?700 |}}
</WRAP>



==== Ajustando os dados de Aveia ====


As variáveis ''plot'' e ''subplot'' não estão nesse conjunto de dados, mas podem ser resgatadas pelo fato de ''plot'' ser o campo onde foi cultivada cada variedade. Só precisaremos criar a variável ''plot'', pois a subdivisão dos campos de cultivo nos diferentes níveis de adubação (que seriam os subplots) não precisa ser indicada no modelo, uma vez que essa divisão representa o menor nível hierárquico e a função ''lmer'' já reconhece esse nível automaticamente. 
 
Para criar ''plot'' vamos usar o menu ''Data'' > ''Manage variable in active data set'' > ''Compute new variable'' e utilize o nome ''P'' (plot) no campo ''New variable name'' e no campo ''Expression to compute'' insira o código abaixo:

<code>
rep(paste("plot", 1:3, sep=""), each=4)
</code>

Esse código irá criar uma nova variável ''P'' com os níveis ''plot1'', ''plot2'' e ''plot3'' para cada campo de cultivo de cada variedade de aveia, nos quais os níveis de adubação estão incluídos (4 níveis). Apesar da variável ''P'' estar no mesmo nível da variável ''V'', fazemos isso para explicitar que ''P'' é diferente de ''V'', sendo que uma representa os campos de cultivo e a outra representa as variedades de aveia cultivadas, respectivamente. Note também que os códigos de ''P''  se repetem entre os blocos, apesar de serem unidades distintas, ou seja, o campo de cultivo ''plot1'' do bloco I é outro campo de cultivo do que o ''plot1'' do bloco II. Por outro lado, as variedades  se repetem entre os blocos, ou seja, a variedade ''Marvellous'' é ''Marvellous'' em qualquer bloco. 

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



<WRAP center round box 80%>
 
==== Atividade ====

  * Construa um modelo cheio tendo como variável resposta o rendimento de cultivo ('' Y '') e como preditoras as variáveis fixas variedade ('' V ''), adubação nitrogenada ('' N '') e a interação entre as duas. Como estrutura aleatória utilize o fator aleatório plot ('' P '') aninhado dentro de bloco ('' B ''), com a sintaxe  ''(1|B/P)''. Note que, como não temos variáveis fixas contínuas, apenas categóricas, não há como modelar a inclinação do modelo, apenas o intercepto.
  * Faça a simplificação da estrutura fixa das preditoras e chegue ao modelo mínimo adequado.
  * Faça a interpretação biológica do resultado baseado nos coeficientes fixos e da variação associada aos fatores aleatórios.


</WRAP>

Se você se interessou pelos modelos mistos e acha que eles se encaixam no seu problema, não deixe de conferir as referências listadas abaixo para se aprofundar. Os modelos mistos são muito flexíveis e uma ferramenta poderosa para comtemplar muitas das dependencias que temos em desenhos experimentais na biologia. 
\\
\\




===== Exercício=====
<WRAP center round help 90%>

**PARA ENTREGAR ANTES DA PRÓXIMA AULA**
Acesse o [[https://forms.gle/T9LpETUubgpq8QN48|formulário]] e responda às questões propostas referentes aos dois exemplos (Praias e Aveia) trabalhados nesse tutorial. 

</WRAP>


/*
===== Glossário =====

<WRAP center round tip 100%>


**Efeitos fixos e aleatórios**

Existe **muita** discussão na literatura sobre como se diferencia efeitos fixos de aleatórios (veja sugestões de leitura no final do roteiro - principalmente McGill 2015). De maneira geral, efeitos fixos são constantes em toda sua amostra, não variam de amostra pra amostra. Além disso, os diferentes níveis existentes no efeitos fixos não variam se você incluir mais amostras (Bates et al 2014). Um exemplo claro é sexo: normalmente sua amostra conterá machos e fêmeas, e aumentar o seu tamanho amostral a quantidade de níveis em seus fatores permanecerá constante. Já os efeitos aleatórios variam de amostra para amostra, por exemplo as praias que foram amostrados no dados que analisamos. Os pesquisadores poderiam ter ido a outras diferentes praias para coletar dados. 


**Máxima verossimilhança**

Máxima verossimilhança (ML) é uma abordagem estatística que estima os parâmetros de um modelo a partir de um dado conjunto de dados. 

Nos modelos mistos normais (que assume distribuição normal dos resíduos), utiliza-se a máxima verossimilhança restrita (REML) para estimar os parâmetros, pois a ML enviesa as estimativas de variância do modelo. 

**Modelo linear**

Modelos lineares descrevem a resposta de uma variável dependente, a que você está interessado em explorar, em função de uma variável preditora, explanatória ou independente.  


Em ecologia, um dos usos mais comuns de modelos lineares é o modelo de regressão, por exemplo. Para fazer um modelo linear no R usamos a função `lm`:

<code>
lm(dependente ~ preditora, data = seus dados)
</code>


**Modelo aninhado**

Um modelo mais específico (mais simples) está aninhado dentro de outro modelo mais geral (mais complexo) de modo que os termos preditores (variáveis preditoras e suas interações) do modelo mais específico formam um subconjunto dos termos do modelo mais geral. 

**Fatores aninhados ou cruzados**

Duas variáveis categóricas (fatores) podem ter efeitos aninhados ou cruzados. No cruzado, todos os níveis de um fator são encontrados nos níveis do outro fator. No caso do Por exemplo, quando temos tratamentos com diferentes adubos (3 níveis) e irrigação (sim ou não) em plantações, o efeito cruzado é ter replicas de todos os 3 tipos de adubação em condição de irrigação e não irrigação. Nos fatores aninhados os níveis de uma variável só são encontrados em um nível do outro fator. Por exemplo, nivel uma medidade crescimento de uma árvore está associada a um fator que é espécie que no caso pode estar aninhada ao fator gênero ou família taxonômica. Não há como ter a espécie replicada em todos os gêneros nesse caso.

</WRAP>

*/

===== Referências e recomendações =====


<WRAP center round box 100%>
Bates, et al. 2014.[[http://arxiv.org/abs/1406.5823|Fitting linear mixed-effects models using lme4. arXiv preprint arXiv:1406.5823.]] (publicação do pacote `lme4`)

Burnham, K. & Anderson, D. 2002.[[http://gen.lib.rus.ec/book/index.php?md5=0572C2F65088CFA05EC3757297DBC173| Model selection and multimodel inference: a practical information-theoretic approach]]. 2nd edn. New York: Springer-Verlag. (Livro sobre a abordagem de seleção de modelos baseada em Teoria da Informação)

McGill, B. 2015.[[https://dynamicecology.wordpress.com/2015/11/04/is-it-a-fixed-or-random-effect/| Is it a fixed or random effect?]] Blog Dynamic Ecology. (Uma boa discussão sobre o que são efeitos fixos e aleatórios)

Winter, B. 2013.[[http://arxiv.org/pdf/1308.5499.pdf| Linear models and linear mixed effects models in R with linguistic applications. arXiv:1308.5499.]] (Nesse excelente roteiro, o autor explica modelos lineares e depois apresenta modelos mistos de uma forma bem didática)

Yates, F. 1935.[[https://www.jstor.org/stable/2983638| Complex experiments. Journal of the Royal Statistical Society Suppl. 2, 181-247]]. (O experimento de aveia que tratamos nesse roteiro é descrito nesse artigo).

Zuur, A., Ieno, E., Walker, N., Saveliev, A. & Smith, G. 2009.[[http://gen.lib.rus.ec/search.php?req=mixed+effect+models+and+extensions+with+r&lg_topic=libgen&open=0&view=simple&res=25&phrase=1&column=def| Mixed effects models and extensions in ecology with R.]] (Livro muito bom e completo sobre modelos mistos e aditivos)

</WRAP>
Laboratório de Ecologia de Florestas Tropicais

Ferramentas do usuário

Ferramentas do site

Ferramentas da página
