#  Exact confidence intervals for functions of parameters in the $k$-sample
multinomial problem

Michael C Sachs∗1, Erin E Gabriel1, Michael P Fay2

( ∗<EMAIL_ADDRESS>
1Section of Biostatistics, University of Copenhagen, Denmark
2Biostatistics Branch, National Institutes of Allergy and Infectious Disease,
National Institutes of Health, USA
)

###### Abstract

When the target of inference is a real-valued function of probability
parameters in the $k$-sample multinomial problem, variance estimation may be
challenging. In small samples, methods like the nonparametric bootstrap or
delta method may perform poorly. We propose a novel general method in this
setting for computing exact p-values and confidence intervals which means that
type I error rates are correctly bounded and confidence intervals have at
least nominal coverage at all sample sizes. Our method is applicable to any
real-valued function of multinomial probabilities, accommodating an arbitrary
number of samples with varying category counts. We describe the method and
provide an implementation of it in R, with some computational optimization to
ensure broad applicability. Simulations demonstrate our method’s ability to
maintain correct coverage rates in settings where the nonparametric bootstrap
fails.
Keywords: exact inference, multinomial, computation.
Word count: 2497

##  1 Introduction

The target of inference in applied settings is sometimes related to a function
of the probability parameters from one or more samples from multinomial random
variables. If the function is complicated, or its variance is otherwise
difficult to calculate, then one may consider using the nonparametric
bootstrap for inference (Efron, 1979). However, the bootstrap may perform
poorly in small samples, e.g., not having the correct coverage or type I error
rate (Bickel and Freedman, 1981) or in large sample sizes when the true
parameter is on the boundary of the range (Bickel and Sakov, 2008). Instead,
we propose a general method in this setting for computing valid p-values and
confidence intervals, that is, tests based on the p-values correctly bound the
type I error rate and confidence intervals that ensure at least nominal
coverage. We present an implementation of the method in R (R Core Team, 2015).
Our approach is general in that it applies to any real-valued function of
multinomial probabilities from an arbitrary number of samples and with
different numbers of categories.

The main focus of previous work in exact methods for multinomial settings has
been testing for independence. Agresti (1992) provides a review of exact
methods for contingency tables from more than 30 years ago. Resin (2023) is a
recent work focusing on fast algorithms to compute exact goodness of fit
tests. Chafai and Concordet (2009) and Malloy et al. (2021) provide exact
confidence regions for a single multinomial parameter vector. Frey (2009)
describes an exact equivalence testing method for one sample multinomial
settings, which is a problem that fits into our general framework. It would be
nice to compare to our method, but no implementation of Frey (2009) could be
found.

In Section 2 we describe the problem set up and the statistical methods for
calculation of p-values and confidence intervals. In Section 3, we describe
some details of the implementation of the statistical methods, including our
approaches for computational optimization, so that the method is more broadly
applicable. We further perform some experiments to assess the accuracy of the
computational approach. In Section 4 we illustrate the method in two specific
settings and assess the performance of the method in comparison to the
nonparametric bootstrap via simulation studies. In Section 5, we describe the
main usage of our implementation of the method in the R package called
xactonomial. Finally, we discuss future directions in Section 6.

##  2 Methods

###  2.1 Setting and notation

We consider the following problem. Let $T_{j}$ be distributed
$\mbox{Multinomial}(n_{j},d_{j},\theta_{j})$ for $j=1,\ldots,k$, where $n_{j}$
is the number of trials, $d_{j}$ the number of mutually exclusive categories,
and $\theta_{j}$ the vector of $d_{j}$ probabilities. Denote
$\boldsymbol{T}=(T_{1},\ldots,T_{k})$ and
$\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{k})$. Suppose one is
interested in the parameter $\psi\equiv\psi(\boldsymbol{\theta})\in\Psi$ where
$\Psi$ is a subset of the real line $\mathbb{R}$. Given a sample of size
$n=\sum_{j=1}^{k}n_{j}$ from $\boldsymbol{T}$, one can estimate
$\boldsymbol{\theta}$ with the sample proportions as
$\hat{\boldsymbol{\theta}}(\boldsymbol{T})$ and then estimate the parameter of
interest as $\psi(\hat{\boldsymbol{\theta}}(\boldsymbol{T}))$. A simple
example for $k=1$ is $\psi(\boldsymbol{\theta})=\theta_{1}$, i.e., the first
element of the probability vector. An example for $k=2$ and $d_{1}=d_{2}=d$ is
the Bhattacharyya coefficient between $T_{1}$ and $T_{2}$ (Bhattacharyya,
1946). In that case we have
$\psi(\boldsymbol{\theta})=\sum_{i=1}^{d}\sqrt{\theta_{1i}\theta_{2i}}$.
Another example is given in Frey (2009), where for $k=1$, and for a fixed
vector of probabilities $\theta_{0}$, it is of interest to assess the distance
between $\theta_{0}$ and the observed data. The distance can be defined in
different ways, but one example is the Euclidean metric such that
$\psi(\boldsymbol{\theta})=\sqrt{\sum_{i=1}^{d}(\theta_{i}-\theta_{0i})^{2}}$.

