Teste t de Poisson

O titulo de post é um nome bem esquisito, e errado também, teste t usa a estatística t, então não da para ter um teste t de Poisson exatamente. Mas a gente acaba associando alguns nomes a situações.
Teste t para avaliar um fator de dois níveis, anova para um fator de vários níveis, ancova para variáveis continuas e variáveis discretas e assim vai.
A questão é que todos esses métodos estão na verdade dentro de uma categoria, existe uma teoria que unifica todos esses métodos, esses e outros, como aqueles que já comentamos usando distribuição binomial, ou mesmo com distribuição de poisson como vemos aqui são todos casos da mesma “estratégia” em estatística aplicada, se meus termos aqui não são ruins. Mas em 1972, John Nelder e Robert Wedderburn unificaram mais ou menos tudo que falamos ali em cima na ideia geral de Modelos gerais linearizados.. Eles mostraram que todas essas técnicas, representadas como tipos diferentes de analises, incluindo chi-quadrados até que a gente ve muito em génetica, podiam todos serem representados como casos especiais de modelos lineares. Quando eu fiz estatística na faculdade, na graduação, eu nunca fui apresentado a essa ideia que tudo é “modelo linear”, que tudo esta interligado e por isso muitas ideias são comuns, teste t e anova eram coisas distintas, técnicas que se você tem a situação a você aplica isso, se tem a situação b, aplica aquilo e se por ventura tem a situação c você se ferrou, porque não tem um teste pronto pra você. E isso gera o problema, que as vezes as pessoas abandonam a ideia, pressupostos da biologia por exemplo, para forçar os dados serem aplicáveis a um teste e temos que lidar com um monte de problemas decorrentes disso. Além é claro de isso muitas vezes diminuir o avanço cientifico, se é que alguém ainda se importa com isso.

Mas então, um modelo geral linear vai ter mais ou menos três componentes.

1. Uma distribuição estatística, usada para descrever a parte aleatória da resposta do modelo, o eixo y dos gráficos que a gente normalmente faz. Essa é a parte estocástica do sistema que estamos tentado descrever.
2. Uma função de ligação, que é aplicada ao valor esperado, um E(y).
3. E um preditor linear, os a+bx, intercepto e inclinação, as variáveis preditoras que na maioria das vezes estamos interessados. se a função de ligação é a função g, então um g(E(y)). Essa é a parte determinística do modelo, veja que aqui entra qualquer função matemática na verdade, mas essa ai é bem comum.

Nesse sentido a regressão linear que a gente faz é o seguinte.

1. Distribuição: y_i \sim Normal(\mu_i , \sigma^2)
2. Função de ligação: Identidade, \mu_i = E(y_i) que é o preditor linear.
3. Preditor linear: \alpha + \beta \cdot x_i

Se você pensa que nunca fez isso, veja que todas as simulações que fazemos aqui a gente meramente vai do item 3 pro item 1, a gente gera os valores determinísticos, transforma pela função de ligação e depois gera a parte aleatória dos dados, e depois roda uma analise para ver se recuperamos os valores do preditor linear.

Mas e um modelo usando Poisson?

Ai teremos o seguinte:

1. Distribuição: C_i \sim Poisson(\lambda_i)
2. Função de ligação: \log{(\lambda_i)} ou \log({E(C_i)})
3. Preditor linear: \alpha + \beta \cdot x_i

Viu, a mesma ideia, só mudamos a função de ligação, que antes era a identidade, que nada mais é que o mesmo valor do número, agora usamos o log, log na base natural, que é o número do euler la, que a gente ja esbarrou por ele aqui.

Bem agora a analise, primeiro vamos entender os dados.

Vamos supor que queremos estudar o número de aranhas que encontramos em arbustos. Então a gente pega um arbusto inteiro, poe ele num sacão e leva para o laboratório, la vasculhamos ele e contamos quantas aranhas conseguimos achar. Particularmente os arbustos são de uma espécie de planta que ocorre tanto na zona rural como na zona urbana, em terrenos baldios. Mas nos achamos que aranhas não conseguem dispersar para dentro da zona urbana, logo esperamos que elas sejam restritas a áreas rurais, certo? Claro que uma ou outra aranha vai acabar chegando na cidade, ainda mais que estamos falando “aranhas” sem considerar as espécies delas. Mas de modo geral esperamos menos aranhas na cidade, ou mais claramente falando, esperamos que em média, arbustos na cidade tenham menos aranhas, ou seja em média as contagens de aranhas serão menor na cidade quando comparado a zona rural.

Então temos a pergunta, existe diferença no número de aranhas entre a área urbana e a área rural?

Bem sabendo o que queremos saber, fica fácil, é so pegar arbustos na cidade e na área rural e contar a aranhas e ver se existe diferença nas contagens de aranhas em cada área. Mas porque não usar um teste t normal? Porque primeiro aranhas so vem inteiras, não existe meia aranha nem um quarto de aranha, como a distribuição normal é continua, ela vai prever uma variável continua, o que é ruim. Mas além disso, a distribuição normal, se eu ver um número muito baixo de aranhas, algo que deixe a média perto de zero, eu vou prever números negativos de aranhas em arbusto, o que não pode existir, se num local existem menos aranhas, eu tenho que ver mais zeros e não números negativos. Esses dois problemas são solucionados usando a distribuição de Poisson, existem mais coisas a se considerar, mas por agora vamos em frente, outras considerações ficam como cenas dos próximos capítulos.

Então fomos la no mato, coletamos os arbustos e contamos a aranhas, e vamos ver o que conseguimos:

02

Olha ai, pegamos 20 arbustos no total, 10 na zona rural e 10 na zona urbana. Como ja suspeitávamos, na cidade não vimos mais que 4 aranhas por arbusto, enquanto no campo olha la, a gente tinha muito mais, até 8 aranhas. Podemos fazer um box-plot.

01

Certo, são os mesmo dados, mas agora acho que visualmente da para ver outra questão, aqui temos uma variação muito maior dos dados na área rural e isso quebra uma das premissas de modelos lineares em geral, como o teste t, que é homogeniedade de variâncias. Mas não é difícil pensar que se a gente sempre vê em média uma aranha, não tem como ter uma variação igual quando a gente tem muitas, pense que não tem como contar algo para menos que zero então um número pequeno implica de certo modo que só podemos ter uma variação pequena, agora se a gente contar em média 100 aranhas nos arbustos, seria fácil pensar por exemplo que erramos e contamos 90, ou 110, e isso faz um desvio muito grande em relação a uma contagem pequena, o ponto é que o desvio tem uma relação com o quanto a gente conta, ou seja com o \lambda, mas na distribuição de poisson é exatamente assim que funciona, o \lambda determina a média e a variação dos dados, no caso quanto maior a contagem média, maior a variação. Temos como mexer um pouco nessa relação, mas isso é outro post, por agora já da para ter uma noção de outra característica do modelo que estamos fazendo, que é em modelos usando distribuição de poisson, a média e o desvio estão relacionados entre si.

