Likelihood ou Verossimilhança.

Essa é uma palavra muito comum nos "outputs" de programas estatísticos, mas que eu demorei um pouco pra entender, se é que compreendo realmente.

Fazendo a disciplina oferecida pelo Brian Caffo no coursera (https://www.coursera.org/course/biostats) ele da um exemplo bem simples e legal.

Então vamos pensar pensar em uma moeda... na verdade não ja que moedas e bolinhas coloridas em caixinhas fazem muito pouco sentido pra mim :)

Voltando aquelas aranhas do post anterior, falamos da teoria de meta-populações, que é a forma que acho mais legal na biologia de "jogar moedas", a chance de um lugar estar ocupado ou não.

Vamos pensar que fomos olhar 4 arbustos, fomos la batemos com um pauzinho nele e em 3 arbustos caiu aranhas. Hummm
Então no primeiro estava ocupado por aranhas, o segundo idem, o terceiro não tinha aranhas e o quarto estava ocupado. Ou seja 3 arbustos com aranhas e 1 sem.

Então se pensarmos que cada arbusto tem uma chance de ter aranhas, tem sua chance de ocupação, qual a taxa de ocupação nesse exemplo?

Podemos perguntar da seguinte forma, se a taxa de ocupação for 50%, poderíamos ver esse resultado?

A resposta é  sim.
Usando o R, podemos simular essa situação usando a função rbinom, que simula amostras tiradas de uma distribuição binomial.
Vamos simular 100 exemplo pra ver se algum fica igual a nossa coleta de dados.

#Vamos somar todos os casos que tivemos 3 arbustos com aranhas
#e 1 arbusto sem aranha, ou seja a amostra tem tamanhos 4, e queremos
#aquelas que em 3 encontramos aranhas
sum(rbinom(100,4,prob=0.5)==3)

E temos:
Temos 22 casos idênticos ao que coletamos.

Mas pera ai, 3 arbustos com aranhas, um sem, talvez a chance de ocupação seja maior que 50%, e se repetirmos essa simulação com 60% de chance de ocupação?

sum(rbinom(100,4,prob=0.6)==3)

E vemos 37,  casos, olha ae muito mais casos ficam igual, estamos melhorando.

Mas podemos subir mais ainda, vamos tentar uma chance de ocupação de 80%

sum(rbinom(100,4,prob=0.8)==3)

Legal, agora vemos quase 50 casos, as vezes mais as vezes menos, claro que o número sempre varia, é uma simulação.

Mas a verosimilhança, ou likelihood, é exatamente isso, na verdade qualquer valor de chance de ocupação pode produzir resultados como os que observamos no campo, no entanto, a diferença é o quão comum veremos resultados idênticos aos que observamos no campo, na vida real para cada chance de ocupação que tentamos.

Ao realizar varias tentativas, vimos que algumas chances de ocupação dificilmente produziriam um resultado semelhante aquele que observamos, se tentarmos chances de ocupação ainda mais extremas, uma chance de ocupação de 20%, de 10%, ainda sim podemos observar algo como o que vimos no campo, mas é muito raro, pouco comum, então algo em torno de 80%  produz mais resultados semelhantes ao que vemos.

A verosimilhança é exatamente isso, essa exploração para saber qual o valor mais provável para um parâmetro baseado num conjunto de dados, nesse caso a chance de aranhas ocuparem os arbustos, e os nossos dados são exatamente os quatro arbustos que observamos.


Então nesse caso temos, a maior verosimilhança ali perto do 80%, o pico, mas a chance de ocupação poderia assumir outros valores também, mas 20% por exemplo seria pouco provável, olhem como da um valor baixo ali no gráfico.
E 1%? Poderia também, mas a chance é quase nenhuma. Ou seja usando o likelihood da pra resumir todo o processo acima a um gráfico.

Bem, eu não vou entrar muito em formulas aqui, até porque eu sou péssimo com elas, mas da pra ter acho que ja da pra ter uma ideia melhor do que diabos é o likelihood, sendo que isso é usado pra todo tipo de parâmetro estimado.

O script para esse grafico é assim:

#

