Ferramentas do usuário

Ferramentas do site


biometria:verossim:06-model

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.


biometria:verossim:06-model [2022/11/24 14:21] (atual) – criada - edição externa 127.0.0.1
Linha 1: Linha 1:
 +
 +
 +^[[start|Modelos Estatísticos: Abordagem da Verossimilhança]]^
 +
 +\\
 +<html>
 +      <font face="Times New Roman" size="5" align="center">
 +       7. Modelos com distribuição binomial e Poisson
 +      </font>
 +</html>
 +
 +\\
 +
 +
 +
 +
 +====== 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").
 +
 +<code>
 +library(bbmle)
 +library(car)
 +library(MASS)
 +</code>
 +
 +Fixamos o gerador pseudo-aleatório de dados para todos obter exatamente o mesmos resultados:
 +<code>
 +set.seed(1234)
 +</code>
 +
 +==Começamos a simulação:==
 +<code>
 +nplan<-1000
 +lambda=5
 +frutos<-rpois(nplan,lambda)
 +</code>
 +
 +Qual é a probabilidade teórica de encontrar uma planta sexualmente madura sem frutos?
 +<code>
 +dpois(0,lambda)
 +</code>
 +Compare com a proporção de plantas sem frutos na amostra simulada:
 +<code>
 +sum(frutos==0)/length(frutos)
 +</code>
 +Para fazer esta avaliação para todos os valores, comparamos a distribuição de probabilidade teórica com as proporções 
 +empíricas e:
 +<code>
 +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)
 +</code>
 +Adicionamos os valores teóricos esperados:
 +<code>
 +prob.tr<-dpois(0:mfrut, lambda)
 +points(0:mfrut,prob.tr, pch=21, col="red")
 +</code>
 +
 +Definimos a função de máxima verossimilhança do parâmetro lambda da distribuição Poisson:
 +
 +<code>
 +x<-frutos
 +poisNLL<-function(lambda){
 +-sum(dpois(x, lambda, log=TRUE))
 +}
 +</code>
 +
 +Aplicamos a função sobre uma série de valores:
 +<code>
 +xvec<-seq(4.85,5.3, length=1000)
 +LLest<-sapply(xvec,poisNLL)
 +</code>
 +
 +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''.
 +
 +<code>
 +minLLest<-min(LLest)
 +lambdaLL.fb<-xvec[ LLest == min(LLest)]
 +lambdaLL.nm<-mle2(poisNLL, start=list(lambda=4))
 +</code>
 +
 +Vamos comparar as estimativas do lambda obtidas com os métodos de força bruta, numéricos e analíticos((já que a média amostral equivale 
 +à solução analítica para o MLE de $$\lambda$$ no caso da distribuição Poisson)):
 +
 +<code>
 +lambdaLL.fb
 +lambdaLL.nm
 +mean(frutos)
 +</code>
 +
 +Plotamos a função de verosimilhança e os pontos obtidos:
 +
 +<code>
 +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")
 +</code>
 +
 +
 +==== 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:
 +
 +<code>
 +set.seed(1234)
 +phos<-runif(100,0,10)
 +a= 1
 +b= 0.3
 +x<-phos
 +</code>
 +E então os valores esperados de frutos em função do nível de fósforo, que seguirá uma relação exponencial:
 +<code>
 +ydet<-exp(a+b*x)
 +plot(x,ydet)
 +</code>
 +Agora para simular um processo Poisson sorteamos amostras Poisson com parâmetro lambda igual a estes valores esperados:
 +<code>
 +fec<-rpois(100,ydet)
 +</code>
 +E plotamos nossos dados, que simulam um processo Poisson com valor esperado que é uma função do nível de fósforo:
 +<code>
 +par(las=1)
 +plot(phos, fec, xlab="Fósforo mg/Kg", ylab="Número de frutos")
 +</code>
 +Definimos a função de verosimilhança para este modelo: 
 +<code>
 +poisglmNLL = function(a,b) {
 +   ypred= exp(a+b*x)
 +  -sum(dpois(fec,lambda=ypred, log=TRUE))
 +}
 +</code>
 +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:
 +<code>
 +mod.pois<-mle2(poisglmNLL, start= list(a=2.5,b= 0.33))
 +</code>
 +Conferimos os parâmetros estimados:
 +<code>
 +summary(mod.pois)
 +</code>
 +
 +E os perfis de verossimilhança das estimativas:
 +<code>
 +mod.pois.prof <- profile(mod.pois)
 +par(mfrow=c(1,2))
 +plot.profmle(mod.pois.prof)
 +par(mfrow=c(1,1))
 +</code>
 +
 +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$$ ((z = raiz quadrada da deviance, detalhes na [[http://finzi.psych.upenn.edu/R/library/stats4/html/profile.mle-class.html|ajuda da função profile]])). 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: 
 +<code>
 +plot(mod.pois.prof)
 +</code>
 +
 +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:
 +
 +<code>
 +confint(mod.pois)
 +</code>
 +
 +
 +Por fim, plotamos as curvas dos valores esperados de frutos por planta com os parâmetros e suas estimativas:
 +<code>
 +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))
 +</code>
 +
 +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).
 +<code>
 +negbinNLL<- function(a,b,k){
 + ypred<-exp(a+b*x)
 + -sum(dnbinom(fec, mu=ypred, size=k, log = TRUE))
 +        }
 +</code>
 +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:
 +<code>
 +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))
 +</code>
 +Conferimos os parâmetros estimados:
 +<code>
 +summary(mod.negbin)
 +</code>
 +Comparamos os modelos com a função AICtab e confira o grau de suporte relativo do modelo Poisson respeito a binomial 
 +negativa:
 +<code>
 +AICtab(mod.pois,mod.negbin, delta=T, sort=T, weights = TRUE)
 +</code>
 +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// (([[http://en.wikipedia.org/wiki/Generalized_linear_model|generalized linear model]])) 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 [[http://finzi.psych.upenn.edu/R/library/stats/html/glm.html|glm]] do R:
 +
 +<code>
 +glm.pois <- glm(fec~phos,family=poisson(link="log"))
 +summary(glm.pois)
 +summary(mod.pois)
 +AIC(mod.pois)
 +AIC(glm.pois)
 +</code> 
 +
 +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:
 +<code>
 +fec1<-fec[fec!=0]
 +num.plant<-length(fec1)
 +</code>
 +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.
 +<code>
 +set.seed(4444)
 +ur<-runif(num.plant, 0.3, 1)
 +</code>
 +
 +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.
 +<code>
 +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)
 +</code>
 +O número de frutos intactos é o complemento dos murchos.
 +<code>
 +frut.sau<-fec1-frut.mur 
 +</code>
 +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:
 +<code>
 +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)
 +</code>
 +
 +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:
 +<code>
 +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))
 + }
 +</code>
 +
 +Otimizamos:
 +<code>
 +mod.bin<-mle2(binomNLL, start=list(a =-3,b =9), data =list(N= fec1, k = frut.mur))
 +</code>
 +
 +Alternativamente podemos usar a interface de fórmula da função ''mle2'' para fazer os dois passos simultaneamente:
 +<code>
 +mod.bin2<-mle2(frut.mur~dbinom(prob = exp(a + b*ur)/(1 + exp(a + b*ur)), size = fec1),
 +start = list(a= -3, b = 9))
 +</code>
 +
 +Avaliamos os perfis de verossimilhança:
 +<code>
 +par(mfrow=c(1,2))
 +plot.profmle(profile(mod.bin))
 +par(mfrow=c(1,1))
 +</code>
 +
 +
 +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.
 +
 +<code>
 +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))
 +</code>
 +
 +Compare os resultados usando a função padrão no R, ''glm'' com distribuição binomial e função de ligação logit:
 +
 +<code>
 +glm.bin <- glm(cbind(frut.mur,frut.sau)~ur, data=datos,family=binomial(link="logit"))
 +</code>
 +
 +====== Exercício =====
 +
 +Baixe o arquivo {{:biometria:verossim:besouro.csv|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 ======
 +
 +{{:biometria:verossim:rems_aula_poisson_e_binomial_negativa.pdf|slides aula}}
 +
 +{{:biometria:verossim:pinecones.r|Códigos do exemplo dos cones de pinheiro}} o conjunto de dados está na [[https://collaborate.nicholas.duke.edu/clark/pinecones.txt/view|página do autor]]
 +
 +{{:biometria:verossim:exercicio_besouros.r|Código do exercício besouro}}
 +
 +
 +
 +
 +====== Pesquisa ======
 +
 +Indique os itens do tutorial e do texto prioritários para a discussão [[http://ecologia.ib.usp.br/bie5782/doku.php?id=modelos:pesquisa1#modelos_com_par%C3%A2metros_que_s%C3%A3o_fun%C3%A7%C3%B5esmodelos_binomial_e_poisson|aqui]]
 +
 +
 +