Introdução

Suponha que \(\theta\) é uma proporção e você deseja construir uma distribuição a priori baseada em estudos anteriores. Como, em geral, os dados de outros estudos não estão disponíveis na íntegra, você pode utilizar estatísticas resumo.

Paper n Cases p CI
Estudo 1 78 30 0.385 [ 0.277 ; 0.493 ]
Estudo 2 60 42 0.700 [ 0.584 ; 0.816 ]
Estudo 3 40 30 0.750 [ 0.616 ; 0.884 ]

Proposta 1

Uma primeira forma de fazer isso seria supor que, antes de observar esses estudos, você não tem informação sobre \(\theta\), e usar a posteriori obtida com os dados do estudos como priori. Assim, seja \(X_i|\theta \sim \text{Bin}(n_i,\theta)\), \(i=1,2,3\), e \(\theta \sim \text{Uniforme}(0,1)\). Temos que \[\theta | x_1 \sim \text{Beta}(1+x_1,1+n_1-x_1)\] \[\theta | x_1,x_2 \sim \text{Beta}(1+x_1+x_2,1+n_1-x_1+n_2-x_2)\] \[\theta | x_1,x_2,x_3 \sim \text{Beta}(1+x_1+x_2+x_3,1+n_1-x_1+n_2-x_2+n3-x_3)\]

Isso é equivalente a juntar todos os estudos e calcular diretamente a posteriori. Seja \(x=x_1+x_2+x_3\) e \(n=n_1+n_2+n_3\). Então, \[\theta \sim \text{Beta}(1,1) \implies \theta | x_1,x_2,x_3 \sim \text{Beta}(1+x,1+n-x)~.\] \(~\)

De fato, se \(X\) e \(Y\) são, condicionalmente independentes dado \(\theta\),

\(f(\theta|x)\) \(\propto f(x|\theta)f(\theta)\)

\(f(\theta|x,y)\) \(\propto f(y|\theta)~f(\theta|x)\) \(\propto f(y|\theta)f(x|\theta)f(\theta)\) \(\propto f(x,y|\theta)f(\theta)~.\)

\(~\)

eixox=c(min(estudos$inf),max(estudos$sup))

PriorsPlot=tibble(theta=seq(eixox[1],eixox[2],length.out=1000),
              beta=dbeta(theta,sum(estudos$Cases)+1,
                            sum(estudos$n)-sum(estudos$Cases)+1)) %>% 
  ggplot() + theme_bw() +
  geom_line(aes(x=theta,y=beta, colour="Beta")) +
  labs(colour="Priori") + xlim(eixox) + 
  xlab("theta") + ylab("Prior")

Meta=estudos %>%
  ggplot(aes(y=Paper))+theme_bw()+
  geom_point(aes(x=p))+
  geom_segment(aes(x=inf,xend=sup,y=Paper,yend=Paper))+
  theme(axis.title.x=element_blank(),axis.title.y=element_blank())+
  xlim(eixox)+geom_vline(xintercept=sum(estudos$Cases)/sum(estudos$n), color='darkgrey', linetype='dashed')

ggpubr::ggarrange(PriorsPlot,Meta,heights=c(2,1),
          ncol = 1, align = "v",common.legend=T,legend="bottom")

Proposta 2

Considere agora que a priori será construída da seguinte forma: para cada estudo será calculada uma “posteriori” supondo que \(\theta \sim \text{Uniforme}(0,1)\) e faremos uma mistura dessas posterioris ponderada pelo tamanho amostral dos estudos. Assim:

\[f(\theta) = \sum_{i=1}^{3} \frac{n_i}{n}~f(\theta~|~a_i=1+x_i~,~b_i=1+n_i-x_i)~,\] em que \(f(\theta~|~a~,~b)\) é a densidade da \(\text{Beta}(a,b)\) e \(n=n_1+n_2+n_3\). No exemplo,

\[f(\theta) = 0.44~f(\theta~|~a_1=31,b_1=49)+0.34~f(\theta~|~a_2=43,b_2=19)+0.22~f(\theta~|~a_3=31,b_3=11)~.\]

As funções de densidade, distribuição e para gerar números aleatórios de misturas de densidades betas podem ser escritas no R como:

dmixbeta=function(theta,w,a,b){
  w = w/sum(w)
  apply(as.matrix(theta),1,function(t){t(w)%*%dbeta(t,a,b)})
}
pmixbeta=function(theta,w,a,b){
  w = w/sum(w)
  apply(as.matrix(theta),1,function(t){t(w)%*%pbeta(t,a,b)})
}
rmixbeta=function(n,w,a,b){
  w = w/sum(w)
  s = rmultinom(n=n,size=1,prob=w)
  apply(t(s),1,function(l){rbeta(1,l%*%a,l%*%b)})
}

Assim, a nova distribuição a priori é apresentada do gráfico a seguir.

Posteriori

Como vimos em aulas anteriores, como a priori é mistura de distribuições conjugadas, a posteriori também será. Suponha então que foi observada uma amostra de tamanho \(n=100\) com \(x=67\) sucessos. A posteriori é

\[f(\theta|x) = 0.0007~f(\theta~|~a_1=98,b_1=82)+0.7007~f(\theta~|~a_2=110,b_2=52)+0.2986~f(\theta~|~a_3=98,b_3=44)~.\]