It is naturally of interest to calculate a p-value for a test of the null
hypothesis $H_{0}:\psi(\boldsymbol{\theta})\leq\psi_{0}$ and construct a
$100(1-\alpha)\%$ confidence interval for $\psi(\boldsymbol{\theta})$.

###  2.2 Computation of p-values

We want to test the null hypothesis, $H_{0}:\psi\leq\psi_{0}$ versus the
alternative hypothesis, $H_{1}:\psi>\psi_{0}.$ Denote the joint probability
mass function
$\textrm{Pr}\\{\boldsymbol{T}=t|\boldsymbol{\theta}\\}=f(t,\boldsymbol{\theta})$.
Given our observed vectors of counts in each sample ${T}_{1},\ldots,{T}_{n}$,
we can compute the proportions in each cell, and then denote the vectors of
sample proportions $\hat{\boldsymbol{\theta}}(\boldsymbol{T})$. Then
$\psi(\hat{\boldsymbol{\theta}}(\boldsymbol{T}))$ is our estimated version of
the parameter of interest.

The parameter space under the null is
$\Theta_{0}(\psi_{0})=\\{\boldsymbol{\theta}:\psi(\boldsymbol{\theta})\leq\psi_{0}\\}.$
Let $\mathcal{S}$ denote the sample space, i.e., the set of possible outcomes
of the $\boldsymbol{T}$ random vector, and let
$\overline{\mathcal{S}}(\boldsymbol{T})=\\{\boldsymbol{t}\in\mathcal{S}:\psi(\hat{\boldsymbol{\theta}}(\boldsymbol{t}))\geq\psi(\hat{\boldsymbol{\theta}}(\boldsymbol{T}))\\}$
be the sub-sample space where the estimated parameter would be equal or more
extreme than the associated estimate using $\boldsymbol{T}$. Here “more
extreme” suggests observations that are less likely under the null hypothesis.
Then a p-value for testing $H_{0}$ can be defined as

$p(\boldsymbol{T},\psi_{0})=\sup_{\boldsymbol{\theta}\in\Theta_{0}(\psi_{0})}\sum_{t\in\overline{\mathcal{S}}(\boldsymbol{T})}f(t,\boldsymbol{\theta}).$

It is difficult to compute $p(\psi_{0})\equiv p(\boldsymbol{T},\psi_{0})$
analytically in general, so we propose the following numeric approximation
procedure.

Data: a large integer $B$, null hypothesis value $\psi_{0}$, subsample space
$\overline{\mathcal{S}}$

$p(\psi_{0})\leftarrow 1/B$;

for  _$i=1,\ldots,B$_ do

$\boldsymbol{\theta}^{(i)}\leftarrow()$ (empty vector);

for  _$j=1,\ldots,k$_ do

$\boldsymbol{\theta}^{(i)}\leftarrow(\boldsymbol{\theta}^{(i)},$ a random draw
from the uniform distribution on the probability simplex with dimension
$d_{j}$$)$;

end for

if  _$\psi(\boldsymbol{\theta}^{(i)}) >\psi_{0}$_ then

next;

else

$p(\psi_{0})\leftarrow\mbox{max}\left\\{p(\psi_{0}),\sum_{t\in\overline{\mathcal{S}}}f(t,\boldsymbol{\theta}^{(i)})\right\\}$;

