Na seção anterior vimos que a função de verossimilhança liga nossos dados a modelos probabilísticos. Nesta unidade veremos como podemos usar a função de verossimilhança para encontrar os parâmetros de uma distribuição de probabilidade que melhor descreve os dados. Em outras palavras, vamos ajustar distribuições de probabilidade a conjuntos de dados. Estes são os modelos estatísticos mais simples.

• Distribuições de Probabilidade como Modelos
• Parâmetros e estimativas de máxima verossimilhança
• Ajuste por soluções analíticas de MLEs
• Ajuste por soluções numéricas de MLEs
• Diagnóstico de ajustes
• Comparação de ajustes

A distribuição geométrica é um modelo simples de ser ajustado a dados porque há uma expressão para o MLE de seu parâmetro:

$$\widehat p \ = \ \frac{N}{N+\sum x_i}$$

Onde $N$ é o total de observações e $\sum x_i$ é a soma dos valores observados.

Como chegamos a isso? Neste caso por meio de manipulações matemáticas. Quando isso nos leva a uma expressão de funções matemáticas conhecidas dizemos que temos uma solução analítica ¹⁾.

Para ter uma ideia mais concreta deste processo, vamos refazer passo a passo a dedução analítica do mle da geométrica, com o auxílio de sistemas de álgebra simbólica. Estes sistemas são programas que guardam muitas regras de manipulação de objetos matemáticos, como a álgebra. Com ele podemos conferir passo a passo a deduções matemáticas, mesmo sem lembrar de todos detalhes necessários em cada passagem da dedução. Você tem duas opções para realizar este tutorial:

1. Em Python, usando um Jupyter notebook, com a interface online Google Colab: clique aqui.

2 . Executar os comandos nesta página mesmo, usando o software Maxima.

Nas caixas abaixo há comandos do Maxima que são executados por um servidor remoto. Para enviar cada comando clique no botão Evaluate ²⁾.

O comando abaixo cria um objeto simbólico da função de densidade da geométrica.

Ao executar você terá de volta apenas a expressão, o que significa que ela pode ser manipulada no Maxima. Demos à expressão o nome P para facilitar a manipulação. Começamos aplicando o logaritmo à função de densidade

e então distribuímos os logaritmos ³⁾:

A log-verossimilhança é a soma da expressão acima aplicada para cada dado $x_i$:

O que é a representação no Maxima para

$$\sum_{i=1}^{N}\ \log p+x_i\,\log (1-p) $$

Como a somatória é apenas em $x_i$, os outros termos são tratados como constantes. Com isso a somatória restringe-se ao termo $\sum x_i$, que chamamos de $S$. O comando a seguir faz esta substituição de termos:

O resultado é a representação no Maxima da subsitutição de $\sum x_i$ por $S$ na expressão acima, que fica então:

$$N \log p\ +\ S\,\log (1-p) $$

O MLE $\widehat p$ é o valor do parâmetro $p$ que minimiza a função de log-verossimilhança. O truque para encontrar isso é achar o valor de $p$ que faz a derivada da função ser igual a zero. Primeiro encontramos a derivada em relação a $p$ da expressão simbólica LL, criada acima:

e finalmente igualamos esta expressão a zero e resolvemos para $p$:

Veja mais notebooks de deduções analíticas de MLEs aqui.

Ajustar uma distribuição a dados consiste em encontrar os valores mais plausíveis dos parâmetros da distribuição escolhida, condicionados aos dados. Ou seja, obter o valor do MLEs. Isso é simples quando há uma solução analítica de seus MLEs. Basta aplicar a expressão desta solução, que é uma função de quantidades obtidas dos dados. Mas além disso temos que analisar o comportamento da função de verossimilhança na vizinhança do MLE.

Vamos fazer isso ajustando a distribuição geométrica a dados de sobrevivência de aves Vanellus vanellus. O conjunto de dados é um data frame com o número de anos que indivíduos anilhados na Grã-Bretanha sobreviveram ⁴⁾

Comece importando os dados para um objeto no R:

(vanellus <- read.csv2("vanellus.csv"))

Note que o data frame tem cada valor de número de anos sobrevividos após a anilhagem e o número de aves que viveu até cada um.

Agora vamos calcular o valor do MLE da geométrica para estes dados, usando sua expressão analítica⁵⁾ :

(n <- sum(vanellus$freq)) #total de observações
(soma.x <- sum(vanellus$tempo*vanellus$freq)) #soma dos tempos
(p <- n / (n+soma.x)) #mle

Para facilitar a construção do gráfico da curva de verossimilhança, vamos criar no R uma função de log-verossimilhança negativa ⁶⁾:

nllgeom <- function(p, n, sx, rel=FALSE, p.mle){
    nll = -n*log(p) - log(1-p) * sx
    if(rel)
        nll = nll - nllgeom(p.mle, n, soma.x)
    nll
}

Com essa função é fácil criar um gráfico com a curva de verossimilhança e traçar o intervalo de verossimilhança (razão de 8). Basta usá-la para calcular a verossimilhança para uma sequência de valores do parâmetro $p$, em torno do valor do MLE:

curve(nllgeom(x, n, soma.x, rel=TRUE, p.mle=p), 0.29, 0.4,
      col="blue", xlab="p morte Vanellus vanellus",
      ylab="Log-Veros. Neg.", cex.lab=1.6, cex.axis=1.4)
abline(v=p, col="red")

A linha vertical marca o MLE. Para delimitar o intervalo de verossimilhança, ou de plausibilidade, marque uma linha na altura de 1/8 da verossimilhança relativa:

abline(h=log(8), col="red", lty=2) 

Agora que você conhece a distribuição geométrica mais plausível para este conjunto de dados você pode comparar o que ela prevê com o que você observou. Isto é um maneira de avaliar a qualidade do ajuste, sobrepondo o número previsto de sobreviventes a cada ano com as frequências observadas:

plot(freq~tempo, data=vanellus, type="h",
     xlab="Tempo sobrevivência", ylab="Frequência")
x <- 0:max(vanellus$tempo)
points(x, dgeom(x, p)*n, col="blue")

Soluções analíticas são muito raras, e não existem para muitos MLEs, Por exemplo, não é possível obter solução analítica para as MLEs dos parâmetros da distribuição Weibull. Nesses casos usamos aproximações numéricas, em geral com ajuda de computadores ⁷⁾

Vamos usar dados de DAP de árvores de floresta de Paragominas, PA: (parago-sobrev.csv). Utilizaremos apenas das árvores com mais de 25cm de DAP (variável “dap” > 25) :

##Leitura dos dados
parag<- read.table("http://cmq.esalq.usp.br/BIE5781/lib/exe/fetch.php?media=04-parametros-constantes:parago-sobrev.csv", header=TRUE)
# Seleciona daps > 25
dap = parag$dap[ parag$dap > 25 ]
par(mfrow=c(1,2))
hist(dap, xlim=c(0,max(dap)))

Para evitarmos o problema do parâmetro de locação da Weibull, vamos trabalhar como se o DAP de 25cm fosse o ponto inicial (zero):

dap0 = dap - 25
hist(dap0, xlim=c(0,max(dap)))
par(mfrow=c(1,1))

Para ajustarmos a dist. Weibull é necessário construir a função de log-verossimilhança negativa no R:

nllweibull = function(escala, forma)
	-sum(dweibull(x =dap0, shape=forma, scale=escala, log=TRUE))

Obtemos o ajuste com a função 'mle2':

parag.wei = mle2(nllweibull, start=list(escala=20, forma=1))

Se quisermos olhar a superfície de log-verossimilhança negativa é necessário vetorizar a função log-verossimilhança negativa ⁸⁾:

nllweibull.V = Vectorize( nllweibull, c("escala","forma") )

Com a função vetorizada podemos construir um gráfico de linhas de contorno:

dx = seq(11,16, length=50)
dy = seq(0.8, 1.1, length=50)
z = outer( dx, dy, "nllweibull.V")
z <- z-min(z)# verossimilhança padronizada
contour(dx, dy, z, xlab="Escala", ylab="Forma", col="blue", cex.lab=1.6, cex.axis=1.4, lwd=2, nlevels=25)
contour(dx, dy, z, levels=2, col="red", add=T)
abline(v=coef(parag.wei)[1], lty=2, col="red")
abline(h=coef(parag.wei)[2], lty=2, col="red")

As curvas de verossimilhança perfilhada são obtidas perfilhando o objeto mle, e utilizando a função plotprofmle do pacote sads:

parag.wei.prof = profile(parag.wei)
par(mfrow=c(1,2))
plotprofmle(parag.wei.prof)
par(mfrow=c(1,1))

Também não há solução analítica para o parâmetro de disperão da binomial negativa. Neste tutorial comparamos os ajustes da binomial negativa e da Poisson⁹⁾ a dados de contagens.

Vamos retomar os tutoriais sobre distribuições discretas em que simulamos uma distribuição agregada de plantas em uma área dividida em quadrículas (parcelas):

set.seed(42)
cx <- runif(20,0,20)
cy <- runif(20,0,20)
px <- rnorm(2000)
py <- rnorm(2000)
x1 <- cx+px
y1 <- cy+py
x2 <- x1[x1>0&x1<20&y1>0&y1<20]
y2 <- y1[x1>0&x1<20&y1>0&y1<20]
x2.parc <- cut(x2,breaks=0:20, labels=1:20)
y2.parc <- cut(y2,breaks=0:20, labels=LETTERS[1:20])

