Data Augmentation para dados de Captura e Recaptura com efeito do tempo como uma variável aleatória

A algum tempo atrás, a gente falou de data augmentation aqui e depois usando data augmentation num modelo com efeito do tempo aqui.

Muito legal, mas particularmente no modelo com efeito do tempo, a gente cogitou (no livro do Marc Kery é cogitado na verdade) a possibilidade de tratar esse efeito como uma variável aleatória, então vamos tentar fazer isso e uma pequena modificação depois para ver o que acontece.

Bem os dados serão os mesmo que para o post de data augumentation com efeito do tempo. Mas vamos mudar o modelo. Na verdade vamos mudar um pouco a geração dos dados, vamos pegar o efeito de tempo de uma distribuição normal, para condizer com nossa analise.

Mas no mais é tudo igual, agora o modelo na linguagem Bugs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()

Basicamente ele continua o mesmo, só mudamos a parte do prior la em cima.

1
2
3
4
5
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo

Veja que agora temos ainda a geração dos efeitos do tempo, de cada tempo de amostragem, mas agora pegamos esses valores de uma distribuição normal, e não de uma distribuição uniforme como da outra vez.

Mas tentamos determinar agora, qual a média e o desvio padrão dessa distribuição.
Para a média partimos de uma distribuição normal com média zero e desvio 0.001, ou seja, um pico no zero, e para o desvio padrão saímos de uma distribuição uniforme de 0 a 5, lembre-se que o desvio padrão tem que ser um número maior que zero, e não pode ser negativo, senão o bugs vai reclamar.

Uma outra coisa da linguagem bugs é que a distribuição normal usa precisão, e não o desvio diretamente na função. Isso não tem nenhum impacto na estatística, tem haver com precisão numérica apenas eu acho. Até hoje não entendi esse conceito perfeitamente, mas isso não muda nada.

1
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)

No fim, tau vai ser 1 dividido pela variância, que é o desvio padrão ao quadrado, veja que o desvio só vai aumentar, e um sobre um valor entre zero e alguma coisa, vai ser cada vez mais um valor pequenininho, e bem é isso.

Fazemos umas cadeias maiores, adicionamos na função de geradores de valores iniciais nossos novos priors

1
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),sigma.tempo = runif(1,0, 5))

E é isso, agora é so rodar o modelo.

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 5000 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 13500 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.825 0.995 99.000 99.000 100.000 100.000 102.000 1.001 7100 p[1] 0.656 0.047 0.562 0.625 0.658 0.689 0.743 1.001 14000 p[2] 0.140 0.034 0.080 0.115 0.138 0.162 0.212 1.002 3000 p[3] 0.070 0.025 0.029 0.051 0.067 0.085 0.124 1.002 2400 p[4] 0.559 0.049 0.464 0.525 0.559 0.593 0.653 1.001 14000 p[5] 0.931 0.027 0.869 0.916 0.935 0.951 0.974 1.002 2000 med.tempo 0.470 0.306 -0.123 0.325 0.469 0.612 1.073 1.001 14000 sigma.tempo 0.579 0.381 0.236 0.359 0.471 0.663 1.608 1.001 14000 omega 0.402 0.031 0.341 0.381 0.402 0.423 0.464 1.001 5100 deviance 436.569 10.169 424.600 428.000 435.400 441.800 461.100 1.001 4200 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = var(deviance)/2) pD = 51.7 and DIC = 488.3 DIC is an estimate of expected predictive error (lower deviance is better). >

E apesar da estimativa do tamanho populacional estar ok, as estimativas do tempo ficam bem ruins, confira os valores que geraram os dados com os que foram estimados.

> variacao_Tempo [1] 1.517343 -1.058690 -1.568697 1.357698 3.792000

Veja que a estimativa do desvio padrão ficou bem abaixo do real usado, que era 2. Ela ficou algo entre zero e 1.6, enquanto o desvio real era 2. A média ficou correta quanto ao usado para gerar os dados, englobando o zero, mas o desvio deixou a desejar.

Se a gente olhar a figura de das estimativas do tamanho da população.

01

Veja que a gente perdeu bastante em variabilidade, já que o tempo variou menos, porque ficou “preso” dentro do nosso desvio padrão subestimado.

Mas ainda existe uma possibilidade. Lembre-se que aqui temos poucas medidas de de tempo, apenas 5. Lembra de alguma distribuição que tentamos tratar do número de amostras para controlar o desvio da média? Se pensou na distribuição t, pensou certo.

