Como já foi visto no tutorial sobre a Função de Verossimihança, a Lei da Verossimilhança pode ser enunciada da seguinte forma:
Dada uma variável aleatória $X$, cujo comportamento pode ser explicado por duas hipóteses: $H_A$ e $H_B$.
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:
$$ \frac{p_A(x)}{p_B(x)}$$.
Tomemos o exemplo de um laboratório que realizou o seguinte experimento: um certo produto foi ministrado a $20$ cobaias para se comparar duas hipóteses:
Um ponto importante é que a observação do número de cobaias mortas é que irá definir qual hipótese é favorecida e qual hipótese é desfavorecida.
Vejamos as probabilidades que a hipótese $H_A$ estabelece para cada uma das observações possíveis (1, 2, …, 20):
pa = dbinom(0:20, 20, prob=0.5) barplot(pa, width=1, space=0.1, col="blue") axis(1, 1, label=0:20, at=0:20*1.1+(1.1/2) )
No caso da hipótese $H_B$ temos:
pb = dbinom(0:20, 20, prob=0.3) barplot(pb, width=1, space=0.1, col="red") axis(1, 1, label=0:20, at=0:20*1.1+(1.1/2) )
A Razão de Verossimilhança para as observações possíveis pode ser facilmente obtida:
raz <- pa/pb barplot(raz, width=1, space=0.1, col="darkgreen") axis(1, 1, label=0:20, at=0:20*1.1+(1.1/2) )
A escala da Razão de Verossimilhança pode facilmente nos confundir. É melhor olharmos a Razão de Verossimilhança em termos da diferença da Log-Verossimilhança Negativa (com a indicação da Razão de 8, como limite definidor):
barplot(-log(raz), width=1, space=0.1, col="darkgreen") axis(1, 1, label=0:20, at=0:20*1.1+(1.1/2) ) abline( h = c(log(8), -log(8)), col="red", lty=2 )
Como a transformação inclui a mudança de sinal, a interpretação é que os valores positivos favorecem a hipótese $H_B$, enquanto que os valores negativos favorecem a hipótese $H_A$.
Resultado:
CONCLUSÕES:
Os dados também indicam, através da verossimilhança, qual o modelo mais plausível numa situação de comparação. O arquivo caxeta-3parcelas.csv tem os dados de DAP (diâmetro à altura do peito) de árvores em três parcelas instaladas em “caixetais” na região do Vale do Ribeira (São Paulo).
Comparemos a distribuição Weibull e a distribuição Gama, como modelos para a distribuição do DAP em cada uma das parcelas.
Primeiramente ler os dados e carregar o pacote “MASS” que possui a função “fitdistr” para ajuste de distribuições por máxima verossimilhança:
cax3p = read.csv("caxeta-3parcelas.csv",header=T,as.is=T,sep=";") library(MASS)
O segundo passo é ajustar os modelos (Weibull e Gama) para cada parcela.
# Ajuste da Dist. Weibull para as três parcelas weib1 = fitdistr( cax3p$dap[ cax3p$parcela==1 ] - 47, "weibull", start=list(shape=1, scale=10) ) weib2 = fitdistr( cax3p$dap[ cax3p$parcela==2 ] - 47, "weibull", start=list(shape=1, scale=10) ) weib3 = fitdistr( cax3p$dap[ cax3p$parcela==3 ] - 47, "weibull", start=list(shape=1, scale=10) ) # Ajuste da Dist. Gamma para as três parcelas gamm1 = fitdistr( cax3p$dap[ cax3p$parcela==1 ] - 47, "gamma", start=list(shape=1, scale=10) ) gamm2 = fitdistr( cax3p$dap[ cax3p$parcela==2 ] - 47, "gamma", start=list(shape=1, scale=10) ) gamm3 = fitdistr( cax3p$dap[ cax3p$parcela==3 ] - 47, "gamma", start=list(shape=1, scale=10) )
Comparação dos modelos nas parcelas uma a uma:
# Comparação Parcela 1 hist( cax3p$dap[ cax3p$parcela==1 ], prob = TRUE ) curve( dweibull(x, shape=weib1$estimate["shape"] , scale=weib1$estimate["scale"]), 40, 250, col="red", add=TRUE) curve( dgamma(x, shape=gamm1$estimate["shape"] , scale=gamm1$estimate["scale"]), 40, 250, col="blue", add=TRUE) AIC(weib1) - AIC(gamm1) # Comparação Parcela 2 hist( cax3p$dap[ cax3p$parcela==2 ], prob = TRUE ) curve( dweibull(x, shape=weib2$estimate["shape"] , scale=weib2$estimate["scale"]), 40, 400, col="red", add=TRUE) curve( dgamma(x, shape=gamm2$estimate["shape"] , scale=gamm2$estimate["scale"]), 40, 400, col="blue", add=TRUE) AIC(weib2) - AIC(gamm2) # Comparação Parcela 3 hist( cax3p$dap[ cax3p$parcela==3 ], prob = TRUE ) curve( dweibull(x, shape=weib3$estimate["shape"] , scale=weib3$estimate["scale"]), 40, 400, col="red", add=TRUE) curve( dgamma(x, shape=gamm3$estimate["shape"] , scale=gamm3$estimate["scale"]), 40, 400, col="blue", add=TRUE) AIC(weib3) - AIC(gamm3)
Questões:
O Princípio da Verossimilhança estabelece que a Função de Verossimilhança contem TODA evidência contida nos dados a respeito de uma dada hipótese. Esta hipótese pode ser um modelo, ou o valor de parâmetro de um modelo. Assim, a Razão de Verossimilhança é uma medida absoluta da força de evidência na comparação de duas hipóteses, independentemente do conjunto de dados utilizado ou das hipóteses comparadas.
Voltemos ao exemplo da aplicação de um produto em cobaias para verificar a taxa de mortalidade. Nesse caso temos dois laboratórios que testaram o produto:
A questão principal agora é saber qual o valor mais plausível para o parâmetro $p$, que indica a probabilidade de morte das cobaias.
Vejamos as curvas de verossimilhança para nos dois laboratórios:
p = seq(0.01, 0.99, by=0.01) lik.binom = dbinom(6, 20, p) # Lab 1: dist. Binomial lik.nbinom = dnbinom(14, 6, p) # Lab 2: dist. Binomial Negativa plot(p, lik.binom, type="l", ylab="Verossimilhança", col="blue") lines(p, lik.nbinom, col="red")
Aparentemente as curvas não são as mesmas. Mas devemos nos lembrar que o Princípio da Verossimilhança generaliza a Razão de Verossimilhança. Portanto devemos apresentar a curva de Verossimilhança Relativa 1), que obtemos dividindo cada valor de verossimilhança pelo máximo:
lik.binom = lik.binom / max(lik.binom) lik.nbinom = lik.nbinom / max(lik.nbinom) plot(p, lik.binom, type="l", ylab="Verossimilhança Relativa", col="blue", lwd=8) lines(p, lik.nbinom, col="red", lwd=2)
CONCLUSÕES:
Consideremos o mesmo exemplo das cobaias, mas comparemos o primeiro laboratório com outros dois laboratórios que possuem mais recursos para o experimento:
Vejamos as curvas de verossimilhança desses 3 laboratórios:
p = seq(0.01, 0.99, by=0.01) lik.binom1 = dbinom(6, 20, p) # Lab 1: dist. Binomial lik.binom2 = dbinom(60, 200, p) # Lab 2: dist. Binomial lik.binom3 = dbinom(600, 2000, p) # Lab 3: dist. Binomial plot(p, lik.binom1, type="l", ylab="Verossimilhança", col="blue") lines(p, lik.binom2, col="red") lines(p, lik.binom3, col="green")
Vejamos as curvas de verossimilhança RELATIVA desses 3 laboratórios:
lik.binom1 = lik.binom1/ max(lik.binom1) lik.binom2 = lik.binom2/ max(lik.binom2) lik.binom3 = lik.binom3/ max(lik.binom3) plot(p, lik.binom1, type="l", ylab="Verossimilhança Relativa", col="blue") lines(p, lik.binom2, col="red") lines(p, lik.binom3, col="green")
Façamos um “zoom” na curva de log-verossimilhança negativa relativa na vizinhaça da MLE:
nlik.binom1 = -log(lik.binom1) nlik.binom2 = -log(lik.binom2) nlik.binom3 = -log(lik.binom3) plot(p, nlik.binom1, type="l", ylab="Log-Veros. Neg. Relativa", col="blue", xlim=c(0.2,0.4)) lines(p, nlik.binom2, col="red") lines(p, nlik.binom3, col="green")
Questões:
A consequência imediata da combinação da Lei e do Princípio da Verossimilhança é que a função de verossimilhança, ou mais especificamente a função de log-verossimilhança negativa é o suporte necessário e suficiente à construção da inferência estatística:
Partindo de um modelo assumido como apropriado, qualquer inferência sobre os parâmetros do modelo, ou funções desses parâmetros2) pode ser realizada tendo a função de log-verossimilhança negativa como suporte.
Voltemos ao exemplo da distribuição de DAP no caxetal (parcela 2):
hist( cax3p$dap[ cax3p$parcela==2 ], prob = TRUE ) curve(dweibull(x, shape=weib2$estimate["shape"], scale=weib2$estimate["scale"]), 40, 400, col="red", add=TRUE)
Vejamos a superfície de log-verossimilhança negativa relativa para inferência sobre os parâmetros:
lweibull = function(forma, escala) -sum(dweibull( cax3p$dap[ cax3p$parcela==2 ]-47, shape=forma, scale=escala, log=TRUE)) vlweibull = Vectorize( lweibull, c("forma","escala") )
forma = seq(0.5, 2.5, by=0.05) escala = seq( 50, 100, by=0.5 )
sup.weibull = outer( forma, escala, vlweibull ) sup.weibull = sup.weibull - min(sup.weibull)
contour(forma, escala, sup.weibull, xlab="Forma", ylab="Escala", col="purple", levels=c(5,10,20,40,80,120))
abline(v=weib2$estimate[1], col="red", lty=2) abline(h=weib2$estimate[2], col="red", lty=2)
contour(forma, escala, sup.weibull, levels=log(8), labels="log(8)", add=T, col="darkgreen" )
Na distribuição Weibull a média (valor esperado) é definido em função dos parâmetros da seguinte forma: $$ \mu = \beta \ \Gamma\left( \frac{\gamma + 1}{\gamma} \right)$$. onde:
Assim podemos construir uma superfície para inferência sobre a Média:
mean.weibull = function(c, b) (b*gamma( (c+1)/c )+47)/10 sup.mean = outer(forma, escala, mean.weibull)
contour(forma, escala, sup.mean, col="darkgreen", xlab="Forma", ylab="Escala",levels=seq(8,18,by=1)) abline(v=weib2$estimate[1], col="red", lty=2) abline(h=weib2$estimate[2], col="red", lty=2)
contour(forma, escala, sup.weibull, levels=log(8), labels="log(8)", add=T, col="orange" ) media.estimada = (weib2$estimate[2] * gamma(1 + (1/weib2$estimate[1])) + 47)/10 contour(forma, escala, sup.mean, levels=media.estimada, col="red", add=T) contour(forma, escala, sup.mean, levels=mean(dap)/10, col="blue", add=T)
Na distribuição Weibull os quantis podem ser determinados a partir da função inversa da função de distribuição: $$d_p = \beta\ \left( \log \frac{ 1}{1-p} \right)^{1/\gamma}$$, onde:
Para construir a superfície para inferência sobre o quantil 95% basta seguir os mesmos passos da construção da superfície sobre a média:
dap95.weibull = function(c, b) (b*( log(1/(1-0.95)) )^(1/c) + 47)/10 sup.dap95 = outer(forma, escala, dap95.weibull) contour(forma, escala, sup.dap95, col="darkseagreen", levels=c(seq(15,25,by=1),seq(30,60,by=10))) abline(v=weib2$estimate[1], col="red", lty=2) abline(h=weib2$estimate[2], col="red", lty=2) contour(forma, escala, sup.weibull, levels=log(8), labels="log(8)", add=T, col="orange" ) dap95.est = dap95.weibull( weib2$estimate[1], weib2$estimate[2]) contour(forma, escala, sup.dap95, levels=dap95.est, labels=round(dap95.est,1), add=T, col="red" ) dap95.amostral = quantile(cax3p$dap[ cax3p$parcela==2 ]/10, 0.95, type=6) contour(forma, escala, sup.dap95, levels=dap95.amostral, labels=round(dap95.amostral,1), add=T, col="blue" )
A inferência sobre modelos consiste na comparação dos modelos dois-a-dois através da razão de verossimilhança.
É comum se utilizar o “Akaike Information Criterion” (AIC) como forma de comparação de modelos, sendo que além de considerar a verossimilhança do modelo ele penaliza o número de parâmetros no modelo.
Elementos que tornam essa abordagem mais simples para inferência sobre modelos quando comparada à abordagem “frequentista” ou “bayesiana” são:
Voltemos ao exemplo das 3 parcelas em caxetais. Foram ajustados as distribuições Weibull e gama para cada uma das parcelas. Assim, podemos usar o AIC para uma comparação parcela-a-parcela
library(bbmle) AICtab(weib1, gamm1, base=TRUE, logLik=TRUE) AICtab(weib2, gamm2, base=TRUE, logLik=TRUE) AICtab(weib3, gamm3, base=TRUE, logLik=TRUE)
Podemos pensar no conjunto dos três ajustes como um só modelo para as três parcelas, com seis parâmetros. Como a log-verossimilhança é aditiva, o AIC para este modelo combinado é a soma do AICs dos modelos componentes:
AIC(weib1) + AIC(weib2) + AIC(weib3) AIC(gamm1) + AIC(gamm2) + AIC(gamm3)
Um modelo mais parcimonioso para as três parcelas seria ajustar um só modelo para os dados agregados.
weib.agr = fitdistr( cax3p$dap - 47, "weibull", start=list(shape=1, scale=10) ) gamm.agr = fitdistr( cax3p$dap - 47, "gamma", start=list(shape=1, scale=10) ) AICtab(weib.agr, gamm.agr, base=TRUE, logLik=TRUE) AIC(weib1) + AIC(weib2) + AIC(weib3) AIC(gamm1) + AIC(gamm2) + AIC(gamm3)
Questões:
Na abordagem frequentista a inferência de intervalo é realizada através do Intervalo de Confiança.
O intervalo de confiança apela para o conceito de probabilidade a longo prazo que implica na repetição ilimitada do procedimento utilizado para gerar os dados, como se os dados fossem uma amostra de uma população infinita de observações possíveis. Dessa forma a construção do intervalo de confiança segue os seguintes passos:
Considere que numa plantação de Eucalyptus grandis foram amostradas aleatoriamente 100 árvores, sendo que dessas 37 se mostraram infectadas por uma certa doença. O objetivo é estimar a taxa de ocorrência da doença com estimativa de intervalo de confiança.
Pela distribuição binomial a MLE da taxa de ocorrência é:
$$ \hat{p} = \frac{37}{100} = 0.37$$.
e o erro padrão dessa estimativa é:
$$ \hat{\sigma} = [ \frac{p (1 - p)}{n} ]^{1/2} = [ \frac{ 0.37 (1 - 0.37) }{100} ]^{1/2} = 0.04828043$$.
Utilizando a aproximação normal para grandes amostras ($n \geq 30$) sabemos que a estatística
$$ \hat{z} = \frac{ \hat{p} - p }{ \hat{\sigma}_p } $$
tem distribuição amostral igual à distribuição Normal padronizada (média zero e desvio padrão um).
Assim, um intervalo com probabilidade 95% para essa estatística é:
$$ P( z_{0.025} \leq \hat{z} \leq z_{0.975} ) = 0.95 $$
$$ P( -1.96 \leq \hat{z} \leq 1.96 ) = 0.95$$
$$ P( -1.96 \leq \frac{ \hat{p} - p }{\hat{\sigma}} \leq 1.96 ) = 0.95 $$
$$ P( \hat{p} -1.96 \hat{\sigma} \leq p \leq \hat{p} + 1.96 \hat{\sigma} ) = 0.95 $$
Assim o intervalo de confiança de 95% para estimativa da taxa de ocorrência de doença $\hat{p}$ é:
$$ \hat{p} \pm 1.96 \sigma = 0.37 \pm (1.96) (0.04828043) = 0.37 \pm 0.09462964 $$.
Intervalo de Verossimilhança
O intervalo de verossimilhança (para razão 8, por exemplo) é obtido inspecionando a vizinhança da MLE $\hat{p}$ na curva de verossimilhança:
p = seq(0.20, 0.50, length=100) lik = dbinom(37, 100, p) lik = lik / max(lik) plot(p, lik, type="l", col="red") abline(h=1/8, lty=2, col="blue") abline(v=37/100, lty=9, col="red")
Suponha agora que a amostra aleatória de árvores de 100 árvores foi obtida, mas nenhuma das árvores se mostrou doente. Isso indica que a ocorrência de doença é rara e, consequentemente, a taxa de infestação é muito pequena, próxima de zero.
Estimativa da taxa: $\hat{p} = 0 / 100 = 0$
Erro padrão da estimativa: $\hat{\sigma} = [ (0 (1 - 0))/ 100 ]^{1/2} = 0$
Como utilizar a aproximação normal nesse caso? Não é possível obter um intervalo de confiança de 95% por essa abordagem. Uma nova abordagem, com outra estatística e outra distribuição amostral, deverá ser concebida.
O que muda no intervalo de verossimilhança? Absolutamente nada:
p = seq(0.0, 0.05, length=100) lik = dbinom(0, 100, p) lik = lik / max(lik) plot(p, lik, type="l", col="red") abline(h=1/8, lty=2, col="blue") abline(v=0, lty=9, col="red")
A forma de teste de hipótese de uso mais geral na estatística frequentista é o o teste de significância.
Essa abordagem consiste em enunciar duas hipóteses:
O teste de significância segue os seguintes passos:
Voltemos ao exemplo da aplicação de um produto em cobaias para verificar a taxa de mortalidade com os dois laboratórios:
A questão agora é testar as seguintes hipóteses:
Laboratório A: o modelo deste experimento é uma distribuição binomial. A probabilidade de obter seis ou menos mortes em 20 tentativas sob a hipótese de que $p=0.5$ é dada pela probabilidade acumulada da binomial:
pbinom(q=6, size=20, prob=0.5)
Laboratório B: o modelo do experimento é uma distribuição binomial negativa. A probabilidade de obter seis mortes em 20 ou mais tentativas é:
1 - pnbinom(q=14, size=6, prob=0.5)
Conclusão: Tomando-se o nivel de significância de 5% (0.05), o laboratório A não rejeitaria a hipótese nula, mas o laboratório B a rejeitaria.
Mesmo nível de significância e mesmos dados, mas conclusões diferentes.
Na inferência por verossimilhança, como vimos, os dois laboratório tem exatamente a mesma curva de verossimilhança e, portanto, a evidência estatística também:
p = seq(0.01, 0.99, by=0.01) lik.binom = dbinom(6, 20, p) # Lab 1: dist. Binomial lik.binom = lik.binom / max(lik.binom) lik.nbinom = dnbinom(14, 6, p) # Lab 2: dist. Binomial Negativa lik.nbinom = lik.nbinom / max(lik.nbinom) plot(p, lik.binom, type="l", ylab="Verossimilhança Relativa", col="blue", lwd=8) lines(p, lik.nbinom, col="red", lwd=2)
Entretanto, não seria apropriado estabelecer as hipóteses na forma
Pois a hipótese A indica um ponto na curva enquanto que a hipótese B indica uma região.