“Fazendo sentido de modelos mixtos”

Bem, entre aspas o titulo, que só estou copiando aqui o que vi em alguns blogs essas semana.
Eu olho sempre um agregador de blogs sobre R, chamado Rbloggers . E esta semana teve duas postagens que achei muito legal. Olhem o original aqui e aqui.

Mas trata-se de um exemplo muito interessante que da sentido aos modelos mixtos.
O exemplo é o seguinte, imagine nove indivíduos medidos sobre cinco tratamentos, podemos pensar em nove aranhas e quanto tempo elas demoram para manipular uma presa em flores de cinco espécies de plantas diferentes e perguntamos “Há uma diferença entre o tempo de manipulação de presa entre espécies de plantas?”.

Então fomos la no campo virtual (código no final do post), medimos o tempo que cada aranha demora para manipular uma presa em cada florzinha diferente.

Primeiro olhamos os dados em um gráfico.

A primeira coisa que vemos pelo gráfico é que parece haver uma diferença clara entre cada espécie de planta. Mas além disso, observem que no gráfico, a gente diferencia perfeitamente cada indivíduo de aranha.

Mas como assim?

Olhem que apesar da diferença dada pela espécie de planta (eixo x), existe uma aranha muito rápida e cada aranha vai ficando mais lenta (esta na ordem do 1 ao 9 o melhor pro pior individuo porque simulamos assim, são as cores, indivíduos estão ligados por uma linha para facilitar a visualização). Mas então esperamos que quando ajustar-mos o modelo, veremos uma grande variação devido ao indivíduo e pouca variação devido a planta.

E ajustando nosso modelo mixto temos:

> summary(m3) Linear mixed model fit by REML Formula: size ~ levs - 1 + (1 | ind) Data: idf AIC BIC logLik deviance REMLdev 86.43 99.08 -36.21 63.82 72.43 Random effects: Groups Name Variance Std.Dev. ind (Intercept) 7.634760 2.76311 Residual 0.079001 0.28107 Number of obs: 45, groups: ind, 9 Fixed effects: Estimate Std. Error t value levsl1 6.1632 0.9258 6.657 levsl2 15.9373 0.9258 17.215 levsl3 2.0436 0.9258 2.207 levsl4 9.8413 0.9258 10.630 levsl5 13.0373 0.9258 14.083 Correlation of Fixed Effects: levsl1 levsl2 levsl3 levsl4 levsl2 0.990 levsl3 0.990 0.990 levsl4 0.990 0.990 0.990 levsl5 0.990 0.990 0.990 0.990

Mas vamos olhar essa saída por parte.
Primeiro o efeito fixo:

Fixed effects: Estimate Std. Error t value levsl1 6.1632 0.9258 6.657 levsl2 15.9373 0.9258 17.215 levsl3 2.0436 0.9258 2.207 levsl4 9.8413 0.9258 10.630 levsl5 13.0373 0.9258 14.083

E olha ae, a estimativa do tempo médio para cada espécie de planta, notem que temos uma diferença, não ter diferença seria as estimativas serem iguais (valores próximos, não idênticos). Lembrando que isto esta assim pois no modelo foi declarado “levs – 1”, esse -1 retirou o intercepto, então ai são as médias, e não um intercepto e as diferenças a partir deste. Mas tudo bem essa é a parte fixa e vemos diferença como no gráfico, senão ele seria uma reta, isso já responde nossa pergunta. Mas agora vamos a parte interessante, o efeito aleatório.

Random effects: Groups Name Variance Std.Dev. ind (Intercept) 7.634760 2.76311 Residual 0.079001 0.28107

Olha ae que coisa linda, temos a maior parte da variação devido aos indivíduos, e so um pouquinho de variação devido a espécie de planta, lembrando que é justamente o valor que usamos na simulação, variação de 8 (na verdade do -4 a 4), e é exatamente isso que vemos no gráfico, existe aranhas muito boas, podemos identificar fácil isso no gráfico, aranhas ruins, e elas variam de acordo com a planta, mas se não houve essa variação devido ao indivíduo, todo mundo seria bem igual, olha a variação residual, que não depende do indivíduo, é quase zero.

Podemos ainda fazer o seguinte, retirar o efeito fixo, ou seja o efeito da planta, e olhar somente a variação entre aranhas.
Como eu disse, se não tivesse a variação devido a planta, esperaríamos uma reta certo, então ao tirar a variação da planta temos…

Agora podemos ainda fazer o seguinte, podemos pergunta como seria se essa variação se ela fosse devido a planta e não ao indivíduo de aranha? Primeira coisa que vem a cabeça é que, não deveríamos conseguir identificar muito bem indivíduos bons ou indivíduos ruins no gráfico certo? Pegamos e misturamos os valores de cada planta, então na planta 1, misturamos os indivíduos, logo cada valor é agora atribuído aleatoriamente para cada individuo.
Vamos olhar o gráfico e temos:

Olha só, agora não da para ver aranhas boas ou aranhas ruins. Então olhando esse gráfico ainda esperamos uma diferença entre as plantas, só olhar a primeira e a segunda, são diferentes, mas não parece que a variação esta associada a aranha.
Então vamos la a ajustamos um modelo mixto e temos:

