Modelos Mistos Generalizados (GLMM) [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 Mistos Generalizados (GLMM) ======

Este roteiro é uma continuação do roteiro de [[cursos:planeco:roteiro:11-lmm|]] (LMM), o qual deve ser feito e entendido bem antes de continuar. Também é importante que você já tenha estudado o roteiro de [[cursos:planeco:roteiro:10-glm|]] (GLMs). Os modelos GLMM são uma extensão dessas duas abordagens reunidas.

Ao final do roteiro de LMMs nós observamos que os resíduos do modelo final não atendiam às premissas de homogeneidade de variâncias e normalidade do resíduos. No próprio roteiro, lembramos que a riqueza de espécies (número de espécies no local) é um dado de contagem (valores sempre inteiros e positivos - discretos) e, portanto, não poderíamos utilizar um modelo que tem como premissa a distribuição normal (valores contínuos). Logo, precisaremos utilizar um modelo com uma distribuição diferente da normal. Neste exemplo, vamos recorrer a um GLM com distribuição de ''Poisson'' - que lida com os dados de contagem, assim como fizemos no roteiro [[cursos:planeco:roteiro:10-glm#glm_contagem| GLM contagem]] .


==== Voltando aos dados da riqueza de espécies das praias ====

Vamos voltar ao exemplo do livro do Zuur et al. (2009), no qual queremos entender se a riqueza de espécies da macro-fauna em 9 praias na costa da Holanda é influenciado pela altura da estação de amostragem em relação à altura média da maré (`NAP`) ((vamos começar com o modelo mais simples)). As unidades amostrais estão aninhadas dentro das praias, ou seja, temos 5 valores por praia que não podem ser consideradas independentes, e por isso incluímos no nosso modelo a praia como um efeito aleatório.

Dados de contagem, ou seja, o número de espécies em cada amostra, são geralmente modelados com a distribuição de **Poisson**, que assume também que a média é igual à variância. Ou seja, se a média aumenta, a variância também. 

No nosso exemplo, essa relação variância-média pode ser visualizada olhando o gráfico de dispersão dos dados:

<WRAP center round box 80%>

{{:cursos:planeco:roteiro:richData.png?600|}}
</WRAP>

Perceba que os pontos vão ficando menos dispersos à medida que os valores de riqueza diminuem (e os de `NAP` aumentam).

==== Modelando GLMMs: contagem ====


Então, mãos à massa! 
==== Lendo os dados no R ====

  * baixe o arquivo {{ :cursos:planeco:roteiro:praia.txt |dados vamos à praia}}
  * no R, mude o diretório de trabalho (pasta) para o diretório onde o dado se encontra
  * leia os dados <wrap hi>praia.txt</wrap> com o seguinte comando ((pode fazer isso pela interface do Rcmdr, mas lembre-se de atribuir o nome **dados** para esse conjunto de dados no Rcmdr))
  * o arquivo é do tipo texto, com colunas separadas por tabulação ("Tabs") e com decimal indicado por ponto (".")
  * crie a variável ''fExposure'' como fator, como fizemos no LMM

<code>

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


</code>

Agora, vamos modificar o modelo do roteiro de LMM apenas em relação à distribuição de probabilidades dos resíduos do modelo. Operacionalmente isso é bastante simples, basta usar a função ''glmer()'' do pacote ''lme4'' com a mesma formulação do modelo usada anteriormente, explicitando a distribuição (''family = "poisson"'').

<code>
library(lme4)
mod.riq <- glmer(Richness ~ NAP +  (1 | Beach), data = dados, 
                family = "poisson")
</code>

Vamos inspecionar o modelo:
<code>
summary(mod.riq)
</code>


==== Predito pelo Modelo ====

Lembrando que, como num GLM, os valores dos coeficientes agora estão na escala da função de ligação, que para a distribuição de <wrap hi>poisson</wrap> é a função <wrap hi>log</wrap>.

Para fazer o gráfico de nosso modelo ajustado, com a <wrap hi>reta média</wrap> (predição) dos efeitos fixos e as <wrap hi>retas do efeito aleatório</wrap>, preditas para cada praia, precisamos primeiro dos valores dos coeficientes do preditor linear. A partir dele podemos estimar qualquer valor predito (esperança) com um certo NAP em uma dada praia. Se tivermos uma sequencia de valores de NAP, em um dado intervalo, podermos construir as curvas preditas e com isso representar o modelo.
O código abaixo cria os seguintes objetos:
 
  * ''randcoef'': guarda os coeficientes da estrutura aleatória 
  * ''fixcoef'': guarda os coeficientes da estrutura fixa do modelo
  *  ''seqNap'': cria uma sequencia de 100 valores no intervalo de -1.5 a 2.5  
  *  beach<wrap hi>XX</wrap>: estima o valore esperado para a praia <wrap hi>XX</wrap> para os valores de ''seqNAP''

 <code>

## coeficientes da estrutura aleatória
randcoef <- coef(mod.riq)$Beach
randcoef
## coeficientes da estrutura fixa do modelo
fixcoef  <- fixef(mod.riq)
fixcoef

## valores da preditora em uma sequencia
seqNAP <- seq(-1.5, 2.5,len=100)

## estimativas para cada praia
beach01 <- randcoef[1,1] + randcoef[1,2] * seqNAP
beach02 <- randcoef[2,1] + randcoef[2,2] * seqNAP
beach03 <- randcoef[3,1] + randcoef[3,2] * seqNAP
beach04 <- randcoef[4,1] + randcoef[4,2] * seqNAP
beach05 <- randcoef[5,1] + randcoef[5,2] * seqNAP
beach06 <- randcoef[6,1] + randcoef[6,2] * seqNAP
beach07 <- randcoef[7,1] + randcoef[7,2] * seqNAP
beach08 <- randcoef[8,1] + randcoef[8,2] * seqNAP
beach09 <- randcoef[9,1] + randcoef[9,2] * seqNAP
beachFix <- fixcoef[1] + fixcoef[2] * seqNAP

</code>

==== Predito na escala da observação ====

Os preditos calculados acima estão na escala do preditor linear, no caso o logaritmo da contagem de espécies. Para colocar os valores na escala de observação precisamos exponenciar. O código abaixo faz isso para cada praia e para o predito pelo estrutura fixa do modelo (a praia média).

<code>

beach01a <- exp(beach01)
beach02a <- exp(beach02)
beach03a <- exp(beach03)
beach04a <- exp(beach04)
beach05a <- exp(beach05) 
beach06a <- exp(beach06)
beach07a <- exp(beach07) 
beach08a <- exp(beach08) 
beach09a <- exp(beach09)
beachFixa <- exp(beachFix)

</code>

==== Gráfico Base ====

Primeiro, vamos fazer o gráfico básico com os dados coletados. O código a seguir  contem as seguintes linhas de códigos:

  * ''cores ...'': cria um objeto com nove cores tiradas de uma paleta do arco-iris com a função ''rainbow(...)''
  * ''par(...)'': função que modifica parâmetros da janela gráfica, no caso modificamos as bordas
  * ''plot(...)'': função que cria o gráfico
      * 1° linha: os dados
      * 2° e 3° linha: modificações de parâmetros do gráfico

<code>

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


</code>

/*

calcular os valores preditos para cada praia de acordo com os valores de ''NAP''. Mas agora temos que nos atentar que queremos fazer a predição das médias de riqueza em cada NAP na escala da variável resposta e não da função de ligação usada no modelo. Operacionalmente isso também é simples, basta colocar o argumento ''type = "response"'' na função ''predict()''((esse valor predito é o que se considera ingênuo (//naive//), quando não estamos incorporando a variabilidade dos efeitos aletórios. Como introdução é válido calcular os valores preditos desta maneira, mas para se aprofundar no tema sugerimos ler Bates et al. (2014) para formas mais apropriados de predição de médias e intervalos de confiança estimados para modelos mistos)) que teremos os valores preditos já na escala da variável resposta (riqueza de espécies) para cada praia. 

*/

==== Predito para cada Praia ====

Agora podemos juntar ao gráfico base as predições para cada praia e o predito pela estrutura fixa do modelo. As primeiras linhas abaixo incluem essas predições em cores diferentes e a última linha insere uma legenda no gráfico.



<code>

lines(beach01a ~ seqNAP, col = cores[1])
lines(beach02a ~ seqNAP, col = cores[2])
lines(beach03a ~ seqNAP, col = cores[3])
lines(beach04a ~ seqNAP, col = cores[4])
lines(beach05a ~ seqNAP, col = cores[5])
lines(beach06a ~ seqNAP, col = cores[6])
lines(beach07a ~ seqNAP, col = cores[7])
lines(beach08a ~ seqNAP, col = cores[8])
lines(beach09a ~ seqNAP, col = cores[9])
lines(beachFixa ~ seqNAP, col = rgb(0,0,0,0.5), lwd=5)
legend("topright", c( "praia média",paste("praia", 1:9)) , lty=1, col= c(1,cores), bty="n", cex=1.5)


</code>

/*

<code>
#predições para cada praia separadamente:
# 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(mod.riq, newdata = novo, type="response")

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


Para construir nossa curva média da predição da riqueza de espécies pelo ''NAP'', use o argumento ''re.form=NA'', assim a função vai ignorar os efeitos aleatórios.

<code>
#predição para todas as praias, ingorando os efeitos aleatórios
novis <- expand.grid(NAP = seq(-1.3,2.2,0.5))
pred2 <- predict(mod.riq, newdata = novis, type="response", re.form=NA)
dados.pred2 <- data.frame(pred = pred2, novis)
</code>

Agora podemos plotar os dados com o ajuste do nosso modelo:

<code>
#usando o pacote ggplot2 e cowplot (estética do gráfico)
#não esqueça de instalar os pacotes antes de carregar 
#Remova o # nas linhas abaixo para instalar:
#install.packages("ggplot2")
#install.packages("cowplot")

library(ggplot2)
library(cowplot)
ggplot(data = dados, aes(x = NAP, y = Richness,   # dados e eixos
                         color = as.factor(Beach))) + 
  geom_point(size = 3, shape = 19) +                   # colocando os pontos
  geom_line(data = dados.pred2, aes(y = pred, x = NAP ), 
            col= "black", size=2) +                    # curva média
  geom_line(data = dados.preditos, aes(y = pred, x = NAP, 
                            col = as.factor(Beach)))   # curva para cada praia
</code>



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

*/

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


Podemos ver nessa figura a predição do nosso modelo em relação aos parâmetros fixos (curva em preto) e as predições para cada praia separadamente.

Compare este gráfico ao gráfico das predições no roteiro de LMM. Qual parece se ajustar melhor aos dados?

Já ajustamos nosso primeiro GLMM! Porém nossas análises têm mais uma variável preditora que poderia entrar no modelo. Então, o próximo passo é incorporar a variável fixa `Exposure` e fazer a escolha dos efeitos aleatórios antes de prosseguir com a seleção de modelos.

===== Efeitos aleatórios =====

Os passos para a escolha dos efeitos aleatórios no GLMM, quando temos mais de uma possibilidade de efeitos aleatórios, são iguais aos apresentados no roteiro de LMM. 


  * 1. Criar o modelo completo, com todos os efeitos fixos a serem utilizados. Aqui incluímos também a variável ''fExposure'' e a interação entre ''fExposure'' e ''NAP''. Nosso modelo é:

//Richness = fExposure * NAP + "efeito(s) aleatório(s)"//


  * 2. Incluir os efeitos aleatórios plausíveis. No nosso exemplo, escolhemos duas possibilidades de efeito aleatório, variações no intercepto entre praias (''(1|Beach)'') e a interação entre praia e ''NAP'', resultando também na variação de inclinação entre praias (''(NAP|Beach)''). Além disso, colocamos para "competir" um modelo glm sem efeitos aleatórios.

<code>
m0 <- glm(Richness ~ fExposure*NAP, data = dados, family= poisson)
m1 <- glmer(Richness ~ fExposure * NAP + (1|Beach), data = dados,
            family = poisson)
m2 <- glmer(Richness ~ fExposure * NAP + (NAP|Beach), data = dados,
            family = poisson)
</code>


Um ponto importante, assim como nos LMM, é que estes modelos serão ajustados utilizando uma forma diferente de estimação, ao invés da máxima verossimilhança (ML), vamos usar a **máxima verossimilhança restrita** (REML). Isso acontece porque a ML é enviesada para as estimativas de variância do modelo e a REML corrige este enviesamento. Na prática, nós não precisamos fazer nada, pois a função ''glmer'' já usa, por padrão, a REML (argumento ''REML = TRUE'').

  * 3. Utilizar o comando ''anova'', como temos feito para tomar a decisão sobre qual modelo aninhado reter, buscando a simplificação do modelo ao mínimo adequado. Aqui utilizamos o argumento ''refit = FALSE'' para garantir que a função irá manter a estimativa do ''REML''.

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

</code>

Aqui o valor de ''p'' é relativamente baixo, indicando que os modelos apresentam diferenças marcantes, sendo assim, devemos reter o modelo mais complexo, ou seja com a estrutura de erro ''(NAP|Beach)'' .

Bom, agora sabemos que a melhor estrutura de efeito aleatório é a variação na relação de inclinação de ''NAP'' entre as praias. Assim, podemos prosseguir com a verificação dos efeitos fixos através de teste de hipóteses e/ou seleção de modelos.





==== Inferência por seleção de modelos pela tabela de ANOVA ====

Agora, vamos avaliar os modelos mistos quanto à sua estrutura fixa através da comparação de modelos por ANOVA, da mesma maneira que foi feito no roteiro de LMM.

Comparando modelo com interação entre ''fExposure'' (lembrando que estamos usando como variável categórica com 3 níveis) e ''NAP'' e sem interação para ver se a interação é importante:

<code>
# modelo com interação entre fExposure e NAP
m3 <- glmer(Richness ~ fExposure * NAP + (NAP|Beach), 
            data = dados, family="poisson")

# modelo sem interação entre fExposure e NAP
m4 <- glmer(Richness ~ fExposure + NAP + (NAP|Beach), 
            data = dados, family="poisson")
</code>

E aplicamos a função ''anova'' nos nossos modelos aninhados:
<code>
anova(m3,m4)
</code>

Como o resultado da comparação entre os modelos não foi significativo, ou seja, o modelo com interação não é significativamente diferente do modelo sem interação, nós ficamos então com o modelo sem interação e continuamos a compará-lo com modelos com apenas uma das variáveis preditoras por vez:

<code>
# removendo a variável NAP
m5 <- glmer(Richness ~ fExposure + (NAP|Beach), 
            data = dados, family="poisson")
anova(m4, m5)
</code>


O resultado acima nos mostra que não parece razoável remover a variável ''NAP'' do modelo, pois há uma diferença significativa entre os modelos com e sem ''NAP''. 

Agora, vamos comparar novamente o modelo com as duas variáveis preditoras com um modelo sem a variável ''fExposure''.

<code>
# removendo a variável `fExposure`
m6<- glmer(Richness ~ NAP + (NAP|Beach), 
            data = dados, family="poisson")
anova(m4, m6)
</code>

Também não parece razoável remover ''fExposure'' e ficar apenas com ''NAP'', então nosso modelo final contém apenas o efeito aditivo das duas variáveis. Para finalizar, vamos comparar este modelo com o nulo.

<code>
# modelo nulo
m7 <- glmer(Richness ~ 1 + (NAP|Beach), data = dados, family="poisson")

anova(m4,m7)
</code>

Então podemos concluir que tanto ''fExposure'' quanto ''NAP'' influenciam na riqueza de espécies, mas de maneira aditiva.

Vamos observar o //output// deste modelo:

<code>
summary(m4)
</code>



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

A forma de se fazer o diagnóstico de um GLMM é ligeiramente diferente de um LMM, isso porque não esperamos normalidade nem homocedasticidade dos dados, e mesmo de um GLM, já que residuos escalonados como o resíduo deviance e de Pearson podem não informar bem se o modelo está mal especificado. Sugiro fortemente olhar [[https://theoreticalecology.wordpress.com/2016/08/28/dharma-an-r-package-for-residual-diagnostics-of-glmms/|este roteiro do pacote `DHARMa`]] que explica bem os problemas de usar técnicas de diagnóstico usuais em LMM e GLM para um GLMM e traz soluções bem interessantes!

Nós vamos seguir o diagnóstico proposto por Florian Hartig usando o pacote `DHARMa`. Essa abordagem é baseada em simulações para criar resíduos escalonados, com valores entre 0 e 1, interpretáveis de um modelo GLMM. Aqui, os detalhes estatísticos e teóricos da forma de realizar estas simulações e criar os resíduos escalonados vai ficar de lado (mas está nas referências do pacote `DHARMa`) para que possamos focar em como olhar os gráficos produzidos e interpretá-los. 

=== Criando os resíduos ===

Os resíduos escalonados (também chamados de resíduos quantílicos) são calculados com a função `simulateResiduals()`, na qual indicamos o nosso modelo GLMM e o número de simulações. Para o número de simulações, o padrão da função é 250 pois pode demora muito, mas para uma análise concreta faça pelo menos 1000. Como os nossos dados são poucos e o modelo simples, vamos utilizar 1000 (mas se se computador demorar, pare o processo e faça com menos simulações).

<code>
#install.packages("DHARMa") #não esquecer de instalar o pacote
library(DHARMa)

residuo <- simulateResiduals(fittedModel = m4, n = 1000)
</code>


=== Plotando os resíduos ===

Vamos fazer a inspeção visual dos nossos resíduos:

<code>
plotSimulatedResiduals(residuo)
</code>

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

À esquerda, temos um gráfico de quantil-quantil, que nos ajuda a ver se os nossos resíduos estão de acordo com o que é esperado dado o nosso modelo. Portanto, nossos pontos deveriam estar todos alinhados em cima da linha vermelha. Se não estiverem, aí temos problemas! A forma desse não-alinhamento nos indica quais são os possíveis problemas de especificação do modelo. No nosso modelo, não parece haver desvios muito fortes dos valores esperados.

À direita, temos um gráfico dos resíduos contra os valores preditos pelo modelo. Este gráfico nos ajuda a detectar desvios da uniformidade na direção y. As linhas vermelhas são regressões quantílicas que mostram os quantis 0.25, 0.50, 0.75. Estas linhas deveriam ser retas horizontais para cada quantil. **Note** que alguns desvios desta reta são esperados ao acaso, mesmo para um modelo perfeito, especialmente quando o tamanho amostral é pequeno. No nosso caso, temos um valor predito muito alto que puxa as retas vermelhas pra baixo. Vamos avaliar melhor a uniformidade dos resíduos a seguir com o teste de uniformidade.


/*

=== Testes de ajustes para os resíduos ===

O teste de uniformidade nos ajuda a interpretar melhor os resultados do gráfico da direita mostrado anteriormente. Este teste roda um teste de Kolmogorov-Smirnov para testar a uniformidade dos resíduos quantílicos.

<code>
testUniformity(residuo)
</code>

Este teste nos mostra que os nossos resíduos simulados podem ser considerados uniforme, mesmo que as retas vermelhas do gráfico dos valores dos resíduos x preditos não sejam linhas horizontais porque estão sendo influenciadas por um único ponto que teve uma predição bem alta.


**OBS:** Como o objetivo deste roteiro é apenas uma breve introdução aos GLMMs, nos limitamos a poucas análises de diagnóstico. Com o pacote `DHARMa`, existe a possibilidade de se explorar:

- A sobredispersão nos resíduos (que indicaria que uma distribuição de Poisson já não é adequada - havendo a possibilidade de usar uma distribuição binomial negativa que lida bem com dados sobredispersos). 

- Presença de dados inflados em zero - que é outro problema para as distribuições Poisson e Binomial.

- A autocorrelação espacial e temporal dos resíduos.

*/

==== Predito pelo Modelo Selecionado ====

Finalmente, vamos plotar um gráfico com os dados e as curvas de predição para o modelo que foi selecionado - efeito aditivo de `Exposure` e `NAP`. Vamos seguir os passos anteriores para calcular os valores preditos para cada praia e para a estrutura fixa do modelo. Note que agora temos duas curvas para cada praia, representando as exposições alta e baixa.

<code>

rand04 <- coef(m04)$Beach
rand04
fix04<- fixef(m04)
fix04
beach01L <- exp(rand04[1,1] + rand04[1,3] * seqNAP)
beach01H <- exp(rand04[1,1] +  rand04[1,2]  + rand04[1,3] * seqNAP)
beach02L <- exp(rand04[2,1] + rand04[2,3] * seqNAP)
beach02H <- exp(rand04[2,1] +  rand04[2,2]  + rand04[2,3] * seqNAP)
beach03L <- exp(rand04[3,1] + rand04[3,3] * seqNAP)
beach03H <- exp(rand04[3,1] +  rand04[3,2]  + rand04[3,3] * seqNAP)
beach04L <- exp(rand04[4,1] + rand04[4,3] * seqNAP)
beach04H <- exp(rand04[4,1] +  rand04[4,2]  + rand04[4,3] * seqNAP)
beach05L <- exp(rand04[5,1] + rand04[5,3] * seqNAP)
beach05H <- exp(rand04[5,1] +  rand04[5,2]  + rand04[5,3] * seqNAP)
beach06L <- exp(rand04[6,1] + rand04[6,3] * seqNAP)
beach06H <- exp(rand04[6,1] +  rand04[6,2]  + rand04[6,3] * seqNAP)
beach07L <- exp(rand04[7,1] + rand04[7,3] * seqNAP)
beach07H <- exp(rand04[7,1] +  rand04[7,2]  + rand04[7,3] * seqNAP)
beach08L <- exp(rand04[8,1] + rand04[8,3] * seqNAP)
beach08H <- exp(rand04[8,1] +  rand04[8,2]  + rand04[8,3] * seqNAP)
beach09L <- exp(rand04[9,1] + rand04[9,3] * seqNAP)
beach09H <- exp(rand04[9,1] +  rand04[9,2]  + rand04[9,3] * seqNAP)
beachFixL <- exp(fix04[1] + fix04[3] * seqNAP)
beachFixH <- exp(fix04[1] + fix04[2]  + fix04[3] * seqNAP)

</code>



==== Gráfico do Modelo ====


Para fazer o gráfico, seguimos os mesmos passos do código do gráfico anterior, com a diferença que agora a exposição alta é representada por uma linha tracejada(''lty = 2'').


<code>

par( mar=c(5,5,2,2))
plot(Richness ~ NAP ,data = dados,
     pch = 19, col = cores[as.factor(Beach)] , las = 1,
     cex=1.5, cex.lab= 1.7, cex.axis = 1.5
     )
lines(beach01L ~ seqNAP, col= cores[1], lty=1)
lines(beach01H ~ seqNAP, col= cores[1], lty=2)
lines(beach02L ~ seqNAP, col= cores[2], lty=1)
lines(beach02H ~ seqNAP, col= cores[2], lty=2)
lines(beach03L ~ seqNAP, col= cores[3], lty=1)
lines(beach03H ~ seqNAP, col= cores[3], lty=2)
lines(beach04L ~ seqNAP, col= cores[4], lty=1)
lines(beach04H ~ seqNAP, col= cores[4], lty=2)
lines(beach05L ~ seqNAP, col= cores[5], lty=1)
lines(beach05H ~ seqNAP, col= cores[5], lty=2)
lines(beach06L ~ seqNAP, col= cores[6], lty=1)
lines(beach06H ~ seqNAP, col= cores[6], lty=2)
lines(beach07L ~ seqNAP, col= cores[7], lty=1)
lines(beach07H ~ seqNAP, col= cores[7], lty=2)
lines(beach08L ~ seqNAP, col= cores[8], lty=1)
lines(beach08H ~ seqNAP, col= cores[8], lty=2)
lines(beach09L ~ seqNAP, col= cores[9], lty=1)
lines(beach09H ~ seqNAP, col= cores[9], lty=2)
lines(beachFixL ~ seqNAP, col= rgb(0,0,0,0.3), lty=1, lwd=5)
lines(beachFixH ~ seqNAP, col= rgb(0,0,0,0.3), lty=2, lwd=5)

legend("topright", c( "praia média",paste("praia", 1:9)) , pch =19,  col= c(rgb(0,0,0,0.5),cores), bty="n", cex=1)

</code>


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


/*
<code>
#predições para cada praia:
# 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))
novo$fExposure <- rep(c("low", "low", "high", "high", "low",
                        "high", "high", "low", "low"),8)

#agora calculamos os valores preditos para esse novo conjunto de dados
preditos <- predict(m4, newdata = novo, type="response")

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

<code>
#predição para todas as praias, ingorando os efeitos aleatórios
novis <- expand.grid(NAP = seq(-1.3,2.2,0.5),
                     fExposure = unique(dados$fExposure))
pred2 <- predict(m4, newdata = novis, type="response", re.form=NA)
dados.pred2 <- data.frame(pred = pred2, novis)
</code>

Agora podemos plotar os dados com o ajuste do nosso modelo:

<code>
ggplot(data = dados, aes(x = NAP, y = Richness,        # dados e eixos
                         color = as.factor(Beach))) + 
  geom_point(size = 3, shape = 19) +                   # colocando os pontos
  geom_line(data = dados.pred2, aes(y = pred, x = NAP, 
                  color = fExposure), show.legend = F, size=1.5) +                    # curva média
  geom_line(data = dados.preditos, aes(y = pred, x = NAP, 
                            col = as.factor(Beach)))  + # curva para cada praia
 theme(legend.position ="none")
</code>

{{:cursos:planeco:planeco:roteiro:glmm-figurafinal.png?600|}}


Repare que agora as curvas para cada praia estão agrupadas ao redor das curvas de cada `fExposure` (ou seja, categorias "low" e "high"), justamente porque cada praia é classificada com uma categoria de `fExposure`. 

*/


==== Referências recomendadas ====

[[http://arxiv.org/abs/1406.5823|Bates, et al. 2014.]].**Fitting linear mixed-effects models using lme4**. arXiv preprint arXiv:1406.5823. (publicação do pacote `lme4`)

Florian Hartig (2016). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level/Mixed) Regression	Models.	R package. version 0.1.2. https://github.com/florianhartig/DHARMa

Zuur, A., Ieno, E., Walker, N., Saveliev, A. & Smith, G. 2009. Mixed effects models and extensions in ecology with R. (Livro muito bom e completo sobre modelos mistos e aditivos)