Podemos modelar a variável aleatória que gera o tempo como uma distribuição t, isso vai adicionar um parâmetro a mais no modelo, mas no entanto, teremos como ter a possibilidade de ter os valores mais longe da média como valores mais comuns para o valor do tempo. Ja que se você ler o post sobre a distribuição t, vera que isso que o novo parâmetro no t faz, o grau de liberdade, quanto mais graus de liberdade, mais parecido com uma distribuição normal, mas quanto menos, deixando a distribuição mais “uniforme”, com uma maior chance de valor mais longe da média, e como temos apenas 5 valores, isso pode ser uma propriedade interessante. É meio maluco pensar que vamos estimar esse parâmetro, já que na distribuição t ele é simplesmente o valor de amostras, mas ok, podemos estimar um valor para ele, para tentar ajustar um modelo mais realista.

De todo aquele modelo, vamos mudar apenar a parte dos prior, mais especificamente onde está a variável aleatória que é o tempo.

1
2
3
4
5
6
7
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
        k.tempo ~ dunif(2, 30) # Hyperprior da variacao do tempo
 
        for (i in 1:tempo){
            p[i] ~ dt(med.tempo, tau.tempo, k.tempo)
        }

Veja que usamos o dt ao invés de dnorm dentro do loop, e temos um parâmetro chamado de k.tempo agora, que vamos tirar de uma distribuição uniforme de 2 a 30, quando ele for 30, temos uma distribuição normal, quanto menor, deixamos tudo mais uniforme. Começamos no 2 porque o valor de k tem que ser no mínimo 2 para a linguagem bugs, veja a definição aqui.
E é isso, agora é só colocar ele no inicializador.

1
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),sigma.tempo = runif(1,0, 5),k.tempo= runif(1,2, 30))

E pronto, podemos rodar. Não esquecendo que temos que colocar ele entre os parâmetros que temos que acompanhar.

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 5000 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 13500 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.866 1.010 99.000 99.000 100.000 100.000 102.000 1.001 11000 med.tempo 0.470 0.294 -0.100 0.322 0.469 0.613 1.065 1.001 14000 sigma.tempo 0.558 0.361 0.217 0.345 0.462 0.648 1.519 1.001 5900 k.tempo 16.565 7.884 2.979 9.950 16.660 23.360 29.310 1.001 5100 omega 0.402 0.031 0.343 0.380 0.402 0.423 0.465 1.001 14000 deviance 437.026 10.301 424.700 428.200 435.800 442.900 461.900 1.001 14000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = var(deviance)/2) pD = 53.1 and DIC = 490.1 DIC is an estimate of expected predictive error (lower deviance is better).

E olhando para a figura temos a mesma coisa

02

E no final das contas, o valor de k ficou entre 2 e 30, ou seja, não conseguimos estimar nada para ele, ja que de 2 a 30 era as possibilidades que abrimos. As estimativas de média para a distribuição de parâmetros continua ok, mas o desvio ainda está subestimado. Então apesar dessa ideia de usar a distribuição t seja legal, ela não contribuiu em nada praticamente, aumentamos um parâmetro e não melhoramos muito o modelo. Veja que o deviance continua a mesma coisa. Bom mas valeu a tentativa.

Eu imagino que dado o pequeno número de tempos, essa ideia de usar ele como uma variável aleatória não seja uma boa ideia, ja que estamos trocando 5 parâmetros por 3, ou 2 no primeiro exemplo. Além disso, precisamos lembrar que essa é uma afirmativa forte, pensar que que o tempo é uma variável aleatória, talvez num span pequeno de tempo, pode estar tudo ok, mas num span de tempo maior, dependendo da espécie que estamos trabalhando, podemos pensar o tanto que falamos em aquecimento global, em degradação do ambiente, em tantas mudanças que falar que nada acontece no tempo é algo bem “forte”. Mas apesar de tudo, talvez para mais medidas e um span pequeno de tempo, essa seja uma possibilidade bem legal.

Bem é isso ai, a gente tem que tentar fazer as coisas para aprender, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Kéry e Shaub 2011 – Bayesian population analysis using WinBUGS. A hierarchical perspective – Academic Press 554pp

John K. Kruschke 2014 Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan.2nd Edition. Academic Press / Elsevier.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
library(R2OpenBUGS)
 
##################################
## Gerando dados
##################################
set.seed(51)
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = rnorm(tempo,0,2)
 
yreal<- matrix(NA, nrow=n,ncol=tempo)
 
for (j in 1:tempo){
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}
 
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)
 
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]
 
##################################
## Modelo com tempo
##################################
 
# Data Augumentation
nz <- 150
yaug <- rbind(yobservado, array(0, dim = c(nz, tempo)))
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(yaug = yaug, amostras = nrow(yaug), tempo = ncol(yaug))
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),
                  sigma.tempo = runif(1,0, 5))
 
# Parametros para saida
params <- c("N", "p","med.tempo","sigma.tempo","omega")
 