thetaVals <- span="span"> seq(0, 1, length = 1000)

maxL <- span="span"> .75 ^ 3 * .25 ^1

lVals <- span="span"> thetaVals ^ 3 * (1 - thetaVals) ^1 / maxL

#postscript("coinLikelihood.eps", width=4, height=4, horizontal = FALSE)

plot(thetaVals, lVals,

     type = "l",

     frame = F,

     xlab = expression(theta),

     ylab = "Likelihood",

     lwd = 2)

lines(range(thetaVals[lVals > 1 / 8]), c(1/8,1/8))

lines(range(thetaVals[lVals > 1 / 16]), c(1/16,1/16))

#

Copiado diretamente da aula do Brian Caffo
Mais material legal dele aqui se você pensa melhor com moedas :)
https://github.com/bcaffo/Caffo-Coursera

Um pequeno exemplo de um modelo com Zeros Inflados

Modelos inflados por zeros podem ser muitos interessantes. Vou tentar fazer uma breve introdução com um exemplo simples. 

Vamos supor a seguinte situação, eu gostaria de saber se o número de aranhas presentes numa especie de arbusto A é diferente do número de aranhas comumente encontrado num arbusto B. Utilizando o número magico de coletas que me foi passado na faculdade, eu procurei trinta arbustos da especie A e trinta arbustos da especie B e contei quantas aranhas tinham neles:

#Simulando os dados

set.seed(13)
população<-factor(rep(c("A","B"),each=30))
contagem<-rpois(60,rep(c(2,4),each=30))

#

 Vemos que no Arbusto B, em média temos 4 aranhas, enquanto no arbusto A temos 2 aranhas.

Podemos testar essa afirmação através de um modelo de regressão de poisson.

#Modelo de Poisson

modelo1<-glm(contagem~população,family="poisson")
summary(modelo1)
exp(coef(modelo1))

#

E obtemos como resultado:

Call:
glm(formula = contagem ~ população, family = "poisson")

Deviance Residuals:
    Min       1Q   Median       3Q      Max 
-1.9833  -0.8255   0.0237   0.5868   2.9973 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   0.6763     0.1302   5.195 2.05e-07 ***
populaçãoB    0.6587     0.1604   4.107 4.01e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 90.678  on 59  degrees of freedom
Residual deviance: 72.885  on 58  degrees of freedom
AIC: 235.15

Number of Fisher Scoring iterations: 5

> exp(coef(modelo1))
(Intercept)  populaçãoB
   1.966667    1.932203 

Vejam só, o modelo nos da varias estrelinhas, nossa afirmação foi na mosca, mas mais importante que isso, temos estimativas do tamanho da população, e em média elas ficam bem próxima de realidade, o intercepto que representa a população nos arbustos A tem 1.966.. aranhas, bem próximo do 2 que usamos pra simular os dados, e no arbusto B temos 1.966 + 1.932.., que da praticamente 4, o numero que usamos na simulação. Muito legal, muitos p significativos e condecorados com variasssss estrelinhas.

Mas nem sempre o mundo é tão simples assim. Como nosso amigo Levins propôs uma vez, nem sempre todos os arbustos vão ter aranhas, eles podem não estar ocupados por aranhas. 

Pra gente contar aranhas em um arbusto, tem que ter mais de uma aranha la , ou a gente vai contar zero aranhas, até ai tudo bem pois alguns arbustos podem ser um péssimo lugar pra se viver e nenhuma aranha fica la, ou seja zero aranhas.

Mas podem existir um outro tipo de zero, o zero de aranhas devido a elas nunca terem chegado la, o arbusto pode ser perfeito pra abrigar aranhas, mas nenhuma chegou la. Agora vamos colocar isso nos dados e ver o que acontece com nosso modelo.

#Simulando ocupação

ocupação<-rbinom(60,size=1,prob=0.60)
contagem.ocup<-contagem*ocupação

#Modelo 2 de poisson

 modelo2<-glm(contagem.ocup~população,family="poisson")
summary(modelo2)
exp(coef(modelo2))

#

Call:
glm(formula = contagem.ocup ~ população, family = "poisson")