end if

end for

Algorithm 1 Algorithm for computation of a p-value.

The approximation gets better for larger values of $B$. The p-value is
initialized at $1/B$ to avoid the possibility of obtaining a p-value equal to
0. The computationally difficult parts include enumeration of the sample space
$\mathcal{S}$, and calculation of $f$ for each element of that sample space.
It will typically be the case that it is only of interest to accurately
compute small values of $p(\psi_{0})$, hence if at some point a $p(\psi_{0})$
is found the exceeds some threshold, the loop can terminate early. We discuss
details of the implementation to be as efficient as possible in Section 3.

###  2.3 Computation of confidence intervals

Given that we can compute the p-value $p(\psi_{0})$, we can construct a
$100(1-\alpha)\%$ confidence interval by inverting the p-value function.
Consider the class of hypotheses indexed by $\psi_{0}$ that was previously
discussed, $H_{0\ell}:\psi(\boldsymbol{\theta})\leq\psi_{0}$ versus
$H_{1\ell}:\psi(\boldsymbol{\theta})>\psi_{0}$, where now we add an index
$\ell$ to denote that this class is associated with a one-sided confidence
interval using only the lower limit. Let $p_{\ell}(\psi_{0})$ the the p-value
associated with those hypotheses. Then the lower one-sided $100(1-\alpha)\%$
confidence limit for $\psi$ is the largest value $\psi_{\ell[\alpha]}$ that
satisfies $p_{\ell}(\psi_{\ell[\alpha]})\leq\alpha$. To compute the upper
confidence limit, we consider the class of hypotheses that switch direction of
the inequalities, $H_{0u}:\psi(\boldsymbol{\theta})\geq\psi_{0}$ versus
$H_{1\ell}:\psi(\boldsymbol{\theta})<\psi_{0}$ and its associated p-value,
$p_{u}(\psi_{0})$. Then the upper one-sided $100(1-\alpha)\%$ confidence limit
for $\psi$ is the smallest value $\psi_{u[\alpha]}$ that satisfies
$p_{u}(\psi_{u[\alpha]})\leq\alpha$. To get a two-sided $100(1-\alpha)\%$
central confidence interval, we simply combine the two one-sided
$100(1-\alpha/2)\%$ limits to get the interval
$(\psi_{\ell[\alpha/2]},\psi_{u[\alpha/2]})$. The central interval ensures
that the coverage error on either side is bounded by $\alpha/2$.

##  3 Implementation

A prerequisite for Algorithm 1 is an enumeration of the subsample space
$\overline{\mathcal{S}}$, which can be computationally expensive. As a first
step, we enumerate the sample space for each $T_{j},j\in 1,\ldots,k$. Let
$S(n,d)$ denote the sample space for a multinomial random variable with $n$
trials and $d$ categories. We compute $S(n,d)$ with a simple recursive
algorithm that we implemented in base R:

$S(d,n)=\begin{cases}\\{(0,n),(1,n-1),\ldots,(n,0)\\}\mbox{ if }d=2\\\
[\texttt{for i}\in(0,\ldots,n):\\{(i,s):s\in S(d-1,n-i)\\}]\mbox{ if
}d>2.\end{cases}$

Then to obtain $\mathcal{S}$ we take all possible combinations of elements,
one from each of the $\mathcal{S}_{1},\ldots,\mathcal{S}_{k}$. This is done in
a loop, at each iteration of which we obtain a $t\in\mathcal{S}$ for which we
also compute $\psi(\hat{\boldsymbol{\theta}}(t))$ and the sum of the log
multinomial coefficients corresponding to $t$ for each $j\in 1,\ldots,k$.

The steps of the p-value computation are nearly identical to those shown in
Algorithm 1. One difference is that the steps of randomly uniformly sampling
$\boldsymbol{\theta}_{i}$ and checking whether
$\psi(\boldsymbol{\theta}_{i})>\psi_{0}$ are vectorized. Thus, a set of
candidate parameters is sampled and checked at a time. The size of the set is
controlled by the argument chunksize, so that chunksize parameters are sampled
and checked at each of maxit iterations, giving a total of $B=$
chunksize$*$maxit iterations. This allows users to balance memory usage, which
increases with larger values of chunksize, with overhead due to starting a new
for loop iteration. We implemented a function to draw $n$ samples from the
probability simplex of dimension $d$ in rust (Matsakis and Klock II, 2014) and
incorporated it into our R package via rextendr (Wilke et al., 2024). The
calculation of the sum of multinomial probabilities over the subsample space
(in the last step of Algorithm 1) is also implemented in rust for speed and
memory safety.

