Ferramentas do usuário

Ferramentas do site


04-parametros-constantes:04-parametros-constantes

4. Modelos com Parâmetros Constantes

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 conjutnos de dados. Estes são os modelos estatísticos mais simples.

Conceitos

  • 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

Tutoriais

Distribuição Geométrica: dedução analítica do mle

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:

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

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

Como sabemos 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 1),

Vamos refazer passo a passo a dedução analítica do mle da geométrica, com o auxílio de um sistema de álgebra simbólica, o Maxima. Nas caixas abaixo há comandos do Maxima que são excutados por um servidor remoto. Para enviar cada comando clique no botão Evaluate 2).

A função de log-verossimilhança

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 logarítmo à função de densidade

e então distribuímos os logarítmos 3):

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

Solução de máxima verossimilhança

O MLE $\hat 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$:

Ajuste da Geométrica

É simples ajustar modelos que têm uma solução analítica de seus MLEs. Basta aplicar a expressão. 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 4)

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.

Obtendo a MLE do parâmetro $p$ da distribuição geométrica:

(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 a 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):

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")
abline(h=log(8), col="red", lty=2)

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

Ajuste numérico da Weibull

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

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

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:

x = dap0
nllweibull = function(escala, forma, x=dap0)
	-sum(dweibull(x, 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 6):

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

Ajuste numérico da Poisson e binomial negativa

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 Poisson7) a dados de contagens.

Vamos retomar os tutoriais sobre o tutorial 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 8). 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 9):

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

Gráficos quantil-quantil

Como avaliar o ajuste de uma distribuição de probabilidade aos nosos 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 10).

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

CODA

Há ainda detalhes téncnicos 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))




Questões motivadoras para discussão

1. Dados agregados

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

a) Na prática

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

ClasseFrequência
0–2 213
2–4 52
4–6 12
6–8 1
8–10 1
10–12 1
b) generalizando

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

2. Dados censurados e truncados

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

3. Binomial com parâmetro p variável

Imagine um experimento de predação similar ao descrito no item 6.2.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?

Exercícios

Os exercícios desta unidade (204.1 e 204.2) estão disponíveis no sistema notaR. Prazo para entrega: até a segunda-feira seguinte à aula desta unidade.

Recursos para Estudo

Leituras

Principal

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

Complementar

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

Na Internet

Resumo da Aula

1)
Isso é bem raro, em geral apenas para casos simples.
2)
Você pode também instalar o Maxima em seu computador e executar estes comandos a partir deste arquivo: mle_geometrica.wxm
3)
lembre-se que $\log(xy)=\log x+\log y$ e que $\log x^a = a \log x$
4)
Haldane, J. B. S. (1953). Some animal life tables. Journal of the Institute of Actuaries, 83-89.
5)
veja o tutorial sobre função de verossimilhança para entender a lógica da otimização computacional em R.
6)
detalhes no help da função Vectorize
7)
cujo mle tem solução analítica
8)
veja o tutorial sobre função de verossimilhança para entender essa otimização computacional.
9)
para isto você vai precisar da função plotprofmle do pacote sads
10)
No dia a dia você pode usar funções do R que já fazem tudo isso automaticamente
11)
ou maior que metade dos dados, se você preferir ;-)
12)
Bolker, B. (2008). Ecological Models and Data in R. Princeton, Princeton University Press.
04-parametros-constantes/04-parametros-constantes.txt · Última modificação: 2016/09/19 20:06 por paulo