%~Mouliné par MaN_auto v.0.28.0 2023-07-31 17:17:47
\documentclass[CRMATH, Unicode, XML]{cedram}

\TopicFR{Statistiques}
\TopicEN{Statistics}

%\usepackage{pstool}
%\usepackage{mathtools}
\usepackage{subfigure}
\usepackage{array}
\usepackage[noadjust]{cite}

\DeclareMathOperator{\Ber}{Ber}
\DeclareMathOperator{\Kum}{Kum}

\newcolumntype{C}[1]{>{\centering\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}

%\newcommand{\mb}{\mbox}

\newcommand{\te}{\theta}
\newcommand{\G}{\Gamma}
\newcommand{\al}{\alpha}
\newcommand{\be}{\beta}

\newcommand*{\dd}{\mathrm{d}}
\newcommand*{\dz}{\dd z}
\newcommand*{\dte}{\dd \theta}

\newcommand{\E}{\mathbb{E}}
\newcommand{\V}{\mathbb{V}\mathrm{ar}}

\newcommand{\foxH}[7]{\,H_{#3\,\,\,#4}^{#1\,\,\,#2}\left[#5\middle|
\begin{array}{c}
#6 \\
#7 \\
\end{array}
\right]}

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}


\title{Exact Posterior distribution of risk ratio in the Kumaraswamy--Binomial model}

\author{\firstname{Jose} \middlename{A. A.} \lastname{Andrade}\IsCorresp}
\address{Department of Statistics and Applied Mathematics, Federal University of Ceara, 60455-670, Fortaleza-Ce, Brazil}
\email{ailton@ufc.br}

\author{\firstname{Pushpa} \lastname{N.} \lastname{Rathie}}
\address{Department of Statistics, University of Brasilia, 70910-900, Brasilia-DF, Brazil}
\email{pushpanrathie@yahoo.com}

\begin{abstract}
In categorical data analysis, the $2\times 2$ contingency tables are commonly used to assess the association between groups and responses, this is achieved by using some measures of association, such as the contingency coefficient, odds ratio, risk relative, etc. In a Bayesian approach, the risk ratio is modeled according to a Beta-Binomial model, which has exact posterior distribution, due to the conjugacy property of the model. In this work, we provide the exact posterior distribution of the relative risk for the non-conjugate Kumaraswamy--Binomial model. The results are based on special functions and we give exact expressions for the posterior density, moments, and cumulative distribution. An example illustrates the theory.
\end{abstract}

%~\keywords{Exact Bayesian computation, Contingency tables, Beta-Binomial Model, Special Functions}
\subjclass{62C10}

\begin{document}
\maketitle

\section{Introduction}

The $2\times 2$ contingency tables is one of the most prominent forms of displaying statistical data, while the tables are very intuitive and easy to interpret, important information also can be extracted from them. In fact, verifying whether an attribute is more frequent in a group rather than another is a question which we often want to investigate. For instance, the association between binary variables is an essential tool for decision-making in many areas, such as classification, clinical trials, cohort studies, etc. In those studies, several measures can be used, to name a few: \emph{the odds ratio}, \emph{contingency coefficient}, \emph{$\phi$ coefficient}, etc. One of the most widely used measure of association is the \emph{risk ratio} (or \emph{relative risk}). In Bayesian analysis, the groups in a $2\times 2$ table are modeled as independent Binomial distributions with unknown probabilities of success, and the corresponding parameters are modeled as Beta distributions with some fixed hyperparameters, conveniently chosen to represent our prior knowledge about the odds of success. Thus, those measures usually are estimated through a Beta-Binomial model, which is a \emph{conjugate} structure, in the sense that the posterior distribution belongs to the same family as the prior distribution. From this setting, the posterior distribution of the risk ratio can be obtained in exact form (Aitchison \& Bacon-Shone \cite{Aitchison91}). The Beta-Binomial model is widely used to obtain posterior estimates for some measures of association, Smith et al. \cite{Smith85} obtain some hypothesis test involving the probabilities of success, Aitkin \& Chadwick \cite{Aitkin03} provides estimates for several measures of association.

