^[[start|Modelos Estatísticos: Abordagem da Verossimilhança]]^ \\ 3. Função de Verossimilhança \\ ====== Conceitos ====== * Verossimilhança * Função de Verossimilhança * Lei e Princípio de Verossimilhança * Log Verossimilhança Negativa * Método da Máxima Verossimilhança * Intervalos e Regiões de Verossimilhança * Parâmetros Inconvenientes (//Nuisance Parameters//) * Verossimilhança Estimada e Perfilhada ====== Tutorial ====== ===== Lei da Verossimilhança ===== Dada uma variável aletória $$X$$, cujo comportamento pode ser explicado por duas hipóteses: $$H_A$$ e $$H_B$$. * A hipótese $$H_A$$ afirma que a observação $$X=x$$ seria observada com probabilidade $$p_A(x)$$. * A hipótese $$H_B$$ afirma que a observação $$X=x$$ seria observada com probabilidade $$p_B(x)$$. A observação $$X=x$$ é uma evidência em favor de $$H_A$$ **vis-a-vis** (face-a-face) $$H_B$$ se, e somente se, $$p_A(x) > p_B(x)$$. A **força de evidência** em favor de $$H_A$$ vis-a-vis $$H_B$$ é dada pela **razão de verossimilhança**: $$ \frac{p_A(x)}{p_b(x)}$$. ==== Exemplo: Ninhada de Cachorro-do-Mato ==== Numa ninhada de 5 filhotes de cachorros-do-mato foi observado apenas um macho. A hipótese H1 estabelece que a probabilidade de nascimento de um filhote do sexo masculino é p=1/2, enquanto que a hipótese H2 estabelece uma probabilidade p=1/3 para filhote macho. A probabilidade de observar 1 macho em 5 é: * Sob H1: > dbinom(1, 5, 1/2) [1] 0.15625 * Sob H2: > dbinom(1, 5, 1/3) [1] 0.3292181 **Questão:** Qual das duas hipóteses é mais **//plausível//**, isto é, mais **//verossímil//**? **Questão:** Qual a **//força de evidência//** em favor da hipótese mais plausível (mais verossímil)? > dbinom(1, 5, 1/3) / dbinom(1, 5, 1/2) [1] 2.106996 **Questão:** É razoável afirmar que a hipótese H2 é **//duas vezes mais plausível//** que a hipótese H1? ===== Função de Verossimilhança ===== A variável $$X$$ é uma variável aleatória com função de densidade $$f_X(x; \theta)$$, onde: * $$x$$ é uma variável matemática de descreve os valores que $$X$$ pode assumir; * $$\theta$$ indica os parâmetros que controlam o comportamento de $$X$$. Para que a função de densidade $$f_X(\cdot)$$ seja utilizada para calcular a **//probabilidade//** da observação $$X=x$$ é necessário que o valor dos parâmetros $$\theta$$ **//sejam conhecidos//**. Na pesquisa científica o valor dos parâmetros nunca são conhecidos. O que a pesquisa gera são as observações a partir das quais se deseja //**estimar**// os parâmetros do modelo. Assim, a função de densidade pode ser utilizada de uma forma diferente, onde o valor da observação $$X=x$$ é conhecido e, portanto, não nos interessa calcular a sua probabilidade, mas os parâmetros terão que ser estimados. Essa é a **//função de verossimilhança//**: $$\mathcal{L} (\theta | X=x) = f_X(\theta | X = x)$$ Assim, na função de verossimilhança o elemento que varia são os parâmetros ($$\theta$$), enquanto que o elemento constante é a observação ($$X=x$$). ==== Exemplo: Ninhada de Cachorro-do-Mato ==== Considere que numa ninhada de cachorro-do-mato nasceu 1 macho em 5 filhotes. Assumindo que o número de machos numa ninhada segue a **//distribuição binomial//**, o que essa observação nos diz sobre a probabilidade de nascimento de filhotes machos ($$p$$)? Vejamos um gráfico da função de verossimilhança para os diferentes valores de $$p$$ que essa observação gera: > curve(dbinom(1, 5, p=x), from=0, to=1) **Questão:** É razoável afirmar que o gráfico descreve uma **//curva de plausibilidade//** para os valores de probabilidade de nascer macho, dada a observação de 1 macho em 5 filhotes? No gráfico construído, a variável foi o parâmetro $$p$$, ou seja a probabilidade de nascimento de filhote macho. Portanto, os valores obtidos pela função não podem ser interpretados como probabilidades. A função nos apresenta os valores de //**verossimilhança**//, sendo razoável assumir que eles medem a //**plausibilidade**// de cada valor do parâmetro $$p$$ numa escala relativa, isto é, em comparação com os demais valores do gráfico. \\ **PONTO de MÁXIMO da FUNÇÃO de VEROSSIMILHANÇA** O valor de **//máxima verossimilhança//** (= máxima plausibilidade) pode ser encontrado e colocado no gráfico: > theta = seq(0,1,length=100) > y = dbinom(1, 5, p=theta) > theta.mle = theta[ y == max(y) ] > abline( v = theta.mle, lty=2, col="orange" ) \\ **VEROSSIMILHANÇA RELATIVA** Podemos, então, re-escalonar o eixo-//y// (ordenadas) do gráfico em relação a esse valor máximo: > yrel = y / max(y) > plot(theta, yrel, type="l", xlab="Probabilidade p", col="blue") > abline( v = theta.mle, lty=2, col="orange" ) **Questão:** Podemos agora chamar esse gráfico de **//verossimilhança relativa//** (plausibilidade relativa) para os valores de probabilidade de nascimento de filhote macho, dada a observação de 1 macho em 5 filhotes. [Observe os valores no eixo-y] \\ **INTERVALO DE VEROSSIMILHANÇA** Imagine que queremos encontrar nesse gráfico um intervalo para todos os valores cujo valor de **//verossimilhança relativa//** é de pelo menos 1/8 em relação ao máximo. Isto é, o valor máximo tem **//verossimilhança relativa//** de até 8 vezes os demais valores do intervalo. > theta.interv = theta[ max(y)/y <= 8 ] > lines(theta.interv, rep(1/8, length(theta.interv)), col="red", lwd=3) ===== Múltiplas Observações ===== Geralmente, os dados científicos são obtidos na forma de **//uma amostra//** da população, sendo que essa amostra é composta de **//várias observações independentes//**. Se foram tomadas 4 observações independentes de uma variável aleatória: $$X$$ ($$\mathbf{X}_4 = \{x_1, x_2, x_3, x_4\}$$) e cada uma delas possui uma função de verossimilhança para o parâmetro ($$\theta$$) do modelo, então a função de verossimilhança da amostra é o **//produto//** das funções de verossimilhança das observações independentes: $$ \mathcal{L} ( \theta | X = \mathbf{X}_4 ) = \mathcal{L} ( \theta | X = x_1 ) \cdot \mathcal{L} ( \theta | X = x_2 ) \cdot \mathcal{L} ( \theta | X = x_3 ) \cdot \mathcal{L} ( \theta | X = x_4 )$$ Generalizando: se uma amostra é composta de $$n$$ observaçoes independentes, então a verossimilhança da amostra é o produto das verossimilhanças das observações individuais: $$ \mathcal{L} ( \theta | X = \mathbf{X}_n ) = \mathcal{L} ( \theta | X = x_1 ) \cdot \mathcal{L} ( \theta | X = x_2 ) \cdot \mathcal{L} ( \theta | X = x_3 ) \cdot \ldots \cdot \mathcal{L} ( \theta | X = x_n )$$ $$ \mathcal{L} ( \theta | X = \mathbf{X}_n ) = \prod_{i=1}^n \mathcal{L} ( \theta | X = x_i )$$ ==== Exemplo: Número de Espécies Arbóreas ==== Num levantamento em floresta ombrófila densa foram contadas o número de espécies por parcelas de 900m2, obtendo-se os valores: 5, 8, 11 e 9. Como se trata dados de contagem, a distribuição de Poisson é candidata à modelagem desses dados. Se construírmos uma **//curva de verossimilhança//** para cada observação tempos: > curve(dpois( 5, lambda=x), 0, 30 , xlab=expression(lambda), ylab="") > curve(dpois( 8, lambda=x), 0, 30 , add=TRUE) > curve(dpois( 9, lambda=x), 0, 30 , add=TRUE) > curve(dpois( 11, lambda=x), 0, 30 , add=TRUE) Entretanto, nossa **amostra** é formada pelas quatro observações **conjuntamente** e precisamos representar uma única curva: > x = seq(0,30, by=0.1) > y1 = dpois(5, lambda=x) > y2 = dpois(8, lambda=x) > y3 = dpois(9, lambda=x) > y4 = dpois(11, lambda=x) > y = y1 * y2 * y3 * y4 > plot(x, y, type="l", col="red") Em um único passo, teríamos: > curve(dpois(5,x)*dpois(8,x)*dpois(9,x)*dpois(11,x), 0, 30 , xlab=expression(lambda), ylab="Verossimilhança", col="red") \\ ===== Função de Log-Verossimilhança Negativa ===== A função de verossimilhança quando aplicada a uma única observação resulta frequentemente (mas não necessariamente) em valores menores do que um. Quando ela é aplicada a uma amostra com muitas observações, o produto das verossimilhanças individuais vai rapidamente em direção ao zero, se tornando um número muito pequeno. Embora se possa trabalhar com a verossimilhança relativa, a função de verossimilhança é geralmente uma função de difícil tratamento matemático. A se aplicacar uma tranformação na função de verossimilhança, que consiste na aplicação do logaritmo e na mudança de sinal, se obtem a **//função de log-verossimilhança negativa//** que é uma função de tratamento matemático mais simple e que resulta em valores numéricos de manipulação mais fácil. ==== Exemplo: Número de Espécies Arbóreas ==== Se compararmos os valores nos vetores ''c(y1,y2,y3,y4)'' e ''y'', veremos que ao fazermos o produto, obtemos números muito pequenos. > range(c(y1,y2,y3,y4)) [1] 0.0000000 0.1754674 > range(y) [1] 0.0000000000 0.0001162465 Em grandes amostras, com 100 ou mais observações, esse produto gera números muito pequenos, sendo problemático o seu processamento no computador. Por isso é conveniente utilizarmos a transformação para a função de log-verossimilhança negativa: > log.y = - log(y) > log.y2= - (log(y1) + log(y2) + log(y3) + log(y4)) > plot(x, log.y, type="l", col="red", xlab=expression(lambda)) > lines(x, log.y2, col="blue") Note que ''log.y'' e ''log.y2'' são identicos. === Como Calcular a Log-Verossimilhança Negativa para uma Amostra === No R, a log-verossimilhança negativa pode ser gerada diretamente das funções de densidade das distribuições utilizando o argumento "log = TRUE" que existe em todas elas. Vejamos uma forma mais prática (e elegante) de se gerar a função de log-verossimilhança negativa da distribuição Poisson para esse conjunto de dados: > library(MASS) # carrega o pacote "MASS" > args(dpois) # verifica os argumentos da função "dpois" > lpois = function(x,lambda) -dpois(x,lambda,log=TRUE) # cria uma função da log-veros. neg. > llikpois = Vectorize( lpois, c("x","lambda") ) # cria uma função "vetorizada" da log-veros. neg. > x = c(5, 8, 9, 11) # "x" agora são as observações > lambda = seq(0, 30, by=0.1) # "lambda" o vetor com os valores do parâmetro > lver.mat = outer( x, lambda, "llikpois" ) # uma matrix a log-veros. neg. de cada observação > lver = apply( lver.mat, 2, sum) # um vetor com a log-veros. neg da amostra > > plot( lambda, lver, type="l", col="red", xlab=expression(lambda) ) \\ Da mesma forma que re-escalonamos a curva de verossimilhança em relação ao seu máximo, podemos re-escalonar essa curva de log-verossimilhança negativa em relação ao seu mínimo: > plot(lambda, lver - min(lver), type="l", col="red", xlab=expression(lambda) ) \\ O ponto que maximiza a verossimilhança é o mesmo que **minimiza a log-verossimilhança negativa**: > lambda.min = lambda[ min(lver) == lver ] > abline(v = lambda.min, lty=2, col="purple") Note que o valor mais verossímil //graficamente// para o parâmetro $$\lambda$$ da Poisson é: > lambda.min [1] 8.3 que é muito próximo a média amostral: > mean(c(5,8,9,11)) [1] 8.25 \\ Também podemos estabelecer graficamente um intervalo de verossimilhança para o parâmetro $$\lambda$$ na curva de log-verossimilhança negativa: > lambda.int = lambda[ lver - min(lver) <= log(8) ] > lines(lambda.int, rep(log(8), length(lambda.int)), col="purple", lwd=3 ) ===== Método da Máxima Verossimilhança ===== O método da Máxima Verossimilhança consiste em realizar o que foi feito nos exemplos acima, isto é, encontrar a estimativa dos parâmetros de um modelo que //**maximiza**// a função de verossimilhança ou, o que é equivalente, **//minimizar//** a função de log-verossimilhança negativa. As estimativas assim geradas são chamadas de //MLE = Maximum Likelihood Estimates//. ==== Minimizando a Função de Log-Verossimilhança no R ==== O R possui duas funções básicas para minimizar expressões matemáticas, ou seja, para buscar os valores dos parâmetros de uma função que resultam no menor valor da função: * [[http://finzi.psych.upenn.edu/R/library/stats/html/optimize.html|optimize]]: otimização unidimensional, isto é, para funções com apenas um parâmetro livre. * [[http://finzi.psych.upenn.edu/R/library/stats/html/optim.html|optim]]: para otimização de funções com mais de um parâmetro livre. A partir destas funções há várias definidas especificamente para minimizar funções de verossimilhança. As duas mais usadas são a [[http://finzi.psych.upenn.edu/R/library/stats4/html/mle.html|mle]] do pacote ''stats4'', e a [[http://finzi.psych.upenn.edu/R/library/bbmle/html/mle2.html|mle2]], do pacote ''bbmle''((veja também o ótimo tutorial deste pacote em [[http://finzi.psych.upenn.edu/R/library/bbmle/doc/mle2.pdf]] )). Como o pacote ''bbmle'' foi feito por um ecólogo, para acompanhar o ótimo livro sobre modelos que usamos como referência básica, usaremos a ''mle2'' em nossos tutoriais. A função que deve ser fornecida para a "mle2" é semelhante às funções apresentadas acima, mas nesse caso ela deve ser uma função que **soma** a log-verossimilhança negativa para cada valor dos parâmetros para **todas observações** na amostra. ==== Exemplo: Número de Plântulas do Palmiteiro Juçara ==== Nesse exemplo utilizaremos os dados de número de plântulas do palmiteiro juçara (//Euterpe edulis//) encontradas em parcelas de regeneração natural. > euplant = c(5, 23, 23, 17, 0, 4, 1, 0, 0, 1, 0, 2, 26, 65, 34, 14, 18, 13, 19, 7) > length(euplant) > mean(euplant) É necessário criar a função de log-verossimilhança negativa da amostra para distribuição Poisson: > nvl = function(lambda) -sum(stats::dpois( euplant, lambda, log=TRUE ) ) Vamos encontrar a estimativa de máxima verossimilhança (MLE) do parâmetro $$\lambda$$: > library(bbmle) # Primeiro carrega o pacote "bbmle" > euplant.mle = mle2(nvl, start=list(lambda=10) ) # gera um objeto "mle" > summary(euplant.mle) # resumo do objeto "mle" > logLik(euplant.mle) # valor da log-verossimilhança negativa do objeto Verifique se a MLE gerada pela função "mle" de fato corresponde à estimativa esperada teoricamente: > coef(euplant.mle) > mean(euplant) Para visualizar a //**curva da log-verossimilhança negativa**// é interessante usar uma função "feita sob encomenda". Grave o seguinte código {{:biometria:verossim:plot-profmle.txt|plot-profmle.txt}} num arquivo texto e utilize o seguinte comando para gerar a função no R: > source("plot-profmle.txt") Agora é possível visualizar facilmente a curva: > plot.profmle( profile(euplant.mle) ) ===== Propriedades dos MLE ===== Para demonstrar as propriedades dos MLE é necessário lançar mão de um expediente de simulação. A simulação consiste nos seguintes passos: * (1) gerar-se um grande número de amostras de uma distribuição conhecida com valores conhecidos dos parâmetros; * (2) obtem-se para cada uma das amostras simuladas as MLE; * (3) constroem-se gráficos da distribuição das MLE para analisar o seu comportamento. ==== Consistência ==== Para verificar a consistência, utilizaremos o MLE do desvio padrão que é sabidamente viciado: $$ \hat{s} = [ \frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2 ]^(1/2) $$ Primeiramente, criaremos uma função para obter esta MLE: > sd.mle = function(x) sqrt( (1/length(x)) * sum( (x - mean(x))^2 ) ) Esse estimador será aplicado a amostras da distribuição Gaussiana com média 100 e desvio padrão 5, de tamanhos crescentes. Gerando a distribuição das MLE para amostras de tamanho 10: > samp10 = sort(rep(1:1000, 10)) # identificador das 1000 amostras de tamanho 10 > vals10 = rnorm(1000*10, mean=10, sd=5) # valores das observaçoes > dens10 = density( tapply(vals10, samp10, sd.mle) ) # obtendo a MLE para cada amostra e a densidade empírica > dens10$y = dens10$y/ max(dens10$y) # convertendo o valor de densidade para valor relativo ao max. > plot( dens10, type="l", xlab="Valor Estimado", ylab="Densidade" ,main="") # gráfico da distribuição das MLE > abline(v = 5, col="red", lwd=2 ) # posição do valor do parâmetro no gráfico Note que a distribuição das MLE das amostras de tamanho 10 não é simétrica e não parece centrada no valor do parâmetro. Gerando a distribuição das MLE para amostras de tamanho 100: > samp100 = sort(rep(1:1000, 100)) > vals100 = rnorm(1000*100, mean=100, sd=5) > dens100 = density( tapply(vals100, samp100, sd.mle) ) > dens100$y = dens100$y/ max(dens100$y) > lines( dens100, col="blue" ) Note que as MLE das amostras de tamanho 100 já tem distribuição simétrica e o ponto central se aproxima bastante do valor do parâmetro. ==== Eficiência Assimptótica ==== A eficiência assimptótica afirma que os MLE terão na menor variância possível à medida que o tamanho de amostra aumenta. Vamos comparar a variâncias dos MLEs do desvio padrão obtidos acima (nossa função "sd.mle") contra os valores do estimador tradicional, calculados para as mesmas amostras (a função ''sd'' do R): > ## Amostras de tamanho 10 > var( tapply(vals10, samp10, sd.mle) ) > var( tapply(vals10, samp10, sd) ) >## Amostras de tamanho 100 > var( tapply(vals100, samp100, sd.mle) ) > var( tapply(vals100, samp100, sd) ) ==== Normalidade Assimptótica ==== Esta propriedade garante que os MLE aproximam-se de uma distribuição normal com o aumento do tamanho da amostra. A melhor forma de investigar a aproximação assimptótica da distribuição Gaussiana é através do gráfico Quantil-Quantil: > par(mfrow=c(2,1)) > qqnorm( tapply(vals10, samp10, sd.mle), main="MLE do Desvio Padrão, N=10" ) > qqline( tapply(vals10, samp10, sd.mle) ) > qqnorm( tapply(vals100, samp100, sd.mle), main="MLE do Desvio Padrão, N=100") > qqline( tapply(vals100, samp100, sd.mle) ) > par(mfrow=c(1,1)) ==== Invariância ==== A propriedade de invariância dos MLE é de difícil de comprovação via simulação, sendo mais fácil exemplificar com estimadores que não tem essa propriedade. O estimador tradicional da variância tem a propiedade de não ser viciado, independentemente do tamanho da amostra. $$ \hat{s^2} = \frac{ \sum_{i=1}^n (x_i - \bar{x})^2 }{ n - 1 }$$ Isso pode ser comprovado via simulação: > var.trad10 = tapply( vals10, samp10, var) # estimativas tradicionais da variância de amostras de tamanho 10 > plot( density(var.trad10), type="l" ) # distribuição das estimativas > abline( v = mean(var.trad10) ) # média das estimativas > abline( v = 5^2, col="red" ) # valor do parâmetro \\ Será que a propriedade de não ser viciado se mantem quando a variância é transformada em desvio padrão ? $$ \hat{s} = [ \hat{s^2} ]^(1/2) = [ \frac{ \sum_{i=1}^n (x_i - \bar{x})^2 }{ n - 1 } ]^(1/2)$$ > sd.trad10 = tapply(vals10, samp10, sd ) > plot(density(sd.trad10), type="l") > abline( v = mean(sd.trad10) ) > abline( v = 5, col="red" ) ===== Superfície e Curvas de Verossimilhança ===== Quando temos mais de parâmetro no modelo estatístico, a função de verossimilhança se torna uma função **multivariada** e, consequentemente, ela só pode ser representada graficamente por **superfícies**. No caso de dois parâmetros, como a distribuição Gaussiana, a função de verossimilhança é **bivariada**, sendo possível construir gráficos de superfície, mas com três ou mais parâmetros a construção de gráficos e sua interpretação se torna complexa. ==== Exemplo: Superfície para uma Variável Gaussiana ==== > dados = rnorm(1000, mean=150, sd=20 ) # gerando os dados > m = seq(125, 175, length=50) # vetor de médias > s = seq(10, 60, length=50) # vetor de desvio padrão > lgauss = function(m, s, dados) -sum(stats::dnorm(dados, mean=m, sd=s, log=T)) # função de log-ver. neg. > llikgauss = Vectorize(lgauss, c("m","s")) # vetorização da função > sup.mat = outer(m, s, llikgauss, dados) # cálculo da superfície > contour(m, s, sup.mat, xlab="Média", ylab="Desvio Padrão", nlevels=20) # gráfico de contorno > abline(v=150, col="red", lty=2) # linha da média > abline(h=20, col="red", lty=2) # linha do desvio padrão > persp(m, s, sup.mat, theta=135, phi=10, col="purple") # gráfico de perspectiva ==== Exemplo: Árvores de Copaifera langsdorffii no Cerradão ==== //Copaifera langsdorffii// é uma espécie árborea que ocorre em várias formações florestais no Estado de São Paulo, sendo muito ambundante no cerradão. O arquivo {{:biometria:verossim:copaifera.csv|copaifera.csv}} apresenta o número de árvores em 64 parcelas na Estação Ecológica de Assis. > copa = read.csv("copaifera.csv",header=T) > mean(copa$narv) > var(copa$narv) A distribuição espacial da //Copaifera// é em geral agregada, de forma que o número de árvores por parcela pode ser descrita com uma distribuição Binomial Negativa com parâmetros: * $$k$$ - parâmetro de agregação e * $$\mu$$ - número médio de árvores por parcela. > > lnbinom.x = function(k, mu, x) -sum(dnbinom(x, size=k, mu=mu, log=T)) # função de log-veros. neg. > lnbinom.vec = Vectorize( lnbinom.x, c("k","mu") ) # função vectorizada > > k.var = seq(0.1, 10, length=50) # vetor do parâmetro k > mu.var = seq(10, 80, length=50) # vetor do parâmetro mu > sup.mat = outer( k.var, mu.var, lnbinom.vec, x=copa$narv ) # matriz com a log-veros. neg. > contour(k.var, mu.var, sup.mat, nlevels=20, xlab="k", ylab=expression(mu)) # gráfico de contorno > persp(k.var, mu.var, sup.mat, theta=-35, phi=20, col="lightblue") # gráfico de perspectiva ===== Verossimilhança Estimada e Verossimilhança Perfilhada ===== Quando o modelo estatístico tem muitos parâmetros, soma-se ao problema de interpretar uma superfície multivariada a questão de alguns dos parâmetros do modelo serem //incovenientes//, isto é, não possuem uma interpretação biológica clara, mas são necessários ao modelo. Há duas formas de contornar esse problema que permitem investigar somente o comportamento dos parâmetros de interesse um-a-um, sem a necessidade de avaliar uma superfície: * a verossimilhança estimada e * a verossimilhança perfilhada. A verossimilhança perfilhada é uma forma mais apropriada de estudar a superfície de verossimilhança, pois representa de modo mais fidedigno a superfície de verossimilhança na vizinhança das MLE. Mas em algumas situações, quando não se tem à disposição os dados primários, a verossimilhança estimada pode ser a única alternativa para se analisar os dados. Vejamos alguns exemplos: ==== Verossimilhança Estimada ==== ==== Exemplo: Produção Média de Florestas ==== Num inventário em floresta plantada de //Pinus// foram medidas 52 parcelas, tendo sido obtida a produção média de 243 $$m^3\, ha^{-1}$$, com desvio padrão de 28 $$m^3\, ha^{-1}$$. Qual o intervalo de verossimilhança (razão de 8) dessas estimativa da produção ((note que a pergunta é sobre o comportamento da estimativa, e não da população de medidas de produtividade de onde ela veio)) ? Pelo Teorema Central do Limite, sabemos que a média amostral segue distribuição Gaussiana, com média igual à média populacional e desvio padrão igual ao desvio padrão populacional dividido pela raiz quadrada do tamanho da amostra. Portanto, podemos assumir que a produção média relatada é uma observação de uma distribuição Gaussiana. Como só dispomos dessa informação poderemos avaliar apenas a variação da média, deixando o desvio padrão constante (igual ao relatado): > prod.var = seq(230, 260, length=100) # intervalo de variação da produção média > lvn.prod = -dnorm(243, mean=prod.var, sd=28/sqrt(52), log=T) # valores de log-veros. negativa > lvn.prod = lvn.prod - min(lvn.prod) # valores de log-veros. negativa RELATIVA > plot(prod.var, lvn.prod, type="l") # gráfico log-veros. neg. relativa > abline(h=log(8), lty=2, col="red") # intervalo de verossimilhança ==== Exemplo: Árvores Bifurcadas ==== Numa floresta planta de eucalipto foram encontradas 7 árvores bifurcadas numa amostra aleatória de 150. Deseja-se a taxa de bifurcação e um intervalo de //verossimilhança//. Neste caso a estimativa é o parâmetro $$p$$ de uma binomial, sendo o outro parâmetro (número de tentativas, $$N$$) conhecido. \\ > p.var = seq(0.01, 0.10, length=100) > lvn.binom = -dbinom(7, 150, prob=p.var, log=T) > lvn.binom = lvn.binom - min(lvn.binom) > plot( p.var, lvn.binom, type="l", col="red", xlab="Prob. Bifurcação") > p.mle = p.var[lvn.binom == min (lvn.binom)] > abline(v=p.mle, lty=2, col="blue") > abline(h=log(8), lty=2, col="blue") ==== Verossimilhança Perfilhada ==== ==== Exemplo: Estrutura de Tamanho em Florestas ==== O arquivo {{:biometria:verossim:diam.csv|diam.csv}} apresenta dados do diâmetro à altura do peito (DAP) de árvores de uma parcela em floresta nativa no município de Bom Jardim, MA. Frequentemente se assume que a distribuição dos DAP das árvores em floresta tropical segue uma **distribuição exponencial**. Existe, entretanto, a **distribuição Weibull**, bem mais geral, da qual a distribuição exponencial é um caso particular. Quando o parâmetro da forma (//shape//) da Weibull é 1, ela se reduz a uma exponencial. Assim, podemos verificar se a distribuição exponencial é de fato um bom modelo para os dados acima, comparando-a com a Weibull. \\ Inicialmente é necessário ler os dados do arquivo > dap = read.csv("diam.csv",header=T)$dap Para ajustar a dist. Weibull aos dados é necessário considerar que só foram medidas as árvores com DAP a partir de 14.3 cm (45cm de circunferência) e, portanto, o valor 14.3 representa um valor mínimo conhecido. > dap.a = dap - 14.3 # 14.3 é o DAP mínimo de medição > lweibull = function(esc, forma){ -sum( dweibull(dap.a, shape=forma, scale=esc, log=T) ) } # função log-ver.neg. da Weibull > dap.mle = mle2( lweibull, start=list(esc=10, forma=1) ) # obtendo as MLE > summary(dap.mle) # resumo das MLE > > dap.mle.p = profile( dap.mle ) # perfilhando as MLE > par(mfrow=c(1,2)) # preparando para 2 gráficos > plot.profmle( dap.mle.p ) # gráficos da veros. perfilhada **Questão:** é razoável assumir que os dados podem ser modelados por uma distribuição exponencial? ==== Exemplo: Árvores de Copaifera langsdorffii ==== Para obtermos a curva de verossimilhança perfilhada para o modelo de Binomial Negativa nos dados de número de árvores de //Copaifera langsdorffii// devemos re-escrever a função de log-verossimilhança negativa e ajustar o modelo com a função ''mle2''. > lnbinom = function(k, mu) -sum(dnbinom(narv, size=k, mu=mu, log=T)) > narv = copa$narv > copa.mle = mle2( lnbinom, start=list(k=2, mu=32) ) > summary(copa.mle) > par(mfrow=c(1,2)) > plot.profmle( profile(copa.mle) ) > > par(mfrow=c(1,1)) > contour(k.var, mu.var, sup.mat, nlevels=20, xlab="k", ylab=expression(mu)) > abline(v=coef(copa.mle)["k"], col="red") > abline(h=coef(copa.mle)["mu"], col="red") > \\ \\ ====== Leituras ====== ==== PRINCIPAL ==== * **Verossimilhança e Máxima Verossimilhança** 2009, Batista, J.L.F. {{:biometria:verossim:verossim.pdf|pdf}} ==== COMPLEMENTARES ==== * **The Concept of Likelihood** 1992. In: Edwards, A.W.F., 1992 Likelihood, cap.2, p.8-24. Baltimore: John Hopkins University Press. ====== Exercícios ====== ==== 1) Ocorrência de Eventos Raros ==== A tabela abaixo apresenta o levantamento de fenômenos relativos a árvores individuais em floresta plantada de //Eucalyptus grandis//: ^ Fenômeno ^ Tamanho da Amostra ^ Árvores Afetadas ^ | Cancro | 250 | 0 | | Árvore Seca | 300 | 10 | Utilizando o modelo binomial, encontre a MLE e os intervalos de verossimilhança (razão de 8) para probabilidade de ocorrência dos fenômenos (parâmetro $$p$$). ==== 2) Cipós e Sobrevivência de Árvores ==== A tabela abaixo apresenta o número de árvores em diferentes classes de incidência de cipó (incidência crescente). Os dados no ano de 1997 se referem a árvores com as quais se iniciou o estudo e os dados no ano de 2000 se referem às árvores que sobreviveram. | ^ Incidência de Cipó ^^^ ^ Ano ^ Baixa ^ Média ^ Alta ^ | 1997 | 1146 | 797 | 1081 | | 2000 | 1044 | 711 | 934 | Use perfis de verossimilhança para avaliar se as taxas de sobrevivência são diferentes entre as classes de incidência de cipó. ==== 3) Cipós e Crescimento de Árvores ==== A tabela abaixo traz a média e desvio padrão do crescimento relativo em diâmetro de árvores sob diferentes incidências de cipó: | | ^ Crescimento em Diâmetro (cm)^^ ^ Incidência de Cipó ^ Tamanho da Amostral ^ Média ^ Desvio Padrão ^ | Baixa | 1044 | 1,717 | 1,965 | | Méda | 711 | 1,175 | 1,319 | | Alta | 934 | 0,696 | 0,963 | Utilizando a função de verossimilhança **estimada** avalie se é realista afirmar que o crescimento médio das árvores diferem entre as classes de incidência de cipó. == DICA == Lembre-se que a média amostral $$\bar x$$ quando obtida de grandes amostras pode ser aproximada por uma variável normal. O desvio-padrão desta variável aleatória normal, também chamado de erro-padrão da média amostral, é $$\sigma_{\bar x} \ = \ frac{\sigma}{\sqrt{n}$$ Onde $$\sigma$$ é o desvio-padrão da população de onde foi tomada a amostra (estimado pelo desvio-padrão amostral), e $$n$$ é o tamanho da amostra. \\ \\ ^[[start|BIE 5187 - Modelos Estatísticos na Pesquisa em Ecologia e Recursos Naturais]]^ ===== Recursos para Estudo ===== ==== Resumo da Aula ==== * {{:biometria:verossim:aula3funcao-verossim.pdf|Verossimilhança e Função de Verossimilhança}} ==== Soluções dos Exercícios ==== * {{:biometria:verossim:exerc_verossimilh.txt|Códigos em R das soluções dos exercícios}} ==== Na Internet ==== * Índice do pacote [[http://finzi.psych.upenn.edu/R/library/bbmle/doc/mle2.pdf|bbmle]] na página do R. * Índice do [[http://finzi.psych.upenn.edu/R/library/emdbook/html/00Index.html|emdbook]] na página do R. Este é um pacote de apoio ao livro [[http://www.zoology.ufl.edu/bolker/emdbook/index.html|Ecological Models and Data in R]]. * [[http://finzi.psych.upenn.edu/R/library/bbmle/doc/mle2.pdf| tutorial]] do pacote ''bbmle''. * [[http://www.cimat.mx/reportes/enlinea/D-99-10.html|"Likelihood"]]: Uma palestra de A. W. F. Edwards, um dos pioneiros na formalização e defesa do princípio da verossimilhança. * [[http://statgen.iop.kcl.ac.uk/bgim/mle/sslike_1.html|Maximum Likelihood Primer]] por S. Purcell. Feito para aplicações em genética, mas uma ótima introdução ao tema. * Scott, C., & Nowak, R. 2004. [[http://cnx.org/content/m11446/1.5/|Maximum Likelihood Estimation.]] Connexions, May 12, 2004. ===== Pesquisa ===== Indique os tutoriais e partes da leitura básica que merecem mais atenção na discussão [[http://ecologia.ib.usp.br/bie5782/doku.php?id=bie5782:verossim:pesquisa1#fun%C3%A7%C3%A3o_de_verossimilhan%C3%A7a|aqui]].