Deviance Residuals:
    Min       1Q   Median       3Q      Max 
-2.0656  -1.8719  -0.3133   0.6582   4.2838 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   0.4906     0.1429   3.434 0.000594 ***
populaçãoB    0.2671     0.1898   1.407 0.159458   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 160.14  on 59  degrees of freedom
Residual deviance: 158.14  on 58  degrees of freedom
AIC: 264.01

Number of Fisher Scoring iterations: 6

> exp(coef(modelo2))
(Intercept)  populaçãoB
   1.633333    1.306122

É muito triste, mas o fato de nem todos os arbustos estarem ocupados por aranhas nos tira muitas muitas estrelinhas, mas pior que isso, agora concluímos que os arbustos A e B tem o mesmo numero de aranhas, isso é mal, muito mal. Olha a estimativa do intercepto 1.63, usamos 2 na simulação, mas os zeros adicionais puxaram a media do numero de aranhas nos arbustos pra baixo.

Mas o que acontece é que nesse caso não estamos levando em conta que temos 2 tipos de zeros nas coletas, existem zeros de o arbusto ter zero aranhas, mas existem zeros devido ao arbusto não estar ocupado, aranhas nunca chegaram la.

O que deveríamos fazer é tentar fazer um modelo onde esses dois tipos de situação são levados em conta.

Temos de declarar ao modelo que existem zeros devido a uma chance de ocupação que cada arbusto tem, ai considerando os arbustos ocupados que deveríamos fazer as contas para o numero de aranhas no arbusto. Essa é a solução de modelos com zero inflados tenta trazer.

#modelo zero inflado usando o pacote pscl

library(pscl)
modelo3<-zeroinfl(contagem.ocup~população|1)
summary(modelo3)
exp(coef(modelo3)[1:2])
plogis(coef(modelo3)[3])

#

Call:
zeroinfl(formula = contagem.ocup ~ população | 1)

Pearson residuals:
    Min      1Q  Median      3Q     Max
-0.9876 -0.8934 -0.2306  0.4150  3.2704

Count model coefficients (poisson with log link):
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   0.8296     0.1540   5.388 7.13e-08 ***
populaçãoB    0.5977     0.2019   2.960  0.00307 **

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4567     0.2829  -1.615    0.106
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 8
Log-likelihood: -104.9 on 3 Df

> exp(coef(modelo3)[1:2])
count_(Intercept)  count_populaçãoB
         2.292305          1.817889

> plogis(coef(modelo3)[3])
zero_(Intercept)
       0.3877701

A primeira boa noticia é que nossas estrelinhas são devolvidas, temos muitos p significativos. Mas vamos olhar as nossas estimativas, agora estimamos para o arbusto A uma média de 2.29 e para o B 2.29 + 1.81, o que é bem próximo de 4, 2 e 4 foram exatamente os valores que usamos.

Olhamos ainda a estimativa de ocupação, é de 38%, tudo bem que usamos 60%, mas é bastante razoável, e temos que levar em consideração o número de amostras, mas voltamos a concluir certo que arbustos A tem em média 2 aranhas, e arbusto B o dobro.

A parte mais interessante, é que apesar de dois processos estarem acontecendo simultaniamente, um que é o ambiente influenciando no número de aranhas que encontramos e o outro é arbustos estarem ocupados ou não, estamos modelando esses dois processos relativamente bem. Resgatando estimativas razoáveis. Se ignoramos um processo, como ignoramos o processo de ocupação no modelo 2, além de péssimas estimativas do tamanho populacional, podemos concluir erroneamente o que esta acontecendo.

Mas ao pensarmos sobre todos os processos que influenciam as nossas observações, podemos selecionar uma estrategia que leve tudo em conta, todos os processos que temos teoria para subsidiar e gerar boas estimativas.

Nesse caso, pensando que o ambiente influencia o número de aranhas que encontramos, mas encontrar aranhas também depende da ocupação, como prevê a teoria de metapopulações, uma regressão de poisson com zero inflados funciona perfeitamente.  E o preço do modelo errado pode ser varias estrelinhas pior, uma conclusão errada do que está acontecendo.