Mas podemos ver as médias das áreas também, para ver como ficaram.

> aggregate(C, by = list(x), FUN = mean) Group.1 x 1 Urbana 2.4 2 Rural 5.3

E a média da área rural é quase o dobro da média na área urbana. Agora podemos tentar definir isso como parâmetros de um modelo linear generalizado usando a distribuição de poisson.

Primeiro temos que definir o modelo na linguagem bugs.

# Definindo o modelo na linguagem Bugs
sink("Poisson.t.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i]
# Ajuste
    residuo[i] <- (C[i] - lambda[i]) / sqrt(lambda[i])
}
}
",fill=TRUE)
sink()

Além do modelo em si, vamos salvar os resíduos.

Agora se pararmos um momento para olhar o likelihood, podemos ver os itens que definimos la em cima.
O modelo é:

# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i]
}

E olha aqui, temos um loop para processar todos as amostras, e dentro dele temos:

1. A distribuição de poisson

C[i] ~ dpois(lambda[i])

2. A função de ligação

log(lambda[i])

3. A parte determinística do modelo.

alpha + beta *x[i]

Esta numa linguagem de computador, mas é exatamente o que definimos la em cima.
Mas agora com o modelo pronto, definimos as configurações do MCMC e rodamos ele.
Com os resultados em mão, antes de fazer qualquer inferência, primeiro de tudo a gente olha se as cadeias convergiram. Existem muitas formas de ver isso, mas primeiro podemos olhar para os valores de Rhat, eles são simples de interpretar, porque é so buscar se estes estão inferiores a 1.1, o valor limite, acima disso temos que pensar bem como prosseguir, mas nesse caso…

03

Todos menores que 1.1.

Podemos fazer um teste rapido, e simplesmente testar se alguém ficou acima de 1.1

which(out$summary[,8] > 1.1)

Certo, outra coisa, podemos olhar agora os resíduos, estes tem que ter uma distribuição normal, se eles não tiverem uma distribuição normal, significaria que alguma coisa muito importante que não estamos levando em conta esta acontecendo, que tem um efeito claramente mudando muito tudo. Se os resíduos estiveram normal, com média 0, quer dizer que estamos errando igualmente tanto para mais quanto para menos.

04

Olha ai, passamos uma linha no zero, e temos mais ou menos a mesma quantidade de pontos pra cima e para baixo, e nenhuma tendência aparente. Podemos olhar também esses resíduos na forma de um histograma.

05

E apesar de umas barras talvez esquisitas, podemos ver também a densidade em vermelho, que é um “histograma” infinito, mais suave, e vemos bem uma distribuição normal, com média perto do zero, marcado em azul na figura.

Existem também testes para testar a normalidade, como o teste de shapiro-wilks.

> shapiro.test(out$mean$residuo) Shapiro-Wilk normality test data: out$mean$residuo W = 0.9374, p-value = 0.2144

Mas a gente tem que tomar cuidado para não entrar num paradoxo, se a gente começar a testar o teste do teste do teste, a gente não vai parar nunca, e nunca voltaremos para a questão que realmente nos interessa.
Mas pelo menos eu ja estou mais que convencido que o modelo teve um ajuste bom, não tem problemas, então poderíamos olhar as conclusões que ele nos traz, quais inferências podemos fazer.

06

E exatamente o \beta é o que nos interessa, ele é a diferença entre a área rural e a área urbana.
Voltando a nossa pergunta inicial, que era se existia diferença no número de aranhas entre a área urbana e a área rural, um valor de \beta=0 significaria que não há diferença. Como em todas as amostras do MCMC, não temos nenhum 0 nesse caso, nem nada menor que zero estamos fortemente propensos a acreditar que existe uma diferença. Podemos calcular o intervalo de confiança de 95%, ou seja, qual intervalo mais comum de valores e falar que este é maior que zero. Poderíamos sei la, calcular um intervalo de 99%, mas aqui o que interessa é que podemos fazer uma afirmação mais ou menos assim.

“Temos 95% de certeza que deve haver uma diferença no número de aranhas presente em arbustos na zona rural e zona urbana.”

Podemos apresentar então quantas aranhas esperamos encontrar em cada local, e a partir daqui é tentar entender melhor porque existe essa diferença.

07

E assim concluímos nosso primeiro glm usando distribuição de Poisson, já que esse post está ficando grande mais. Mas da para ver como os glms abrem uma grande gama de perguntas que podemos abordar. Espero que não tenha falado muita bobeira no caminho, mas modelos lineares usando distribuição normal, binomial e de poisson cobrem a grande maioria das coisas que lemos comumente na literatura científica relacionada a biologia.

Bem o script está disponível no repositório recologia no github, além de aqui em baixo, e acho que é isso, até o próximo post.

Referencias:
Marc Kéry (2010) Introduction to WinBUGS for Ecologists. Academic Press, Burlington.

set.seed(123)
### Gerando os dados
#Gerando o preditor linear, o item 3 do post.
n.arbustos <- 10
x <- gl(n = 2, k = n.arbustos, labels = c("Urbana", "Rural"))
n <- 2*n.arbustos
 
#Aplicando a função de ligação, exponencial é o contrario de logaritimo,
#estamos no caminho contrario aqui, esse é o item 2 do post
lambda <- exp(0.69 + 0.92*(as.numeric(x)-1))
 
#Adicionando a parte estocastica, o item 1 do post
C <- rpois(n = n, lambda = lambda)
 
#Vamos olhar o que criamos
aggregate(C, by = list(x), FUN = mean)
 
boxplot(C ~ x, col = "grey", xlab = "Tipo de ambiente", ylab = "Contagem de aranhas", las = 1,frame=F)
 
stripchart(C ~ x,xlab = "Tipo de ambiente", ylab = "Contagem de aranhas", las = 1,frame=F,
           method="stack",vertical=T,pch=19,col=1,at=c(1.4,2.1),offset = 0.8)
legend("topleft",pch=19,legend="Um arbusto")
 
 
### Analise Bayesiana
# Definindo o modelo na linguagem Bugs
sink("Poisson.t.txt")
cat("
model {
 
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i]
 
# Ajuste
    residuo[i] <- (C[i] - lambda[i]) / sqrt(lambda[i])
 
}
 
}
",fill=TRUE)
sink()
 
# Juntando os dados para mandar para o Openbugs
dados.bugs<- list(C = C, x = as.numeric(x)-1, n = length(x))
 
