The contour probability for $H0:θ0=0.25$ is obtained with
two computational approaches solving the same problem: measuring the
area under the beta posterior density curve with limits defined such
that $p(θ|y) > p(θ0|y)$. When θ0 is less than the
mode of the posterior density, this area corresponds to $P(θ0
< θ < δ)$, with density $p(δ | y)=p(θ0)$. In this
scenario, the complement $pb$ of the probability $P(θ0 <
θ < δ)$ can be interpreted as the counterpart of a two-sided
p-value. This is the contour probability.
In the first approach, we simply look for $δ$ such that
$(P(θ_0|y) - P(δ|y))^2=0$. This is accomplished with the
numerical optimization R program optimize. Again assuming
$θ0$ to be the lower limit of integration, the complement of
$pb$ is then obtained from:
\begin{align}
\begin{split}
P(θ0 < θ < δ) =& ∫θ_{0}δ
P(θ | y) dθ
=&∫θ_{0}δ\frac{1}
{\operatorname{Beta}(\bar{α}, \bar{β})} θ\bar{α-1}
(1 - θ)\bar{β - 1} \
=& F(δ | y) - F(θ0 | y)
\end{split}
\end{align}
with the R digital table for the beta cumulative distribution
$F$. In the second approach, we first obtain the mode $\bar{θ}$
of the posterior beta density with
$\frac{\bar{α}-1}{\bar{α} + \bar{β}-2}$. Then, for the
case where $θ0$ is less than $\bar{θ}$, we search for
$δ$ in the interval $[\bar{θ}, 1]$ such that $P(θ0 |
y)=P(δ | y)$. This is accomplished with the R program
uniroot. Once $δ$ is obtained, the complement of the contour
probability $pb$ is measured by following Equation eq_q3int. Both
approaches give return the same contour probabilities.
As in Questions 2 and 3, we give dose 0 an informative conjugate prior
with $α=y-1=29$ and $β=n-y-1=193$. The other doses are given
vague beta priors with $α=1$ and $β=1$ which correspond to a
uniform distribution with support from 0 to 1. In Table tab_contours,
the contour probabilities are presented for each dose. The probability
at dose 500 couldn’t converge because 1 is the maximum of the beta
distribution and the posterior mode as a value of 1. Convergence could
have be obtained with an informative prior. For every dose, the
probability of the true $θ$ being 0.25 is less than 0.05. Thus,
we conclude that the $H_0$ is implausible and reject it. Remarkably
with dose 0, we would have kept $H0$ had the prior been
non-informative (i.e. $p_b=0.68$).
| Dose |
$θ_0$ |
$δ$ |
p_b |
1 - p_b |
Obj. |
| 0 |
0.25 |
0.137 |
0.001 |
0.999 |
0 |
| 62.5 |
0.25 |
0.140 |
0.035 |
0.965 |
0 |
| 125 |
0.25 |
0.947 |
0.000 |
1.000 |
0 |
| 250 |
0.25 |
1.000 |
0.000 |
1.000 |
0 |
| 500 |
0.25 |
1.000 |
0.000 |
1.000 |
NA |
In Figure fig_contourplots, the complement of the contour probabilities are
illustrated for each dosage. These are coloured in grey and the
plotted curves are the posterior beta distributions.
Our goal here is to find the dosage associated with an increase of
0.01 in excess risk $q$ of foetus malformation. It is useful to
reformulate the question in order to clarify the problem and its
solution.
The excess risk $q$ is associated with a probability of malformation
$P(d*)$ by:
\begin{equation}
P(d*)=q(1-P(0)) + P(0)
\end{equation}
where $P(0)$ is the prevalence of foetus malformation at a dose 0 of
diglyme. The logit of $P(d)$, the probability of malformation at a
given dose, can be estimated with the linear model $α + β*
\mathrm{dose}$. As such, the dose we are looking for can be deduced
with simple algebra from
\begin{equation}
\mathrm{dose}=\mathrm{BMD}=\frac{\operatorname{logit}(P(d*))
\end{equation}
provided we have estimates for $α$ and $β$. This dose is
called the benchmark dose (BMD). From Question 3, we posess samples
from the posterior distributions of $β$ and $α$. Thus, we can
not only furnish a point estimate of the BMD, but also express its
uncertainty with a certitude interval. Only, as the dose variable is
standardized, Equation eq_q10dose becomes
\begin{equation}
\mathrm{BMD}=\left(\frac{\operatorname{logit}(P(d*)) - α}
{β}\right) s\mathrm{dose} + \operatorname{avg}(\mathrm{dose})
\end{equation}
where $s$ and $\operatorname{avg}$ respectively designate the standard
deviation and mean of the dose covariate.
Assuming that the prevalence at dose 0, $P(0)$, is known and obtained
from the proportion in the provided table (i.e. 0.238), we get a mean
BMD of 32.17 with a $95\%$ equal-tails certitude interval of 20.53 to
42.79.
