Ferramentas do usuário

Ferramentas do site


03-funcao-veros:03-funcao-veros

3. Função de Verossimilhança

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

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:

{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/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?

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

Portanto: É 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: cria-se com ela uma função dos parâmetros, mantidos os dados constantes.

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)

Portanto: É 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" )

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.

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

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

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 use a função 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 tem várias propriedades de um 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:

  1. geramos um grande número de amostras de uma distribuição conhecida com valores conhecidos dos parâmetros;
  2. obtemos para cada uma das amostras simuladas as MLE;
  3. construímos gráficos da distribuição das MLE para analisar o seu comportamento.

Consistência

A 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 viciado:

$$ \hat{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.

Eficiência Assimptótica

A eficiência assimptó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 Assimptótica

Esta é a propriedade de distribuição assimptótica que garante que a distribuição dos MLE 3) 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 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 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 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 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


Será que a propriedade de não ser viciado se mantém quando a variância é transformada em desvio padrão ?

$$ \hat{s} = [ \hat{s^2} ]^{(1/2)} = \left[ \frac{ \sum_{i=1}^n (x_i - \bar{x})^2 }{ n - 1 } \right]^{(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 abundante 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 4) ?

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 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 notaR. Prazo: até a segunda-feira seguinte ao dia da aula deste tópico.

Questões motivadoras para a dicussão

Proponha questões para discutirmos em nosso subfórum. Criamos um subfórum para isso aqui.

Questões escolhidas

  1. 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?
  2. 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 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?
  3. 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

Complementares

  • 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 Apêndice A é uma excelente introdução).

Na Internet

Resumos de Aula

1)
veja também o ótimo tutorial deste pacote em http://finzi.psych.upenn.edu/R/library/bbmle/doc/mle2.pdf
2)
taxa média de ocorrência de plantas por parcelas, no caso.
3)
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.
4)
note que a pergunta é sobre o comportamento da estimativa, e não da população de medidas de produtividade de onde ela veio
03-funcao-veros/03-funcao-veros.txt · Última modificação: 2016/09/19 11:24 por paulo