Ferramentas do usuário

Ferramentas do site


05-binomial-poisson:05-binomial-poisson

5. Modelos Binomial e Poisson


Conceitos

  • Modelos não-gaussianos
  • Modelos lineares generalizados
  • Função de ligação
  • Função logística, probito e logito
  • Sobredispersão

Tutoriais

Modelos Poisson: Fecundidade de Solanum sp.

Neste tutorial vamos usar modelos Poisson para descrever a distribuição de número de frutos em plantas adultas de Solanum sp.. Mas ao invés de dados reais, criaremos os valores de uma distribuição Poisson. A vantagem da simulação é que podemos avaliar o desempenho dos modelos, pois conhecemos os valores dos parâmetros.

Modelo Poisson sem variáveis preditoras

Passos iniciais

Iniciamos o R e carregamos os pacote bbmle, car e sads1). O pacote bbmle foi criado por Ben Bolker e tem algoritmos de optimização para obter a máxima verosimilhança dos parâmetros. Se você ainda não tem algum dos pacotes instalado use o comando install.packages (“pacote”).

library(bbmle)
library(car)
library(sads)

Fixamos o gerador pseudo-aleatório de dados para todos obtermos exatamente o mesmos resultados:

set.seed(1234)
Começamos a simulação:

Vamos sortear 1000 valores de uma distribuição Poisson com parâmetro $\lambda=5$, o que simula uma população de mil plantas com uma média de cinco frutos por indivíduo. Ao usarmos uma Poisson simulamos que os frutos ocorrem a uma taxa constante por planta, e que a cada fruto é um evento independente dos demais.

nplan<-1000
lambda=5
frutos<-rpois(nplan,lambda)
Análise exploratória

Qual é a probabilidade teórica de encontrar uma planta sexualmente madura sem frutos?

dpois(0,lambda)

Compare com a proporção de plantas sem frutos na amostra simulada:

sum(frutos==0)/length(frutos)

Para fazer esta avaliação para todos os valores, comparamos a distribuição de probabilidade teórica com as proporções empíricas. Comece com um gráfico de barras das proporções:

mfrut<-max(frutos)
fa<-factor(frutos, levels=0:mfrut)
prob.obs<-table(fa)/nplan
par(las=1)
plot(0:mfrut,prob.obs, xlab="Numero de frutos", 
     ylab="Probabilidade", type="h", lwd=5)

E adicione os pontos com os valores teóricos esperados:

prob.tr<-dpois(0:mfrut, lambda)
points(0:mfrut,prob.tr, pch=21, col="red")
Ajuste do modelo

Definimos a função de máxima verossimilhança do parâmetro lambda da distribuição Poisson:

x<-frutos
poisNLL<-function(lambda){
-sum(dpois(x, lambda, log=TRUE))
}

Há um teorema que garante que essa função de verossimilhança é unimodal e tem um único mínimo. Por isso podemos buscar a estimativa de máxima verossimilhança de $\lambda$ simplesmente experimentando diversos valores até encontrar um mínimo. Para isso criamos uma sequência de valores na vizinhança da média amostral e aplicamos a cada um a função de log-verossimilhança:

xvec<-seq(4.85,5.3, length=1000)
LLest<-sapply(xvec,poisNLL)

Em seguida buscamos o valor mínimo:

minLLest<-min(LLest)
lambdaLL.fb<-xvec[ LLest == min(LLest)]

Vamos comparar as estimativas do lambda obtidas com os métodos de força bruta, numéricos e analíticos2).

Faça a minimização numérica com a função mle2:

lambdaLL.nm<-mle2(poisNLL, start=list(lambda=4))

E agora compare os três valores:

lambdaLL.fb
lambdaLL.nm
mean(frutos)

E vamos fazer a comparação gráfica, plotando a função de verosimilhança e os pontos obtidos:

mfrutos<-mean(frutos)
LLest2<-LLest-min(LLest)
plot(xvec,LLest2, typ="l", xlab="frutos", ylab="loglik")
abline(v=mfrutos, col="blue", lwd = 3)
abline(v=coef(lambdaLL.nm),col ="darkgray")
abline(v=lambdaLL.fb, col="red")

