relative risk simulation
relative risk simulation
10 September 2018
Illustration of simulation to evaluate the posterior for the relative risk, given data from two Binomial experiments. \[
\begin{aligned}
Y_{A} |N_{A}, \theta_{A} & \sim \mathrm{Binomial}(\theta_{A},10) \\
Y_{B} |N_{B}, \theta_{B} & \sim \mathrm{Binomial}(\theta_{B},10)
\end{aligned}
\] \[
\begin{aligned}
p(\theta_{A},\theta_{B}) & = p(\theta_{A})p(\theta_{B}) \\
\theta_{A} & \sim ~ \mathrm{Beta}(1,1) \\
\theta_{B} & \sim \mathrm{Beta}(1,1)
\end{aligned}
\] We are interested in the relative risk, \(\mathrm{rr}= \theta_{A}/\theta_{B}\)
###set-up the data
YA <- 7 #successes
Na <- 10
YB <- 3
Nb <- 10
Under our model assumptions \[
\begin{aligned}
p(\theta_{A},\theta_{B}|\mathrm{data}) & = p(\theta_{A}|Y_{A},N_{A})p(\theta_{B}|Y_{B},N_{B}) \\& = \mathrm{Beta}(1+Y_{A},N_{A}-Y_{A}+1)\mathrm{Beta}(1+Y_{B},N_{B}-Y_{B}+1)
\end{aligned}
\] Monte Carlo evaluation for relative risk is therefore
nsim <- 10000 ## fix simulation sample size
postA <- rbeta(n=nsim,shape1=1+YA,shape2=Na-YA+1)
postB <- rbeta(n=nsim,shape1=1+YB,shape2=Nb-YB+1)
str(postA)
## num [1:10000] 0.655 0.713 0.604 0.742 0.559 ...
str(postB)
## num [1:10000] 0.216 0.359 0.496 0.125 0.303 ...
RR <- postA / postB
plot(density(RR), main="posterior density for RR")
summary(RR)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3209 1.5439 2.0574 2.4521 2.8899 25.7447
##Get useful summary quantiles
qRR <- quantile(RR,probs=c(0.025,0.5,0.975))
qRR
## 2.5% 50% 97.5%
## 0.9070335 2.0573733 6.3425039