# Configurações do MCMC
ni <- 5000
nt <- 2
nb <- 500
nc <- 3
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura1
jpeg("01.jpg")
hist(out_tempo$sims.list$N,breaks=seq(90,110,by=1), col = "gray", main = "", xlab = "Tamanho Populacional", las = 1,xlim=c(95,105))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
dev.off()
#####################################
## Efeito aleatorio com distribuiçao t
#####################################
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
 
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
        k.tempo ~ dunif(2, 30) # Hyperprior da variacao do tempo
 
        for (i in 1:tempo){
            p[i] ~ dt(med.tempo, tau.tempo, k.tempo)
        }
 
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),
                         sigma.tempo = runif(1,0, 5),k.tempo= runif(1,2, 30))
inits()
# Parametros para saida
params <- c("N","med.tempo","sigma.tempo","k.tempo","omega")
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura2
jpeg("02.jpg")
hist(out_tempo$sims.list$N,breaks=seq(90,110,by=1), col = "gray", main = "", xlab = "Tamanho Populacional", las = 1,xlim=c(95,105))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
dev.off()

Anova com efeito aleatório Bayesiana

Opa, ja vimos em alguns momentos os modelos mixtos, mas eles são um pouco complicados as vezes. Mas talvez tentar ver parte a parte ajude.

Então vamos voltar um pouco no exemplo anterior quanto medimos plantinhas. Medimos várias plantinhas controle e vamos focar somente nelas. O que temos que lembrar é que caso não tenhamos medidos todas as plantinhas do mesmo lugar, vindas da mesma população, teremos uma variação devido a espécie e outra devido a população.

Vamos lembrar do conceito de herdabilidade. Se tinha um lugar for colonizado por uma plantinha grande, esperamos que os descendentes sejam grandes também, logo toda aquela população será de plantinhas grandes. Agora se um outro lugar for colonizado por uma planta pequenininha, todos os descendentes serão pequenos, só lembrar que se nossos pais são altos, seus filhos são altos também, mas não exatamente a mesma altura, existe uma variabilidade dentro da família.

Aqui o problema começa é que se você fizer medidas em uma população e em outra população, devem haver diferenças e qualquer outro experimento, por exemplo o que fizemos no post anterior, a gente vai ter problemas se pegar uma população para um nível do tratamento e outra população para outro nível do tratamento, porque a gente vai se perder no que é diferença entre populações e o que é diferença no tratamento.

Agora outra coisa seria fazer os níveis do tratamento nas duas populações, mas ai a gente vai inflar a variação, porque vão ter diferenças devido a população e o tratamento.

A ultima idéia é estimar um parâmetro para cada população, a média de tamanho dela, mas ai quanto mais populações, mais parâmetros, e o pior é que são mais parâmetros que não temos interesse.
Mas aqui a gente pode tentar separar o que é a variação da população e o que é a espécie. Bem se tudo esta ficando confuso, vamos fazer umas coletas “virtuais” e olhar um gráfico como ficam.

01

Bem essa figura é o seguinte, nos medimos 10 plantinhas em 15 populações, aqui pensando em população como um lugar. Então fomos em 15 lugares e cada lugar mediamos 10 plantinhas diferentes, as 10 plantinhas estão resumidas em um boxplot. Aqui a gente pode ver que existem algumas populações como a 1, 2, 14 e 15 que são bem extremas, ou muito de plantas muito grandes ou muito pequenas, mas em média todo mundo é bem parecido, vide o as populações ali no meio do gráfico.

Então podemos supor que a população 1 pode estar assim porque um indivíduo grandão colonizou esse lugar, mas os tamanhos médios são mais abundantes em todas as populações, logo eles colonizarão mais lugares e assim eles são mais comuns.

Agora vamos pensar isso em termos do modelo. Primeiro tentamos ilustrar com uma figura.

ranova

Então vamos debaixo para cima, temos uma plantinha medida, essa plantinha veio de uma distribuição com uma média e um desvio, notem que até aqui a figura é idêntica a figura do post anterior, agora adicionamos mais um andar, um hyperparametro, que é a distribuição que faz as medidas da população, e essa distribuição, é a variação da espécies.

Então fazemos um modelo para passar para o BUGS assim:

for (i in 1:npop){
pop.med[i] ~ dnorm(med.sp, tau.sp) # Prior para media de cada populacao
effe[i] <- pop.med[i] – med.sp # Efeito da populacao (efeito derivado) (para comparar com lme4)
}
med.sp ~ dnorm(0,0.001) # Hyperprior da media geral
sigma.sp ~ dunif(0, 10) # Hyperprior da variacao da sp
sigma.res ~ dunif(0, 10) # Prior da variacao da pop
 