O numero de plantas por parcela é obtido com:

cont2 <- data.frame(table(x2.parc,y2.parc))$Freq

A avaliação visual feita no tutorial anterior indica que a variável binomial negativa é uma descrição mais adequada destes dados, comparada com a Poisson. Vamos ajustar os dois modelos e usar a função de verossimilhança para confirmar isso.

Primeiro definimos funções de log-verossimilhança negativa para cada modelo:

LL.pois <- function(lam){
  -sum(dpois(cont2,lambda=lam,log=T))
}
 
LL.nbin <- function(media,k){
  -sum(dnbinom(cont2,mu=media,size=k,log=T))
}

Em seguida minimizamos numericamente essas funções com o R ¹⁰⁾. Para isto carregue o pacote bbmle e use a função mle2. É preciso fornecer valores iniciais razoáveis, no argumento start:

library(bbmle) # basta uma vez por seção 
mod1 <- mle2(LL.pois,start=list(lam=mean(cont2)))
mod2 <- mle2(LL.nbin,start=(list(media=mean(cont2),k=0.1)))

Como esperado, a binomial negativa é um modelo muito mais plausível:

> logLik(mod1)
'log Lik.' -1786.502 (df=1)
> logLik(mod2)
'log Lik.' -1014.969 (df=2)

Você pode obter um resumo dos modelo com o comando summary:

summary(mod1)
summary(mod2)

E podemos fazer um gráfico dos valores previsto pelos dois modelos com:

##MLEs de cada modelo
(cf1 <- coef(mod1))
(cf2 <- coef(mod2))
 
##grafico
cont2.f <- factor(cont2, levels=0:max(cont2))
plot(table(cont2.f)/400, xlab="N de indivíduos na parcela", ylab="Proporção das parcelas")
points(x=0:max(cont2),y=dpois(0:max(cont2),lambda=cf1), type="b", col="blue", lty=2)
points(x=0:max(cont2),y=dnbinom(0:max(cont2),mu=cf2[1],size=cf2[2]), type="b", col="red", lty=2)

Por fim, avalie o perfil de verossimilhança dos dois parâmetros da binomial negativa ¹¹⁾:

library(sads)
p.mod2 <- profile(mod2)
par(mfrow=c(1,2))
plotprofmle(p.mod2)
par(mfrow=c(1,1))

Como avaliar o ajuste de uma distribuição de probabilidade aos nossos dados? O gráfico quantil-quantil (ou qqplot) oferece um critério simples: os quantis empíricos e previstos pelo modelo de distribuição devem ser iguais.

O padrão esperado são pontos ao longo da linha de equivalência (1:1) em um gráfico de dispersão dos quantis observados em função dos teóricos. Como somos capazes de perceber desvios pequenos em relação a retas, este critério visual é mais sensível do que muitos testes estatísticos de aderência.

Apesar da regra de decisão simples, a maneira de construir o qqplot é menos intuitiva. Neste tutorial vamos refazê-la passo a passo, para fins didáticos ¹²⁾.

Para começar vamos perguntar ao R quais são os quantis empíricos, que são os valores que delimitam uma certa proporção dos dados. O mais conhecido deles é a mediana, o valor que é menor do que metade dos dados ¹³⁾. A mediana é o quantil de 50%.

Vamos usar os objetos criados no tutorial anterior de ajuste de contagens à Poisson e binomial negativa. Começamos calculando os decis dos dados, que são os quantis que delimitam cada décimo do total de valores. A função do R que calcula quantis empíricos é a quantile. Vamos usá-la para obter os quantis das contagens que estão no objeto cont2.

(qqs <- data.frame(decis =seq(.1,.9,by=.1)))
## Quantis observados
qqs$q.obs <- quantile(cont2, qqs$decis)
qqs

Agora adicionamos os quantis esperados para os mesmos decis, pelos modelos Poisson e binomial ajustados:

qqs$q.pois <- qpois(qqs$decis, lambda=cf1)
qqs$q.bn <- qnbinom(qqs$decis, mu=cf2[1], size=cf2[2])
qqs

Qual dos dois modelos faz uma melhor previsão dos decis de 10% a 90%? O gráfico quantil-quantil confronta estes valores:

par(mfrow=c(1,2))
plot(q.obs~q.pois, data=qqs,
     xlab="Decis previstos",
     ylab="Decis observados",
     main=paste("Poisson, lambda =",round(cf1,2)))
abline(0,1, col="red")
plot(q.obs~q.bn, data=qqs,
     xlab="Decis previstos",
     ylab="Decis observados",
     main=paste("Bin. Neg., mu =",round(cf2[1],2),
         ", k =", round(cf2[2],2)))
