====== 3. Função de Verossimilhança ====== A função de verossimilhança liga nossos dados a modelos probabilísticos. Enquanto distribuições de probabilidades atribuem probabilidades a cada observação possível, a função de verossimilhança expressa a plausibilidade de diferentes modelos, dada uma mesma observação ou conjunto de dados. Os tutoriais a seguir mostram como construir funções de verossimilhança em linguagem computacional, e também como encontrar o valor máximo destas funções com otimização computacional. Vamos fazer também algumas simulações para demonstrar propriedades das funções e dos estimadores de máxima verossimilhança, que são o resultado dessas otimizações. Como sempre, o objetivo é usar recursos computacionais para tornar mais concretos conceitos matemáticos e estatísticos. ====== Conceitos ====== * Verossimilhança * Função de Verossimilhança * Log Verossimilhança Negativa * Método da Máxima Verossimilhança * Propriedades dos MLE * Curva e Superfície de Verossimilhança * Parâmetros Inconvenientes (//Nuisance Parameters//) * Verossimilhança Estimada e Perfilhada * Intervalos e Regiões de Verossimilhança ====== Tutoriais ====== ===== 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**: | ${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/3, enquanto que a hipótese H2 estabelece uma probabilidade p=1/2 para filhote macho. A probabilidade de observar 1 macho em 5 é: * Sob H1: > dbinom(1, 5, 1/3) [1] 0.3292181 * Sob H2: > dbinom(1, 5, 1/2) [1] 0.15625 **Questão:** Qual das duas hipóteses é mais **//plausível//**, isto é, mais **//verossímil//**? **R:** A hipótese H1, pois atribui mais probabilidade à observação de um macho em uma ninhada de 5 do que a hipótese H2 ($ p_H1(x) > p_H2(x)$ **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 **R:** Usando a razão de verossimilhanças como medida de força de evidência dos dados, a hipótese H1 é **//duas vezes mais plausível//** que a hipótese H2. Ou seja, a observação é uma evidência mais forte a favor da hipótese H1 do que da hipótese H2. ===== 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: cria-se com ela uma função dos parâmetros, mantidos os dados constantes. Essa é a **//função de verossimilhança//**, que é a função de densidade com parâmetros desconhecidos e condicionada aos dados, que são conhecidos: $$\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) **Portanto:** 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$ da distribuição binomial. Apesar, deste caso específico o parâmetro ser interpretável como a probabilidade de um sucesso, isso não é a probabilidade atribuída à observação, é apenas um parâmetro da distribuição. 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" ) **Portanto:** 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 temos: 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. Se aplicamos o logaritmo e mudança de sinal à função de verossimilhança relativa, obtemos a **//função de log-verossimilhança negativa//** que é uma função de tratamento matemático mais simples 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, a solução analítica para o estimador de máxima verossimilhança $\lambda$ da Poisson: mean(c(5,8,9,11)) #[1] 8.25 A diferença que encontramos é apenas por causa da precisão desta solução gráfica. Que é didática mas pouco precisa. \\ 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 de uma maneira mais precisa 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. Vamos ajustar a distribuição Poisson ao dados, o que siginifica buscar a estimativa de máxima verossimilhança do único parâmetro da distribuição Poisson ((taxa média de ocorrência de plantas por parcelas, no caso.)). 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( 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**// use a função [[http://finzi.psych.upenn.edu/R/library/sads/html/plotprofmle.html|plotprofmle]], do pacote sads. Instale e carregue o pacote e então será possível visualizar facilmente a curva: library(sads) plotprofmle( profile(euplant.mle) ) ===== Propriedades dos MLE ===== Os estimadores de máxima verossimilhança têm várias propriedades de um [[http://en.wikipedia.org/wiki/Estimator|estimador]] ideal. Para demonstrar estas propriedades dos MLEs é necessário lançar mão de um expediente de simulação. A simulação consiste nos seguintes passos: - geramos um grande número de amostras de uma distribuição conhecida com valores conhecidos dos parâmetros; - obtemos para cada uma das amostras simuladas as MLE; - construímos gráficos da distribuição das MLE para analisar o seu comportamento. ==== Consistência ==== A [[http://en.wikipedia.org/wiki/Consistent_estimator|consistência]] é a propriedade de convergir em probabilidade para o valor do parâmetro com o aumento da amostra. Para verificar a consistência, utilizaremos o MLE do desvio padrão da Gaussiana, que é sabidamente enviesado: $$ \widehat{s} = \left[ \frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2 \right] ^{(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. Mais precisamente, a probabilidade de que probabilidade das estimativas estarem a uma certa distância do valor real fica cada vez menor com o aumento da amostra. ==== Eficiência Assintótica ==== A [[http://en.wikipedia.org/wiki/Efficient_estimator|eficiência assintótica]] afirma que os MLE terão a 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 Assintótica ==== Esta é a propriedade de [[http://en.wikipedia.org/wiki/Asymptotic_normality|distribuição assintótica]] que garante que a distribuição dos MLE ((Esta é a distribuição teórica dos MLEs que seriam obtidos se os dados fossem tomados e o modelo ajustado a eles várias vezes.)) aproximam-se de uma distribuição normal com o aumento do tamanho da amostra. Pode ser entendida como uma extensão de Teorema Central do Limite para todos os MLEs. A melhor forma de investigar a aproximação assintó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 por simulação, sendo mais fácil exemplificar com estimadores que não têm essa propriedade. O estimador tradicional da variância tem a propriedade de não ser viciado, independentemente do tamanho da amostra. $$ \widehat{s^2} = \frac{ \sum_{i=1}^n (x_i - \bar{x})^2 }{ n - 1 }$$ Isso pode ser comprovado com uma 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 \\ No entanto, este não é um estimador de máxima verossimilhança. Por isso a propriedade de não ser viciado não se mantém quando a variância é transformada em desvio padrão , que é uma transformação monotônica: $$ \widehat{s} = \sqrt{ \widehat{s^2}} \ = \ \sqrt{ \frac{ \sum_{i=1}^n (x_i - \bar{x})^2 }{ n - 1 } }$$ Verifique isto repetindo a avaliação de viés feita antes para a variância: 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 abundante no cerradão. O arquivo {{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 {{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 distribuição 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 plotprofmle( 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)) plotprofmle( 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") \\ \\ ====== Exercícios ====== Faça os exercícios 203.1 a 203.3 no sistema [[http://notar.ib.usp.br|notaR]]. ====== Questões motivadoras para a dicussão ====== Proponha questões para discutirmos em nosso subfórum. Criamos um subfórum para isso [[http://bie5781.64486.x6.nabble.com/Questoes-para-discussao-de-verossimilhanca-f146.html|aqui]]. ===Questões escolhidas=== - Na inferência por verossimilhança o espaço amostral é irrelevante, uma vez que não se considera o que poderia se amostrar, mas sim o que foi amostrado. Quais as implicações dessa característica na definição do delineamento amostral? Seria diferente do delineamento amostral quando se usa Inferência Clássica? - Na análise de dados usando a Lei da Verossimilhança não se compara uma determinada hipótese a uma hipótese nula, e sim comparam-se hipóteses entre si (podendo haver mais de duas). Podemos dizer que esse tipo de análise depende mais fortemente da habilidade do pesquisador para formular hipóteses e criar modelos que explicarão melhor seus dados (já que as hipóteses não estão prontas, elas devem ser formuladas)? Seria a estatística frequentista mais adequada a pesquisadores menos experientes, menos hábeis a formularem modelos mais coerentes com a realidade? Seria possível uma transposição da ideia de hipótese nula para a abordagem de verossimilhança? - O Princípio da Verossimilhança afirma que a função de verossimilhança contém toda informação que um conjunto de dados tem sobre um dado modelo. Quais as vantagens de se realizar a inferência com base apenas na função de verossimilhança? Quais as desvantagens? ====== Recursos para Estudo ====== ===== Leituras ===== ==== Principal ==== * {{:03-funcao-veros:verossim.pdf|Batista, J.L.F. 2009 Verossimilhança e Máxima Verossimilhança.}} ==== Complementares ==== * {{:leituras:jlfb-tese-livre-docencia.pdf|Batista, J.L.F. 2014 Biometria Florestal segundo o Axioma da Verossimilhança.}} ESALQ-USP, Tese de Livre Docência, 391p. * Likelihood and all that. capítulo 6 de: Bolker, B.M. 2008 Ecological Models and Data in R Princeton: Princeton University Press. * The Concept of Likelihood , cap.2 de: Edwards, A.W.F., 1992 Likelihood. Baltimore: John Hopkins University Press. * Hobbs, N.T. & Hilborn, R. 2006. Alternatives to statistical hypothesis testing in ecology: A guide to self-teaching. Ecological Applications: 16(1): 5-19. (o [[http://www.esapubs.org/archive/appl/A016/001/appendix-A.htm|Apêndice A]] é uma excelente introdução). ===== 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''. Ótimo roteiro para aplicação do método de máxima verossimilhança em R. * [[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. * [[http://elaestaemtudo.ime.usp.br/?portfolio=new-mobile-app-3-2-2-5-2-2-2|Renata e a máxima verossimilhança]], explicação super didática, do projeto [[http://elaestaemtudo.ime.usp.br/|Ela está em tudo (Sim, Matemática é coisa de menina)]] , do IME -USP e UFABC. * Scott, C., & Nowak, R. 2004. [[http://cnx.org/content/m11446/1.5/|Maximum Likelihood Estimation.]] Connexions, May 12, 2004. ====Notebooks de deduções analíticas de MLEs==== Nossos notebooks interativos em Python demonstrando os passos para a dedução de alguns estimadores de máxima verossimilhança (MLEs). * Você pode executar os notebooks online no servidor [[https://mybinder.org/v2/gh/piLaboratory/bie5782/master?filepath=jupyter%2F|Binder]]. Demora um pouco para carregar. As instruções passo a passo estão na seção [[04-parametros-constantes:04-parametros-constantes#extras|"Extras"]] do tutorial sobre dedução analíticas da MLE da exponencial. * Ou você pode ver as páginas estáticas dos mesmos códigos, já com os resultados (mais rápido): * [[https://nbviewer.jupyter.org/github/piLaboratory/bie5782/blob/master/jupyter/.ipynb_checkpoints/MLE%20exponencial-checkpoint.ipynb|MLE da distribuição exponencial]] * [[https://nbviewer.jupyter.org/github/piLaboratory/bie5782/blob/master/jupyter/MLE%20Poisson.ipynb| MLE da Poisson]] * [[https://nbviewer.jupyter.org/github/piLaboratory/bie5782/blob/master/jupyter/MLE%20geom%C3%A9trica.ipynb| MLE da distribuição geométrica]] * [[https://github.com/piLaboratory/bie5782|Repositório dos códigos no GitHub]] (( Contibuições são bem-vindas. Envie-nos um Pull request ou informe um [[https://github.com/piLaboratory/bie5782/issues| problema ou sugestão]] ))