# Função geradora de valores iniciais
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))}
inits()
 
# Parametros a estimar
params <- c("alpha", "beta","residuo")
 
# Configurações do MCMC
nc <- 3
ni <- 3000
nb <- 1000
nt <- 2
 
# Iniciando o Gibss Sampler
library(R2OpenBUGS)
out <- bugs(data=dados.bugs, inits=inits, parameters.to.save=params,model.file="Poisson.t.txt",
            n.thin=nt, n.chains=nc, n.burnin=nb,n.iter=ni)
 
 
### Checando convergencia e ajuste do modelo
 
hist(out$summary[,8],col="grey",main="Valores de Rhat")
which(out$summary[,8] > 1.1)
 
plot(out$mean$residuo, las = 1,frame=F,ylab="Resíduos",xlab="Amostras",xlim=c(0,20))
abline(h = 0)
 
hist(out$mean$residuo,xlab="Resíduos",freq = F)
abline(v = 0,lwd=2,col="blue",lty=3)
points(density(out$mean$residuo),type="l",lty=2,lwd=2,col="red")
 
shapiro.test(out$mean$residuo)
 
### Inferencia baseada no modelo que criamos
print(out, dig = 3)
 
#Função do Livro Doing Bayesian Analysis
HDIofMCMC <- function( sampleVec , credMass=0.95 ) {
    sortedPts = sort( sampleVec )
    ciIdxInc = floor( credMass * length( sortedPts ) )
    nCIs = length( sortedPts ) - ciIdxInc
    ciWidth = rep( 0 , nCIs )
    for ( i in 1:nCIs ) {
        ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
    }
    HDImin = sortedPts[ which.min( ciWidth ) ]
    HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
    HDIlim = c( HDImin , HDImax )
    return( HDIlim )
}
 
int95<-HDIofMCMC(out$sims.list$beta)
hist(out$sims.list$beta, col = "grey", las = 1, xlab = "Coeficiente para area rural", main = "",freq=F)
points(density(out$sims.list$beta),type="l",lty=2,lwd=3,col="red")
arrows(int95[1],0.25,int95[2],0.25,col="blue",angle =90,lwd=3,lty=3,code=3)
legend("topleft",legend="Intervalo de confiança de 95%",lwd=3,lty=3,col="blue",bty="n")
 
par(mfrow = c(2,1))
int95<-HDIofMCMC(exp(out$sims.list$alpha))
hist(exp(out$sims.list$alpha), main = "Área urbana",col = "grey", xlab = "", xlim = c(0,10), breaks = 20,freq=F)
points(density(exp(out$sims.list$alpha)),type="l",lty=2,lwd=3,col="red")
arrows(int95[1],0.1,int95[2],0.1,col="blue",angle =90,lwd=3,lty=3,code=3)
legend("topright",legend="Intervalo de confiança de 95%",lwd=3,lty=3,col="blue",bty="n")
 
int95<-HDIofMCMC(exp(out$sims.list$alpha + out$sims.list$beta))
hist(exp(out$sims.list$alpha + out$sims.list$beta), main = "Área rural",freq=F,
     col = "grey", xlab = "Número de aranhas esperado", xlim = c(0,10), breaks = 20)
points(density(exp(out$sims.list$alpha + out$sims.list$beta)),type="l",lty=2,lwd=3,col="red")
arrows(int95[1],0.1,int95[2],0.1,col="blue",angle =90,lwd=3,lty=3,code=3)
 
### Analise Frequentista
poisson.t <- glm(C ~ x, family = poisson)
summary(poisson.t)
anova(poisson.t, test = "Chisq")

Teste t Binomial – GLM com prevalências.

A gente ja usou um bocado estatística bayesiana usando Openbugs. Os modelos feitos no Openbugs são facilmente replicados no Winbugs, Jags e com um pouquinho de dedicação no Stan. Mas a maioria dos posts anteriores era apenas usando a distribuição normal para descrever os dados.

Mas as vezes não é possível descrever o mundo assim. Um exemplo é quando a gente trabalha com algo como a prevalência de um parasita.
Vamos pegar um exemplo aqui, passarinhos sendo parasitados por Haemoproteus columbae.

Haemoproteus é um genero de protozoa que parasita passarinhos, repteis e anfíbios. O seu nome é derivado do grego, “Haima” = Sangue e
“Proteus” = Deus dos mares que tem o poder de assumir varias formas diferentes. Dai da para sacar que que ele vive no sangue e tem mais de uma forma.

vol102(3)_foto_capa2

Então temos passarinhos, e eles podem ou não estar contaminado. Mas como descobrimos?

Bem existem dois jeitos, ou a gente faz uma lâmina do sangue, prepara ela e olha no microscópio e fica procurando o Haemoproteus ou a gente prepara o sangue, faz uma pcr colocando primers que replicam o DNA do haemoproteus mas não do passarinho hospedeiro (devem vir varias células, logo vários DNAS diferentes num pouquinho de sangue) e se houver bastante DNA amplificado no final, bingo, o parasita tem que estar ali.

Mas no final das contas o que vamos falar é o seguinte, vamos ter sangue de vários hospedeiros, e nos dois casos vamos apenas falar se o indivíduo, o passarinho esta ou não parasitado pelo Haemoproteus.

Individuo rep.dados.Sp..each...15. Parasita 1 1 Columbina talpacoti Sim 2 2 Columbina talpacoti Sim 3 3 Columbina talpacoti Sim 4 4 Columbina talpacoti Sim 5 5 Columbina talpacoti Não . . . 25 10 Scardafella squammata Sim 26 11 Scardafella squammata Não 27 12 Scardafella squammata Não 28 13 Scardafella squammata Sim 29 14 Scardafella squammata Não 30 15 Scardafella squammata Não . . .

Acho que da para entender como estão os dados, temos uma identidade para cada indivíduo, cada amostra, ou seja um passarinho é uma amostra. Cada amostra pertence a uma espécie, e sabemos se ele esta ou não parasitada, temos algo como Sim ou Não, parasitado ou não parasitado, como quiser chamar.

Bem como são sim ou não, normalmente as pessoas resumem essa informação em contagens de quantos passarinhos estão parasitados e quantos não estão parasitados.

> tabela Parasita Sp Sim Não Columbina talpacoti 32 30 Scardafella squammata 10 47

Aqui mantemos a informação de quantos passarinhos foram amostrados (soma de todos as células da matriz), quantos de cada espécie, quantos parasitados ou não por espécie, notem que com essa matriz podemos reconstruir a planilha la em cima, isso se não nos importamos com a identidade de cada passarinho para outras informações, mas temos aqui a mesma informação que la em cima.