> summary(m4) Linear mixed model fit by REML Formula: size ~ levs - 1 + (1 | ind) Data: idf_rand AIC BIC logLik deviance REMLdev 220.2 232.9 -103.1 214.3 206.2 Random effects: Groups Name Variance Std.Dev. ind (Intercept) 0.0000 0.0000 Residual 7.7138 2.7774 Number of obs: 45, groups: ind, 9 Fixed effects: Estimate Std. Error t value levsl1 6.1632 0.9258 6.657 levsl2 15.9373 0.9258 17.215 levsl3 2.0436 0.9258 2.207 levsl4 9.8413 0.9258 10.630 levsl5 13.0373 0.9258 14.082 Correlation of Fixed Effects: levsl1 levsl2 levsl3 levsl4 levsl2 0.000 levsl3 0.000 0.000 levsl4 0.000 0.000 0.000 levsl5 0.000 0.000 0.000 0.000

Como fizemos anteriormente, primeiro olhamos os fatores fixos e temos:

Fixed effects: Estimate Std. Error t value levsl1 6.1632 0.9258 6.657 levsl2 15.9373 0.9258 17.215 levsl3 2.0436 0.9258 2.207 levsl4 9.8413 0.9258 10.630 levsl5 13.0373 0.9258 14.082

Ainda exatamente a mesma coisa, diferença entre plantas. Mas agora ao olhar o fator aleatório, ou a variação devido a indivíduos vemos:

Random effects: Groups Name Variance Std.Dev. ind (Intercept) 0.0000 0.0000 Residual 7.7138 2.7774

Olha ai que lindo, a variação devido a indivíduos é zero, mas agora temos uma grande variação residual, ou seja a variação vem toda da planta, e é isso que esperamos certo, porque agora no gráfico você vê uma variação na planta, mas ela não é devido ao individuo, você não identifica direito indivíduos diferentes no gráfico. E a variação que era de 8 continua toda ae.

Eu não tenho muita certeza quanto as contas, ou como os ajustes são feitos, olhando o código da função eu acho bem complexo, mas o modelo parece que recupera perfeitamente tudo que foi simulado, e funciona lindamente. E acompanhando exemplos como esse ajudam a ver como modelos mixtos funcionam além de melhorar nossa capacidade de interpretar os nossos resultados.

#########################
#Modelo mixto comparação#
########################
library(lme4)
library(ggplot2)
set.seed(13)
 
#Criando os niveis do fator
levs <- as.factor(c("l1","l2","l3","l4","l5"))
 
#Definindo as médias por fator
#Compare esses valores com os valores fixos estimados
f_means <- c(6,16,2,10,13)
 
#Colocando individuos como um fator
 
ind <- as.factor(paste("i",1:9,sep=""))
 
#Efeito do individuo
#note a escolha da variação de -4 a 4 e compare esse
#intervalo com o fator aleatorio depois
 
i_eff <- seq(-4,4,length=9)
 
#Simulando medidas repetidas para cada individuo
idf <- data.frame(matrix(0,ncol=3,nrow=45))
colnames(idf) <- c("size","ind","levs")
counter <- 1
for(i in 1:length(levs)){
for(j in 1:length(ind)){
idf$size[counter] <- rnorm(1,f_means[i]+i_eff[j],.3)
idf$ind[counter] <- ind[j]
idf$levs[counter] <- levs[i]
counter <- counter + 1
}
}
idf$ind <- rep(ind,5)
idf$levs <- sort(rep(levs,9))
 
#olhando o que foi criado
ggplot(idf,aes(x=levs,y=size,group=ind,colour=ind))+geom_point()+geom_path()
 
#ajustando o modelo
m3 <-lmer(size~levs – 1 +(1|ind), data=idf)
summary(m3)
 
#Vendo efeitos fixos e aleatorios
fixef(m3)
ranef(m3)$ind
 
#retirando o efeito fixo
idf$adjusted <- idf$size – rep(fixef(m3), each = 9)
 
#olhando somente a variação
ggplot(idf, aes(x = levs, y = adjusted, group = ind, colour = ind)) +
geom_point() + geom_path()
 
#variação devido ao individuo
var(unlist(ranef(m3)))
 
#Agora aleatorizando os individuos,
#a variação não esta mais associada ao indiduo
 
idf_rand <- idf
for(i in 1:5){
idf_rand$ind[idf_rand$levs==levs[i]] <- sample(idf$ind[idf$levs==levs[i]],
9,replace=F)
}
 
#Olhando o resultado
 
ggplot(idf_rand,aes(x=levs,y=size,group=ind,colour=ind))+geom_point(
)+geom_path()
 
#ajustando o modelo
m4 <-lmer(size~levs – 1 +(1|ind), data=idf_rand)
summary(m4)
 
#efeitos fixos e aleatorio
fixef(m4)
ranef(m4)
#note como não tem efeito para nenhum individuo

Leave a Reply

Your email address will not be published. Required fields are marked *