The confidence interval is computed by finding the roots of
$p_{\ell}(\psi_{0})-\alpha/2$ and $p_{u}(\psi_{0})-\alpha/2$. To minimize the
number of p-value function evaluations during this process, we use the
interpolate, truncate, project (ITP) root finding algorithm (Oliveira and
Takahashi, 2020), which converges with fewer function evaluations compared to
the Brent method that is implemented in the R function uniroot. We implemented
the ITP algorithm in base R, though an add on package called itp also exists
(Northrop, 2023). To further reduce computation time, we have an optional
parameter that terminates the inner loop of the p-value computation early if
the current p-value is larger than the parameter, which takes the default
value of $\alpha/2+0.001$.

##  4 Simulations

We examine the numeric performance of our method of constructing confidence
intervals in comparison to the nonparametric bootstrap. We do this using two
examples, the Bhattacharyya coefficient, and symbolic nonparametric bounds. In
all cases we used 500 bootstrap replicates, and used a maximum of 2500
iterations of our method. Code to reproduce the numerical experiments is
available as Supplementary Material.

###  4.1 Bhattacharyya coefficient

Suppose $d_{j}=d$ and $n_{j}=n$ for all $j=1,\ldots,k$, and that we observe
the multinomial random vectors $T_{1},\ldots,T_{k}$ which have parameters
$\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{k}$. Let

$\psi(\boldsymbol{\theta})=\sum_{i=1}^{d}\left({\prod_{j=1}^{k}\theta_{ji}}\right)^{1/k},$

which is the extension of the Bhattacharyya coefficient to more than 2
distributions, as proposed in Kang and Wildes (2015). We examine the
performance of our proposed method in comparison to the nonparametric
bootstrap in terms of empirical 95% confidence interval coverage, for
different values of $d,n$ and $\boldsymbol{\theta}$. The results are shown in
Table 1. The nonparametric bootstrap has very poor coverage at small to
moderate sample sizes, while our proposed method appears to have at least
nominal coverage in all cases (although in many cases the coverage is
conservative).

Table 1: Bhattacharyya coefficient simulation: Empirical coverage (percent) of nominal 95% confidence intervals using the proposed method in comparison to the nonparametric bootstrap for the Bhattacharyya coefficient under different sample sizes, dimensions, and true values.  |  |  |  | coverage
---|---|---|---|---
sample size (n) | samples (k) | dimension (d) | true value | xactonomial | bootstrap
5 | 4 | 3 | 0.85 | 99.2 | 8.9
5 | 5 | 3 | 0.97 | 95.0 | 0.1
10 | 2 | 3 | 0.73 | 97.5 | 64.8
10 | 2 | 3 | 0.77 | 96.4 | 55.4
10 | 2 | 4 | 0.98 | 96.5 | 43.9
10 | 2 | 4 | 0.99 | 97.9 | 67.0
10 | 3 | 3 | 0.81 | 97.9 | 40.5
20 | 2 | 3 | 0.73 | 95.2 | 76.9
20 | 2 | 3 | 0.77 | 94.8 | 72.9
20 | 2 | 4 | 0.98 | 98.4 | 71.2
20 | 2 | 4 | 0.99 | 99.1 | 83.8

###  4.2 Stability of the approximation

Here we examine the degree to which the calculated p-value changes as the
number of iterations increases. Using the Bhattycharya coefficient example
with $d=10$, $k=2$, $n=10$, and a fixed observation of data, we computed
p-values at different values of $\psi_{0}$ corresponding to different true
p-values, and show how the estimate changes over a sequence of iterations of
the algorithm. Figure 1 shows the sequence of p-values. Since we are taking
the maximum at each iteration, the lines are nondecreasing. In all of these
scenarios, the computation stabilizes after around 2500 iterations. The
stabilization appears to happen even earlier for the smaller p-values, where
accuracy is more important. This and a few other experiments (not reported)
were used to guide the default arguments in our package xactonomial, but in
critical settings it is recommended to do a similar assessment for the
function $\psi$ of interest.

