Ferramentas do usuário

Ferramentas do site


Barra lateral

CMQ
Centro de Métodos Quantitativos


USP ESALQ
Depto. de Ciências Florestais
ESALQ
UNIVERSIDADE de SÃO PAULO
Av. Pádua Dias, 11
Caixa Postal 09
13418-900 - Piracicaba - SP
BRASIL
biometria:verossim:06-model


7. Modelos com distribuição binomial e Poisson


Conceitos

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

Tutoriais

Modelos Poisson: Fecundidade de Solanum sp.

A planta Solanum sp. é de grande importância econômica. Neste tutorial vamos simular que o número de frutos nas plantas adultas é variável Poisson, com parâmetro $$\lambda$$ fixo ou dependente de co-variáveis.

Em seguida ajustaremos modelos Poisson a estes dados simulados.

Simulando e estimando o lambda de um modelo Poisson sem variáveis preditoras

Uma forma de conferir se a estimação dos parâmetros está correta é simular os dados e verificar se os parâmetros utilizados na simulação concordam com os parâmetros obtidos mediante máxima verosimilhança.

Passos iniciais

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

library(bbmle)
library(car)
library(MASS)

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

set.seed(1234)
Começamos a simulação:
nplan<-1000
lambda=5
frutos<-rpois(nplan,lambda)

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

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)

Adicionamos os valores teóricos esperados:

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

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

Aplicamos a função sobre uma série de valores:

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

Determinamos o valor de máxima verossimilhança ou equivalentemente os valores mínimos para função de log-verossimilhança negativa por força bruta (procurando o valor minimo da função de verossimilhança negativa em uma serie de valores) e por métodos de otimização numérica usando a função mle2.

minLLest<-min(LLest)
lambdaLL.fb<-xvec[ LLest == min(LLest)]
lambdaLL.nm<-mle2(poisNLL, start=list(lambda=4))

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

lambdaLL.fb
lambdaLL.nm
mean(frutos)

Plotamos 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

Em um cenário mais realista, a fecundidade das plantas poderia incrementar-se com um aumento da concentração de nutrientes no solo, por exemplo o fósforo.

Pensemos que as plantas se encontram em solos com diferentes graus de concentração de fósforo (mg/kg). Neste caso, lambda deixará de ser um valor fixo e será determinado em função pela variável preditora fósforo, mediante a função exponencial,

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

Onde $$X_i$$ são as variáveis preditoras. A vantagem desta transformação é que independentemente dos valores da variável preditora, a variável de resposta será sempre positiva, sendo compatível com os valores que o parâmetro $$\lambda$$ da distribuição Poisson pode ter.

Estimar os parâmetros usando esta transformação exponencial do parêmetro é equivalente a analisar os dados com um modelo linear generalizado de distribuição Poisson e função de ligação log.

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, que seguirá uma relação exponencial:

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, que simulam um processo Poisson com valor esperado que é uma função do nível de fósforo:

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

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.

Pergunta: Você tem alguma ideia de como obter uma estimativa inicial dos parâmetros?

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

mod.pois<-mle2(poisglmNLL, start= list(a=2.5,b= 0.33))

Conferimos os parâmetros estimados:

summary(mod.pois)

E os perfis de verossimilhança das estimativas:

mod.pois.prof <- profile(mod.pois)
par(mfrow=c(1,2))
plot.profmle(mod.pois.prof)
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$$ 2). 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 é:

$$σ^2= μ + μ^2/k$$

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

med<-mean(x)
vari<-var(x)
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.

Função glm no R

Uma regressão de uma variável Poisson como a que fizemos equivale a um glm 3) 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 a variável dependente.

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.

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 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 estimar o parâmetro $$p$$ nos modelos de distribuição binomial.

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

Otimizamos a função de verossimilhança e obtemos os intervalos de confiança dos parâmetros. Como estimamos os valores iniciais dos parâmetros? Uma pista: tem a ver com plotar a razão dos frutos murchos/total dos frutos na escala logit. A função logit lineariza a curva logística. A equação que transforma as proporções (frutos murchos/total frutos) para escala logit é ln(p/(1-p).

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

Otimizamos:

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

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))
plot.profmle(profile(mod.bin))
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))

Compare os resultados usando a função padrão no R, glm com distribuição binomial e função de ligação logit:

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

Exercício

Baixe o arquivo besouro.csv. Os dados descrevem o efeito de diferentes concentrações de insecticida na mortalidade dos besouros expostos à substância. Estime a probabilidade da mortalidade em função do nível do inseticida, plote a função logistica prevista e estime os intervalos de verossimilhança dos parâmetros.

Resumos da aula

Pesquisa

Indique os itens do tutorial e do texto prioritários para a discussão aqui

1)
já que a média amostral equivale à solução analítica para o MLE de $$\lambda$$ no caso da distribuição Poisson
2)
z = raiz quadrada da deviance, detalhes na ajuda da função profile
biometria/verossim/06-model.txt · Última modificação: 2015/08/10 20:48 (edição externa)