# Likelihood
for (i in 1:n) {
y[i] ~ dnorm(mean[i], tau.res)
mean[i] <- pop.med[x[i]]
}
 
# Quantidades derivadas
tau.sp <- 1 / (sigma.sp * sigma.sp)
tau.res <- 1 / (sigma.res * sigma.res)

Notem que adicionamos essa parte:

for (i in 1:npop){
pop.med[i] ~ dnorm(med.sp, tau.sp) # Prior para media de cada populacao
}
med.sp ~ dnorm(0,0.001) # Hyperprior da media geral
sigma.sp ~ dunif(0, 10) # Hyperprior da variacao da sp

Que é, fazemos um loop para criar uma média para cada população, esses valores vem de uma média e um desvio que são a média e os desvio da espécie, ai para cada população vemos os valores que cada população gera.

Dai ajeitamos os dados para mandar para o Openbugs, fazemos uma função para gerar os valores iniciais, mandamos tudo para o openbugs e ficamos com…

Inference for Bugs model at "re.anova.txt", Current: 3 chains, each with 1200 iterations (first 200 discarded), n.thin = 2 Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff med.sp 2.328 0.288 1.768 2.147 2.326 2.511 2.922 1.001 2600 pop.med[1] 3.732 0.148 3.435 3.633 3.733 3.830 4.022 1.001 3000 pop.med[2] 1.166 0.147 0.871 1.066 1.168 1.263 1.451 1.001 3000 pop.med[3] 1.694 0.150 1.403 1.595 1.690 1.793 1.992 1.001 3000 pop.med[4] 1.856 0.146 1.563 1.761 1.860 1.950 2.137 1.001 3000 pop.med[5] 2.112 0.147 1.826 2.012 2.112 2.210 2.402 1.001 2200 pop.med[6] 2.069 0.149 1.771 1.971 2.072 2.166 2.356 1.003 1000 pop.med[7] 2.859 0.145 2.580 2.762 2.855 2.954 3.143 1.003 1600 pop.med[8] 1.892 0.147 1.599 1.793 1.891 1.988 2.183 1.001 3000 pop.med[9] 2.410 0.149 2.122 2.311 2.407 2.506 2.701 1.000 3000 pop.med[10] 1.571 0.149 1.282 1.470 1.573 1.670 1.864 1.003 990 pop.med[11] 2.249 0.146 1.959 2.154 2.248 2.347 2.539 1.001 3000 pop.med[12] 3.739 0.148 3.448 3.639 3.739 3.838 4.033 1.002 1200 pop.med[13] 2.224 0.144 1.946 2.125 2.224 2.322 2.506 1.001 3000 pop.med[14] 4.484 0.149 4.205 4.381 4.479 4.586 4.784 1.001 3000 pop.med[15] 0.962 0.146 0.677 0.862 0.963 1.059 1.258 1.001 3000 sigma.sp 1.093 0.241 0.740 0.924 1.057 1.214 1.680 1.001 3000 sigma.res 0.468 0.029 0.416 0.448 0.468 0.486 0.531 1.001 3000 deviance 196.697 6.064 187.100 192.300 196.000 200.300 210.500 1.001 3000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 15.7 and DIC = 212.4 DIC is an estimate of expected predictive error (lower deviance is better).

Agora olhem so, estimamos a média para cada espécie, pop.med[1]… até o 15, são a média das medidas de cada população, so olhar para o gráfico e comparar com o valor aqui, assumimos que todas as populações tem o mesmo desvio padrão, que é o sigma.res, e temos as medidas interessantes que são o med.sp e sigma.sp, que são a media e o desvio que geram as medias das populações. vamos tentar ver novamente o gráfico.

02

Bem legal, passamos um risco vermelho na média geral e vemos que a maioria das populações estão perto desse risco, elas tem um desvio igual ao risco azul que são duas vezes o sigma.sp. O risco azul cobre quase (95%) todas as possibilidades de como as populações devem ser.

E agora uma pequena comparação:

Original B.med (Intercept) pop.med[1] 3.705 3.732 3.728 pop.med[2] 1.288 1.166 1.172 pop.med[3] 1.722 1.694 1.697 pop.med[4] 1.880 1.856 1.856 pop.med[5] 1.876 2.112 2.113 pop.med[6] 2.268 2.069 2.069 pop.med[7] 2.727 2.859 2.859 pop.med[8] 2.233 1.892 1.894 pop.med[9] 2.339 2.410 2.410 pop.med[10] 1.448 1.571 1.573 pop.med[11] 2.348 2.249 2.245 pop.med[12] 3.485 3.739 3.741 pop.med[13] 2.188 2.224 2.225 pop.med[14] 4.443 4.484 4.484 pop.med[15] 0.847 0.962 0.961

