Modelos Lineares Mistos [Laboratório de Ecologia de Florestas Tropicais]
cursos:planeco:roteiro:11-lmm
Esse roteiro é baseado em material desenvolvido por alunos de nosso programa 1)Melina Leite, Marília Gaiarsa e Lucas Medeiros e vem sendo modificado, ao longo dos anos, para adaptá-lo ao curso.
Para acompanhar o roteiro é importante ter uma compreensão básica de análises estatísticas partimos da premissa que você já sabe interpretar o resumo de um modelo linear usando as funções como lm() e summary(). Como os modelos mistos ainda não foram incorporados na interface do Rcdmr, será preciso rodar no R os códigos que aparecem em caixas de textos nesse roteiro.
cat("Isso é um código do R")
Para isso, cole as linhas de texto na janela superior do Rcmdr e clique no botão Submite (ou Submeter) para executar a linha em que o cursor do mouse se encontra, ou selecione um conjunto de linhas para executá-las juntas. É possível rodar os códigos a partir de um console básico do R também. Caso queira usar um escript direto no R veja o tutorial Introdução ao R
O objetivo desse roteiro e que você seja capaz de:
Para auxiliar a entender terminologias estatísticas específicas, ou caso queira relembrar alguns termos, ao final do roteiro encontra-se um pequeno glossário.
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 anterior feito seja pedido, sugerimos fortemente que instale os seguintes pacotes antes de abrir o Rcmdr, ou logo no inicio da sessão. Para isso, copie e cole as seguintes linhas de comando:
pkgs <- c("lme4", "RVAideMemoire")
pkgsNot <- pkgs[!(pkgs %in% rownames(installed.packages()))]
install.packages(pkgsNot)
Para explicar o que são modelos mistos e sua importância em ecologia, usaremos como exemplo um conjunto de dados presente no capítulo 5 de Zuur et al. (2009). Esses dados contêm a riqueza de espécies da macro-fauna 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 variável NAP influencia a riqueza de espécies nessas praias, deixando de lado a variável de exposição por hora. Podemos, por exemplo, construir um modelo linear da seguinte forma:
Aqui estamos modelando a riqueza de cada ponto de amostra em função da variável NAP. O coeficiente $\alpha$ é o intercepto do modelo, $\beta$ é o coeficiente angular (inclinação) de NAP e $\varepsilon$ é o erro do nosso modelo (resíduos), isto é, a variação da riqueza que não conseguimos explicar com nossa variável preditora NAP. No entanto, esse modelo tem um problema sério. Estamos violando uma premissa fundamental de modelos lineares, a de que os dados são independentes uns dos outros (Winter, 2013). Os dados obtidos em uma mesma praia não são independentes entre si (dependência espacial). Podemos imaginar diversas características de cada praia que podem influenciar a riqueza de espécies, como o tipo de grão de areia ou a força das ondas, que não estão contempladas em nosso modelo e podem fazer com que as amostras de uma mesma praia sejam mais parecidas, independente das nossas variáveis de interesse ( no caso apenas NAP). Veja o gráfico abaixo, em que cada cor representa uma praia, para se convencer de que a relação entre riqueza e NAP varia entre praias:
Nesta figura, temos as retas do modelo entre Riqueza e NAP aplicado para cada praia (usando a função lm()). Entretanto, nossa pergunta inicial não diz respeito a cada praia separadamente, nós queremos modelar esta relação para todas as praias amostradas, sem especificar a relação para cada praia. A inclusão de unidades amostrais que apresentam dependência entre elas é chamada de pseudo-replicação ou falta de independência das observações. Queremos considerar essa dependência sem falar especificamente de nenhuma praia e, ao mesmo tempo, ser possível ter alguma ideia da variação que existem entre as praias, inclusive as que não foram amostradas.
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 nestas 9 praias, mas poderiam ter feito o mesmo estudo escolhendo outras praias. O efeito aleatório praia organiza parte da variação nos nossos dados que não conseguimos explicar, presente no erro $\varepsilon$ do modelo simples descrito no início desse roteiro (Winter, 2013). Por outro lado, as variáveis preditoras que estamos acostumados a encontrar em modelos lineares são chamadas de aqui de efeitos fixos. No exemplo das praias, a variável NAP é um efeito fixo e estamos diretamente interessados em seu efeito sobre a variável resposta. O nome misto refere-se ao fato de que existirem efeitos fixos e aleatórios no modelo. Esses modelos são chamados também de hierárquico ou multinível.
Na seção a seguir veremos como explorar modelos mistos no R. Veremos que 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).
Nosso modelo será:
$$b_i ~ N(0, \sigma_{praia})$$
$$ \varepsilon ~ N(0, \sigma_{res})$$
A riqueza da praia i é explicada pelo efeito fixo $\beta*NAP$ e pelo efeito aleatório $b_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 $b_i$ e, portanto, o modelo permite que cada uma das praias apresente um intercepto diferente, mas mesmas inclinações $\beta$. O $\beta$ faz parte do efeito fixo, e não possuí nenhum índice i atrelado a ele, portanto é o mesmo para todas as praias. Finalmente, $\varepsilon$ representa o resíduo associado ao desvio do valor observados em relação ao predito pela reta da praia relativa à observação.
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. A seguir, vamos colocar a mão na massa e analisar os dados provenientes do capítulo 5 de Zuur et al. (2009) que exploramos na seção anterior.
Antes de tudo, precisamos baixar e instalar o pacote lme4:
##caso não tenha instalado o pacote tire o comentário da linha seguinte
Os dados estão disponíveis no site do livro do Zuur et al. (2009)2), mas podem ser baixados no link abaixo.
Aqui vamos apenas transformar a variável fExp em fator com dois níveis (low, high), nesta ordem, criando uma nova variável fExposure
Agora, vamos construir nosso modelo. Relembrando, estamos interessados em ver se existe um efeito de NAP na riqueza de espécies. Nossos dados contam com nove diferentes praias e cinco diferentes amostras para cada uma delas. Dessa forma, temos como efeito fixo a variável NAP, e como efeito aleatório a variável praia (Beach'), que significa que cada praia terá um intercepto diferente. Assim, nosso modelo é:
A função summary aplicada a um objeto de modelos, retorna as informações mais relevantes do mesmo.
Vamos dar uma olhada na parte dos efeitos aleatórios (Random effects). O desvio padrão (Std.Dev.) é uma medida do quanto a variabilidade da nossa variável dependente - riqueza - é devida aos dois efeitos aleatórios que estamos analisando (os interceptos das praias e os resíduos). Podemos ver o desvio padrão associado às diferenças de intercepto entre praias (o $b_i$ do modelo). A última linha nos dá o resíduo, que indica o quanto da variabilidade não é prevista pela variável aleatória praia nem pelo NAP, que nada mais é do que o $\varepsilon$ acima explicado.
Os efeitos fixos indicam os coeficientes estimado pra cada um dos fatores que estamos considerando como fixos. No caso, temos o intercepto (medio de todas as praias) que é a riqueza quando o NAP é zero (6.5819), e o coeficiente angular de NAP (-2.5684), o que indica que há uma diminuição de riqueza da ordem de mais de 2.5 espécies a cada unidade de NAP acrescida.
fixLMM <- fixef(modelo.riqueza)
Além desses coeficientes, os modelos mistos também fazem a estimativa aproximada para cada um dos níveis da variável aleatória e que pode ser acessado da saguinte forma:
randLMM <- coef(modelo.riqueza)
Note que a inclinação é a mesma para todas as praias e o intercepto médio dos coeficientes da estrutura aleatória do modelo é o intercepto da estrutura fixa.
mean(randMM[,1])
Para fazer o gráfico de nosso 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 primeiro calcular os valores preditos para cada praia.
O predito pelo modelo pode ser calculado a partir dos coeficientes. Como a intenção é construir uma curva, utilizamos uma sequência de valores da preditora (variável no eixo x, no caso NAP) e estimamos valores em intervalos curtos para podermos construir a curva do modelo.
beach01 <- randLMM[1,1] + randLMM[1,2] * seqNAP
beach02 <- randLMM[2,1] + randLMM[2,2] * seqNAP
beach03 <- randLMM[3,1] + randLMM[3,2] * seqNAP
beach04 <- randLMM[4,1] + randLMM[4,2] * seqNAP
beach05 <- randLMM[5,1] + randLMM[5,2] * seqNAP
beach06 <- randLMM[6,1] + randLMM[6,2] * seqNAP
beach07 <- randLMM[7,1] + randLMM[7,2] * seqNAP
beach08 <- randLMM[8,1] + randLMM[8,2] * seqNAP
beach09 <- randLMM[9,1] + randLMM[9,2] * seqNAP
beachFix <- fixLMM[1] + fixLMM[2] * seqNAP
Vamos primeiro construir um gráfico base, com os valores observados. Para tornar o gráfico mais informativo vamos colocar cores para identificar de que praia a observação foi coletada. No código abaixo temos as seguintes sequencia de eventos:
cores…: cria um objeto com nove cores tiradas de uma paleta das cores do arco-iris
par(…) : modifica parâmetros da janela gráfica, no caso modificamos as bordas
Vamos usar os valores que calculamos para representar as esperanças do modelo e criar as curvas para cada praia e a curva média. No cṕdigo abaixo as primeiras linhas incluem as curvas associadas a cada uma das praias amostradas, a penultima linha inclui o predito pela estrutura fixa do modelo e a última a legenda.
Existem, entretanto, diversas maneiras de criar seu modelo, escolher os parâmetros importantes, e decidir quais são fixos e quais são aleatórios. Nesse exemplo, poderíamos colocar a variável `Exposure` também como um efeito fixo. Outra complicação que poderíamos inserir no nosso modelo é inserir mais um efeito aleatório de interação entre a praia (efeito aleatório) e NAP (efeito fixo), fazendo com que cada praia possa também ter sua relação com a NAP (inclinações de reta diferentes).
Existem modelos e, portanto, perguntas e delineamentos amostrais, que requerem apenas um efeito aleatório para indicar o agrupamento dos dados. Entretanto, como colocado no final da 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, 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 aos nosso dados.
1. A primeira coisa a ser feita é ajustar um modelo com todos os efeitos fixos que estão sendo testados (modelo "completo"). No nosso exemplo tínhamos apenas NAP, mas vamos incluir Exposure (como fator fixo) para exemplificar melhor, e vamos incluir a interação entre NAP e Exposure. No nosso modelo "completo" a estrutura fixa é:
Richness ~ fExposure + NAP + fExposure:NAP + "efeito(s) aleatório(s)"
Podemos usar o símbolo (*) para simplificar a sintaxe, pois esse símbolo significa incluir o efeito sozinho e também suas interações. Escrevendo o mesmo modelo de forma simplificada (usando o asterisco):
Richness ~ fExposure * NAP + "efeito(s) aleatório(s)"
2. Depois que escolhemos o modelo "completo", adicionamos 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)). Podemos 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 num melhor ajuste do modelo:
Um ponto importante é 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 é enviesado 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 `lmer` já usa, por padrão, a REML (argumento REML = TRUE).
3. Vamos utilizar o comando anova como temos feito para fazer a seleção dos modelos aninhados, buscando sempre 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.
Aqui o valor de p é alto indicando que os modelos não apresentam diferenças marcantes, sendo assim, devemos reter o modelo mais simple e seguimos simplificando.
anova(m1, m0, refit= FALSE)
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.
Nessa parte, depois de já escolhermos a estrutura aleatória do nosso modelo, podemos averiguar qual a real influência dos efeitos fixos na riqueza de espécies, utilizando o procedimento que temos adotado nessa disciplina, o Teste de Hipóteses através da comparação de modelos pela tabela de ANOVA. E, depois de selecionado 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.
Nesse caso, precisamos ajustar nosso modelo por máxima verossimilhança (ML), pois é o indicado para compararmos modelos com diferentes efeitos fixos, mas com mesmo efeito aleatório. Então colocamos o argumento REML = FALSE no nosso modelo:
m1 <- lmer(Richness ~ fExposure + NAP + fExposure:NAP + (1|Beach), data = dados, REML = F)
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):
anova(m1, m6)
Sim! 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.
O primeiro diagnóstico do modelo é verificar se os resíduos são normalmente distribuídos, para isso geralmente usamos um qqplot (gráfico de quantil-quantil da distribuição normal). Além disso, podemos também checar visualmente a homogeneidade de variância dos resíduos (homocedasticidade), ou seja, se a variabilidade dos resíduos se mantém constante em relação aos valores ajustados.
OPS! Olhando o gráfico da direita, vemos que os dados não são tão homocedásticos como gostaríamos, vemos algo parecido com um funil se abrindo da esquerda para a direita, o que indica que estamos violando esta premissa. Olhando o gráfico da esquerda (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).
Bates, Douglas; Martin Mächler, Ben Bolker, Steve Walker. 2015. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software. DOI 10.18637/jss.v067.i01. https://www.jstatsoft.org/ article/view/v067i01/.
Bolker, Ben. 215. Linear and Generalized Mixed Models. Chapter 13 in Fox et al. Statistical Ecology. Oxford University Press.
A variável `Exposure` nos dados originais tem 3 níveis, mas porque o valor mais baixo foi observado apenas em uma praia, é necessário reclassificar esta praia para o segundo valor mais baixo (Zuur et al. 2009). Os dados do nosso site já tem isso ajustado