Another perspective may find the Beta distribution rather restrictive, since other prior distributions can better represent our knowledge about the parameters, however using other prior distributions may require approximate methods such as MCMC, numerical integration, etc., since there might not have a conjugate structure. This difficulty makes researchers to choose the facility of the Beta-Binomial model. Nevertheless, the well-established concept, that only conjugate prior distributions (or some specific reference priors) could yield exact posterior distributions was broken by Andrade \& Rathie \cite{Andrade15}, where it was showed that a great variety of nonconjugate models can still produce exact posterior distributions, in the sense that the posterior distribution and its quantities can be explicitly written in a computable form. In this new field, Andrade \& Rathie \cite{Andrade15} and Andrade et al.~\cite{Andrade18} propose very wide classes of exact posterior nonconjugate models for a scale and location parameters, respectively. A generalization for location-scale structures is provided by Andrade \& Rathie~\cite{Andrade17}. The main tool used in these works are the special functions, in particular the Fox's \emph{H-function}~\cite{Fox61}, which generalizes most of the distributions in Statistics.

In the context of Binomial models, Andrade \cite{Andrade20} proposes an alternative to the Beta-Binomial model, the \emph{Kumaraswamy-Binomial} model, in which the posterior quantities are also obtained in exact form. The theory provides exact expressions for the posterior moments, cumulative posterior and predictive distributions. Due to some friendly mathematical properties (e.g. cdf, moments, median, etc), the Kumaraswamy distribution (Kumaraswamy \cite{Kumaraswamy80}) has been widely used as an alternative to Beta distribution (see for instance, Andrade et al. \cite{Andrade18} and Fletcher \& Ponnambalam~\cite{Fletcher96}). In this work, we use the results of Andrade~\cite{Andrade20}, concerning the nonconjugate Kumaraswamy--Binomial model, in order to obtain the exact posterior distribution and its quantities of the risk ratio in $2\times 2$ contingency tables.

The text is organized as the following. Section~\ref{sec:hfunc} provides some preliminary definitions and results useful to the theory. In Section~\ref{sec:main} the main results are presented, in which the explicit posterior distribution and the posterior quantities are provided. In Section~\ref{sec:example} an illustrative example given, in which the exact results are compared with those obtained by the MCMC method. In Section~\ref{sec:conc} some general comments consider possible generalizations of the theory.


\section{Preliminaries}\label{sec:hfunc}

Andrade \cite{Andrade20} gives the exact posterior quantities for the Kumaraswamy--Binomial model. Firstly, the results are expressed through \emph{unormalized moments}, that is, if a random variable $X\sim f(x)$ such that $f(x)=c\times h(x)$, where $c$ is the normalizing constant and $h(x)$ is the kernel of the distribution, the unnormalized moment of order $r\in \mathbb{Z}$ is defined as $I(r)=\int_{\mathcal{X}} x^r h(x)\dd x$. Once we obtain an expression for $I(r)$ we can access most of the posterior quantities, such as the normalizing constant itself, the posterior moments, etc.

Consider the Kumaraswamy--Binomial model
\begin{equation}\label{eq:mod0}
\begin{cases}
y_1,\dots,y_n|\te\sim \Ber(\te)\,\,\text{iid} \\
\te\sim \Kum(\al,\be),
\end{cases}
\end{equation}
where $0<\te<1$ is the parameter of interest, $\al,\be$ are hyperparameters and $\Kum(\al,\be)$ stands for the Kumaraswamy distribution $p(\te)=\al\be \te^{\al-1}(1-\te^{\al})^{\be-1}$. Andrade \cite{Andrade20} gives an exact expression for the unnormalized posterior moment of order $r$, given by
\begin{equation}\label{eq:exp1}
I_y(r)=\al\be \G(\be)\G(n-y+1) \sum_{h=0}^\infty \frac{(-1)^h\G(r+y+\al(h+1))}{h!\G(r+n+\al(h+1)+1)\G(\be-h)},
\end{equation}
where $\al>1$ and $y=\sum_i y_i$.

An important result (Luke \cite{Luke79}) uses the fact that, in general, for $0<b<1$ and $c>0$, using the binomial theorem, we have
\begin{equation}\label{eq:r1}
\int_0^bz^{c-1}(1-z)^d \dz=\sum_{h=0}^{\infty} \frac{(-d)_h}{h!}\int_0^b z^{c+r-1}\dz=\sum_{h=0}^{\infty} \frac{(-d)_h}{h!}\frac{b^{c+h}}{c+h},
\end{equation}
where $(-d)_h$ is the Pochhammer notation, $(-d)_h=\frac{(-d)(-d-1)\dots,(-d-h+1)}{h!}=\frac{\G(h-d)}{\G(-d)}$.