Podemos ver os valores originais das populações, que usamos para fazer a simulação e comparar com os valores que estimamos através da estatística bayesiana e estatística frequentista. Vemos que ambas estimam relativamente bem os valores.

Bem aqui demos uma olhada na idéia de variáveis aleatória, e como podemos ao invés de sofrer com o fato de que cada população tem características próprias, podemos ajustar o nosso modelo ao que achamos, mas isso sozinho não muda muita coisa. O próximo passo é somar essa possibilidade de modelagem com a análise anterior de anova e ver a diferença entre níveis de tratamento levando em conta a variação das populações.

Referência:
Marc Kéry (2010) Introduction to WinBUGS for Ecologists.

set.seed(321)
 
#Gerando dados de exemplo
#quantidades amostras
npop <- 15
namostra <- 10
n <- npop * namostra
 
#Caracteristicas
med.geral.sp <- 2 #Média geral da especies
sp.sd <- 1 #Variação de tamanhos médio da especie
pop.med <- rnorm(n=npop,mean=med.geral.sp,sd=sp.sd)
sigma <- 0.5 #Variação observada em uma população
eps <- rnorm(n, 0, sigma)
 
x <- rep(1:npop, rep(namostra, npop))
X <- as.matrix(model.matrix(~ as.factor(x)-1))
y <- as.numeric(X %*% as.matrix(pop.med) + eps) # as.numeric is ESSENTIAL
 
#Figura 1
boxplot(y ~ x, col = "grey", xlab = "População", ylab = "Medidas",
main = "",frame.plot=F,xlim=c(0,15),ylim=c(0,6))
 