Figure 1: Illustration of the stability of the p-value computation for
increasing numbers of iterations.

###  4.3 Causal bounds on the risk difference

The setting of interest for this example is depicted in the directed acyclic
graph (DAG) of Figure 2, where $Z,X,Y$ are all binary taking values 0 or 1,
and $U$ is an unmeasured variable of arbitrary dimension and range. Assume we
have an independent and identically distributed sample
$(Z_{i},X_{i},Y_{i}),i=1,\ldots,n$ from a data generating mechanism that
coincides with the DAG.

$Z$$X$$Y$$U$ Figure 2: DAG depicting the setting of interest. The observed
variables $Z$, $X$, and $Y$ are all binary, while the unmeasured confounders
$U$ are arbitrary in dimension and range.

The DAG translates to the following structural equation model:

$\displaystyle z$ $\displaystyle=$ $\displaystyle g_{Z}(\epsilon_{z})$ (1)
$\displaystyle x$ $\displaystyle=$ $\displaystyle
g_{X}(\boldsymbol{u},z,\epsilon_{x})$ $\displaystyle y$ $\displaystyle=$
$\displaystyle g_{Y}(\boldsymbol{u},x,\epsilon_{y}),$

where the error terms $(\epsilon_{z},\epsilon_{x},\epsilon_{y})$ are
independent of each other, and of $U$.

We are interested in the parameter $\beta=p\\{Y(X=1)=1\\}-p\\{Y(X=0)=1\\}$,
where the notation $Y(X=x)$ denotes potential outcome which represents the
random variable $Y$ if everyone in the population were intervened upon to have
the treatment $X=x$ (Rubin, 1974). Under the DAG, and without any other
assumptions, $\beta$ is not identified, however tight bounds have been derived
for this setting in terms of observable conditional probabilities of the form
$p_{xy\cdot z}=p\\{X=x,Y=y|Z=z\\}$ (Balke, 1995; Balke and Pearl, 1997;
Swanson et al., 2018). They are

$\displaystyle\beta$ $\displaystyle\geq$
$\displaystyle\mbox{max}\left.\begin{cases}-1+p_{00\cdot 1}+p_{11\cdot 1},\\\
-1+p_{00\cdot 1}+p_{11\cdot 0}\\\ -1+p_{00\cdot 0}+p_{11\cdot 0},\\\
-1+p_{00\cdot 0}+p_{11\cdot 1},\\\ -2+2p_{00\cdot 0}+p_{01\cdot 1}+p_{11\cdot
0}+p_{11\cdot 1},\\\ -2+p_{00\cdot 0}+p_{00\cdot 1}+p_{10\cdot 0}+2p_{11\cdot
1},\\\ -2+2p_{00\cdot 1}+p_{01\cdot 0}+p_{11\cdot 0}+p_{11\cdot 1},\\\
-2+p_{00\cdot 0}+p_{00\cdot 1}+p_{10\cdot 1}+2p_{11\cdot
0},\end{cases}\right\\}$ (2)

$\displaystyle\beta$ $\displaystyle\leq$
$\displaystyle\mbox{min}\left.\begin{cases}1-p_{10\cdot 1}-p_{01\cdot 0},\\\
1-p_{10\cdot 1}-p_{01\cdot 1},\\\ 1-p_{10\cdot 0}-p_{01\cdot 0},\\\
1-p_{10\cdot 0}-p_{01\cdot 1},\\\ 2-2p_{10\cdot 1}-p_{01\cdot 0}-p_{01\cdot
1}-p_{11\cdot 0},\\\ 2-p_{00\cdot 1}-p_{10\cdot 0}-p_{10\cdot 1}-2p_{01\cdot
0},\\\ 2-2p_{10\cdot 0}-p_{01\cdot 0}-p_{01\cdot 1}-p_{11\cdot 1},\\\
2-p_{00\cdot 0}-p_{10\cdot 0}-p_{10\cdot 1}-2p_{01\cdot
1}.\end{cases}\right\\}$

