Assista ao vídeo abaixo para um breve introdução sobre os modelos lineares mistos e siga acompanhando pelo tutorial e outras videoaulas ao longo do roteiro.

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

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²⁾. 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.

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

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 nesse curso. Porém, lembre-se que estamos utilizando o pacote lme4 para a construção destes modelos. Então, precisamos carregá-lo quando o Rcmdr é aberto. Para isso, vá em Tools → load packages e escolha lme4⁴⁾.

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⁵⁾:

Richness ~ NAP +  (1 | Beach)

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

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

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⁶⁾ 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⁷⁾. 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 ⁸⁾.

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 ⁹⁾.

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

fixLMM <- fixef(lmm.riq)
fixLMM

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

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

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}$.

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:

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
     )

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:

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])

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”:

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)

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

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.

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¹⁰⁾ é 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.

Os dados que estamos usando estão disponíveis no site 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.

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.

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.

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:

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

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

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

• Compare os modelos anteriores com o comando anova buscando a simplificação do modelo ao mínimo adequado¹¹⁾, usando o menu do Rcmdr Models > Hipothesis tests > Compare two models… e selecionando os modelos lmmint e lmmfull.

anova(lmmint, lmmfull)

anova(lmmint, lmmfull, refit= FALSE)

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.

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.

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:

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

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

Richness ~ fExposure + NAP +  (1 | Beach)

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 ¹³⁾.

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

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

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

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

Vamos nos debruçar um pouco mais nos modelos mistos olhando outros exemplos, buscando sempre entender o conceito de variáveis fixas e aleatórias e treinando a interpretação do resultado dos modelos

O interesse deste estudo foi entender se a deciduidade (perda de folha síncrona) de árvores está relacionada ao crescimento do indivíduo. Os pesquisadores não estavam interessados se espécies distintas tem crescimentos diferentes, entretanto a perda de folha síncrona é característica da espécie e eles precisaram incluir essa variável no desenho experimental. No estudo em questão amostraram 5 espécies decíduas e outras 5 perenes (não perdem as folhas no mesmo período do ano) e, para cada espécie, anotaram o crescimento de '10 árvores. Como indivíduos de uma mesma espécie são mais parecidos entre sí do que indivíduos de outras espécies, essa observações não são independentes e precisam ser contempladas no modelo. Um outro complicador é que árvores de diferentes tamanhos ou idades tem crescimentos diferentes. Por isso inclui-se no estudo o tamanho da árvore como uma preditora fixa e com isso os pesquisadores condicionaram o crescimento à variável tamanho.

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¹⁵⁾ é bastante complexo e contempla diferentes níveis de dependência das observações que precisam ser contempladas 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:

Para entender os dados utilize o menu Help > Help on active data set (if available). Leia da documentação que irá se abrir onde há uma descrição dos dados:

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:

rep(paste("plot", 1:3, sep=""), each=4)

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.

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.