### Modelo Bayesiano
sink("re.anova.txt")
cat("
model {
 
# Priors and some derived things
for (i in 1:npop){
pop.med[i] ~ dnorm(med.sp, tau.sp) # Prior para media de cada populacao
effe[i] <- pop.med[i] – med.sp # Efeito da populacao (efeito derivado)
}
med.sp ~ dnorm(0,0.001) # Hyperprior da media geral
sigma.sp ~ dunif(0, 10) # Hyperprior da variacao da sp
sigma.res ~ dunif(0, 10) # Prior da variacao da pop
 
# Likelihood
for (i in 1:n) {
y[i] ~ dnorm(mean[i], tau.res)
mean[i] <- pop.med[x[i]]
}
 
# Quantidades derivadas
tau.sp <- 1 / (sigma.sp * sigma.sp)
tau.res <- 1 / (sigma.res * sigma.res)
}
",fill=TRUE)
sink()
 
# Junte os dados numa lista para mandar para o bugs
dados.bugs <- list(y=y, x=x, npop = npop, n = n)
 
# Função para iniciar a cadeia
inits <- function(){
list(med.sp=runif(1,0,10),sigma.sp=runif(1,0,2),sigma.res=runif(1,0,2))
}
 
# Parametros a estimar
parameters <- c("med.sp", "pop.med", "effe", "sigma.sp", "sigma.res")
 
# MCMC settings
ni <- 1200
nb <- 200
nt <- 2
nc <- 3
 
# Start WinBUGS
library(R2OpenBUGS)
out<-bugs(dados.bugs, inits, parameters, "re.anova.txt", n.thin=nt,
n.chains=nc, n.burnin=nb, n.iter=ni)
 
# Inspect estimates
print(out, dig = 3)
 
boxplot(y ~ x, col = "grey", xlab = "População", ylab = "Medidas",
main = "",frame.plot=F,xlim=c(0,15),ylim=c(0,6))
abline(h = med.geral.sp,lwd=2,lty=2,col="red")
points(x=c(0,0),y=c(2+1.96,2-1.96),pch=19,col="blue",type="b",lwd=2,lty=2,
cex=2)
legend("topleft",legend=c("Média Geral","Variação da espécie"),bty="n",
col=c("red","blue"),lty=2,lwd=2,pch=c(NA,19))
 
#Analise lme4
library(lme4)
pop <- as.factor(x)
lme.mod <- lmer(y ~ 1 + 1 | pop, REML = TRUE)
lme.mod
 
round(cbind(Original=pop.med,B.med=out$summary[2:16,1],
lme4=coef(lme.mod)$pop),digits=3)

“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

Introdução aos modelos mixtos e variáveis aleatórias.

Modelos mixtos são bastante comuns hoje em dia, na verdade esse nome me confunde, pois um modelo mixto tem que ter variáveis aleatórias e variáveis fixas, mas eu não sei como é o nome de um modelo somente com variáveis aleatórias. Mas como regressão linear, anova e ancova são modelos lineares, eu vou chamar o exemplo que vou colocar aqui de modelo mixto e espero não estar muito errado.

Mas vamos fazer um exemplo para facilitar o entendimento. Vamos supor que queremos avaliar a relação entre o tamanho corporal de cobras e sua massa. A ideia é que cobras maiores tem uma massa maior. Mas coletamos cobras em vários locais, cada local em que coletamos vamos considerar que são uma população distinta. Então coletamos em 56 populações, 10 cobrinhas em cada população.

Como o interesse principal é a relação entre o tamanho das cobras e sua massa, uma possibilidade seria ignorar a população de origem das cobras medidas.

Mas isso, além dos problemas estatísticos, não faz muito sentido biologicamente falando, porque esperamos que cada população tenha características próprias. Algumas populações podem vir de locais pobres em alimento ou onde a competição é acirrada, o que faz com que todas as cobras sejam pequenas e mesmo quando crescem sejam magrinhas. Um outro extremo seria um local onde o alimento é farto, pouca competição, o que leva a cobras grandes e que crescem e engordam facilmente. Mas esses são os extremos, a maioria das populações teriam um crescimento comum, e um tamanho médio nas cobrinhas. Tudo isso dito, o ponto principal é que cada população tem características próprias, características que não deveriam ser ignoradas se queremos entender o que esta acontecendo. Mas isso leva a um problema que é a limitação das nossas conclusões. Então temos que fazer um modelo que leva em consideração as características próprias de cada população, mas ai qualquer conclusão que tirarmos desses dados deveria se limitar a essas populações. Mas vamos olhar por um momento o que esta acontecendo.

Olhando a figura, parece que a proposta é fazer 56 regressões, isso será um mundo de regressão, mas existe uma solução melhor.
Primeiro existe lugares muito bons para se viver, lugares muitos ruins, e lugares médios. Mas lugares muitos bons ou muito ruins são raros, normalmente os lugares serão meio iguais, e isso se reflete nos parâmetros que estimamos. Se uma regressão linear tem a formula:

massa = intercepto + inclinação * tamanho

Podemos dizer que cada população pode ter seu próprio intercepto e sua própria inclinação. Pensando somente em interceptos, são 56. Mas se se viéssemos a coletar a população 57? Mais uma? Bem o que da para ser feito é imaginar que os 56 interceptos saíram de uma distribuição e que a próxima população também sairá dessa distribuição. Ou seja podemos pensar que essas 56 populações tem uma média e um desvio de valores de intercepto, ou seja dificilmente veremos valores extremos, muitos altos ou muitos baixos no intercepto. A mesma ideia vale para as inclinações.

massa = intercepto[pop. de origem] + inclinação[pop. de origem] * tamanho

Olhando a figura 2, agora podemos tentar ver isso, primeiro pensando no intercepto. A população 9 tem um intercepto bem baixo enquanto a população 1 tem um intercepto alto, todo mundo é mais gordinho. Mas talvez mais fácil olhar as inclinações. Olhem como a população um tem uma inclinação decrescente enquanto a 21 tem uma inclinação muito forte, mas esses são casos isolados, a maioria das inclinações parece igual. Isso nos remete a uma distribuição normal, ou seja esses parâmetros intercepto e inclinação devem estar saindo de uma distribuição normal.

Pensar nisso traz um ganho muito interessante. Que é poder fazer um modelo para as populações que eu estou vendo, e poder até pensar sobre as populações que eu não vi. Sempre somos limitados, quanto a tempo disponível ou mesmo recurso de o quanto podemos coletar, mas usar essa ideia faz com que não tenhamos que coletar todas as populações existentes para concluir algo para essa espécie, se coletarmos algumas populações, podemos pensar como seria as populações que não coletamos.

Mas como ajustamos esse tipo de modelo? No R existem 2 pacotes mais comuns para isso, o nlme que é mantido pelo Core Team do R e o lme4, o lme4 é legal pois ele esta na versão 0.999…, desde que foi programado pelo autor, o Douglas Bates nunca achou que o pacote esta bom o suficiente para virar a verão 1, mas computadores tem um limite de casas depois da virgula que salvam os números, então cedo ou tarde vai faltar casa depois da virgula.
Mas esses pacotes ajustam modelos mixtos e vamos usar eles para ajustar nossos exemplos.

Apesar de ser difícil de entender, ajustar eles com o R é fácil, no caso do lme4 é so descrever a formula, colocando de uma forma diferente as variáveis aleatórias:

summary(mod1.lme4) Linear mixed model fit by REML Formula: massa ~ tamanho + (1 | pop) AIC BIC logLik deviance REMLdev 5934 5951 -2963 5933 5926 Random effects: Groups Name Variance Std.Dev. pop (Intercept) 406.98 20.174 Residual 2108.30 45.916 Number of obs: 560, groups: pop, 56 Fixed effects: Estimate Std. Error t value (Intercept) 234.177 3.321 70.51 tamanho 64.293 2.014 31.92 Correlation of Fixed Effects: (Intr) tamanho 0.000

No nlme, ja é um pouco diferente, a gente descreve a formula normalmente, mas as variáveis aleatórias são descritas num argumento diferente.

summary(mod1.nlme) Linear mixed-effects model fit by REML Data: NULL AIC BIC logLik 5933.979 5951.277 -2962.99 Random effects: Formula: ~1 | pop (Intercept) Residual StdDev: 20.1739 45.91627 Fixed effects: massa ~ tamanho Value Std.Error DF t-value p-value (Intercept) 234.17738 3.321511 503 70.50326 0 tamanho 64.29322 2.014287 503 31.91859 0 Correlation: (Intr) tamanho 0 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.49417909 -0.61376308 -0.01319703 0.56456887 3.15519916 Number of Observations: 560 Number of Groups: 56

Esses 2 modelos são os conhecidos como “random intercepts”, ou seja, somente os interceptos são aleatórios, a inclinação é fixa, é o mesmo valor para todas as populações, so bater o olho no nosso gráfico que reconhecemos não é isso que queremos. Logos precisamos tanto de interceptos como inclinações aleatórias, então escrevemos o modelo assim:

summary(mod2.lme4) Linear mixed model fit by REML Formula: massa ~ tamanho + (1 | pop) + (0 + tamanho | pop) AIC BIC logLik deviance REMLdev 5603 5625 -2797 5602 5593 Random effects: Groups Name Variance Std.Dev. pop (Intercept) 502.33 22.413 pop tamanho 1254.65 35.421 Residual 818.39 28.608 Number of obs: 560, groups: pop, 56 Fixed effects: Estimate Std. Error t value (Intercept) 231.356 3.258 71.01 tamanho 65.210 4.916 13.27 Correlation of Fixed Effects: (Intr) tamanho 0.000

e no nlme:

summary(mod2.nlme) Linear mixed-effects model fit by REML Data: NULL AIC BIC logLik 5602.506 5628.452 -2795.253 Random effects: Formula: ~1 + tamanho | pop Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 22.39236 (Intr) tamanho 35.35674 0.241 Residual 28.61556 Fixed effects: massa ~ tamanho Value Std.Error DF t-value p-value (Intercept) 231.37580 3.255510 503 71.07206 0 tamanho 65.16188 4.907942 503 13.27683 0 Correlation: (Intr) tamanho 0.213 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.4933849 -0.6276726 -0.0029472 0.6140966 2.7334334 Number of Observations: 560 Number of Groups: 56

A primeira coisa que eu devo destacar é o fato de que o ajuste do lme4 não tem valores p, como ja falei no post anterior, as estimativas são muito mais informativas que valores p de modo geral, mas para salvação de todos, o nlme fornece valores p.
O fato de o lme4 não fornecer os valores p ja deu muitas discussões acaloradas, não sei dizer se isso inclusive tem haver com o fato de existirem esses 2 pacotes para o mesmo fim praticamente, o lme4 e o nlme, mas o Bates explicando sobre valores p é bem legal.

Mas vamos pegar o resultado do lme4 e entender melhor.

Primeiro temos uma tabelinha parecida com a de uma regressão.

Fixed effects: Estimate Std. Error t value (Intercept) 231.356 3.258 71.01 tamanho 65.210 4.916 13.27

Nessa parte estamos vendo os valores médios de intercepto e inclinação, os valores médios para populações observadas e o que esperamos para outras populações que coletaremos.

Random effects: Groups Name Variance Std.Dev. pop (Intercept) 502.33 22.413 pop tamanho 1254.65 35.421

Aqui, nos vemos a variação, o desvio da média para cada um dos parâmetros dados as diferenças entre populações, ou seja o intercepto tem em média o valor apresentado em cima e esse desvio visto aqui, o mesmo vale para a inclinação. E se coletamos outra população deve ser algo por ai. Podemos comparar esse valor com o que usamos para simular os dados e vemos que, apesar de provavelmente entender muito pouco da conta como um todo, ela funciona bem, resgata bem os valores que usamos na simulação.

Mas e se quisermos discutir população por população? Talvez não todas, mas queremos saber parâmetros caso a caso, com a função coef() vemos os parâmetros estimados para cada população, e podemos inclusive comparar com o que usamos na simulação para ver se estamos tendo boas estimativas.

data.frame(coef(mod2.lme4)$pop,intercepitos,inclinações) X.Intercept. tamanho intercepitos inclinações 1 227.3287 -9.604714 219.7503 -1.993714 2 272.4481 85.263805 276.7211 86.336342 3 240.9874 87.270267 243.0946 102.992448 4 228.8367 45.596422 227.2121 62.895739 5 234.3060 -2.264402 219.9204 -9.558904 6 216.0660 54.134587 231.1457 46.185633 ... (omitindo um pouco dos resultados) 55 241.5680 9.347568 238.6680 5.602856 56 256.0827 68.437325 260.5831 69.414471

E por ultimo, comparamos o modelo somente com interceptos aleatórios com o modelo com interceptos e inclinações aleatórios e vemos que o com interceptos e inclinações aleatórios descreve melhor os dados, afinal fizemos ele assim 🙂

anova(mod1.lme4,mod2.lme4) Data: Models: mod1.lme4: massa ~ tamanho + (1 | pop) mod2.lme4: massa ~ tamanho + (1 | pop) + (0 + tamanho | pop) Df AIC BIC logLik Chisq Chi Df Pr(>Chisq) mod1.lme4 4 5941.5 5958.8 -2966.7 mod2.lme4 5 5612.3 5633.9 -2801.1 331.17 1 < 2.2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Essa ideia de modelos mixto é extensível para muitas possibilidades, aqui vemos um exemplo simples, de varias populações e a relação entre tamanho e massa, mas, qualquer coisa como regressão múltipla , ancova e outros modelos podem ser pensado dessa forma, só precisamos de um cuidado especial para imaginar o que faz sentido ser efeito aleatório e o que faz sentido ser efeito fixo, senão acaba pensando tudo como aleatório. E é bom pensar que para algo ser aleatório, temos que esperar que outras categorias existam e que não foi possível coleta-las, mas queremos inferir algo sobre elas.

Esse exemplo eu tirei do excelente livro do Marc Kéry . Espero que esse exemplo sirva para ter uma ideia do que se trata a ideia de modelos mixtos, espero que a próxima vez que falar de modelo mixto, o post fique menor.

####################################################
# Modelo linear Mixto #
####################################################
 
set.seed(13)
### Gerando os dados
n.grupos <- 56 # Número de populações
n.amostras <- 10 # Número de Cobras por amostras
n <- n.grupos * n.amostras # Total de amostras (Total de cobras)
pop <- gl(n = n.grupos, k = n.amostras) # Indicador da população
 
# Tamanho do corpo (cm)
tamanho.original <- runif(n, 45, 70)
mn <- mean(tamanho.original)
sd <- sd(tamanho.original)
 
#Mesmo que usar o comando scale(), isso faz a média ser zero
tamanho <- (tamanho.original – mn) / sd
 
Xmat <- model.matrix(~pop*tamanho-1-tamanho) #Aqui estamos fazendo o
#contrario que o comando formula faz, para simular uai
 
intercepito.mean <- 230	 # intercepto médio
intercepito.sd <- 20 # desvio dos interceptos
inclinação.mean <- 60 # inclinação média
inclinação.sd <- 30 # desvio das inclinações
 
intercepitos<-rnorm(n = n.grupos, mean = intercepito.mean, sd = intercepito.sd)
inclinações <- rnorm(n = n.grupos, mean = inclinação.mean, sd = inclinação.sd)
todos.os.efeitos <- c(intercepitos,inclinações) # Juntando tudo
 
lin.pred <- Xmat[,] %*% todos.os.efeitos # Preditores lineares
eps <- rnorm(n = n, mean = 0, sd = 30) # residuos
# resposta = preditor linear + residuo
massa <- lin.pred + eps
 
#pegando cores mais bonitas
library(RColorBrewer)
#escolhendo na palheta
display.brewer.all()
cols <- brewer.pal(n.grupos, "Set1")
 
#figura 1
plot(massa~tamanho,type="p",frame.plot=F,col=cols,pch=19)
abline(lm(massa~tamanho),lwd=4)
 
#figura 2
library("lattice")
xyplot(massa ~ tamanho | pop,
panel = function(x, y) {
panel.grid(h = -1, v = 2)
panel.xyplot(x, y,pch=19,col="black")
panel.lmline(x,y,col="red",lwd=1.5)
})
 
### analise com interceptos aleatorios
library(lme4)
mod1.lme4 <-lmer(massa ~ tamanho + (1 | pop))
summary(mod1.lme4)
 
library(nlme)
mod1.nlme<-lme(massa ~ tamanho,random= ~ 1 | pop)
summary(mod1.nlme)
 
### Analise com interceptos e inclinações aleatorios
 
mod2.lme4<- lmer(massa ~ tamanho + (1 | pop) + ( 0+ tamanho | pop))
summary(mod2.lme4)
data.frame(coef(mod2.lme4)$pop,intercepitos,inclinações)
 
mod2.nlme<-lme(massa ~ tamanho,random= ~ 1 + tamanho | pop)
summary(mod2.nlme)
 
#comparando modelos
anova(mod1.lme4,mod2.lme4)