Let

$T_{0}=\left(\begin{array}[]{c}X=0,Y=0|Z=0\\\ X=1,Y=0|Z=0\\\ X=0,Y=1|Z=0\\\
X=1,Y=1|Z=0\end{array}\right)\mbox{ and
}T_{1}=\left(\begin{array}[]{c}X=0,Y=0|Z=1\\\ X=1,Y=0|Z=1\\\ X=0,Y=1|Z=1\\\
X=1,Y=1|Z=1\end{array}\right).$

$T_{0}$ and $T_{1}$ can be viewed as independent multinomial random variables
with probabilities $(p_{00\cdot 0},p_{10\cdot 0},p_{01\cdot 0},p_{11\cdot 0})$
and $(p_{00\cdot 1},p_{10\cdot 1},p_{01\cdot 1},p_{11\cdot 1})$, respectively.
The expressions for the lower and upper bounds are thus real-valued functions
of these multinomial probabilties.

Focusing on the lower bound, let $\boldsymbol{p}=(p_{00\cdot 0},p_{10\cdot
0},p_{01\cdot 0},p_{11\cdot 0},p_{00\cdot 1},p_{10\cdot 1},p_{01\cdot
1},p_{11\cdot 1})$ denote the vector of probabilities, and let
$\psi(\boldsymbol{p})$ be the lower bound (the right-hand-side of expression
2). We consider different values for the vector of parameters $\boldsymbol{p}$
that give a range of true lower bounds and several sample sizes, over 2000
simulation replicates and present the empirical confidence interval coverage
of our method in comparison to the nonparametric bootstrap. The results are
shown in Table 2. The nonparametric bootstrap has less than nominal coverage
at small to moderate sample sizes, whereas our method has somewhat
conservative coverage in these settings.

Table 2: Causal Bounds Simulation: Empirical coverage (percent) of the nominal 95% confidence intervals of the proposed method in comparison to the nonparametric bootstrap for the lower bound of the causal risk difference under different sample sizes and true values.  |  | coverage
---|---|---
Total sample size (n) | True value | xactonomial | bootstrap
5 | -0.15 | 99.2 | 88.9
5 | -0.05 | 100.0 | 15.8
5 | 0.40 | 99.9 | 72.7
5 | 0.55 | 99.6 | 85.4
10 | -0.15 | 98.8 | 92.5
10 | -0.05 | 100.0 | 39.1
10 | 0.40 | 98.6 | 84.2
10 | 0.55 | 99.2 | 89.7
15 | -0.15 | 98.1 | 91.0
15 | -0.05 | 100.0 | 62.1
15 | 0.40 | 98.6 | 87.6
15 | 0.55 | 98.2 | 90.0
20 | -0.15 | 98.2 | 92.0
20 | -0.05 | 100.0 | 74.7
20 | 0.40 | 98.9 | 91.0
20 | 0.55 | 98.2 | 91.7

##  5 Description of the package

This method is implemented in the R package xactonomial, available at
https://sachsmc.github.io/xactonomial, whose primary feature is its eponymous
function. Users must provide three required arguments to the xactonomial
function: the function $\psi$ which maps the vector of multinomial parameters
to the real line, the data which must be provided as a list of $k$ vectors,
each of which represents the vector of observed counts in each cell for each
of the $k$ samples, and a vector giving the lower and upper limits of the
possible range of $\psi$. Optionally, the user can change the default value of
$\alpha$, which is 0.05, and the maximum number of iterations for the p-value
computation ($B$) is controlled by two arguments, the maximum number of
iterations, and the chunk size, that is, the size of the sample of random
draws from the $d_{j}$ unit simplex at each iteration.

It also only requires a few lines of R code to run:

    true_theta <- c(.45, .15, .3, .1, .05, .15, .4, .4)

    sample_data <- function(n) {

      T1 <- rmultinom(1, n[1], prob = true_theta[1:4])
      T2 <- rmultinom(1, n[2], prob = true_theta[5:8])

      list(T1 = c(T1), T2 = c(T2))

    }

    psi <- function(theta) {

      theta1 <- theta[1:4]
      theta2 <- theta[5:8]
      sum(sqrt(theta1 * theta2))

    }

    psi(true_theta)

    # [1] 0.8464102

    set.seed(2024)
    data <- sample_data(n = c(10, 10))
    xactonomial::xactonomial(psi, data, alpha = .05, psi0 = 0.75, psi_limits = c(0,1))

    # $estimate
    # [1] 0.8343818
    #
    # $conf.int
    # [1] 0.6225000 0.9970739
    #
    # $p.value
    # [1] 0.1051631