\section{Posterior distribution of the risk ratio}\label{sec:main}

We consider an experiment where $n_1$ and $n_2$ subjects are randomly allocated into the groups 1 and 2, say treatment 1 and 2, the responses to the treatments are marked as ``Success/Failure'' (1/0), thus a general $2\times 2$ contingency table is of the form
\[
\begin{tabular}{l|cc}
\hline Group & \multicolumn{2}{c}{Response} \\ \cline{2-3} & 1 & 0 \\ \hline Treat. 1 & $\te_1$ & $1-\te_1$ \\
Treat. 2 & $\te_2$ & $1-\te_2$ \\
\hline
\end{tabular}
\]

The probabilities of success within the two groups are $\te_1$ and $\te_2$ ($0<\te_1,\te_2<1$). Thus, we let $P(X_i=1|\te_1)=\te_1$ ($i=1,\dots,n_1$) and $P(Y_j=1|\te_2)=\te_2$ ($j=1,\dots,n_2$). Note that $X_i$ and $X_j$ are independent, as well as $\te_1$ and $\te_2$, hence we consider the Binomial model
\begin{equation}\label{eq:mod1}
\begin{cases}
x_1,\dots,x_{n_1}|\te_1\sim \Ber(\te_1)\,\,\text{iid} \\
y_1,\dots,y_{n_2}|\te_2\sim \Ber(\te_2)\,\,\text{iid} \\
\te_j\sim \Kum(\al_j,\be_j),\,\,j=1,2.
\end{cases}
\end{equation}
where $\al_j,\be_j$ are hyperparameters for the Kumaraswamy distribution $p_j(\te_j)=\al_j\be_j \te^{\al_j-1}(1-\te^{\al_j})^{\be_j-1}$ ($j=1,2$). The kernel of the posterior distribution is given by
\begin{align}
p(\te_1,\te_2|x,y) & \propto L(\te_1,\te_2) p_1(\te_1) p_2(\te_2) \nonumber \\
& \propto \te_1^{x}(1-\te_1)^{n_1-x}p_1(\te_1)\times \te_2^{y}(1-\te_1)^{n_2-y} p_2(\te_2), \label{eq:mod3}
\end{align}
where $x=\sum_i x_i$, $y=\sum_i y_i$ and $L(\te_1,\te_2)$ is the likelihood function. Clearly, the marginal posterior distributions will not be in the Kumaraswamy family, that is, the likelihood and prior distribution will not combine their kernels, thus in some regular Bayesian view, some approximate method would be required.

In our approach, we first obtain the unnormalized moment of order $r$ and $s$, that is
\begin{align*}
I_{x,y}(r,s) & =\int_0^1\int_0^1\te_1^r\te_2^s p(\te_1,\te_2|x,y)\dte_1\dte_2 \\
& = \int_0^1\te_1^{x+r}(1-\te_1)^{n_1-x} p_1(\te_1)\dte_1\times\int_0^1 \te_2^{y+s}(1-\te_1)^{n_2-y} p_2(\te_2)\dte_2=:I_x(r)\times I_y(s), 
\end{align*}
where $I_x(r)$ and $I_y(s)$ are given by~\eqref{eq:exp1} and the computable representations~\eqref{eq:exp1}, replacing conveniently $\al$ and $\be$ by $\al_j, $ and $\be_j$ ($j=1,2$), respectively. The posterior normalizing constant is obtained by evaluating the series up to a required precision, the normalizing constant is $I_{x,y}(0,0)=I_x(0)\times I_y(0)$, hence the exact posterior distribution is obtained dividing~\eqref{eq:mod3} by $I_{x,y}(0,0)$, which can be written in terms of the \emph{binomial theorem},
\begin{multline}
p(\te_1,\te_2|x,y) =\frac{\al_1\be_1\al_2\be_2}{I_x(0) I_y(0)} \sum_{j=0}^{n_1-x}\sum_{\ell=0}^{n_2-y} \binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell} \te_1^{x+\al_1+j-1}\te_2^{y+\al_2+\ell-1} \\
 \times (1-\te_1^{\al_1})^{\be_1-1}(1-\te_2^{\al_2})^{\be_2-1}\label{eq:mod4}