Modelos Poisson com variáveis preditoras

Vamos supor que a fecundidade das plantas poderia incrementar-se com um aumento da concentração de fósforo no solo. Simulamos esta situação fazendo da taxa de frutos por planta uma função de variáveis preditoras, como o fósforo:

$\lambda \ = \ e^{a_0+a_iX_i}$

O que é o mesmo que

$$\ln(\lambda) \ = \ a_0 + a_iX_i$$

Ou seja, o logarítmo de $\lambda$ é uma função linear das variáveis preditoras $X_i$3). Essa é a funcão de ligação logarítmica. Sua vantagem é preservar a relação linear com as preditoras, e garantir valores positivos de $\lambda$, como deve ser na distribuição Poisson.

Então vamos simular a dependência de $\lambda$ ao fósforo! Primeiro definimos valores dos parâmetros:

set.seed(1234)
phos<-runif(100,0,10)
a= 1
b= 0.3
x<-phos

E então os valores esperados de frutos em função do nível de fósforo, conforme a função de ligação:

ydet<-exp(a+b*x)
plot(x,ydet)

Agora para simular um processo Poisson sorteamos amostras Poisson com parâmetro lambda igual a estes valores esperados:

fec<-rpois(100,ydet)

E plotamos nossos dados:

par(las=1)
plot(phos, fec, xlab="Fósforo mg/Kg", ylab="Número de frutos")

Definimos a função de verosimilhança para este modelo:

poisglmNLL = function(a,b) {
   ypred= exp(a+b*x)
  -sum(dpois(fec,lambda=ypred, log=TRUE))
}

Agora que simulamos um processo Poisson com valor esperado que é uma função do nível de fósforo, vamos ver o o ajuste de um modelo que descrevesse exatamente isso. Para estimar a máxima verosimilhança temos que fornecer na função de otimização um valor inicial para os parâmetros a e b.

chute <- list(a=2.5, b= 0.33)

Pergunta

Definir bons pontos de partida (starting values) para o otimizador é crucial, especialmente para funções de verossimilhança mais complexas, com muitas preditoras. Aqui o caso é simples e pouco sensível aos valores iniciais. Por isso “chutamos” um valor razoável. Mas mesmo neste caso há maneiras de obter bons palpites e ajudar o otimizador. Você tem alguma ideia de como obter uma estimativa inicial dos parâmetros?

E então otimizamos a função de verossimilhança, com a função mle2:

mod.pois<-mle2(poisglmNLL, start= chute)

Conferimos o resumo do modelo:

summary(mod.pois)

E os perfis de verossimilhança das estimativas:

mod.pois.prof <- profile(mod.pois)
par(mfrow=c(1,2))
plotprofmle(mod.pois.prof) ## carregue o pacote sads para ter esta função
par(mfrow=c(1,1))

Numa abordagem frequentista, os perfis são usados para estimar intervalos de confiança dos parâmetros, transformando a log-verossimilhança negativa em uma variável normal-padrão $z$ 4). Esta aproximação à normal é válida se os perfis são parábolas, como parece ser o caso.

Aplicando a função plot ao objeto da classe perfil você obtem perfis na escala desta transformação:

plot(mod.pois.prof)

As linhas vermelhas neste gráfico marcam no eixo $x$ os quantis correspondentes a diferentes probabilidades acumuladas da normal. Estes quantis são usados como os limites dos intervalos de confiança para estas probabilidades. Compare os valores na figura para probabilidade de 95% com os obtidos com o comando:

confint(mod.pois)

Por fim, plotamos as curvas dos valores esperados de frutos por planta com os parâmetros e suas estimativas:

par(las=1)
plot(phos,fec, xlab="Fósforo mg/Kg", ylab="Número de frutos" )
a.est<-coef(mod.pois)[1]
b.est<-coef(mod.pois)[2]
curve(exp(a+b*x),add=TRUE, col="red")
curve(exp(a.est +b.est*x), add=TRUE, col="blue", lty=2) #estimada
legend("topleft", c("Parâmetro","Estimativa"),col=c("red","blue"), lty=c(1,2))

O que aconteceria se usássemos outro modelo para descrever estes dados? Neste caso sabemos qual é o modelo correto, mas vamos simular esta situação, imaginando que o pesquisador está experimentando diferentes modelos. Usaremos a função de verosimilhança de uma binomial negativa, que permite agregação. Vamos calcular a máxima verosimilhança e comparar os modelos com o critério de informação de Akaike (AIC).

negbinNLL<- function(a,b,k){
	ypred<-exp(a+b*x)
	-sum(dnbinom(fec, mu=ypred, size=k, log = TRUE))
        }

Tentamos um valor inicial de k usando o método dos momentos. Como sabemos que a variância da binomial negativa é:

$$\sigma^2= \mu + \frac{\mu^2}{k}$$

Então podemos estimar um valor aproximado de k com:

med<-mean(fec)
vari<-var(fec)
k.init <-med^2/(vari-med)
k.init

mod.negbin<- mle2(negbinNLL, start=list(a=2.5, b= 0.33, k=k.init))

Conferimos os parâmetros estimados:

summary(mod.negbin)

Comparamos os modelos com a função AICtab e confira o grau de suporte relativo do modelo Poisson respeito a binomial negativa:

AICtab(mod.pois,mod.negbin, delta=T, sort=T, weights = TRUE)

O modelo com melhor suporte é o modelo com distribuição Poisson, porém a diferença nos valores do AIC é menor que 2.

Modelos com distribuição binomial: Infestação por Fusarium sp.

Agora pensemos que os frutos das mesmas plantas solanáceas são atacados pelo fungo Fusarium sp., murchando-os. Pensemos que esse fungo habita no solo e portanto só ataca os frutos quando eles amadurecem e caem.

A probabilidade de infestação dos frutos pode ser constante ou depender de algum fator como por exemplo o grau de umidade relativa do micro-habitat onde se encontra a planta. O número de “tentativas” (parâmetro $N$ da binomial) será variável entre as plantas, pois corresponde ao número de frutos que ela tenha. Vamos simular esta situação e depois verificar como um modelo binomial se sai em descrevê-la.

Primeiro descartamos as plantas sem frutos:

fec1<-fec[fec!=0]
num.plant<-length(fec1)

E simulamos o número de frutos atacados por planta:

Vamos fazer o grau de umidade do solo na vizinhança da planta variar entre 30 a 100%, com probabilidade uniforme.

set.seed(4444)
ur<-runif(num.plant, 0.3, 1)

Usamos a função de ligação logística para determinar e posteriormente simular o grau de infestação do fungo de acordo a umidade relativa do ambiente. Esta função estabelece seguinte relação do parâmetro $p$ com as variáveis preditoras $X_i$

$$p \ = \frac{e^{a_0+a_iX_i}}{1+e^{a_0+a_iX_i}}$$

A logística restringe os resultados esperados ao intervalo entre 0 a 1, mesmo que a variável preditora tenha valores muito grandes ou pequenos, por isso é adequada para descrever o parâmetro $p$ nos modelos de distribuição binomial. Assim como a função de ligação log para processos Poisson, a logística “preserva” uma relação linear das preditoras com a resposta. Mais formalmente, a ligação faz com que uma transformação do parâmetro seja uma função linear das variáveis preditoras. No caso do modelo binomial, aplicar a logística acima é o mesmo que transformar $p$ para seu logito:

$$\ln \frac{p}{1-p} \ = \ a_0+a_iX_i $$

E novamente encontramos uma transformação do parâmetro -o logito $\ln(p/(1-p))$- que pode ser expressa como uma função linear das preditoras sem violar pressupostos da distribuição usada no modelo, no caso a binomial5).

Vamos então usar tudo isso para construir a simulação de efeito das preditoras sobre a proporção de frutos murchos:

a<- -4
b<- 7
ydet<-exp(a+b*ur)/(1 + exp(a +b*ur))
plot(ur,ydet)
frut.mur<-rbinom(num.plant,size=fec1,prob= ydet)

O número de frutos intactos é o complemento dos murchos.

frut.sau<-fec1-frut.mur 

Para uma avaliação visual, plotamos os fruto sãos e murchos por planta e arranjamos as plantas em ordem crescente exposição à umidade relativa:

datos<-data.frame(frut.sau,frut.mur,ur)
datos.ur<-datos[order(datos$ur),]
conjun<-datos.ur[,1:2]
conjun.m<-as.matrix(conjun)

barplot(t(conjun.m),beside=F, col=c("darkgray","darkblue"), ylim=c(0,80), ylab= "Número de Frutos",
xlab= "Plantas", cex.names=0.4)
legend("topright", c("frutos intactos", "frutos murchos"), col = c("darkgray","darkblue"), pch=c(19,19))
arrows(15.8,62,100,62)
text(60,64, "Umidade", cex=1.8)

Definimos a função de verossimilhança:

binomNLL<- function(a,b){
	prob.det=exp(a+b*ur)/(1 +exp(a+b*ur))
	-sum(dbinom(k,size=N,prob=prob.det, log=TRUE))
	}

E agora vamos encontrar o mínimo da função de log-verossimilhança negativa com otimização computacional com a função mle2.

Otimizamos:

mod.bin<-mle2(binomNLL, start=list(a =-3,b =9), data =list(N= fec1, k = frut.mur))

Pergunta

Como no tutorial anterior, temos que definir valores dos parâmetros para o otimizador iniciar. O comando acima usou um “chute informado”. Descubra como chegamos a ele. Uma pista: tem a ver com plotar a razão dos frutos murchos/total dos frutos na escala de logitos. Lembre-se que a função logito lineariza a curva logística.

Alternativamente podemos usar a interface de fórmula da função mle2 para fazer os dois passos simultaneamente:

mod.bin2<-mle2(frut.mur~dbinom(prob = exp(a + b*ur)/(1 + exp(a + b*ur)), size = fec1),
start = list(a= -3, b = 9))

Avaliamos os perfis de verossimilhança:

par(mfrow=c(1,2))
plotprofmle(profile(mod.bin)) ## carregue o pacote sads para ter esta função
par(mfrow=c(1,1))

Finalmente comparamos a função logística simulada com a estimada. No gráfico, o tamanho do ponto e proporcional ao número de frutos que a planta tinha.

prop.mur<-frut.mur/fec1 #Proporção de plantas murchas
par(las=1)
plot(ur, prop.mur,cex=fec1/25, xlab="Umidade relativa", ylab="Proporção Frutos Infectados")
a.esti<-coef(mod.bin)[1]
b.esti<-coef(mod.bin)[2]
curve(exp(a+b*x)/(1+exp(a+b*x)),add=TRUE, col="red")
curve(exp(a.esti+b.esti*x)/(1+ exp(a.esti +b.esti*x)), add=TRUE, col="blue", lty=2) #estimada
legend("topleft", c("Parâmetro","Estimativa"),col=c("red","blue"), lty=c(1,2))

GLMs

Para modelos Poisson e binomiais como os desses tutorial há uma teoria específica e um método de ajuste mais eficiente. São os glm (generalized linear models).

Compare os resultados usando a função padrão no R para esses modelos, glm com os resultados obtidos. Os principais argumentos são family, que define a distribuição, e neste o argumento 'link', para definir a função de ligação. O padrão de link para Poisson é a logarítmica e a para binomial a logística, como usamos em ajustes numéricos. Para saber mais consulte a ajuda da função e as leituras indicadas.

glm Poisson

Uma regressão de uma variável Poisson como a que fizemos equivale a um glm com distribuição Poisson e função de ligação log. Compare os resultados do modelo Poisson que ajustamos com o mle2 com os resultados obtidos com a função glm do R:

glm.pois <- glm(fec~phos,family=poisson(link="log"))
summary(glm.pois)
summary(mod.pois)
AIC(mod.pois)
AIC(glm.pois)

Cada vez que fazemos um glm com distribuição Poisson e função de ligação log, compramos a premissa de uma relação exponencial entre a variável preditora e o valor esperado da variável dependente. Isso é o mesmo que supor uma relação linear entre a variável preditora e o logarítmo do valor esperado da dependente.

glm Binomial

Uma regressão de uma variável binomial como a que fizemos equivale a um glm com distribuição binomial e função de ligação logística. Compare os resultados do modelo binomial que ajustamos com o mle2 com os resultados obtidos com a função glm do R. Note que para esta função a variável-resposta, do lado esquerdo da fórmula do modelo, é uma matriz com números de observações poisitivas e negativas.

glm.bin <- glm(cbind(frut.mur,frut.sau)~ur, 
               data=datos,family=binomial(link="logit"))
summary(mod.bin)
summary(glm.bin)               




Questões para discussão

Predador com saciação

Em um experimento, girinos em diferentes quantidades ficaram expostos a um predador por um certo tempo. Os dados são o número inicial de girinos expostos e o número predado. Use modelos de regressão binomial para testar se as observações dão suporte à hipótese de que o número final predado é afetado pelo número inicial.

DICAS:

  1. Comece inspecionando gráficos do número de (a) girinos predados pelo número inicial e (b) da proporção de predados em relação ao número inicial.
  2. Quando está abaixo da saciação o predador não ncessariamente captura 100% das preasas disponíveis. Pode ser que ele capture uma proporção fixa das presas. A função de ligação logística pode ser generalizada para qualquer assíntota $M$:

$$p \ = \ \frac{M e^{a + bX}}{1+ e^{a + bX}}$$

Dados

Número de parasitas por hospedeiro

Tomou-se uma amostra de 120 peixes, dos quais contaram-se todos os ectoparasitas. De cada peixe tomou-se o comprimento e registrou-se o sexo. Use as funções glm e glm.nb para investigar a relação entre número de parasitas e atributos dos peixes.

DICA:

Comece inspecionando os gráficos das relações entre a variável resposta e cada uma das preditoras.

Dados

Modelando média e dispersão

Use o conjunto de dados a seguir para construir e comparar modelos de regressão Poisson e binomial negativa. Avalie todas as hipóteses de efeito linear da variável preditora com função de ligação log.

DICA:

Comece enumerando os modelos possíveis. Para a regressão binomial negativa o efeito pode ser sobre a média e/ou o parâmetro de dispersão. Para cada distribuição há os modelos de ausência de efeitos.

Dados

Exercícios

Os exercícios desta unidade estão disponíveis no sistema notaR.

Recursos para Estudo

Leituras

Principais

  • Bolker, B.M. 2008 Ecological Models and Data in R Princeton: Princeton University Press:
    • Cap.6 - Likelihood and all that, seção 6.3.1
    • Cap.9 - Standard Statistics Revisited, seção 9.4

Complementares

  • Crawley, M.J. 2007. The R Book. (caps.13,14, e 16)
Para saber mais sobre glms
  • Dobson, A.J. 1990. An Introduction to Generalized Linear Models. London: Chapman and Hall.
  • McCullagh P. & Nelder, J.A. 1989. Generalized Linear Models. London: Chapman and Hall.

Resumo da Aula

1)
para a funcao plotprofmle
2)
já que a média amostral equivale à solução analítica para o MLE de $\lambda$ no caso da distribuição Poisson
3)
no caso temos uma variável preditora, mas isso pode ser generalizado para quantas você quiser
4)
z = raiz quadrada da deviance, detalhes na ajuda da função profile
5)
funções de ligação são convenientes por outras razões. Se esses detalhes interessarem, veja as leituras recomendadas.
05-binomial-poisson/05-binomial-poisson.txt · Última modificação: 2016/09/23 12:31 por paulo