The xactonomial function returns a list that contains the estimate
$\psi(\hat{\boldsymbol{\theta}})$, the vector of the central $100(1-\alpha)\%$
confidence interval (which is the intersection of the two one-sided confidence
intervals, each at level $100(1-\alpha/2)$), and a p-value for the test at the
given null value. If psi0 is NULL, then only a confidence interval will be
calculated, and if not null and conf.int = FALSE, then only a p-value will be
calculated.

As mentioned above, some of the underlying functions are implemented in rust,
and incorporated into our package with the rextendr package (Wilke et al.,
2024). This led to dramatic gains in speed and memory usage, allowing us to do
the computations in a matter of minutes for total sample sizes up to about 60.
This is illustrated in Figure 3, which shows the timings of the p-value
computation after having enumerated the subsample space with the base R
implementation and the rust implementation. These experiments were done with
the Bhattacharyya coefficient setting with $k=2$, $d=4$, and equal sample
sizes in the two groups. The calculations were done on a laptop running Ubuntu
22 with an AMD Ryzen 7 and 32 GB of RAM. Further speed improvements could be
achieved by optimizing and/or memoising (i.e., pre-computing and storing) the
enumeration of the sample space.

Figure 3: Comparison of computation speed of p-values between the base R and
rust implementations.

##  6 Discussion

Our primary motivation for developing this method was to perform valid
inference on causal bounds in small sample settings or in other cases where
the bootstrap may not work. Symbolic nonparametric causal bounds have recently
been developed for a variety of complex settings (Sachs et al., 2023; Gabriel
et al., 2022, 2023b, 2023a). The expressions for the bounds are similar to our
example, and given estimates of those probabilities obtained from a sample,
one can obtain estimated bounds. Inference, however, has been largely ignored
in the bounds/partial identification literature, with most papers suggesting
the nonparametric bootstrap for obtaining confidence regions on the upper and
lower limits of the bounds. However, in small samples, bootstrap inference may
perform quite poorly. Our alternative method is much more accurate, has little
additional computational cost, and is broadly applicable to these settings.

When the parameter space is unbounded, the method still applies. However, in
practice, one must sample candidate parameters from the parameter space, and
the ITP algorithm requires upper and lower bounds for the root search. Hence
it may be advantageous to perform a transformation to a closed interval, e.g.,
log or logit, prior to calculating p-values and/or confidence intervals.

## Supplementary Material

Supplementary material includes the R code used to run the simulation study
and numerical examples. The R package implementing the method is available at
https://sachsmc.github.io/xactonomial

## Acknowledgements