\end{multline}

In order to make inferences about $\rho$, we need to obtain the posterior distribution of $\rho$, thus we consider the posterior distribution~\eqref{eq:mod4} to obtain the following exact posterior quantities.


\begin{theo}[Cumulative posterior distribution of the risk ratio] \label{cdfrr1}
Let $\rho=\te_1/\te_2$ be the risk ratio. Considering~\eqref{eq:mod4}, the cumulative posterior distribution of $\rho$ is given by
\begin{enumerate}\romanenumi
\item \label{1i} For $0<\rho_0<1$:
\begin{multline}
P(\rho<\rho_0|x,y) =\frac{\be_1\be_2}{I_x(0)I_y(0)}\frac{\G(\be_2)}{\G(1-\be_1)}\sum_{j=0}^{n_1-x}\sum_{\ell=0}^{n_2-y}\binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell}\rho_0^{x+j+\al_1}\\
\times\sum_{h=0}^{\infty}\frac{\rho_0^{\al_1 h}\G(1-\be_1+h)\G\left[\frac{y+\ell+x+j+\al_1+\al_2}{\al_2}+\frac{\al_1h}{\al_2}\right]}{h!\left[\frac{x+j+\al_1}{\al_1}+h\right]\G\left[\frac{y+\ell+x+j+\al_1+\al_2+\al_2\be_2}{\al_2}+\frac{\al_1h}{\al_2}\right]}. \label{eq:posrr1}
\end{multline}

\item \label{1ii} For $\rho_0>1$:
\begin{multline}
P(\rho<\rho_0|x,y) =1-\frac{\al_2\be_1\be_2\G(\be_1)}{I_x(0)I_y(0)\G(1-\be_2)}\sum_{j=0}^{n_1-x}\sum_{\ell=0}^{n_2-y}\binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell}\rho_0^{-(y+\ell+\al_2)}\\
\times\sum_{h=0}^{\infty}\frac{\rho_0^{-\al_2 h}\G(1-\be_2+h)\G[\frac{y+\ell+x+j+\al_1+\al_2+\al_1 h}{\al_1}]}{h!(y+\ell+\al_2+\al_2h)\G\left[\frac{y+\ell+x+j+\al_1+\al_2+\al_2 h+\be_1\al_1}{\al_1}\right]}. \label{eq:posrr2}
\end{multline}
\end{enumerate}
\end{theo}

\begin{proof}

\begin{proof}[(\ref{1i})]
For $0<\rho_0<1$: Let $\rho=\te_1/\te_2$ and $\te_2=\te_2$, using~\eqref{eq:mod4},
\begin{multline}
\begin{aligned}
P(\rho<\rho_0|x,y) & = P(\te_1<\rho_0\te_2|x,y) = \int_0^1\int_0^{\rho_0 \te_2} p(\te_1,\te_2|x,y)\dte_1\dte_2 \\
&=\frac{\al_1\al_2\be_1\be_2}{I_x(0)I_y(0)}\sum_{j=0}^{n_1-x}\sum_{\ell=0}^{n_2-y} \binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell}\times
\end{aligned} \\
\times \int_0^1 \te_2^{y+\al_2+\ell-1}(1-\te_2^{\al_2})^{\be_2-1}\int_0^{\rho_0\te_2} \te_1^{x+\al_1+j-1}(1-\te_1^{\al_1})^{\be_1-1}\dte_1 \dte_2. \label{eq:post001}
\end{multline}

Let $\te=\te_1^{\al_1}$, then
\begin{align*}
A(\te_2) & = \int_0^{\rho_0\te_2} \te_1^{x+\al_1+j-1}(1-\te_1^{\al_1})^{\be_1-1}\dte_1 \\
& = \frac{1}{\al_1}\int_0^{(\rho_0\te_2)^\al_1} \te^{\frac{x+j}{\al_1}}(1-\te)^{\be_1-1}\dte \\
& = \frac{1}{\al_1}\sum_{h=0}^{\infty} \frac{(1-\be_1)_h (\rho_0\te_2)^{\al_1(\frac{x+j}{\al_1}+h+1)}}{h!(\frac{x+j}{\al_1}+h+1)},
\end{align*}
where the latest identity is due to~\eqref{eq:r1}.

