Ferramentas do usuário

Ferramentas do site


biometria:verossim:02-veros


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:

  • optimize: otimização unidimensional, isto é, para funções com apenas um parâmetro livre.
  • 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 mle do pacote stats4, e a mle2, do pacote bbmle1).

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 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 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 2) ?

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 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. 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.



Recursos para Estudo

Resumo da Aula

Soluções dos Exercícios

Na Internet

Pesquisa

Indique os tutoriais e partes da leitura básica que merecem mais atenção na discussão aqui.

1)
veja também o ótimo tutorial deste pacote em http://finzi.psych.upenn.edu/R/library/bbmle/doc/mle2.pdf
2)
note que a pergunta é sobre o comportamento da estimativa, e não da população de medidas de produtividade de onde ela veio
biometria/verossim/02-veros.txt · Última modificação: 2022/11/24 14:21 por 127.0.0.1