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: | ||
+ | |||
+ | \\ | ||
+ | < | ||
+ | <font face=" | ||
+ | 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 (" | ||
+ | |||
+ | < | ||
+ | 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< | ||
+ | lambda=5 | ||
+ | frutos< | ||
+ | </ | ||
+ | |||
+ | Qual é a probabilidade teórica de encontrar uma planta sexualmente madura sem frutos? | ||
+ | < | ||
+ | dpois(0, | ||
+ | </ | ||
+ | Compare com a proporção de plantas sem frutos na amostra simulada: | ||
+ | < | ||
+ | sum(frutos==0)/ | ||
+ | </ | ||
+ | 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< | ||
+ | fa< | ||
+ | prob.obs< | ||
+ | par(las=1) | ||
+ | plot(0: | ||
+ | </ | ||
+ | Adicionamos os valores teóricos esperados: | ||
+ | < | ||
+ | prob.tr< | ||
+ | points(0: | ||
+ | </ | ||
+ | |||
+ | Definimos a função de máxima verossimilhança do parâmetro lambda da distribuição Poisson: | ||
+ | |||
+ | < | ||
+ | x< | ||
+ | poisNLL< | ||
+ | -sum(dpois(x, | ||
+ | } | ||
+ | </ | ||
+ | |||
+ | Aplicamos a função sobre uma série de valores: | ||
+ | < | ||
+ | xvec< | ||
+ | LLest< | ||
+ | </ | ||
+ | |||
+ | 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 '' | ||
+ | |||
+ | < | ||
+ | minLLest< | ||
+ | lambdaLL.fb< | ||
+ | lambdaLL.nm< | ||
+ | </ | ||
+ | |||
+ | 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)): | ||
+ | |||
+ | < | ||
+ | lambdaLL.fb | ||
+ | lambdaLL.nm | ||
+ | mean(frutos) | ||
+ | </ | ||
+ | |||
+ | Plotamos a função de verosimilhança e os pontos obtidos: | ||
+ | |||
+ | < | ||
+ | mfrutos< | ||
+ | LLest2< | ||
+ | plot(xvec, | ||
+ | abline(v=mfrutos, | ||
+ | abline(v=coef(lambdaLL.nm), | ||
+ | abline(v=lambdaLL.fb, | ||
+ | </ | ||
+ | |||
+ | |||
+ | ==== 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< | ||
+ | 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< | ||
+ | plot(x, | ||
+ | </ | ||
+ | Agora para simular um processo Poisson sorteamos amostras Poisson com parâmetro lambda igual a estes valores esperados: | ||
+ | < | ||
+ | fec< | ||
+ | </ | ||
+ | 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=" | ||
+ | </ | ||
+ | Definimos a função de verosimilhança para este modelo: | ||
+ | < | ||
+ | poisglmNLL = function(a, | ||
+ | | ||
+ | -sum(dpois(fec, | ||
+ | } | ||
+ | </ | ||
+ | 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**: | ||
+ | |||
+ | E então otimizamos, com a função mle2: | ||
+ | < | ||
+ | mod.pois< | ||
+ | </ | ||
+ | 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, | ||
+ | plot.profmle(mod.pois.prof) | ||
+ | par(mfrow=c(1, | ||
+ | </ | ||
+ | |||
+ | Numa abordagem frequentista, | ||
+ | |||
+ | Aplicando a funçã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, | ||
+ | a.est< | ||
+ | b.est< | ||
+ | curve(exp(a+b*x), | ||
+ | curve(exp(a.est +b.est*x), add=TRUE, col=" | ||
+ | legend(" | ||
+ | </ | ||
+ | |||
+ | 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< | ||
+ | ypred< | ||
+ | -sum(dnbinom(fec, | ||
+ | } | ||
+ | </ | ||
+ | Tentamos | ||
+ | |||
+ | $$σ^2= μ + μ^2/k$$ | ||
+ | |||
+ | Então podemos estimar um valor aproximado de k com: | ||
+ | < | ||
+ | med< | ||
+ | vari< | ||
+ | k.init < | ||
+ | k.init | ||
+ | |||
+ | mod.negbin< | ||
+ | </ | ||
+ | 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, | ||
+ | </ | ||
+ | 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:// | ||
+ | Compare os resultados do modelo Poisson que ajustamos com o '' | ||
+ | |||
+ | < | ||
+ | glm.pois <- glm(fec~phos, | ||
+ | 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 " | ||
+ | pois corresponde ao número de frutos que ela tenha. | ||
+ | |||
+ | Primeiro descartamos as plantas sem frutos: | ||
+ | < | ||
+ | fec1< | ||
+ | num.plant< | ||
+ | </ | ||
+ | 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< | ||
+ | </ | ||
+ | |||
+ | 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< | ||
+ | plot(ur, | ||
+ | frut.mur< | ||
+ | </ | ||
+ | O número de frutos intactos é o complemento dos murchos. | ||
+ | < | ||
+ | frut.sau< | ||
+ | </ | ||
+ | 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< | ||
+ | datos.ur< | ||
+ | conjun< | ||
+ | conjun.m< | ||
+ | |||
+ | barplot(t(conjun.m), | ||
+ | xlab= " | ||
+ | legend(" | ||
+ | arrows(15.8, | ||
+ | text(60,64, " | ||
+ | </ | ||
+ | |||
+ | 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? | ||
+ | frutos murchos/ | ||
+ | transforma as proporções (frutos murchos/ | ||
+ | |||
+ | Definimos a função de verossimilhança: | ||
+ | < | ||
+ | binomNLL< | ||
+ | prob.det=exp(a+b*ur)/ | ||
+ | -sum(dbinom(k, | ||
+ | } | ||
+ | </ | ||
+ | |||
+ | Otimizamos: | ||
+ | < | ||
+ | mod.bin< | ||
+ | </ | ||
+ | |||
+ | Alternativamente podemos usar a interface de fórmula da função '' | ||
+ | < | ||
+ | mod.bin2< | ||
+ | start = list(a= -3, b = 9)) | ||
+ | </ | ||
+ | |||
+ | Avaliamos os perfis de verossimilhança: | ||
+ | < | ||
+ | par(mfrow=c(1, | ||
+ | plot.profmle(profile(mod.bin)) | ||
+ | par(mfrow=c(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< | ||
+ | par(las=1) | ||
+ | plot(ur, prop.mur, | ||
+ | a.esti< | ||
+ | b.esti< | ||
+ | curve(exp(a+b*x)/ | ||
+ | curve(exp(a.esti+b.esti*x)/ | ||
+ | legend(" | ||
+ | </ | ||
+ | |||
+ | Compare os resultados usando a função padrão no R, '' | ||
+ | |||
+ | < | ||
+ | glm.bin <- glm(cbind(frut.mur, | ||
+ | </ | ||
+ | |||
+ | ====== Exercício ===== | ||
+ | |||
+ | Baixe o arquivo {{: | ||
+ | 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 [[http:// | ||
+ | |||
+ | |||
+ | |||