This work utilized the computational resources of the NIH HPC Biowulf cluster.
(https://hpc.nih.gov), and thanks to Vivian Callier for help with interfacing
with it.

## Funding

EEG was partially supported by a grant from Novo Nordisk fonden
NNF22OC0076595.

## Declaration of Interest

The authors report there are no competing interests to declare.

## References

  * Agresti (1992) Alan Agresti.  A survey of exact inference for contingency tables.  _Statistical science_ , 7(1):131–153, 1992.
  * Balke and Pearl (1997) A. Balke and J. Pearl.  Bounds on treatment effects from studies with imperfect compliance.  _Journal of the American Statistical Association_ , 92:1171–1176, 1997.
  * Balke (1995) Alexander A Balke.  _Probabilistic Counterfactuals: Semantics, Computation, and Applications_.  PhD thesis, Department of Computer Science, University of California, Los Angeles, 1995.
  * Bhattacharyya (1946) A. Bhattacharyya.  On a measure of divergence between two multinomial populations.  _Sankhyā: The Indian Journal of Statistics (1933-1960)_ , 7(4):401–406, 1946.  ISSN 00364452.  URL http://www.jstor.org/stable/25047882.
  * Bickel and Freedman (1981) Peter J Bickel and David A Freedman.  Some asymptotic theory for the bootstrap.  _The annals of statistics_ , 9(6):1196–1217, 1981\.
  * Bickel and Sakov (2008) Peter J Bickel and Anat Sakov.  On the choice of m in the m out of n bootstrap and confidence bounds for extrema.  _Statistica Sinica_ , 18:967–985, 2008.
  * Chafai and Concordet (2009) Djalil Chafai and Didier Concordet.  Confidence regions for the multinomial parameter with small sample size.  _Journal of the American Statistical Association_ , 104(487):1071–1079, 2009.
  * Efron (1979) B Efron.  Bootstrap methods: Another look at the jackknife.  _The Annals of Statistics_ , 7(1):1–26, 1979\.
  * Frey (2009) Jesse Frey.  An exact multinomial test for equivalence.  _Canadian Journal of Statistics_ , 37(1):47–59, 2009.
  * Gabriel et al. (2022) Erin E Gabriel, Michael C Sachs, and Arvid Sjölander.  Causal bounds for outcome-dependent sampling in observational studies.  _Journal of the American Statistical Association_ , 117(538):939–950, 2022.
  * Gabriel et al. (2023a) Erin E Gabriel, Michael C Sachs, and Arvid Sjölander.  Sharp nonparametric bounds for decomposition effects with two binary mediators.  _Journal of the American Statistical Association_ , 118(544):2446–2453, 2023a.
  * Gabriel et al. (2023b) Erin E Gabriel, Arvid Sjölander, and Michael C Sachs.  Nonparametric bounds for causal effects in imperfect randomized experiments.  _Journal of the American Statistical Association_ , 118(541):684–692, 2023b.
  * Kang and Wildes (2015) Soo Min Kang and Richard P Wildes.  The n-distribution bhattacharyya coefficient.  _York University_ , 2015.
  * Malloy et al. (2021) Matthew L Malloy, Ardhendu Tripathy, and Robert D Nowak.  Optimal confidence sets for the multinomial parameter.  In _2021 IEEE International Symposium on Information Theory (ISIT)_ , pages 2173–2178. IEEE, 2021.
  * Matsakis and Klock II (2014) Nicholas D Matsakis and Felix S Klock II.  The rust language.  In _ACM SIGAda Ada Letters_ , volume 34, pages 103–104. ACM, 2014\.
  * Northrop (2023) Paul J. Northrop.  _itp: The Interpolate, Truncate, Project (ITP) Root-Finding Algorithm_ , 2023.  https://paulnorthrop.github.io/itp/, https://github.com/paulnorthrop/itp.
  * Oliveira and Takahashi (2020) Ivo FD Oliveira and Ricardo HC Takahashi.  An enhancement of the bisection method average performance preserving minmax optimality.  _ACM Transactions on Mathematical Software (TOMS)_ , 47(1):1–24, 2020.
  * R Core Team (2015) R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2015.  URL http://www.R-project.org/.
  * Resin (2023) Johannes Resin.  A simple algorithm for exact multinomial tests.  _Journal of Computational and Graphical Statistics_ , 32(2):539–550, 2023.
  * Rubin (1974) D.B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of Educational Psychology_ , 66(5):688–701, 1974.
  * Sachs et al. (2023) Michael C Sachs, Gustav Jonzon, Arvid Sjölander, and Erin E Gabriel.  A general method for deriving tight symbolic bounds on causal effects.  _Journal of Computational and Graphical Statistics_ , 32(2):567–576, 2023.
  * Swanson et al. (2018) Sonja A Swanson, Miguel A Hernán, Matthew Miller, James M Robins, and Thomas S Richardson.  Partial identification of the average treatment effect using instrumental variables: review of methods for binary instruments, treatments, and outcomes.  _Journal of the American Statistical Association_ , 113(522):933–947, 2018.
  * Wilke et al. (2024) Claus O. Wilke, Andy Thomason, Mossa M. Reimert, Ilia Kosenkov, and Malcolm Barrett.  _rextendr: Call Rust Code from R using the ’extendr’ Crate_ , 2024\.  URL https://extendr.github.io/rextendr/.  R package version 0.3.1.9000.