abline(0,1, col="red")
par(mfrow=c(1,1))

Os dados deste tutorial têm 400 observações, mas fizemos um gráfico com apenas nove valores. Para usar toda a informação contida nos dados precisamos calcular o quantil esperado para cada um. Começamos ordenando os valores e depois usamos a função ppoints para atribuir um percetil a cada valor:

## dados ordenados:
cont3 <- sort(cont2)
## percentis
cont3.p <- ppoints(cont3)

Em seguida calculamos os quantis esperados por cada modelo

## quantis da poisson associados a cada percentil
cont3.qp <- qpois(cont3.p, lambda=cf1)
## quantis da bn associados a cada percentil
cont3.qbn <- qnbinom(cont3.p, mu=cf2[1], size=cf2[2])

E finalmente fazemos os gráficos

par(mfrow=c(1,2))
plot(cont3.qp, cont3, xlab="Qantis previstos",
     ylab="Quantis observados",
     main=paste("Poisson, lambda =",round(cf1,2)))
abline(0,1, col="red")
plot(cont3.qbn, cont3, xlab="Qantis previstos",
     ylab="Quantis observados",
     main=paste("Bin. Neg., mu =",round(cf2[1],2), ", k =", round(cf2[2],2)))
abline(0,1, col="red")
par(mfrow=c(1,1))

Há ainda detalhes técnicos como a maneira de definir os percentis empíricos ou traçar a linha de equivalência. Há muitas funções no R com variações do qqplot. Uma bastante flexível é a qqPlot, do excelente pacote car:

library(car)
par(mfrow=c(1,2))
qqPlot(cont2, distribution="pois",
       lambda=cf1, line="robust")
qqPlot(cont2, distribution="nbinom",
       mu=cf2[1], size=cf2[2], line="robust")
par(mfrow=c(1,1))

É possível definir uma função de verossimilhança para dados agregados em intervalos de classe?

Proponha um código em R para fazer o ajuste de máxima verossimilhança a uma log-normal dos dados abaixo:

Definindo:

• $F_i$ : a frequência de observações na classe $i$
• $a_i$, $b_i$: os limites inferior e superior do intervalo de classe $i$
• $N$: número de classes
• $f(x)$: uma função de densidade probabilística

Proponha uma expressão genérica para a função de log-verossmilhança de dados agregados em classes.

==Um exemplo em R== * Código mostrado em sala

• Explique e exemplifique o que são dados truncados e censurados.
• Como ajustar um modelo de distribuição a dados truncados e censurados com o método de máxima verossimilhança?

• Explique e exemplifique o que são dados de contagens com excesso de zeros
• Como ajustar um modelo de distribuição a dados com excesso de zeros com o método de máxima verossimilhança?

Imagine um experimento de predação similar ao descrito no item 6.3.1.1. de Bolker (2008), mas que tenha duas quantidades de girinos: 10 e 100 por tanque. Como ajustar uma distribuição binomial ao número de girinos predados permitindo que a probabilidade de predação dependa do número inicial de girinos no tanque? Outros detalhes sobre o sistema de estudo estão na seção 2.5.2 do livro.

Faça os exercícios desta unidade (204.1 e 204.2) no sistema notaR.

• Likelihood and All That, seções 6.1 a 6.2.1.2 do Capítulo 6 de: Bolker, B.M. 2008 Ecological Models and Data in R Princeton: Princeton University Press.

• Anderson, D. R. (2008). Model based inference in the life sciences: a primer on evidence. New York, Springer: capítulos 2 e 3.

• Um roteiro de ajuste de modelos na página do Laboratório de Ecologia Teórica do IB-USP.
• Outro roteiro, bem mais resumido, da disciplina de introdução à modelagem com verossimilhança de Chad Brassil.
• Excelentes exercícios de simulação e ajustes de distribuições no site de apoio de Bolker (2008). ¹⁴⁾

Nossos notebooks interativos em Python demonstrando os passos para a dedução de alguns estimadores de máxima verossimilhança (MLEs). Você pode executá-los online pelo Google Colab:

• MLE da distribuição exponencial
• MLE da Poisson
• MLE da Geométrica.

Ou você pode ver as páginas estáticas dos mesmos códigos, já com os resultados:

• MLE da distribuição exponencial
• MLE da Poisson
• MLE da distribuição geométrica

Todos estes notebooks estão no Repositório da disciplina no GitHub ¹⁵⁾

• Nadarajah, S., & Kotz, S. (2006). R Programs for Truncated Distributions. Journal of Statistical Software, 16(Code Snippet 2), 1 - 8
• Frederick Novomestky and Saralees Nadarajah (2016). truncdist: Truncated Random Variables. R package https://CRAN.R-project.org/package=truncdist