Thus, replacing $A(\te_2)$ in~\eqref{eq:post001}, we have
\begin{multline*}
P(\rho<\rho_0|x,y) = \frac{\al_1\al_2\be_1\be_2}{I_x(0)I_y(0)}\sum_{j=0}^{n_1-x}\sum_{\ell=0}^{n_2-y} \binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell}\times \nonumber \\
\times \sum_{h=0}^{\infty} \frac{(1-\be_1)_h}{h!}\frac{\rho_0^{x+j+\al_1+\al_1 h}}{x+j+\al_1+\al_1 h}\int_0^1 \te_2^{y+\al_2+\ell+x+j+\al_1+\al_1h-1}(1-\te_2^{\al_2})^{\be_2-1}\dte_2
\end{multline*}
The result follows by evaluating the integral on $\te_2$ as before.
\let\qed\relax
\end{proof}

\begin{proof}[(\ref{1ii})]
For $\rho_0>1$: Note that $P(\rho<\rho_0|x,y) = 1- P(\rho>1/\rho_1|x,y)$, where $\rho_1=1/\rho_0$, hence $0<\rho_1<1$, then we consider
\begin{multline}
\begin{aligned}
P(\rho>1/\rho_1|x,y) & = P(\te_2<\rho_1\te_1|x,y) = \int_0^1\int_0^{\rho_1 \te_1} p(\te_1,\te_2|x,y)\dte_2\dte_1 \\
&=\frac{\al_1\al_2\be_1\be_2}{I_x(0)I_y(0)}\sum_{j=0}^{n_1-x}\sum_{\ell=0}^{n_2-y} \binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell}
\end{aligned}\\
\times \int_0^1 \te_1^{x+\al_1+j-1}(1-\te_1^{\al_1})^{\be_1-1}\int_0^{\rho_1\te_1} \te_2^{y+\al_2+\ell-1}(1-\te_2^{\al_2})^{\be_2-1}\dte_2 \dte_1. \label{eq:post002}
\end{multline}
The result follows by the same strategy as before, that is, replacing $\rho_1=1/\rho_0$ and using $P(\rho<\rho_0|x,y) = 1- P(\rho>\rho_0|x,y)$.
\end{proof}
\let\qed\relax
\end{proof}

\begin{coro}[Posterior density of the risk ratio] \label{rrpos1}
Considering~\eqref{cdfrr1}, the posterior distribution of the risk ratio is given by
\begin{enumerate}\romanenumi
\item \label{2i} For $0<\rho_0<1$:
\begin{multline}
p(\rho|x,y) =\frac{\al_1\be_1\be_2}{I_x(0)I_y(0)}\sum_{j=0}^{n_1-x}\sum_{\ell=0}^{n_2-y}\binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell}\frac{\G(\be_2)}{\G(1-\be_1)}\\
\times\sum_{h=0}^{\infty} \frac{\rho^{x+j+\al_1+\al_1h-1}\G(1-\be_1+h)\G\left[\frac{y+\ell+x+j+\al_1+\al_2}{\al_2}+\frac{\al_1h}{\al_2}\right]}{h!\G\left[\frac{y+\ell+x+j+\al_1+\al_2+\al_2\be_2}{\al_2}+\frac{\al_1h}{\al_2}\right]}. \label{eq:posrr11b}
\end{multline}

\item \label{2ii} For $\rho_0>1$:
\begin{multline}
P(\rho|x,y) =\frac{\al_2\be_1\be_2\G(\be_1)}{I_x(0)I_y(0)\G(1-\be_2)}\sum_{j=0}^{n_1-x}\sum_{\ell=0}^{n_2-y}\binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell}\\
\times\sum_{h=0}^\infty \frac{\rho^{-(y+\ell+\al_2+\al_2 h+1)} \G(1-\be_2+h)\G[\frac{y+\ell+x+j+\al_2+\al_1+\al_2 h}{\al_1}]}{h!\G\left[\frac{y+\ell+x+j+\al_2+\al_2 h+\al_1+\be_1\al_1}{\al_1}\right]}.\label{eq:posrr22b}
\end{multline}
\end{enumerate}
\end{coro}