Virou uma tabelinha assim, a primeira coisa que vem a cabeça é teste de chi-quadrado para ver se temos um desvio nas proporções. Como assim? Uai se o parasita ocorre numa frequência constante, independente da espécie, então esperaríamos uma proporção igual entre as linhas e colunas, agora se algo acontece, uma espécie é mais parasitada ou, sei la, a prevalência é baixinha nas espécies, ou a gente pode observar se o desvio das contagens é devido ao acaso ou há razões para desconfiarmos que algo a mais acontece.

Então realizamos o teste de chi-quadrado:

> teste.chi <- chisq.test(tabela) > teste.chi Pearson's Chi-squared test with Yates' continuity correction data: tabela X-squared = 13.6387, df = 1, p-value = 0.0002216 > teste.chi$observed # Contagens observadas Parasita Sp Sim Não Columbina talpacoti 32 30 Scardafella squammata 10 47 > teste.chi$expected # Contagens esperada sobre a hipotese nula Parasita Sp Sim Não Columbina talpacoti 21.88235 40.11765 Scardafella squammata 20.11765 36.88235

E temos o clássico p muito menor que 0.05, e mais um monte de valor. Bem, como não fornecemos as probabilidades para o esperado, ele calcula automaticamente, mas podemos ver o que esperávamos no resultado olhando o $expected do modelo ajustado. Bem eu ja fiquei perdido nesses números. Quando a gente olha por exemplo o resultado de cruzamentos de Drosophila melanogaster em relação a algumas mutações, é fácil estabelecer mentalmente o que é o esperado, 1:1 em caso do sexo, 3:1 no caso de uma mutação dominante e assim vai, mas essa pergunta pode ser mais complicada para alguns casos. Então a pergunta, “o que você espera de resultado” é sempre uma pergunta interessante. Mas de qualquer forma aqui a gente vê que algo acontece, agora vamos olhar mais a frente.

Outra forma de olhar os dados é calcular a prevalência, que é, para uma espécie nesse exemplo, o quanto de casos de parasitados tivemos do total de dados.

> dados Sp N Prevalence simul 1 Columbina talpacoti 62 0.516 32 2 Scardafella squammata 57 0.193 10

E temos a prevalência, que é mais ou menos a chance de um indivíduo estar parasitado. E Para nossa surpresa, pelo menos pra minha, dividir o número de indivíduos parasitados pelo total de indivíduos amostrados da a máxima verosimilhança (Maximum likelihood) da distribuição Binomial, mas no próximo post vamos olhar isso.
Mas antes de continuar, vamos ver o que estamos querendo falar ao dizer que 0.516, ou 52% é a prevalência do Haemoproteus na espécie Columbina talpacoti

01

Bem, são dados de contagem aqui, mas temos que ter em mente que essa contagem tem um limite, que no primeiro caso é 62 e no segundo 57. Deixamos riscado no gráfico esse limite.
Agora para o gráfico 1, onde esta a bolinha mais alta? Oras, no 32, o maximum likelihood, mas existe um erro inerente a qualquer amostragem, mas para essa mesma prevalência, ainda seria razoável pensar nela como o valor real de prevalência, a real chance de um passarinho ser contaminado mesmo com 31 passarinhos doentes. Agora pense no valor de 30 ou 29, ja estamos ficando com metade da altura que tem o maximum likelihood, aqui talvez não seja muito legal pensar em 51.6% como prevalência, e pensar que a prevalência é um pouco mais baixa. O mesmo vale para outra especie.

Mas então temos crenças diferentes para diferentes valores de prevalência para cada uma das especies. Mas além de estabelecer quanto é a prevalência de cada espécie, podemos tentar estimar a diferença entre essas prevalência. Que seria o que? Algo como o quanto deve ser a diferença entre elas.

So que aqui tem um esquema, meramente falando em diferença, tudo fica meio confuso, a gente não esta meramente diminuindo 0.516 de 0.193. A gente quer saber a distribuição, assim como essas que estimamos para cada especie, a distribuição das diferenças, porque não vamos simplismente falar um valor, ele pode estar errado, a gente tem que ver todas as possibilidades de diferenças.

E a gente pode fazer isso como similar a forma como fizemos o teste t aqui, so que para usar a distribuição binomial, precisamos de uma função de ligação, que nesse caso é a função logistica.

Vamos ver primeiro como é a função na linguagem Bugs:

sink("Binomial.t.test.txt")
cat("
model {
 
# Priors
 alpha ~ dnorm(0,0.01)
 beta ~ dnorm(0,0.01)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dbin(p[i], N[i])
    logit(p[i]) <- alpha + beta *species[i]
 }
# Derived quantities
 prev.C_talpacoti <- exp(alpha) / (1 + exp(alpha))
 prev.S_squammata <- exp(alpha + beta) / (1 + exp(alpha + beta))
 prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste
}
",fill=TRUE)
sink()

A gente ja fez quase isso aqui, sem usar MCMC. Mas agora a diferenças que a gente tem, ao invés de apenas 1 parâmetro, um valor de probabilidade (que aqui a gente le prevalência), a gente tem dois, uma para uma espécie e outro para outra espécie.

Então fazemos aqui temos que fazer isso para as duas espécies.

Então o likelihood tem um loop que vai de 1:2, aqui o n é dois quando a gente preenche os dados.

# Likelihood
 for (i in 1:n) {
    "Meu modelo aqui"
 }

Dentro desse loop a gente vai por nosso modelo.
E nosso modelo é como?
Tem uma parte estocástica e uma parte determinística.

Como é a parte estocástica?

C[i] ~ dbin(p[i], N[i])

Oras, a parte estocástica é uma probabilidade, que depende do N, o número de amostras, é uma distribuição binomial. Mas esse p é determinado pela nossa função determinística, e como é ela?

logit(p[i]) <- alpha + beta *species[i]

Igual foi em todo modelo linear, a+bx, intercepto + inclinação * espécie. Vezes a espécie porque eu vou dar um valor para cada espécie. Aqui que ta o macete, eu vou atribuir a probabilidade de uma espécie como alpha, e o da outra como alpha main beta. so que beta é multiplicado por um valor, que é a hora que eu mando junto dos dados o 0 ou 1, ou seja, vamos supor um exemplo simples, eu tenho espécies A e B, A tem prevalência de 20 e B tem prevalência de 30
Nesse caso, para esse modelo, o intercepto é 20, a inclinação é 10 e o beta da especie A é 0 e o beta da especie B é 1.
Assim ficamos com o seguinte:

Para a espécie A A = interceopto + inclinação * especie A = 20 + 10 * 0 A = 20

Como qualquer número multiplicado por zero é zero, esse 10 some para a especie 1, e ela é meramente o valor do intercepto.
Ja pra a espécie B, a espécie 1, temos o seguinte:

Para a espécie B B = interceopto + inclinação * especie B = 20 + 10 * 1 B = 30

Oras, 10 vezes 1 é 10 mesmo, que somado com o 20 da 30, a prevalência da segunda espécie.
Isso é só pra facilitar a visualização do que ta dentro do likelihood ali, mas observe que o valor de p(de prevalência) esta dentro da função logit, a função logística, usado para generalizar para linearizar a relação e assim usar a distribuição binomial. Por isso o a+bx da equação vai sair com um número esquisito de olhar, que não é porcentagem, não é um número entre 0 e 1, entre 0% e 100%. Para olhar ele depois temos que destransformar, e fazemos isso usando as quantidades derivadas no Openbugs.

# Derived quantities
 prev.C_talpacoti <- exp(alpha) / (1 + exp(alpha))
 prev.S_squammata <- exp(alpha + beta) / (1 + exp(alpha + beta))
 prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste

Agora olha essa equação ai e lembra da formula da equação logística, poxa aqui a gente so ta destransformando a continha que falamos acima, veja a formula no wikipedia, é exatamente isso a formula da função logística, so lembrar que o inverso de log é exponencial. So que além de acessar os valores de prevalência das duas espécies, olhamos também a diferença entre eles, e ja teremos uma saída direto.

Agora vamos rodar o modelo e ver o que acontece.

> out <- bugs(data=bugs.dados, inits=inits, parameters.to.save=params,model.file="Binomial.t.test.txt", + n.thin=nt, n.chains=nc, n.burnin=nb,n.iter=ni) > print(out, dig = 3) Inference for Bugs model at "Binomial.t.test.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 alpha 0.065 0.253 -0.432 -0.113 0.068 0.243 0.539 1.001 3000 beta -1.636 0.424 -2.452 -1.929 -1.637 -1.342 -0.807 1.001 2400 prev.C_talpacoti 0.516 0.062 0.394 0.472 0.517 0.561 0.632 1.001 3000 prev.S_squammata 0.177 0.049 0.095 0.142 0.174 0.208 0.285 1.001 3000 prev.Diff 0.339 0.079 0.180 0.284 0.342 0.395 0.481 1.003 2800 deviance 10.488 1.855 8.609 9.127 9.914 11.250 15.331 1.001 2100 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 = 1.9 and DIC = 12.4 DIC is an estimate of expected predictive error (lower deviance is better).

Primeiro olhamos o Rhat para ver a convergência dos valores, senão nem adianta olhar o resultado, como todos estão ok a gente procede para olhar os resultados, se você não lembra o que é o Rhat, pode dar uma conferida aqui.

Mas então o alpha e beta, que é o a+bx da equação da reta ta ae, mas são números que estão transformados pela função logística. Vamos olhar direto a prevalência das espécies, que ja destransformamos os dados diretamente no bugs, e la temos, dado que simulamos dados com essas mesmas características, recuperamos muito bem os valores de prevalência, como recuperamos bem, vamos olhar a diferença, e a diferença de prevalências tem média de 0.339 e desvio de 0.079. Ta mas mais legal que isso vamos olhar qual o intervalo de confiança de 95% dele, é algo que vai de 0.180 a 0.481, ou seja dentro desse intervalo não tem nenhum valor de 0, ou negativo, isso quer dizer que no minimo existe uma diferença de 0.18 e no máximo 0.481 para as prevalência, então se pegarmos por exemplo a média de C. talpacoti, que é 0.516 de prevalência, nesse caso o S. squammata tem no minimo 0.516-0.180 e no maximo 0.516-0.481.
Essa conta é esse teste aqui:

prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste

Mas fazemos ele para cada valor do MCMC, por isso temos a cadeiazinha dele e assim uma distribuição.
Vamos olhar esses três valores das quantidades derivadas num gráfico.

02

O primeiro e o segundo, são as distribuições de qual a prevalência, e nossa confiança no valor. O risco vermelho é a media da distribuições, bem próximo das distribuições que tinhamos originalmente certo?

So que aqui também calculamos a distribuição da diferença, e vemos que ela é quase que certamente maior que zero, a gente não vê a distribuição nem encostando no zero, e a diferença ser zero é o mesmo que não haver diferença, se você é zero cm mais alto que eu, eu posso afirmar que temos a mesma altura, mas se você é no minimo três centímetros mais alto que eu, eu não posso afirmar que temos a mesma altura, simples assim.

Então nesse caso podemos afirmar que “Existe uma diferença na prevalência de Haemoproteus nessas duas espécies”, a mais ainda, podemos afirmar que “Existe uma diferença, e ela deve ser algo entre 18% e 34% menos prevalente em S. squammata do que em C. talpacoti“, agora porque isso ocorre, precisamos ler mais do que eu li para escrever esse post.

Bem, para constar, podemos fazer o mesmo teste usando estatística frequentista assim.

> modelo.glm<-glm(cbind(dados$simul, dados$N - dados$simul) ~ dados$Sp, family = binomial) > summary(modelo.glm) Call: glm(formula = cbind(dados$simul, dados$N - dados$simul) ~ dados$Sp, family = binomial) Deviance Residuals: [1] 0 0 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.06454 0.25413 0.254 0.799529 dados$SpScardafella squammata -1.61210 0.43111 -3.739 0.000184 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1.5693e+01 on 1 degrees of freedom Residual deviance: -3.5527e-15 on 0 degrees of freedom AIC: 12.551 Number of Fisher Scoring iterations: 3

Notem que temos os mesmos valores para alpha e beta praticamente, não esperamos nada diferente, já que o prior que usamos foi pouco informativo, mas como fizemos aqui a analise de três formas diferentes, chi-quadrado, glm num framework bayesiano e glm num framework frequentistca, e tudo convergiu para a mesma resposta, existe uma diferença entre a espécies quanto a sua prevalência de Haemoproteus, devemos estar no caminho certo.

Bem acho que acabei falando um monte de coisa aqui, mas acho que agora com o tempo da para ir destrinchando o glm melhor nos próximos posts. Mas glms são muito legais, ja que da para usar, nesse mesmo esquema, muitos tipos de distribuições, cobrindo uma gama muito maior de tipos de dados que com os modelo lineares mas usando todas as estrategias que ja tinhamos para cobrir os modelos lineares, aqueles em que os dados tem uma distribuição normal.

Referencias
Adriano E.A. & Cordeiro N.S. 2001 – Prevalence and Intensity of Haemoproteus columbae in Three Species of Wild Doves from Brazil – Mem. Inst. Oswaldo Cruz 96(2) : 175-178

Marc Kéry 2010 – Introduction to WinBUGS for Ecologists.

