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.

• 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

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,

A força de evidência em favor de $H_A$ vis-a-vis $H_B$ é dada pela razão de verossimilhança:

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.

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

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

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)

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 )$$

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")

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.

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.

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 )

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.

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 bbmle²⁾.

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.

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

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

Os estimadores de máxima verossimilhança têm 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.

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

A 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) )

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

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" )

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.

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 

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

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:

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 ⁵⁾ ?

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

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")

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?

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")

Faça os exercícios 203.1 a 203.3 no sistema notaR.

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

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

• Batista, J.L.F. 2009 Verossimilhança e Máxima Verossimilhança.

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

• Índice do pacote bbmle na página do R.
• Índice do emdbook na página do R. Este é um pacote de apoio ao livro Ecological Models and Data in R.
• tutorial do pacote bbmle. Ótimo roteiro para aplicação do método de máxima verossimilhança em R.
• "Likelihood": Uma palestra de A. W. F. Edwards, um dos pioneiros na formalização e defesa do princípio da verossimilhança.
• Maximum Likelihood Primer por S. Purcell. Feito para aplicações em genética, mas uma ótima introdução ao tema.
• Renata e a máxima verossimilhança, explicação super didática, do projeto Ela está em tudo (Sim, Matemática é coisa de menina) , do IME -USP e UFABC.
• Scott, C., & Nowak, R. 2004. Maximum Likelihood Estimation. Connexions, May 12, 2004.

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 Binder. Demora um pouco para carregar. As instruções passo a passo estão na seção "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):
  • MLE da distribuição exponencial
  • MLE da Poisson
  • MLE da distribuição geométrica
• Repositório dos códigos no GitHub ⁶⁾