\begin{proof}
Straightforward by differentiating the~\eqref{eq:posrr1} and~\eqref{eq:posrr2} with respect to $\rho_0$.
\end{proof}


\begin{coro}[Posterior moments of the risk ratio] \label{mrr2}
Considering Corollary~\eqref{rrpos1}, the posterior moment of order $r\in \mathbb{Z}$ of $\rho$ is given by
\begin{equation}\label{eq:cdfrr1}
\E[\rho^r|x,y] = \sum_{j=0}^{n_1-x} \sum_{\ell=0}^{n_2-y}\sum_{h=0}^{\infty}\left[\frac{A_{j\ell(h)}}{x+j+\al_1 h+r+\al_1}+\frac{B_{j\ell}(h)}{y+\ell+\al_2r-r+\al_2}\right],
\end{equation}
where
\begin{align}
A_{j\ell} & =\frac{\al_1\be_1\be_2\G(\be_2)}{I_x(0)I_y(0)\G(1-\be_1)}\binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell} \frac{\G(1-\be_1+h)\G\left[\frac{y+\ell+x+j+\al_1+\al_2}{\al_2}+\frac{\al_1h}{\al_2}\right]}{h!\G\left[\frac{y+\ell+x+j+\al_1+\al_2+\al_2\be_2}{\al_2}+\frac{\al_1h}{\al_2}\right]} \label{eq:posrr11c}
\end{align}
and
\begin{align}
B_{j\ell}(h) & =\frac{\al_2\be_1\be_2\G(\be_1)}{I_x(0)I_y(0)\G(1-\be_2)}\binom{n_1-x}{j}\binom{n_2-y}{\ell}(-1)^{j+\ell} \frac{ \G(1-\be_2+h)\G[\frac{y+\ell+x+j+\al_2+\al_1+\al_2 h}{\al_1}]}{h!\G\left[\frac{y+\ell+x+j+\al_2+\al_2 h+\al_1+\be_1\al_1}{\al_1}\right]}.\label{eq:posrr22c}
\end{align}
\end{coro}

\begin{proof}
Straightforward by integrating
\begin{multline*}
\E[\rho^r|x,y] = \int_0^1 \sum_{j=0}^{n_1-x} \sum_{\ell=0}^{n_2-y}\sum_{h=0}^{\infty} A_{j\ell}(h)\rho^{r+x+j+\al_1+\al_1h-1} \dd\rho\\
 +\int_1^\infty \sum_{j=0}^{n_1-x} \sum_{\ell=0}^{n_2-y}\sum_{h=0}^{\infty} B_{j\ell}(h) \rho^{r-(y+\ell+\al_2+\al_2 h+1)} \dd\rho
\end{multline*}
\end{proof}


\section{Example}\label{sec:example}

In this section we illustrate the theory through an illustrative example, consider the following $2\times 2$ contingency table, considering either the data or the probabilities, that is
\[
\begin{tabular}{l|cc}
\hline Group & \multicolumn{2}{c}{Response} \\ \cline{2-3} & Yes & No \\ \hline Treat. 1 & 2 & 18 \\
Treat. 2 & 3 & 8 \\
\hline
\end{tabular}\qquad\Longrightarrow\qquad
\begin{tabular}{l|cc}
\hline Group & \multicolumn{2}{c}{Response} \\ \cline{2-3} & Yes & No \\ \hline Treat. 1 & $\te_1$ & $1-\te_1$ \\
Treat. 2 & $\te_2$ & $1-\te_2$ \\
\hline
\end{tabular}
\]
Here $\te_1=P(\mathrm{Response}=\mathrm{Yes}|\mathrm{Treat. 1})$ and $\te_2=P(\mathrm{Response}=\mathrm{Yes}|\mathrm{Treat. 2})$ are the parameters that will be used to estimate the risk ratio. We apply Model~\eqref{eq:mod1}, noting that $n_1=20$, $x=2$, $n_2=11$ and $y=3$.