set.seed(123)
### Gerando dados
dados<-data.frame(Sp=c("Columbina talpacoti","Scardafella squammata"),N=c(62,57),
                  Prevalence=c(0.516,0.193))
 
dados$simul<-0
 
for(i in 1:2) {
    dados[i,"simul"]<-rbinom(1, dados[i,"N"], prob = dados[i,"Prevalence"])
    }
 
layout(matrix(1:2,ncol=1,nrow=2))
for(i in 1:2) {
    plot(dbinom(x=0:dados$N[i],size=dados$N[i],prob=dados$Prevalence[i])~c(0:dados$N[i]),frame=F,pch=19,
         main=paste(dados$Sp[i]," N=",dados$N[i],sep=""),ylab="Frequência",xlim=c(0,70),xlab="",ylim=c(0,0.15))
    legend("topright",legend=paste("Prevalência =",dados$Prevalence[i]*100,"%"),bty="n")
    lines(c(dados$N[i],dados$N[i]),c(0.05,0),col="blue",lty=2,lwd=3)
    text(dados$N[i]+5,0.05,paste("Limite N =",dados$N[i]),cex=0.7)
    }
 
tabela<- as.table(cbind(dados$simul, dados$N - dados$simul))
dimnames(tabela) <- list(Sp = c("Columbina talpacoti","Scardafella squammata"),Parasita = c("Sim","Não"))
tabela
 
teste.chi <- chisq.test(tabela)
 
str(teste.chi)
teste.chi$observed   # Contagens observadas
teste.chi$expected   # Contagens esperada sobre a hipotese nula
teste.chi$residuals  # Residuo de Pearson
teste.chi$stdres     # Residos "standardized"
 
### Análise usando Bugs
# Definindo Modelo
 
sink("Binomial.t.test.txt")
cat("
model {
 
# Priors
 alpha ~ dnorm(0,0.01)
 beta ~ dnorm(0,0.01)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dbin(p[i], N[i]) 
    logit(p[i]) <- alpha + beta *species[i]
 }
# Derived quantities
 prev.C_talpacoti <- exp(alpha) / (1 + exp(alpha))
 prev.S_squammata <- exp(alpha + beta) / (1 + exp(alpha + beta))
 prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste
}
",fill=TRUE)
sink()
 
# Junte os dados numa lista
bugs.dados <- list(C = dados$simul, N = dados$N, species = c(0,1), n = length(dados$Sp))
 
# Função para iniciar os dados
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))}
 
# Parametros para acompanhar
params <- c("alpha", "beta", "prev.C_talpacoti", "prev.S_squammata", "prev.Diff")
 
# Configurações do MCMC
nc <- 3
ni <- 1200
nb <- 200
nt <- 2
 
# Iniciando o Gibbs sampling
library(R2OpenBUGS)
out <- bugs(data=bugs.dados, inits=inits, parameters.to.save=params,model.file="Binomial.t.test.txt",
            n.thin=nt, n.chains=nc, n.burnin=nb,n.iter=ni)
 
print(out, dig = 3)
 
out$summary
 
#Figura 2
str(out)
par(mfrow = c(3,1))
hist(out$sims.list$prev.C_talpacoti, col = "grey", xlim = c(0,1), xlab="",
	main = "Prevalência de C. talpacoti", breaks = 30,ylim=c(0,300))
abline(v = out$mean$prev.C_talpacoti, lwd = 3, col = "red")
hist(out$sims.list$prev.S_squammata, col = "grey", xlim = c(0,1), xlab= "",
	main = "Prevalência de S. squammata", breaks = 30,ylim=c(0,300))
abline(v = out$mean$prev.S_squammata, lwd = 3, col = "red")
hist(out$sims.list$prev.Diff, col = "grey", xlim = c(0,1), xlab = "",
	main = "Diferença nas Prevalências", breaks = 30,ylim=c(0,300))
abline(v = 0, lwd = 3, col = "red")
 
### Analise usando GLM
modelo.glm<-glm(cbind(dados$simul, dados$N - dados$simul) ~ dados$Sp, family = binomial)
summary(modelo.glm)

Levando em conta a correlação filogenética em uma análise de dados.

A cada ano mais e mais filogenias são publicadas, esses dados estão cada vez mais fáceis de serem obtidos, mas ainda talvez não tão utilizados como deveriam.
No R existe o pacote ape, que traz várias ferramentas para trabalhar com filogenia.
Mas uma em particular eu achei bem legal. E vamos tentar usar ela em um exemplo.

Bem pensar num exemplo, vamos comparar o tamanho de 10 espécies, 5 de uma ilha, 5 de outra ilha e ver se existe diferença entre o tamanho delas.
Então fomos em uma ilha, medimos o tamanho médio de 5 passarinhos e depois fomos a outra ilha e medimos o tamanho de 5 outras espécies de passarinho.

Com os dados em mão vamos tentar demonstrar que o tamanho dessas espécies são diferentes entre as ilhas.
Primeiro fazemos um gráfico dos dados:

01

Hummm, tudo parece bem diferente, ai fazemos o famoso teste T entre dois grupos e lindo, as ilhas são diferentes, é estrelinha de significativo para todo lado.

> summary(lm(tam~pop)) Call: lm(formula = tam ~ pop) Residuals: Min 1Q Median 3Q Max -0.82754 -0.18734 0.07906 0.29716 0.55995 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.3190 0.1994 11.631 2.72e-06 *** popB 1.5201 0.2820 5.391 0.000653 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.4458 on 8 degrees of freedom Multiple R-squared: 0.7842, Adjusted R-squared: 0.7572 F-statistic: 29.06 on 1 and 8 DF, p-value: 0.0006529

Mas nos dias de hoje facilmente temos acesso a filogenia, talvez não tão facilmente assim, mas com um pouco de perseverança sim, então vamos olhar a filogenia dessas espécies em relação ao seus tamanhos em dois exemplos de filogenias diferentes.

02

So para entender melhor essa figura, eu representei a filogenia e o tamanho das bolinhas é relativo ao tamanho da espécie. Eu também coloquei duas cores, uma cor para cada ilha.
A questão é que aqui a gente vê a diferença entre as espécies(tamanho das bolinhas) e vê que nos vários ramos existem tanto bolinhas grandes como pequenas, assim como nos vários ramos existem tanto bolinhas verdes como vermelha. Ou seja muitas oportunidades devem ter existido de colonizar essas ilhas, e apesar de aberta a todas as espécies, alguma coisa deve ter feito somente espécies grandes se sobressair em uma ilha, enquanto na outra ilha somente espécies pequenas. Ja vimos uma diferença entre os tamanhos, mas com a filogenia conseguimos pensar um pouco além, que todo mundo pode ter chegado a ilha, mas a ilha talvez tenha selecionado os tamanhos. Podemos então colocar essa filogenia na analise, e controlar o quanto acreditamos na diferenças entre espécies vem por causa do “ambiente” e o retirando o parentesco (efeito filogenético).
Então levando a filogenia em conta temos o seguinte resultado:

> compar.gee(tam~pop, phy =phylo1) Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27 running glm to get initial regression estimate (Intercept) popB 2.319003 1.520127 Call: compar.gee(formula = tam ~ pop, phy = phylo1) Number of observations: 10 Model: Link: identity Variance to Mean Relation: gaussian QIC: 10.59292 Summary of Residuals: Min 1Q Median 3Q Max -0.7536925 -0.1102418 0.1550745 0.3753477 0.6337977 Coefficients: Estimate S.E. t Pr(T > |t|) (Intercept) 2.245158 0.2963257 7.576656 0.01027012 popB 1.515783 0.3384849 4.478140 0.03350552 Estimated Scale Parameter: 0.2059916 "Phylogenetic" df (dfP): 4.373501

Bem, chegamos a mesma conclusão (p=0.034), existe uma diferença entre o tamanho, não esperaríamos nada diferente disso. Mas agora sabemos que deve ter algo nas ilhas levando a essa diferença.
Agora vamos pensar em outro exemplo:

03

Bem agora comparando esse caso com o primeiro (salvo que minha árvore ta meia mal feita), a questão é que todas as espécies pequenas estão em um ramo, enquanto todas as espécies grandes estão em outro ramo, e essa ainda por cima é exatamente a divisão das ilhas. Bem não é nada sobrenatural, podemos pensar que 2 ancestrais colonizaram aquelas ilhas, um grande foi para uma ilha, um pequeno para outra. O problema é que se simplesmente fazemos um teste T e falamos que existe uma distribuição normal para descrever uma ilha, isso não vai dar certo, porque deve haver uma dependência dessa ancestralidade. Não deve ter como mudar muito em poucas gerações, e agora podemos pensar que essa diferença não deve ter haver com algo na ilha, mas a ancestralidade. Então existe uma forte correlção entre o tamanho das espécies, e deveríamos levar isso em conta no teste T. Podemos fazer a analise seguindo a ideia do Paradis e temos:

> compar.gee(tam~pop, phy =phylo2) Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27 running glm to get initial regression estimate (Intercept) popB 2.319003 1.520127 Call: compar.gee(formula = tam ~ pop, phy = phylo2) Number of observations: 10 Model: Link: identity Variance to Mean Relation: gaussian QIC: 14.66343 Summary of Residuals: Min 1Q Median 3Q Max -0.87844244 -0.19787730 0.05506749 0.29373377 0.50904782 Coefficients: Estimate S.E. t Pr(T > |t|) (Intercept) 2.369908 0.321294 7.376135 0.008810888 popB 1.466297 0.509500 2.877914 0.077381531 Estimated Scale Parameter: 0.2003873 "Phylogenetic" df (dfP): 4.545455

Opa, agora não acreditamos mais que existe uma diferença significativa (p=0.077) entre as ilha, oras é só olhar a filogenia que nossa crença sobre o que esta acontecendo, quais as forças que impulsionam as espécies a serem desse jeito muda totalmente no exemplo 2.
A ancestralidade é muito forte e ela tem que ser levada em conta nos nossos intervalos de confiança. Esse esquema que o Paradis fez foi levar em conta a filogenia, como uma matriz de correlação entre as espécies junto com gee. Um tipo de correção de auto-correlação nos graus de liberdade. Salvo a amostra pequena (n=10 espécies) e a árvore que fiz que ficou meio feia no exemplo 2, a do exemplo 1 foi feita ao acaso, logo não mudou a conclusão, da para entender a idéia.

Essa abordagem pode ajudar muito na interpretação dos dados e enriquecer muito a discussão. E os ingredientes necessários são os dados e uma filogenia, nada mais. Como os dados a gente vai ter, é so procurar ou fazer uma árvore filogenética. Talvez seja difícil fazer uma árvore, mas existem muitas por ai que podem ser usadas. Existem muitas abordagens para incluir filogenia nas análises, mas essa é de fácil acesso (so ter o pacote ape e os ingredientes certo) e pode ser aplicada em qualquer caso que se usaria um glm. As contas podem ser difíceis de imaginar, mas a interpretação pode ser mais fácil que quando não levado em conta a filogenia como deve ter sido nesse exemplo.

Para entender melhor isso tudo da para ler o artigo do Paradis que explica essa analise, além de ele ser uma pessoa bastante acessível e legal na lista de filogenia do R.

Paradis, E. and Claude J. (2002) Analysis of comparative data using generalized estimating equations. Journal of theoretical Biology, 218, 175–185.

########################
#Gerando dados
library(ape)
set.seed(21)
pop<-factor(rep(c("A","B"),each=5))
tam<-c(rnorm(5,2,0.4),rnorm(5,4,0.4))
ind<-1:10
phylo1<-rtree(10,tip.label=ind)
cat("(((1:1,2:1):1,(3:1,4:1):1,5:2):1,((6:1,7:1):1,(8:1,9:1):1,10:1):4);",
file = "ex.tre", sep = "n")
phylo2 <- read.tree("ex.tre")
 
#Grafico e Teste T
boxplot(tam~pop,frame.plot=F,xlab="Locais",ylab="Tamanho")
 
summary(lm(tam~pop))
 
#Grafico e Analise com Exemplo 1
plot.phylo(phylo1,direction="downwards",use.edge.length=F,show.tip.label=F)
points(phylo1$tip.label,rep(0,10),cex=tam,pch=19,col=as.numeric(pop)+1)
title("Exemplo 1")
 
compar.gee(tam~pop, phy =phylo1)
 
#Grafico e analise com Exemplo 2
plot.phylo(phylo2,direction="downwards",use.edge.length=F,show.tip.label=F)
points(phylo2$tip.label,rep(0,10),cex=tam,pch=19,col=as.numeric(pop)+1)
title("Exemplo 2")
 
compar.gee(tam~pop, phy =phylo2)
##################################

Um pequeno exemplo de um modelo com Zeros Inflados

Modelos inflados por zeros podem ser muitos interessantes. Vou tentar fazer uma breve introdução com um exemplo simples.
Vamos supor a seguinte situação, eu gostaria de saber se o número de aranhas presentes numa especie de arbusto A é diferente do número de aranhas comumente encontrado num arbusto B. Utilizando o número magico de coletas que me foi passado na faculdade, eu procurei trinta arbustos da especie A e trinta arbustos da especie B e contei quantas aranhas tinham neles:

#Simulando os dados
set.seed(13)
população<-factor(rep(c("A","B"),each=30))
contagem<-rpois(60,rep(c(2,4),each=30))
#

Vemos que no Arbusto B, em média temos 4 aranhas, enquanto no arbusto A temos 2 aranhas.
Podemos testar essa afirmação através de um modelo de regressão de poisson.

#Modelo de Poisson
modelo1<-glm(contagem~população,family="poisson")
summary(modelo1)
exp(coef(modelo1))
#

E obtemos como resultado:

Call: glm(formula = contagem ~ população, family = "poisson") Deviance Residuals: Min 1Q Median 3Q Max -1.9833 -0.8255 0.0237 0.5868 2.9973 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.6763 0.1302 5.195 2.05e-07 *** populaçãoB 0.6587 0.1604 4.107 4.01e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 90.678 on 59 degrees of freedom Residual deviance: 72.885 on 58 degrees of freedom AIC: 235.15 Number of Fisher Scoring iterations: 5 > exp(coef(modelo1)) (Intercept) populaçãoB 1.966667 1.932203

Vejam só, o modelo nos da varias estrelinhas, nossa afirmação foi na mosca, mas mais importante que isso, temos estimativas do tamanho da população, e em média elas ficam bem próxima de realidade, o intercepto que representa a população nos arbustos A tem 1.966.. aranhas, bem próximo do 2 que usamos pra simular os dados, e no arbusto B temos 1.966 + 1.932.., que da praticamente 4, o numero que usamos na simulação. Muito legal, muitos p significativos e condecorados com variasssss estrelinhas.
Mas nem sempre o mundo é tão simples assim. Como nosso amigo Levins propôs uma vez, nem sempre todos os arbustos vão ter aranhas, eles podem não estar ocupados por aranhas.
Pra gente contar aranhas em um arbusto, tem que ter mais de uma aranha la , ou a gente vai contar zero aranhas, até ai tudo bem pois alguns arbustos podem ser um péssimo lugar pra se viver e nenhuma aranha fica la, ou seja zero aranhas.
Mas podem existir um outro tipo de zero, o zero de aranhas devido a elas nunca terem chegado la, o arbusto pode ser perfeito pra abrigar aranhas, mas nenhuma chegou la. Agora vamos colocar isso nos dados e ver o que acontece com nosso modelo.

#Simulando ocupação
ocupação<-rbinom(60,size=1,prob=0.60)
contagem.ocup<-contagem*ocupação
 
#Modelo 2 de poisson
modelo2<-glm(contagem.ocup~população,family="poisson")
summary(modelo2)
exp(coef(modelo2))
#

Call: glm(formula = contagem.ocup ~ população, family = "poisson") Deviance Residuals: Min 1Q Median 3Q Max -2.0656 -1.8719 -0.3133 0.6582 4.2838 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.4906 0.1429 3.434 0.000594 *** populaçãoB 0.2671 0.1898 1.407 0.159458 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 160.14 on 59 degrees of freedom Residual deviance: 158.14 on 58 degrees of freedom AIC: 264.01 Number of Fisher Scoring iterations: 6 > exp(coef(modelo2)) (Intercept) populaçãoB 1.633333 1.306122

É muito triste, mas o fato de nem todos os arbustos estarem ocupados por aranhas nos tira muitas muitas estrelinhas, mas pior que isso, agora concluímos que os arbustos A e B tem o mesmo numero de aranhas, isso é mal, muito mal. Olha a estimativa do intercepto 1.63, usamos 2 na simulação, mas os zeros adicionais puxaram a media do numero de aranhas nos arbustos pra baixo.
Mas o que acontece é que nesse caso não estamos levando em conta que temos 2 tipos de zeros nas coletas, existem zeros de o arbusto ter zero aranhas, mas existem zeros devido ao arbusto não estar ocupado, aranhas nunca chegaram la.
O que deveríamos fazer é tentar fazer um modelo onde esses dois tipos de situação são levados em conta.
Temos de declarar ao modelo que existem zeros devido a uma chance de ocupação que cada arbusto tem, ai considerando os arbustos ocupados que deveríamos fazer as contas para o numero de aranhas no arbusto. Essa é a solução de modelos com zero inflados tenta trazer.

#modelo zero inflado usando o pacote pscl
library(pscl)
modelo3<-zeroinfl(contagem.ocup~população|1)
summary(modelo3)
exp(coef(modelo3)[1:2])
plogis(coef(modelo3)[3])
#
Call: zeroinfl(formula = contagem.ocup ~ população | 1) Pearson residuals: Min 1Q Median 3Q Max -0.9876 -0.8934 -0.2306 0.4150 3.2704 Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.8296 0.1540 5.388 7.13e-08 *** populaçãoB 0.5977 0.2019 2.960 0.00307 ** Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.4567 0.2829 -1.615 0.106 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Number of iterations in BFGS optimization: 8 Log-likelihood: -104.9 on 3 Df > exp(coef(modelo3)[1:2]) count_(Intercept) count_populaçãoB 2.292305 1.817889 > plogis(coef(modelo3)[3]) zero_(Intercept) 0.3877701

A primeira boa noticia é que nossas estrelinhas são devolvidas, temos muitos p significativos. Mas vamos olhar as nossas estimativas, agora estimamos para o arbusto A uma média de 2.29 e para o B 2.29 + 1.81, o que é bem próximo de 4, 2 e 4 foram exatamente os valores que usamos.
Olhamos ainda a estimativa de ocupação, é de 38%, tudo bem que usamos 60%, mas é bastante razoável, e temos que levar em consideração o número de amostras, mas voltamos a concluir certo que arbustos A tem em média 2 aranhas, e arbusto B o dobro.
A parte mais interessante, é que apesar de dois processos estarem acontecendo simultaniamente, um que é o ambiente influenciando no número de aranhas que encontramos e o outro é arbustos estarem ocupados ou não, estamos modelando esses dois processos relativamente bem. Resgatando estimativas razoáveis. Se ignoramos um processo, como ignoramos o processo de ocupação no modelo 2, além de péssimas estimativas do tamanho populacional, podemos concluir erroneamente o que esta acontecendo.
Mas ao pensarmos sobre todos os processos que influenciam as nossas observações, podemos selecionar uma estratégia que leve tudo em conta, todos os processos que temos teoria para subsidiar e gerar boas estimativas.
Nesse caso, pensando que o ambiente influencia o número de aranhas que encontramos, mas encontrar aranhas também depende da ocupação, como prevê a teoria de metapopulações, uma regressão de poisson com zero inflados funciona perfeitamente. E o preço do modelo errado pode ser varias estrelinhas pior, uma conclusão errada do que está acontecendo.