It follows that the unnormalized marginal posterior moment of order $r$ of $\te_i$ ($i=1,2$) is given~by
\begin{equation}\label{eq:mome1}
I_{D_i}(r)=\al\be \G(\be)\G(n_i-D_i+1) \foxH{1}{1}{2}{2}{1}{(1-(r+D_i+\al),\al),(\be,1)}{(0,1),(-(r+n_i+\al),\al)},
\end{equation}
and the corresponding computable expression is given by
\begin{equation}\label{eq:expe1}
I_{D_i}(r)=\al\be \G(\be)\G(n-D_i+1) \sum_{h=0}^\infty \frac{(-1)^h\G(r+D_i+\al(h+1))}{h!\G(r+n+\al(h+1)+1)\G(\be-h)},
\end{equation}
where $D_1=x$, $D_2=y$ and $\al>1$ and $\be>0$ are the hyperparameters of the Kumaraswamy distribution. We define the same prior distribution $\te_i\sim \mathrm{Kum}(\al=1.1,\be=1.1)$ ($i=1,2$), which are flat (nearly uniform) prior distributions.

As the results, we provide comparisons of the posterior quantities obtained from MCMC algorithm ({\tt OpenBugs} package with 10,000 cycles and 8,001 burned in) with the exact form equivalents. In Table~\ref{tab:results}, one can notice that the posterior quantities of the posterior expectation, variance, cumulative distribution are similar in the two methods. Notice that the cumulative posterior distribution can be used to obtain credible intervals. According to~\eqref{rrpos1}, the posterior density of $\rho$ assumes different expression for $0<\rho\leq 1$ and $\rho>1$, thus we separate the MCMC simulated values according to these intervals. Note that the comparison between exact and approximated densities in Figure~\ref{fig:fig1} shows some similarities, however our result is exact.

\begin{table}[ht]
\caption{Comparison of the exact posterior estimates with their equivalent obtained by MCMC (10,000, discarding 8,000)}\label{tab:results}
\begin{tabular}{l|C{2cm}|C{2cm}}
\hline\hline Posterior & \multicolumn{2}{c}{Estimates} \\ \cline{2-3} Quantity & MCMC & Exact \\ \hline $\E[\rho|x,y]$ & 0.5526 & 0.4881 \\
$\V[\rho|x,y]$ & 0.1526 & 0.0739 \\
$P(\rho\leq 0.8|x,y)$ & 0.8120 & 0.8167 \\
\hline\hline
\end{tabular}
\end{table}

\begin{figure}[tp]
\subfigure[For values of $0<\rho\leq 1$]{
\centering\label{fig:fig1a}
%\psfrag{rho}[c][l][1.2][0]{$\rho$}
%\psfrag{MCMC}[c][l][1.2][0]{\hspace{0.8cm} MCMC}
%\psfrag{Exact}[c][l][1.2][0]{\hspace{0.8cm} Exact}
\includegraphics[width=0.49\textwidth]{conttable1.pdf}}
\subfigure[For values of $\rho> 1$]{
\centering\label{fig:fig1b}
%\psfrag{rho}[c][l][1.2][0]{$\rho$}
%\psfrag{MCMC}[c][l][1.2][0]{\hspace{0.8cm} MCMC}
%\psfrag{Exact}[c][l][1.2][0]{\hspace{0.8cm} Exact}
\includegraphics[width=0.49\textwidth]{conttable2}}
\caption{Exact and approximated posterior density of the risk ratio $\rho$.}\label{fig:fig1}
\end{figure}

\section{Concluding remarks}\label{sec:conc}

In categorical data analysis there are many measures of association to handle $2\times 2$ tables, the most prominent of them are the odds ratio and risk ratio. In Bayesian approaches, the risk ratio is commonly modeled through Beta-Binomial model, which is rather straightforward to apply, due to the conjugancy of the model. In this work, we establish results based on special functions that bring the same facility of the conjugate models, but with the Kumaraswamy distribution.

Future work will search for some exact models involving the odds ratio, which will increase the complexity of the model. Note that, once we have set the Binomial and prior distributions for the proportion of favorable answers, we will have to obtain a posterior distribution for the odds ratio, we will need to use special functions to assess the distribution of the odds ratio, this may lead to some very complex distribution.

\bibliographystyle{crplain}
\bibliography{CRMATH_Andrade_20221108}
\end{document}