# How to understand the $X(2900)$?

Bo Wang<EMAIL_ADDRESS>School of Physical Science and Technology, Hebei
University, Baoding 071002, China  Key Laboratory of High-precision
Computation and Application of Quantum Field Theory of Hebei Province, Baoding
071002, China  Shi-Lin Zhu<EMAIL_ADDRESS>School of Physics and Center of
High Energy Physics, Peking University, Beijing 100871, China

###### Abstract

In this work, the $S$\- and $P$-wave $\bar{D}^{\ast}K^{\ast}$ interactions are
studied in a coupled-channel formalism to understand the recently observed
$X_{0}(2900)$ and $X_{1}(2900)$ at LHCb. The experimental event distributions
can be well described, and two states with $I(J^{P})=0(0^{+})$ and $0(1^{-})$
are yielded in a unified framework. The masses of the $0^{+}$ and $1^{-}$
states are consistent with the experimental data, but the width of the $0^{+}$
state is larger than that of the $1^{-}$ one. The $X_{1}(2900)$ can be
interpreted as the $P$-wave excitation of the ground-state $X_{0}(2900)$ in
the hadronic molecular picture. The $S$\- and $P$-wave multiplets in the
$\bar{D}^{\ast}K^{\ast}$ system have many members, so the present peak in the
$D^{-}K^{+}$ invariant mass distributions might contain multi substructures.

##  I Introduction

Recently, the LHCb Collaboration observed a clear peak in the $D^{-}K^{+}$
invariant mass spectrum of the $B^{+}\to D^{+}D^{-}K^{+}$ decay Aaij  _et al._
(2020a, b), where the helicity angle distribution shows evident $P$-wave
behavior. The peak is fitted with spin-$0$ and spin-$1$ states. Their
resonance parameters are determined to be

$\displaystyle X_{0}(2900):M$ $\displaystyle=$ $\displaystyle 2866\pm 7\pm
2~{}\text{MeV},$ $\displaystyle\Gamma$ $\displaystyle=$ $\displaystyle 57\pm
12\pm 4~{}\text{MeV};$ $\displaystyle X_{1}(2900):M$ $\displaystyle=$
$\displaystyle 2904\pm 5\pm 1~{}\text{MeV},$ $\displaystyle\Gamma$
$\displaystyle=$ $\displaystyle 110\pm 11\pm 4~{}\text{MeV}.$

The $D^{-}K^{+}$ decay channel implies their quark components should be
$\bar{c}\bar{s}ud$, which means that they are the fully open flavor exotic
hadrons.

Many theoretical interpretations have been proposed to understand the inner
structures of $X(2900)$, such as the molecular states from the
$\bar{D}^{\ast}K^{\ast}$ and $\bar{D}_{1}K$ interactions Chen  _et al._
(2020); He and Chen (2021); Liu  _et al._ (2020a); Hu  _et al._ (2021); Agaev
_et al._ (2020), the compact tetraquarks $\bar{c}\bar{s}ud$ Chen  _et al._
(2020); Karliner and Rosner (2020); He  _et al._ (2020); Wang (2020); Zhang
(2021); Wang  _et al._ (2021), and kinematic effects from the triangle
singularities Liu  _et al._ (2020b); Burns and Swanson (2021a). The production
and decay properties were investigated in Refs. Huang  _et al._ (2020); Chen
_et al._ (2021a); Burns and Swanson (2021b); Xiao  _et al._ (2021). One can
also consult Refs. Albuquerque  _et al._ (2021); Lü  _et al._ (2020); Mutuk
(2021); Tan and Ping (2020); Abreu (2021); Qi  _et al._ (2021); Chen (2021);
Hsiao and Yu (2021); Duan  _et al._ (2021); Kong  _et al._ (2021); Dong and
Zou (2021); Bondar and Milstein (2020); Chen  _et al._ (2021b) for other
pertinent works.

In Ref. Karliner and Rosner (2020), the $X_{0}(2900)$ was interpreted as the
$S$-wave compact tetraquark in a string-junction picture, while the subsequent
calculations in Ref. He  _et al._ (2020) supported it to be the radial excited
tetraquark with $J^{P}=0^{+}$. Meanwhile, the obtained masses with the refined
quark model calculations in Refs. Wang  _et al._ (2021); Lü  _et al._ (2020);
Tan and Ping (2020) are much lower than the measured mass of $X_{0}(2900)$. In
other words, the isosinglet $S$-wave compact tetraquark $\bar{c}\bar{s}ud$ can
not reconcile with the $X_{0}(2900)$. A typical feature of $X_{0}(2900)$ is
its mass below the $\bar{D}^{\ast}K^{\ast}$ threshold about $30$ MeV. The QCD
sum rule calculations Chen  _et al._ (2020); Agaev  _et al._ (2020) and the
one-boson-exchange inspired models He and Chen (2021); Liu  _et al._ (2020a)
all supported the $S$-wave $\bar{D}^{\ast}K^{\ast}$ molecule for
$X_{0}(2900)$. Besides, Refs. Huang  _et al._ (2020); Xiao  _et al._ (2021)
investigated the decays of $X_{0}(2900)$ and indicated large
$\bar{D}^{\ast}K^{\ast}$ component in its wave function. Another hint comes
from the width of $X_{0}(2900)$, which is close to the $K^{\ast}$ width, i.e.,
one may infer that the $X_{0}(2900)$ has an intrinsic $K^{\ast}$ width, and
the residual part arises from its decays. Thus, it seems the $X_{0}(2900)$ is
more likely to be the $\bar{D}^{\ast}K^{\ast}$ molecular state.

Most of the previous works cannot give a unified description for $X_{0}(2900)$
and $X_{1}(2900)$ if they are the genuine states. For example, Ref. Karliner
and Rosner (2020) relegated the $X_{1}(2900)$ to the $\bar{D}^{\ast}K^{\ast}$
rescattering effect. Ref. He and Chen (2021) gave a virtual state explanation
for $X_{1}(2900)$ that was generated from the $\bar{D}_{1}K$ interaction and
the bound solution in the $\bar{D}_{1}K$ channel needs an unnaturally large
cutoff. In Ref. Liu  _et al._ (2020a), the $X_{1}(2900)$ cannot be interpreted
as the $P$-wave $\bar{D}^{\ast}K^{\ast}$ molecule. If the $X_{1}(2900)$ is
indeed the $P$-wave compact tetraquark as suggested in Refs. Chen  _et al._
(2020); He  _et al._ (2020), then where is the $S$-wave ground state? It
definitely cannot be the $X_{0}(2900)$ as the results shown in Refs. He  _et
al._ (2020); Wang  _et al._ (2021); Lü  _et al._ (2020); Tan and Ping (2020).
In addition, we also need to answer the different dynamics in $b$-decays for
the formations of $X_{0}(2900)$ and $X_{1}(2900)$. In order to eschew the
dilemmas as mentioned above, this work is devoted to describing these two
states in a unified framework, and understanding their internal configurations
in a more natural way.

With the discoveries of more and more near-threshold exotic states Chen  _et
al._ (2016); Guo  _et al._ (2018); Liu  _et al._ (2019); Lebed  _et al._
(2017); Esposito  _et al._ (2017); Brambilla  _et al._ (2020); Chen  _et al._
(2021c, 2022); Meng  _et al._ (2022), we find the connections between hadronic
physics and nuclear physics become closer and closer. The concepts of
$\chi$EFT Weinberg (1990, 1991), which have been successfully substantialized
to describe the nuclear forces Bernard  _et al._ (1995); Epelbaum  _et al._
(2009); Machleidt and Entem (2011); Meißner (2016); Hammer  _et al._ (2020);
Rodriguez Entem  _et al._ (2020), can also be generalized to depict the
interactions between heavy hadrons Meng  _et al._ (2022); AlFiky  _et al._
(2006); Fleming  _et al._ (2007); Baru  _et al._ (2011); Valderrama (2012);
Braaten (2015); Schmidt  _et al._ (2018); Wang  _et al._ (2019a); Meng  _et
al._ (2019); Wang  _et al._ (2019b, 2020); Meng  _et al._ (2021); Chen  _et
al._ (2021d). This work dedicates to model the $\bar{D}^{\ast}K^{\ast}$
interactions with $S$\- and $P$-waves in a coupled-channel formalism via
mimicking the manipulations in $\chi$EFT. By fitting the $D^{-}K^{+}$ event
distributions in experiments, we analyze the pole positions in the scattering
$T$-matrix to see whether the $X_{0}(2900)$ and $X_{1}(2900)$ can be assigned
as the $S$-wave and $P$-wave excited $\bar{D}^{\ast}K^{\ast}$ molecules. This
can definitely resolve the problem we are facing now, and give a uniform
description of the $X_{0}(2900)$ and $X_{1}(2900)$.

This paper is organized as follows. In Sec. II, we establish the effective
potentials of $\bar{D}^{\ast}K^{\ast}$ and iterate them into the coupled-
channel Lippmann-Schwinger equations. In Sec. III, we present our numerical
results and discussions. In Sec. IV, we conclude this work with a short
summary.

##  II $\bar{D}^{\ast}K^{\ast}$ rescattering in coupled-channel formalism

In this work, we assume the $X_{0}(2900)$ and $X_{1}(2900)$ are the isosinglet
$S$\- and $P$-wave hadronic molecules which are generated from the
$\bar{D}^{\ast}K^{\ast}$ interactions, respectively. Thus the flavor wave
functions in $\bar{D}^{\ast}K^{\ast}$ and $\bar{D}K$ channels can be written
as

$\displaystyle|\bar{D}^{\ast}K^{\ast},I=0\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}[\bar{D}^{\ast 0}K^{\ast
0}-D^{\ast-}K^{\ast+}],$ $\displaystyle|\bar{D}K,I=0\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}[\bar{D}^{0}K^{0}-D^{-}K^{+}].$ (1)

The mass of $K^{\ast}$ is close to the nucleon mass, thus $K^{\ast}$ can be
regarded as a heavy matter field. We can generalize the approaches of
$\chi$EFT as that in the $NN$ system Bernard  _et al._ (1995); Epelbaum  _et
al._ (2009); Machleidt and Entem (2011); Meißner (2016); Hammer  _et al._
(2020); Rodriguez Entem  _et al._ (2020) to the $\bar{D}^{\ast}K^{\ast}$ case.
The effective potentials in the $\bar{D}^{\ast}K^{\ast}$ system can be
parameterized as follows,

$\displaystyle\mathcal{V}=\sum_{i}V_{i}(\bm{p}^{\prime},\bm{p})\mathcal{O}_{i}(\bm{p}^{\prime},\bm{p},\bm{\varepsilon},\bm{\varepsilon}^{\dagger},\bm{\varepsilon}^{\prime},\bm{\varepsilon}^{\prime\dagger}),$
(2)

where $\bm{p}$ and $\bm{p}^{\prime}$ represent the three momenta of the
initial and final states in the center of mass system (c.m.s) of
$\bar{D}^{\ast}K^{\ast}$, respectively. $\bm{\varepsilon}^{(\prime)}$ and
$\bm{\varepsilon}^{(\prime)^{\dagger}}$ denote the polarization vectors of
initial and final $\bar{D}^{\ast}$ ($K^{\ast}$), respectively. $V_{i}$ are the
scalar functions of $\bm{p}$ and $\bm{p}^{\prime}$. $\mathcal{O}_{i}$ are the
pertinent operators that can be constructed from the scalar products among
unit operator $\bm{1}$, vectors $\bm{p}^{(\prime)}$,
$\bm{\varepsilon}^{(\prime)}$ and $\bm{\varepsilon}^{(\prime)^{\dagger}}$,
which read

$\displaystyle\mathcal{O}_{1}$ $\displaystyle=$
$\displaystyle(\bm{\varepsilon}^{\dagger}\cdot\bm{\varepsilon})(\bm{\varepsilon}^{\prime\dagger}\cdot\bm{\varepsilon}^{\prime}),$
$\displaystyle\mathcal{O}_{2}$ $\displaystyle=$
$\displaystyle(\bm{\varepsilon}^{\prime\dagger}\cdot\bm{\varepsilon})(\bm{\varepsilon}^{\dagger}\cdot\bm{\varepsilon}^{\prime})-(\bm{\varepsilon}^{\prime\dagger}\cdot\bm{\varepsilon}^{\dagger})(\bm{\varepsilon}\cdot\bm{\varepsilon}^{\prime}),$
$\displaystyle\mathcal{O}_{3}$ $\displaystyle=$
$\displaystyle(\bm{q}\cdot\bm{\varepsilon}^{\prime\dagger})(\bm{q}\cdot\bm{\varepsilon})(\bm{\varepsilon}^{\dagger}\cdot\bm{\varepsilon}^{\prime})+(\bm{q}\cdot\bm{\varepsilon}^{\dagger})(\bm{q}\cdot\bm{\varepsilon}^{\prime})(\bm{\varepsilon}^{\prime\dagger}\cdot\bm{\varepsilon}),$
$\displaystyle\mathcal{O}_{4}$ $\displaystyle=$
$\displaystyle(\bm{q}\cdot\bm{\varepsilon}^{\dagger})(\bm{q}\cdot\bm{\varepsilon}^{\prime\dagger})(\bm{\varepsilon}^{\prime}\cdot\bm{\varepsilon})+(\bm{q}\cdot\bm{\varepsilon})(\bm{q}\cdot\bm{\varepsilon}^{\prime})(\bm{\varepsilon}^{\prime\dagger}\cdot\bm{\varepsilon}^{\dagger}),$
$\displaystyle\mathcal{O}_{5}$ $\displaystyle=$
$\displaystyle(\bm{k}\cdot\bm{\varepsilon}^{\prime\dagger})(\bm{k}\cdot\bm{\varepsilon})(\bm{\varepsilon}^{\dagger}\cdot\bm{\varepsilon}^{\prime})+(\bm{k}\cdot\bm{\varepsilon}^{\dagger})(\bm{k}\cdot\bm{\varepsilon}^{\prime})(\bm{\varepsilon}^{\prime\dagger}\cdot\bm{\varepsilon}),$
$\displaystyle\mathcal{O}_{6}$ $\displaystyle=$
$\displaystyle(\bm{k}\cdot\bm{\varepsilon}^{\dagger})(\bm{k}\cdot\bm{\varepsilon}^{\prime\dagger})(\bm{\varepsilon}^{\prime}\cdot\bm{\varepsilon})+(\bm{k}\cdot\bm{\varepsilon})(\bm{k}\cdot\bm{\varepsilon}^{\prime})(\bm{\varepsilon}^{\prime\dagger}\cdot\bm{\varepsilon}^{\dagger}),$
$\displaystyle\mathcal{O}_{7}$ $\displaystyle=$
$\displaystyle(\bm{\varepsilon}^{\prime\dagger}\cdot\bm{\varepsilon}^{\prime})(\bm{\varepsilon}^{\dagger}\times\bm{\varepsilon})\cdot(\bm{k}\times\bm{q})$
(3)
$\displaystyle+(\bm{\varepsilon}^{\dagger}\cdot\bm{\varepsilon})(\bm{\varepsilon}^{\prime\dagger}\times\bm{\varepsilon}^{\prime})\cdot(\bm{k}\times\bm{q}),\dots,$

with $\bm{q}=\bm{p}-\bm{p}^{\prime}$ the transferred momentum and
$\bm{k}=(\bm{p}^{\prime}+\bm{p})/2$ the average momentum. The ellipsis denotes
the higher terms with more $\bm{q}$ and $\bm{k}$. The operators
$\mathcal{O}_{1}$, $\mathcal{O}_{2}$, $\mathcal{O}_{3,\dots,6}$ and
$\mathcal{O}_{7}$ account for the central force, the spin-spin interaction,
the tensor force, and the spin-orbital ($SL$) coupling, respectively.

In order to describe the transitions between $\bar{D}^{\ast}K^{\ast}$ and
$\bar{D}K$, we define the corresponding transition operators,

$\displaystyle\mathcal{O}_{1}^{\prime}$ $\displaystyle=$
$\displaystyle\bm{\varepsilon}\cdot\bm{\varepsilon}^{\prime}\text{~{}or~{}}\bm{\varepsilon}^{\dagger}\cdot\bm{\varepsilon}^{\prime\dagger},$
$\displaystyle\mathcal{O}_{2}^{\prime}$ $\displaystyle=$
$\displaystyle(\bm{q}\cdot\bm{\varepsilon})(\bm{q}\cdot\bm{\varepsilon}^{\prime})\text{~{}or~{}}(\bm{q}\cdot\bm{\varepsilon}^{\dagger})(\bm{q}\cdot\bm{\varepsilon}^{\prime\dagger}),$
$\displaystyle\mathcal{O}_{3}^{\prime}$ $\displaystyle=$
$\displaystyle(\bm{k}\cdot\bm{\varepsilon})(\bm{k}\cdot\bm{\varepsilon}^{\prime})\text{~{}or~{}}(\bm{k}\cdot\bm{\varepsilon}^{\dagger})(\bm{k}\cdot\bm{\varepsilon}^{\prime\dagger}).$
(4)

For performing the partial wave decompositions in the following, the
polarization vectors in above equations can be transformed into the spin
transition operator $\bm{S}_{t}$ with the approach given in the appendix C of
Ref. Wang  _et al._ (2019b).

The effective potentials of the elastic
($\bar{D}^{\ast}K^{\ast}\to\bar{D}^{\ast}K^{\ast}$) and the inelastic
($\bar{D}^{\ast}K^{\ast}\to\bar{D}K$, or in reverse) channels have both
contributions from the contact and one-pion-exchange (OPE) interactions. The
short-range contact terms for the elastic (el) and inelastic (in) scatterings
are given as, respectively,

$\displaystyle\mathcal{V}_{\text{ct}}^{\text{el}}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{7}C_{i}\mathcal{O}_{i}+(C_{8}\bm{q}^{2}+C_{9}\bm{k}^{2})\mathcal{O}_{1}$
(5) $\displaystyle+(C_{10}\bm{q}^{2}+C_{11}\bm{k}^{2})\mathcal{O}_{2},$
$\displaystyle\mathcal{V}_{\text{ct}}^{\text{in}}$ $\displaystyle=$
$\displaystyle
C_{2}^{\prime}\mathcal{O}_{2}^{\prime}+C_{3}^{\prime}\mathcal{O}_{3}^{\prime}+(C_{4}^{\prime}\bm{q}^{2}+C_{5}^{\prime}\bm{k}^{2})\mathcal{O}_{1}^{\prime},$
(6)

where $C_{i}^{(\prime)}$ are the low energy constants (LECs).

The long-distance interaction for $\bar{D}^{\ast}K^{\ast}$ is depicted by the
OPE contribution. The OPE effective potential can be easily obtained from the
chiral Lagrangians, in which the anticharmed meson and pion coupling reads
Wise (1992); Manohar and Wise (2000)

$\displaystyle\mathcal{L}_{\tilde{\mathcal{H}}\varphi}=g_{\varphi}\langle\bar{\tilde{\mathcal{H}}}\gamma^{\mu}\gamma_{5}u_{\mu}\tilde{\mathcal{H}}\rangle,$
(7)

where $g_{\varphi}\simeq 0.57$ is the axial coupling. $\tilde{\mathcal{H}}$
and $u_{\mu}$ are the superfield for anticharmed mesons and axial-vector
current for pion, respectively, their expressions can be found in Refs. Meng
_et al._ (2019); Wang  _et al._ (2019b). In SU(2) case, if we treat the
$\bar{s}$ quark as a relatively heavy quark, then the $K^{\ast}$ and pion
coupling can be formulated as the same form as that in Eq. (7). Its axial
coupling constant $g_{\varphi}\simeq 1.12$ is determined from the decay width
of $K^{\ast}\to K\pi$ Zyla  _et al._ (2020).

Figure 1: The OPE diagrams in time-ordered perturbation theory. The black
solid and dashed lines denote the $\bar{D}^{(\ast)}/K^{(\ast)}$ and pion,
respectively. The red-dashed (horizontal) line indicates the time at which the
intermediate state is evaluated. In diagram (a), the particle $B$ feels the
effect of a pion that from the source particle $A$ emitted at an earlier time,
thus the $A$ and $B$ particles interact in this case through a retarded
propagator. In diagram (b), from the $B$'s point of view, the effect is felt
before the source $A$ emitted the pion, thus the pion propagator in this case
is called an advanced propagator.

The OPE effective potentials for elastic and inelastic scatterings read

$\displaystyle\mathcal{V}_{\text{OPE}}^{\text{(a)}}$ $\displaystyle=$
$\displaystyle-\frac{3g_{\varphi}^{2}}{4f_{\varphi}^{2}}\frac{\mathcal{O}^{\text{ch}}}{2E_{\pi}[E_{\pi}+E_{A^{\prime}}+E_{B}-E]},$
(8) $\displaystyle\mathcal{V}_{\text{OPE}}^{\text{(b)}}$ $\displaystyle=$
$\displaystyle-\frac{3g_{\varphi}^{2}}{4f_{\varphi}^{2}}\frac{\mathcal{O}^{\text{ch}}}{2E_{\pi}(E_{\pi}+E_{A}+E_{B^{\prime}}-E)},$
(9)

where the $\mathcal{V}_{\text{OPE}}^{\text{(a)}}$ and
$\mathcal{V}_{\text{OPE}}^{\text{(b)}}$ represent the contributions from Figs.
1(a) and 1(b), respectively. The
$\mathcal{O}^{\text{ch}}=-\bm{q}^{2}\mathcal{O}_{2}+\mathcal{O}_{3}-\mathcal{O}_{4}$
for the elastic channel, and
$\mathcal{O}^{\text{ch}}=\mathcal{O}_{2}^{\prime}$ for the inelastic channels,
respectively. $g_{\varphi}=\sqrt{0.57\times 1.12}$ is the redefined coupling
constant, $f_{\varphi}=92.4$ MeV is the pion decay constant, and
$m_{\pi}\simeq 140$ MeV is the pion mass. Note that the non-static effect (the
energy dependence) for the OPE potentials is considered with the time-ordered
perturbation theory (see Refs. Baru  _et al._ (2011); Du  _et al._ (2020);
Wang  _et al._ (2018); Baru  _et al._ (2019)). The energies appearing in Eqs.
(8) and (9) are given as

$\displaystyle E_{\pi}$ $\displaystyle=$
$\displaystyle\sqrt{\bm{q}^{2}+m_{\pi}^{2}},$ (10) $\displaystyle
E_{i^{(\prime)}}$ $\displaystyle=$ $\displaystyle
m_{i^{(\prime)}}+\frac{p^{(\prime)2}}{2m_{i^{(\prime)}}},\quad i=A,B,$ (11)

and $E$ stands for the total energy of the system.

The $X_{0}(2900)$ and $X_{1}(2900)$ are observed in the $b$-decay process
$B^{+}\to D^{+}D^{-}K^{+}$, so we need to simulate this production process.
The corresponding Feynman diagrams are shown in Fig. 2. The reactions in Fig.
2 can be expressed with the coupled-channel Lippmann-Schwinger equations
(LSEs),

$\displaystyle\mathcal{U}_{\alpha}^{j}(E,\bm{p})=$
$\displaystyle\mathcal{M}_{\alpha}^{j}+\sum_{\beta}\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\mathcal{V}_{\alpha\beta}^{j}(E,\bm{p},\bm{q})\mathcal{G}_{\beta}(E,\bm{q})\mathcal{U}_{\beta}^{j}(E,\bm{q}),$
(12)

where the subscript $\alpha(\beta)=1,2$ denotes the corresponding channels
[the $\bar{D}K$ and $\bar{D}^{\ast}K^{\ast}$ channels are labeled as $1$ and
$2$, respectively, e.g., see Fig. 2], while the superscript $j$ represents the
fixed total angular momentum. $\mathcal{M}$ is the direct production
amplitude. $E$ stands for the invariant mass of the $D^{-}K^{+}$ system.
$\mathcal{G}$ is the Green's function of the intermediate states, which read

$\displaystyle\mathcal{G}_{\beta}(E,\bm{q})=\frac{2\mu_{\beta}}{\bm{p}_{\beta}^{2}-\bm{q}^{2}+i\epsilon},|\bm{p}_{\beta}|=\sqrt{2\mu_{\beta}(E-m_{\text{th}}^{\beta})},$
(13)

with $\mu_{\beta}$ and $m_{\text{th}}^{\beta}$ the reduced mass and the
threshold of the $\beta$-th channel. The width of $K^{\ast}$ is $50.8$ MeV
Zyla  _et al._ (2020), which is considered in $m_{\text{th}}^{\beta}$ by using
a complex mass $m-i\Gamma/2$ for the $K^{\ast}$ Du  _et al._ (2020).

Figure 2: Diagrams (a) and (b) describe the direct production and rescattering
contribution, respectively. The gray circle with cross represents the
effective $B^{+}\to D^{+}\bar{D}K$ and $B^{+}\to D^{+}\bar{D}^{\ast}K^{\ast}$
coupling, while the gray box in diagram (b) signifies the rescattering
$T$-matrix of the $\bar{D}^{\ast}K^{\ast}$ system. The dashed line denotes the
spectator $D^{+}$, while the solid lines stand for the involved $\bar{D}K$
($\bar{D}^{\ast}K^{\ast}$) in the rescattering.

In order to solve the LSEs in Eq. (II), we need to make the partial wave
decomposition. The effective potentials in Eq. (2) are given in the plane wave
helicity basis, which can be transformed to the partial wave ($\ell sj$) basis
via the approach in Ref. Golak  _et al._ (2010). If the $X_{0}(2900)$ and
$X_{1}(2900)$ are the $S$\- and $P$-wave $\bar{D}^{\ast}K^{\ast}$ molecules
with $J^{P}=0^{+}$ and $J^{P}=1^{-}$, respectively. Then the corresponding
effective potentials in partial wave bases with the coupled-channels

$\displaystyle
j=0:|\bar{D}K\rangle_{{}^{1}S_{0}},|\bar{D}^{\ast}K^{\ast}\rangle_{{}^{1}S_{0}},|\bar{D}^{\ast}K^{\ast}\rangle_{{}^{5}D_{0}}$
(14)

and

$\displaystyle j$ $\displaystyle=$ $\displaystyle
1:|\bar{D}K\rangle_{{}^{1}P_{1}},|\bar{D}^{\ast}K^{\ast}\rangle_{{}^{1}P_{1}},|\bar{D}^{\ast}K^{\ast}\rangle_{{}^{3}P_{1}},|\bar{D}^{\ast}K^{\ast}\rangle_{{}^{5}P_{1}}$

can be written as the $3\times 3$ and $4\times 4$ matrices, respectively.

The contact interactions through the partial wave decompositions in $S$-wave
with the coupled-channels read

$[\mathcal{V}_{\text{ct}}]_{\ell,\ell^{\prime}}^{S}=\left[\begin{array}[]{ccc}0&\mathcal{C}_{12}(p^{2}+p^{\prime
2})&\mathcal{C}_{13}p^{2}\\\ \mathcal{C}_{12}(p^{2}+p^{\prime
2})&\mathcal{C}_{22}+\mathcal{C}_{22}^{\prime}(p^{2}+p^{\prime
2})&\mathcal{C}_{23}p^{2}\\\ \mathcal{C}_{13}p^{\prime
2}&\mathcal{C}_{23}p^{\prime 2}&0\end{array}\right],$ (16)

where

$\displaystyle\mathcal{C}_{12}$ $\displaystyle=$
$\displaystyle-\frac{\pi}{\sqrt{3}}(4C_{2}^{\prime}+C_{3}^{\prime}+12C_{4}^{\prime}+3C_{5}^{\prime}),$
$\displaystyle\mathcal{C}_{13}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}}{\sqrt{3}}\pi(4C_{2}^{\prime}+C_{3}^{\prime}),$
$\displaystyle\mathcal{C}_{22}$ $\displaystyle=$ $\displaystyle
4\pi(C_{1}-2C_{2}),$ $\displaystyle\mathcal{C}_{22}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\pi}{3}(8C_{3}+24C_{4}+2C_{5}+6C_{6}+12C_{8}+3C_{9}$
$\displaystyle-24C_{10}-6C_{11}),$ $\displaystyle\mathcal{C}_{23}$
$\displaystyle=$
$\displaystyle-\frac{\sqrt{2}\pi}{3}(8C_{3}+12C_{4}+2C_{5}+3C_{6}),$ (17)

while the ones for the $P$-wave are given as

$[\mathcal{V}_{\text{ct}}]_{\ell,\ell^{\prime}}^{P}=\left[\begin{array}[]{cccc}0&\mathbb{C}_{12}pp^{\prime}&0&\mathbb{C}_{14}pp^{\prime}\\\
\mathbb{C}_{12}pp^{\prime}&\mathbb{C}_{22}pp^{\prime}&0&\mathbb{C}_{24}pp^{\prime}\\\
0&0&\mathbb{C}_{33}pp^{\prime}&0\\\
\mathbb{C}_{14}pp^{\prime}&\mathbb{C}_{24}pp^{\prime}&0&\mathbb{C}_{44}pp^{\prime}\end{array}\right],$
(18)

where

$\displaystyle\mathbb{C}_{12}$ $\displaystyle=$
$\displaystyle\frac{2\pi}{3\sqrt{3}}(4C_{2}^{\prime}-C_{3}^{\prime}+12C_{4}^{\prime}-3C_{5}^{\prime}),$
$\displaystyle\mathbb{C}_{14}$ $\displaystyle=$
$\displaystyle\frac{2\sqrt{5}\pi}{3\sqrt{3}}(4C_{2}^{\prime}-C_{3}^{\prime}),$
$\displaystyle\mathbb{C}_{22}$ $\displaystyle=$
$\displaystyle\frac{2\pi}{9}(-8C_{3}-24C_{4}+2C_{5}+6C_{6}-12C_{8}+3C_{9}$
$\displaystyle+24C_{10}-6C_{11}),$ $\displaystyle\mathbb{C}_{24}$
$\displaystyle=$
$\displaystyle\frac{2\pi}{9\sqrt{5}}(-40C_{3}-60C_{4}+10C_{5}+15C_{6}+12C_{7}),$
$\displaystyle\mathbb{C}_{33}$ $\displaystyle=$
$\displaystyle\frac{\pi}{3}(12C_{3}-3C_{5}-8C_{7}-8C_{8}+2C_{9}$
$\displaystyle+8C_{10}-2C_{11}),$ $\displaystyle\mathbb{C}_{44}$
$\displaystyle=$
$\displaystyle-\frac{\pi}{45}(220C_{3}-55C_{5}+144C_{7}+120C_{8}$ (19)
$\displaystyle-30C_{9}+120C_{10}-30C_{11}).$

In above potentials, we switch off the $\bar{D}K$ interaction for reducing the
free parameters, which may be described by the chiral perturbation theory Liu
_et al._ (2009). However, the $\bar{D}K$ lies far below the
$\bar{D}^{\ast}K^{\ast}$ threshold, thus their interactions should be
irrelevant to the physics around the $\bar{D}^{\ast}K^{\ast}$ threshold.

The LSEs in Eq. (II) is regulated with a form factor, in which the
$\mathcal{V}_{\ell,\ell^{\prime}}$ is multiplied by the Gaussian regulator,

$\displaystyle\mathcal{V}_{\ell,\ell^{\prime}}\to\mathcal{V}_{\ell,\ell^{\prime}}\exp\left[-\frac{p^{2}}{\Lambda^{2}}-\frac{p^{\prime
2}}{\Lambda^{\prime 2}}\right].$ (20)

At the $X(2900)$ energy, the center of mass momentum of the $\bar{D}K$ channel
is $0.74$ GeV, and hence we use a relatively hard cutoff $\Lambda=1.0$ GeV for
this channel. For the elastic channel scattering, we vary the cutoff
$\Lambda^{\prime}$ [in this case, the $\Lambda$ and $\Lambda^{\prime}$ in Eq.
(20) have the equal value] over a wide range, i.e., from a soft scale $0.3$
GeV to a hard scale $1.0$ GeV.

Additionally, we also need to mimic the production vertex. We assume the
$D^{+}D^{-}K^{+}$ and $D^{+}D^{\ast-}K^{\ast+}$ are produced from a point-like
source. Then the $S$\- and $P$-wave production amplitudes can be parametrized
as Back  _et al._ (2018)

$\displaystyle\mathcal{M}^{S}_{\\{B^{+}\to D^{+}D^{-}K^{+},B^{+}\to
D^{+}D^{\ast-}K^{\ast+}\\}_{\alpha}}$
$\displaystyle~{}\quad=\\{\tilde{g}_{s},\tilde{g}_{s}^{\prime}\bm{\varepsilon}^{\dagger}\cdot\bm{\varepsilon}^{\prime\dagger}\\},$
(21) $\displaystyle\mathcal{M}^{P}_{\\{B^{+}\to D^{+}D^{-}K^{+},B^{+}\to
D^{+}D^{\ast-}K^{\ast+}\\}_{\alpha}}$
$\displaystyle~{}\quad=\\{\tilde{g}_{p}\bm{k}_{1}\cdot\bm{k}_{2},\tilde{g}_{p}^{\prime}(\bm{k}_{1}\cdot\bm{k}_{2}^{\prime})(\bm{\varepsilon}^{\dagger}\cdot\bm{\varepsilon}^{\prime\dagger})$
$\displaystyle~{}\quad+\tilde{g}_{p}^{\prime\prime}(\bm{k}_{1}\cdot\bm{\varepsilon}^{\dagger})(\bm{k}_{2}^{\prime}\cdot\bm{\varepsilon}^{\prime\dagger})+\tilde{g}_{p}^{\prime\prime\prime}(\bm{k}_{1}\cdot\bm{\varepsilon}^{\prime\dagger})(\bm{k}_{2}^{\prime}\cdot\bm{\varepsilon}^{\dagger})\\},$
(22)

where $\tilde{g}_{s(p)}^{(\prime)}$ represents the production strengths in
corresponding channels. $\bm{k}_{1}$ and $\bm{k}_{2}^{(\prime)}$ denote the
three momentum of $D^{+}$ and $D^{-}$ ($D^{\ast-}$) in the c.m.s of
$D^{-}K^{+}$ ($D^{\ast-}K^{\ast+}$).

Projecting the production amplitudes to the bases in Eqs. (14) and (II),
respectively, one easily obtains

$\displaystyle\mathcal{M}^{j=0}_{\\{B^{+}\to D^{+}D^{-}K^{+},B^{+}\to
D^{+}D^{\ast-}K^{\ast+}\\}_{i}}$
$\displaystyle~{}\quad=\mathcal{N}_{s}\\{1,\mathcal{R}_{s},0\\},$ (23)
$\displaystyle\mathcal{M}^{j=1}_{\\{B^{+}\to D^{+}D^{-}K^{+},B^{+}\to
D^{+}D^{\ast-}K^{\ast+}\\}_{i}}$
$\displaystyle~{}\quad=\mathcal{N}_{p}\\{k_{1}k_{2},\mathcal{R}_{p}k_{1}k_{2}^{\prime},\mathcal{R}_{p}^{\prime}k_{1}k_{2}^{\prime},\mathcal{R}_{p}^{\prime\prime}k_{1}k_{2}^{\prime}\\},$
(24)

where the irrelevant factors are absorbed into the $\mathcal{N}_{s(p)}$, while
the $\mathcal{R}_{s(p)}$ describes the relative strengths between
$|\bar{D}K\rangle_{{}^{1}S_{0}}$ and
$|\bar{D}^{\ast}K^{\ast}\rangle_{{}^{1}S_{0}}$
($|\bar{D}K\rangle_{{}^{1}P_{1}}$ and
$|\bar{D}^{\ast}K^{\ast}\rangle_{{}^{1}P_{1}}$) channels, and similar meanings
for the $\mathcal{R}_{p}^{\prime,\prime\prime}$.

With the above preparations, the differential decay width for $B^{+}\to
D^{+}D^{-}K^{+}$ reads

$\displaystyle\frac{d\Gamma}{dE}=\mathcal{N}\frac{|\tilde{\bm{k}}_{1}||\bm{k}_{2}^{\ast}|}{4(2\pi)^{3}m_{B}^{2}}\left[\mathcal{R}_{sp}\left|\mathcal{U}_{1}^{j=0}\right|^{2}+\left|\mathcal{U}_{1}^{j=1}\right|^{2}\right],$
(25)

where an overall normalization factor $\mathcal{N}$ is used to match the event
distributions in experiments. $\tilde{\bm{k}}_{1}$ and $\bm{k}_{2}^{\ast}$ are
the three momentum of the spectator $D^{+}$ in the c.m.s of $B^{+}$ and the
three momentum of $D^{-}$ in the c.m.s of $D^{-}K^{+}$, respectively.
$\mathcal{R}_{sp}$ is introduced to designate the ratio
$\mathcal{N}_{s}^{2}/\mathcal{N}_{p}^{2}$.

##  III Numerical results and discussions

Unlike the $D_{s0}(2317)$ and $D_{s1}(2460)$ Zyla  _et al._ (2020), there is
no $\bar{c}\bar{s}$ bare core in the hadron spectrum, so the $X(2900)$
provides us a relatively clean environment to study the coupled-channel
effects between open-charm and open-strange mesons.

The values of partial wave LECs are constrained by fitting to the $D^{-}K^{+}$
event candidates that measured by the LHCb Aaij  _et al._ (2020b), which are
given in Fig. 3 for the $S$\- and $P$-waves. One can find that some LECs are
sensitive to the variations of cutoff, such as $\mathcal{C}_{23}$,
$\mathbb{C}_{22}$, $\mathbb{C}_{33}$, $\mathbb{C}_{44}$ and $\mathbb{C}_{24}$
[in effective field theory, the LECs are cutoff dependent, $C_{i}\equiv
C_{i}(\Lambda)$, i.e., the cutoff dependence is absorbed by the LECs, which
makes the observables cutoff independent]. The distribution behaviors of these
LECs with the cutoff is similar to that of the spin singlet case of $NN$
scattering Nogga  _et al._ (2005), while the LECs of the spin triplet $NN$
(the channel for deuteron) scattering show limit-cycle-like behavior (an
unsmooth dependence on the cutoff). Meanwhile, the other LECs are not very
sensitive to cutoffs, such as the $\mathcal{C}_{12}$, $\mathcal{C}_{13}$,
$\mathbb{C}_{12}$ and $\mathbb{C}_{14}$. A more general analysis on the
nonrelativistic two-body scattering of the single-channel case with
renormalization-group invariance constraints was provided by Birse et al Birse
_et al._ (1999).

Figure 3: Changes of the $S$\- and $P$-wave LECs when the elastic channel's
cutoff is varied from $0.3$ GeV to $1.0$ GeV and the inelastic channel's
cutoff is fixed at $1.0$ GeV.

The fitted line shape of the $S$\- and $P$-wave components, as well as their
total contributions are plotted in Fig. 4. We find the results are also
insensitive to cutoffs, the $\chi^{2}/{\text{d.o.f}}$ is around $1.3$ when
$\Lambda^{\prime}\in[0.3,1.0]$ GeV. The fit with $\Lambda=1.0$ GeV,
$\Lambda^{\prime}=0.5$ GeV gives the least $\chi^{2}$ ($\simeq 1.25$), so we
adopt the results in this fit as our outputs. The corresponding values of LECs
with the errors are given in Table. 1.

Table 1: The values of LECs $(\mathcal{C}_{22},\mathcal{C}_{22}^{\prime},\mathcal{C}_{23},\mathcal{C}_{12},\mathcal{C}_{13},\mathbb{C}_{22},\mathbb{C}_{33},\mathbb{C}_{44},\mathbb{C}_{24},\mathbb{C}_{12},\mathbb{C}_{14})/10^{4}$ with errors that fitted at $\Lambda=1.0$ GeV and $\Lambda^{\prime}=0.5$ GeV. $\mathcal{C}_{22}$ (GeV-2) |  $\mathcal{C}_{22}^{\prime}$ (GeV-4) |  $\mathcal{C}_{23}$ (GeV-4) |  $\mathcal{C}_{12}$ (GeV-4) |  $\mathcal{C}_{13}$ (GeV-4) |
---|---|---|---|---|---
$0.428\pm 0.0823$ | $-0.0466\pm 0.0151$ | $-1.18\pm 0.324$ | $-0.102\pm 0.0435$ | $0.0336\pm 0.0117$ |
$\mathbb{C}_{22}$ (GeV-4) |  $\mathbb{C}_{33}$ (GeV-4) |  $\mathbb{C}_{44}$ (GeV-4) |  $\mathbb{C}_{24}$ (GeV-4) |  $\mathbb{C}_{12}$ (GeV-4) |  $\mathbb{C}_{14}$ (GeV-4)
$-0.172\pm 0.0216$ | $-1.05\pm 0.251$ | $-0.548\pm 0.134$ | $-0.734\pm 0.188$ | $0.175\pm 0.0534$ | $-0.0352\pm 0.0163$

Two pronounced peaks around the $\bar{D}^{\ast}K^{\ast}$ threshold are
obtained. One can notice that the signal of $P$-wave is much stronger than
that of the $S$-wave, this can naturally explain why the angular distribution
of $D^{-}K^{+}$ in experiments is overwhelmed by the $P$-wave structure Aaij
_et al._ (2020a, b).

Figure 4: The fitted $D^{-}K^{+}$ invariant mass distributions for $B^{+}\to
D^{+}D^{-}K^{+}$ decay. The experimental data are extracted from Ref. Aaij
_et al._ (2020b), where the reflection contributions from the charmonia are
subtracted. The dashed, dot-dashed and solid lines denote the $S$-, $P$\- and
$S+P$-wave contributions, respectively. The lineshapes are obtained with the
cutoffs $\Lambda=1.0$ GeV, and $\Lambda^{\prime}=0.5$ GeV.

We then study the pole structures of the production $\mathcal{U}$-matrix in
different Riemann sheets (the poles of the $\mathcal{U}$-matrix are the same
as the $T$-matrix in this sense), which can be achieved through the analytical
continuation of the Green's functions in Eq. (13):

$\displaystyle\mathcal{G}_{1}(E)$ $\displaystyle\to$
$\displaystyle\mathcal{G}_{1}(E)+\zeta_{1}\frac{i\mu_{1}p_{1}}{4\pi^{2}},$
$\displaystyle\mathcal{G}_{2}(E)$ $\displaystyle\to$
$\displaystyle\mathcal{G}_{2}(E)+\zeta_{2}\frac{i\mu_{2}p_{2}}{4\pi^{2}}.$
(26)

Two channel coupling generally introduces four Riemann sheets Frazer and
Hendry (1964), then the four Riemann sheets $(\zeta_{1},\zeta_{2})$ in the
complex energy plane are defined as

$\displaystyle{\rm Sheet~{}I}:(\zeta_{1},\zeta_{2})$ $\displaystyle=$
$\displaystyle(0,0),$ $\displaystyle{\rm Sheet~{}II}:(\zeta_{1},\zeta_{2})$
$\displaystyle=$ $\displaystyle(1,0),$ $\displaystyle{\rm
Sheet~{}III}:(\zeta_{1},\zeta_{2})$ $\displaystyle=$ $\displaystyle(1,1),$
$\displaystyle{\rm Sheet~{}IV}:(\zeta_{1},\zeta_{2})$ $\displaystyle=$
$\displaystyle(0,1),$ (27)

such as the physical sheet is denoted as $(0,0)$ in this definition. This is
analogous to the commonly used notations in complex momentum plane Frazer and
Hendry (1964). Considering a lower channel-$1$ with particles $A$ and $B$, and
a higher channel-$2$ with particles $C$ and $D$, a pole in Sheet II near
$m_{C}+m_{D}$ corresponds to a bound state of $CD$ (channel-$2$), such a pole
is manifested as a resonance in channel-$1$ ($AB$) scattering process. For
more general discussions on the pole distributions and classifications in
different Riemann sheets, we refer to Frazer and Hendry (1964); Eden and
Taylor (1964); Badalian  _et al._ (1982).

The $S$\- and $P$-wave peaks correspond to two poles of the
$\mathcal{U}$-matrix in the Sheet II. The masses and widths of the $S$-wave
state with $I(J^{P})=0(0^{+})$ and $P$-wave state with $0(1^{-})$ are
extracted from the pole positions $E=m-i\Gamma/2$. Their values are
insensitive to different cutoffs. We obtain

$\displaystyle[m,\Gamma]_{0^{+}}$ $\displaystyle\simeq$
$\displaystyle[2873.6^{+11.3}_{-13.5},72.1^{+10.5}_{-9.8}]~{}\rm{MeV},$
$\displaystyle{[m,\Gamma]_{1^{-}}}$ $\displaystyle\simeq$
$\displaystyle[2892.7^{+6.8}_{-10.2},54.2^{+8.2}_{-1.5}]~{}\rm{MeV}.$ (28)

Therefore, they correspond to the two bound states of $\bar{D}^{\ast}K^{\ast}$
in $S$\- and $P$-waves, respectively. We find the masses of the $0^{+}$ and
$1^{-}$ states are very consistent with the $X_{0}(2900)$ and $X_{1}(2900)$
Aaij  _et al._ (2020a, b), respectively, while their widths deviate a lot. The
analyses from the experimental work give $\Gamma_{X_{1}}>\Gamma_{X_{0}}$ Aaij
_et al._ (2020a, b), whereas this relation is reversed in our results. The
result in Ref. Huang  _et al._ (2020) agrees with ours.

Inspecting the $S$\- and $P$-wave peaks in Fig. 4, one may intuitively infer
that the width of $P$-wave state is larger than that of the $S$-wave. But as
we mentioned above, we do have $\Gamma_{0^{+}}>\Gamma_{1^{-}}$ from the pole
analyses. This is because the $P$-wave production vertex is momentum
dependent, which broadens the line shape, but the pole position cannot be
changed.

The experimental data may consist of multi substructures, similar to the
stories of the $P_{c}$ states [the $P_{c}(4450)$ observed in 2015 by the LHCb
Aaij  _et al._ (2015) was proved to contain two distinct states $P_{c}(4440)$
and $P_{c}(4457)$ in 2019 Aaij  _et al._ (2019)]. The $S$\- and $P$-wave
multiplets for the $\bar{D}^{\ast}K^{\ast}$ system are $[0^{+},1^{+},2^{+}]$
and $\\{1^{-},[0^{-},1^{-},2^{-}],[1^{-},2^{-},3^{-}]\\}$, respectively. The
present experimental data can be well reproduced with the inclusion of the
$0^{+}$ and $1^{-}$ states, but the other states in the $S$\- and $P$-wave
multiplets might also exist. In 2010, Molina et al predicted the existence of
the $S$-wave multiplets in $\bar{D}^{\ast}K^{\ast}$ system Molina  _et al._
(2010), in which the $0^{+}$ state mass is roughly consistent with our result.
The present resolution cannot distinguish the spin-spin and spin-orbital
interactions, thus more experimental data are still needed. However, one
should also note that these possibly existed substructures might be hard to be
discerned in experiments given that their widthes were similar to those of the
$X_{0,1}(2900)$ states.

##  IV Summary

We have investigated the $D^{-}K^{+}$ invariant mass distributions for the
decay process $B^{+}\to D^{+}D^{-}K^{+}$, where the $X_{0}(2900)$ and
$X_{1}(2900)$ were observed by the LHCb Collaboration Aaij  _et al._ (2020a,
b). We study the $S$\- and $P$-wave $\bar{D}^{\ast}K^{\ast}$ interactions in a
coupled-channel formalism via mimicking the $\chi$EFT. The short- and long-
distance forces are both incorporated in our calculations.

The event distributions can be well described, and two sharp peaks that
correspond to the $\bar{D}^{\ast}K^{\ast}$ molecular states with $0(0^{+})$
and $0(1^{-})$ are obtained. Their masses are in good agreement with the
$X_{0}(2900)$ and $X_{1}(2900)$, respectively, but the width of $0^{+}$ state
is larger than that of the $1^{-}$. We can simultaneously reproduce the
$X_{0}(2900)$ and $X_{1}(2900)$ in an unified framework. Our calculations
support the molecular interpretations for these two states. In other words,
the $X_{1}(2900)$ may be the $P$-wave excitation of the ground state
$X_{0}(2900)$. If this conclusion is verified in future, it shall be the first
time that a ground-state hadronic molecule and its orbital excitation are
synchronously observed in experiments.

However, the present experimental data cannot pin down the fine structures
behind the $\bar{D}^{\ast}K^{\ast}$ interactions, i.e., we still cannot
ascertain whether the other members in the $S$\- and $P$-wave multiplets
exist. The single whole peak now might contain multi substructures. In order
to resolve the structures hidden in the peak, some refined measurements in
this decay channel are still needed.

## Acknowledgments

B. Wang is very grateful to L. Meng for helpful discussions and carefully
reading the manuscript. This work is supported by the National Natural Science
Foundation of China under Grants No. 12105072, No. 11975033 and No.
12070131001. B.W is also supported by the Youth Funds of Hebei Province (No.
A2021201027) and the Start-up Funds for Young Talents of Hebei University (No.
521100221021).

## References

  * Aaij  _et al._ (2020a) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 125, 242001 (2020a), arXiv:2009.00025 [hep-ex] .
  * Aaij  _et al._ (2020b) R. Aaij _et al._ (LHCb), Phys. Rev. D 102, 112003 (2020b), arXiv:2009.00026 [hep-ex] .
  * Chen  _et al._ (2020) H.-X. Chen, W. Chen, R.-R. Dong, and N. Su, Chin. Phys. Lett. 37, 101201 (2020), arXiv:2008.07516 [hep-ph] .
  * He and Chen (2021) J. He and D.-Y. Chen, Chin. Phys. C 45, 063102 (2021), arXiv:2008.07782 [hep-ph] .
  * Liu  _et al._ (2020a) M.-Z. Liu, J.-J. Xie, and L.-S. Geng, Phys. Rev. D 102, 091502 (2020a), arXiv:2008.07389 [hep-ph] .
  * Hu  _et al._ (2021) M.-W. Hu, X.-Y. Lao, P. Ling, and Q. Wang, Chin. Phys. C 45, 021003 (2021), arXiv:2008.06894 [hep-ph] .
  * Agaev  _et al._ (2020) S. S. Agaev, K. Azizi, and H. Sundu, (2020), arXiv:2008.13027 [hep-ph] .
  * Karliner and Rosner (2020) M. Karliner and J. L. Rosner, Phys. Rev. D 102, 094016 (2020), arXiv:2008.05993 [hep-ph] .
  * He  _et al._ (2020) X.-G. He, W. Wang, and R. Zhu, Eur. Phys. J. C 80, 1026 (2020), arXiv:2008.07145 [hep-ph] .
  * Wang (2020) Z.-G. Wang, Int. J. Mod. Phys. A 35, 2050187 (2020), arXiv:2008.07833 [hep-ph] .
  * Zhang (2021) J.-R. Zhang, Phys. Rev. D 103, 054019 (2021), arXiv:2008.07295 [hep-ph] .
  * Wang  _et al._ (2021) G.-J. Wang, L. Meng, L.-Y. Xiao, M. Oka, and S.-L. Zhu, Eur. Phys. J. C 81, 188 (2021), arXiv:2010.09395 [hep-ph] .
  * Liu  _et al._ (2020b) X.-H. Liu, M.-J. Yan, H.-W. Ke, G. Li, and J.-J. Xie, Eur. Phys. J. C 80, 1178 (2020b), arXiv:2008.07190 [hep-ph] .
  * Burns and Swanson (2021a) T. J. Burns and E. S. Swanson, Phys. Lett. B 813, 136057 (2021a), arXiv:2008.12838 [hep-ph] .
  * Huang  _et al._ (2020) Y. Huang, J.-X. Lu, J.-J. Xie, and L.-S. Geng, Eur. Phys. J. C 80, 973 (2020), arXiv:2008.07959 [hep-ph] .
  * Chen  _et al._ (2021a) Y.-K. Chen, J.-J. Han, Q.-F. Lü, J.-P. Wang, and F.-S. Yu, Eur. Phys. J. C 81, 71 (2021a), arXiv:2009.01182 [hep-ph] .
  * Burns and Swanson (2021b) T. J. Burns and E. S. Swanson, Phys. Rev. D 103, 014004 (2021b), arXiv:2009.05352 [hep-ph] .
  * Xiao  _et al._ (2021) C.-J. Xiao, D.-Y. Chen, Y.-B. Dong, and G.-W. Meng, Phys. Rev. D 103, 034004 (2021), arXiv:2009.14538 [hep-ph] .
  * Albuquerque  _et al._ (2021) R. M. Albuquerque, S. Narison, D. Rabetiarivony, and G. Randriamanatrika, Nucl. Phys. A 1007, 122113 (2021), arXiv:2008.13463 [hep-ph] .
  * Lü  _et al._ (2020) Q.-F. Lü, D.-Y. Chen, and Y.-B. Dong, Phys. Rev. D 102, 074021 (2020), arXiv:2008.07340 [hep-ph] .
  * Mutuk (2021) H. Mutuk, J. Phys. G 48, 055007 (2021), arXiv:2009.02492 [hep-ph] .
  * Tan and Ping (2020) Y. Tan and J. Ping, (2020), arXiv:2010.04045 [hep-ph] .
  * Abreu (2021) L. M. Abreu, Phys. Rev. D 103, 036013 (2021), arXiv:2010.14955 [hep-ph] .
  * Qi  _et al._ (2021) J.-J. Qi, Z.-Y. Wang, Z.-F. Zhang, and X.-H. Guo, (2021), arXiv:2101.06688 [hep-ph] .
  * Chen (2021) H.-X. Chen, (2021), arXiv:2103.08586 [hep-ph] .
  * Hsiao and Yu (2021) Y.-K. Hsiao and Y. Yu, (2021), arXiv:2104.01296 [hep-ph] .
  * Duan  _et al._ (2021) M.-X. Duan, J.-Z. Wang, Y.-S. Li, and X. Liu, (2021), arXiv:2104.09132 [hep-ph] .
  * Kong  _et al._ (2021) S.-Y. Kong, J.-T. Zhu, D. Song, and J. He, (2021), arXiv:2106.07272 [hep-ph] .
  * Dong and Zou (2021) X.-K. Dong and B.-S. Zou, Eur. Phys. J. A 57, 139 (2021), arXiv:2009.11619 [hep-ph] .
  * Bondar and Milstein (2020) A. E. Bondar and A. I. Milstein, JHEP 12, 015 (2020), arXiv:2008.13337 [hep-ph] .
  * Chen  _et al._ (2021b) H. Chen, H.-R. Qi, and H.-Q. Zheng, Eur. Phys. J. C 81, 812 (2021b), arXiv:2108.02387 [hep-ph] .
  * Chen  _et al._ (2016) H.-X. Chen, W. Chen, X. Liu, and S.-L. Zhu, Phys. Rept. 639, 1 (2016), arXiv:1601.02092 [hep-ph] .
  * Guo  _et al._ (2018) F.-K. Guo, C. Hanhart, U.-G. Meißner, Q. Wang, Q. Zhao, and B.-S. Zou, Rev. Mod. Phys. 90, 015004 (2018), arXiv:1705.00141 [hep-ph] .
  * Liu  _et al._ (2019) Y.-R. Liu, H.-X. Chen, W. Chen, X. Liu, and S.-L. Zhu, Prog. Part. Nucl. Phys. 107, 237 (2019), arXiv:1903.11976 [hep-ph] .
  * Lebed  _et al._ (2017) R. F. Lebed, R. E. Mitchell, and E. S. Swanson, Prog. Part. Nucl. Phys. 93, 143 (2017), arXiv:1610.04528 [hep-ph] .
  * Esposito  _et al._ (2017) A. Esposito, A. Pilloni, and A. D. Polosa, Phys. Rept. 668, 1 (2017), arXiv:1611.07920 [hep-ph] .
  * Brambilla  _et al._ (2020) N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C.-P. Shen, C. E. Thomas, A. Vairo, and C.-Z. Yuan, Phys. Rept. 873, 1 (2020), arXiv:1907.07583 [hep-ex] .
  * Chen  _et al._ (2021c) S. Chen, Y. Li, W. Qian, Y. Xie, Z. Yang, L. Zhang, and Y. Zhang, (2021c), arXiv:2111.14360 [hep-ex] .
  * Chen  _et al._ (2022) H.-X. Chen, W. Chen, X. Liu, Y.-R. Liu, and S.-L. Zhu, (2022), arXiv:2204.02649 [hep-ph] .
  * Meng  _et al._ (2022) L. Meng, B. Wang, G.-J. Wang, and S.-L. Zhu, (2022), arXiv:2204.08716 [hep-ph] .
  * Weinberg (1990) S. Weinberg, Phys. Lett. B 251, 288 (1990).
  * Weinberg (1991) S. Weinberg, Nucl. Phys. B 363, 3 (1991).
  * Bernard  _et al._ (1995) V. Bernard, N. Kaiser, and U.-G. Meissner, Int. J. Mod. Phys. E 4, 193 (1995), arXiv:hep-ph/9501384 .
  * Epelbaum  _et al._ (2009) E. Epelbaum, H.-W. Hammer, and U.-G. Meissner, Rev. Mod. Phys. 81, 1773 (2009), arXiv:0811.1338 [nucl-th] .
  * Machleidt and Entem (2011) R. Machleidt and D. R. Entem, Phys. Rept. 503, 1 (2011), arXiv:1105.2919 [nucl-th] .
  * Meißner (2016) U.-G. Meißner, Phys. Scripta 91, 033005 (2016), arXiv:1510.03230 [nucl-th] .
  * Hammer  _et al._ (2020) H. W. Hammer, S. König, and U. van Kolck, Rev. Mod. Phys. 92, 025004 (2020), arXiv:1906.12122 [nucl-th] .
  * Rodriguez Entem  _et al._ (2020) D. Rodriguez Entem, R. Machleidt, and Y. Nosyk, Front. in Phys. 8, 57 (2020).
  * AlFiky  _et al._ (2006) M. T. AlFiky, F. Gabbiani, and A. A. Petrov, Phys. Lett. B 640, 238 (2006), arXiv:hep-ph/0506141 .
  * Fleming  _et al._ (2007) S. Fleming, M. Kusunoki, T. Mehen, and U. van Kolck, Phys. Rev. D 76, 034006 (2007), arXiv:hep-ph/0703168 .
  * Baru  _et al._ (2011) V. Baru, A. A. Filin, C. Hanhart, Y. S. Kalashnikova, A. E. Kudryavtsev, and A. V. Nefediev, Phys. Rev. D 84, 074029 (2011), arXiv:1108.5644 [hep-ph] .
  * Valderrama (2012) M. P. Valderrama, Phys. Rev. D 85, 114037 (2012), arXiv:1204.2400 [hep-ph] .
  * Braaten (2015) E. Braaten, Phys. Rev. D 91, 114007 (2015), arXiv:1503.04791 [hep-ph] .
  * Schmidt  _et al._ (2018) M. Schmidt, M. Jansen, and H. W. Hammer, Phys. Rev. D 98, 014032 (2018), arXiv:1804.00375 [hep-ph] .
  * Wang  _et al._ (2019a) B. Wang, Z.-W. Liu, and X. Liu, Phys. Rev. D 99, 036007 (2019a), arXiv:1812.04457 [hep-ph] .
  * Meng  _et al._ (2019) L. Meng, B. Wang, G.-J. Wang, and S.-L. Zhu, Phys. Rev. D 100, 014031 (2019), arXiv:1905.04113 [hep-ph] .
  * Wang  _et al._ (2019b) B. Wang, L. Meng, and S.-L. Zhu, JHEP 11, 108 (2019b), arXiv:1909.13054 [hep-ph] .
  * Wang  _et al._ (2020) B. Wang, L. Meng, and S.-L. Zhu, Phys. Rev. D 102, 114019 (2020), arXiv:2009.01980 [hep-ph] .
  * Meng  _et al._ (2021) L. Meng, B. Wang, and S.-L. Zhu, Sci. Bull. 66, 1413 (2021), arXiv:2012.09813 [hep-ph] .
  * Chen  _et al._ (2021d) K. Chen, B. Wang, and S.-L. Zhu, Phys. Rev. D 103, 116017 (2021d), arXiv:2102.05868 [hep-ph] .
  * Wise (1992) M. B. Wise, Phys. Rev. D 45, R2188 (1992).
  * Manohar and Wise (2000) A. V. Manohar and M. B. Wise,  _Heavy quark physics_ , Vol. 10 (2000).
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Du  _et al._ (2020) M.-L. Du, V. Baru, F.-K. Guo, C. Hanhart, U.-G. Meißner, J. A. Oller, and Q. Wang, Phys. Rev. Lett. 124, 072001 (2020), arXiv:1910.11846 [hep-ph] .
  * Wang  _et al._ (2018) Q. Wang, V. Baru, A. A. Filin, C. Hanhart, A. V. Nefediev, and J. L. Wynen, Phys. Rev. D 98, 074023 (2018), arXiv:1805.07453 [hep-ph] .
  * Baru  _et al._ (2019) V. Baru, E. Epelbaum, A. A. Filin, C. Hanhart, A. V. Nefediev, and Q. Wang, Phys. Rev. D 99, 094013 (2019), arXiv:1901.10319 [hep-ph] .
  * Golak  _et al._ (2010) J. Golak _et al._ , Eur. Phys. J. A 43, 241 (2010), arXiv:0911.4173 [nucl-th] .
  * Liu  _et al._ (2009) Y.-R. Liu, X. Liu, and S.-L. Zhu, Phys. Rev. D 79, 094026 (2009), arXiv:0904.1770 [hep-ph] .
  * Back  _et al._ (2018) J. Back _et al._ , Comput. Phys. Commun. 231, 198 (2018), arXiv:1711.09854 [hep-ex] .
  * Nogga  _et al._ (2005) A. Nogga, R. G. E. Timmermans, and U. van Kolck, Phys. Rev. C 72, 054006 (2005), arXiv:nucl-th/0506005 .
  * Birse  _et al._ (1999) M. C. Birse, J. A. McGovern, and K. G. Richardson, Phys. Lett. B 464, 169 (1999), arXiv:hep-ph/9807302 .
  * Frazer and Hendry (1964) W. R. Frazer and A. W. Hendry, Phys. Rev. 134, B1307 (1964).
  * Eden and Taylor (1964) R. J. Eden and J. R. Taylor, Phys. Rev. 133, B1575 (1964).
  * Badalian  _et al._ (1982) A. M. Badalian, L. P. Kok, M. I. Polikarpov, and Y. A. Simonov, Phys. Rept. 82, 31 (1982).
  * Aaij  _et al._ (2015) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 115, 072001 (2015), arXiv:1507.03414 [hep-ex] .
  * Aaij  _et al._ (2019) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 122, 222001 (2019), arXiv:1904.03947 [hep-ex] .
  * Molina  _et al._ (2010) R. Molina, T. Branz, and E. Oset, Phys. Rev. D 82, 014010 (2010), arXiv:1005.0335 [hep-ph] .


# Low incidence rate of COVID-19 undermines confidence in estimation of the
vaccine efficacy

Yasin Memari<EMAIL_ADDRESS>MRC Cancer Unit, University of Cambridge,
Cambridge CB2 0XZ, UK

## Abstract

Knowing the true effect size of clinical interventions in randomised clinical
trials is key to informing the public health policies. Vaccine efficacy is
defined in terms of the relative risk or the ratio of two disease risks.
However, only approximate methods are available for estimating the variance of
the relative risk. In this article, we show using a probabilistic model that
uncertainty in the efficacy rate could be underestimated when the disease risk
is low. Factoring in the baseline rate of the disease, we estimate broader
confidence intervals for the efficacy rates of the vaccines recently developed
for COVID-19. We propose new confidence intervals for the relative risk. We
further show that sample sizes required for phase 3 efficacy trials are
routinely underestimated and propose a new method for sample size calculation
where the efficacy is of interest. We also discuss the deleterious effects of
classification bias which is particularly relevant at low disease prevalence.

## Introduction

Vaccines are seen as the best control measure for the coronavirus pandemic. In
this context, understanding the true efficacy of the vaccines and clinical
interventions is crucial. Randomised clinical trials are conducted to
systematically study the safety and the efficacy of an intervention in a
subset of the population before it is widely used in the general population.
In placebo-controlled vaccine trials, participants are randomised into
vaccinated and unvaccinated groups where cases of the disease or infection are
allowed to accrue over time. In planning a clinical trial, advance sample size
calculation determines the size of the trial population needed to detect a
minimal clinically relevant difference between the two groups if such a
difference exists. The indicator for effectiveness of a vaccine is usually
reduction of the cases in the vaccinated group relative to the control group.
However, it is sometimes naively assumed that the trial participants who do
not experience the event provide no information. Consequently, the event rate
or the incidence rate of the disease receives inadequate attention. For rare
diseases, it is often simply accepted that the accrual of the cases takes
longer. Human clinical trials are also an area where theory and practice are
seldom consistent, as experiments in human populations are hardly fully
controlled experiments, not least due to unrealistic assumptions, loss to
follow-up, noncompliance, heterogeneity of treatment effect and the trial
population, etc. [1] Therefore, it is not uncommon that, by the time an
interim analysis declares a significant finding, the original assumptions used
to define the statistical power of the study and the sample size are
neglected.

In this article our interest is on evaluating the impact of the event rate,
insofar as it could affect the estimation of the efficacy rate. We show that
low incidence rate of the disease could lead to overestimation of confidence
in the estimated efficacy rates. We propose a new method for posterior
probability of the vaccine efficacy that has a more subtle relationship with
the event rate. Using our approach, we obtain broader confidence intervals for
the efficacy of the vaccines recently developed for COVID-19. Based on our
findings, we propose new confidence intervals for the relative risk. A new
method for sample size calculation in controlled efficacy trials is proposed
which is more robust at low disease prevalence. Also highlighted is the impact
of classification bias which could have large consequences when the disease
risk is low.

## Methods

Vaccine efficacy is defined as the proportionate reduction in the risk of
disease or infection in a vaccinated group compared to an unvaccinated group.
It is defined as (1-RR)$\times$100% in terms of the relative risk or the risk
ratio, $\textrm{RR}=\pi_{v}/\pi_{c}$, where $\pi$ are the incidence of the
disease among those exposed in the vaccinated and control groups. Throughout
this paper we interchangeably use the terms, incidence rate, disease risk,
prevalence and event rate.

It is important to remember that the variables $\pi_{v}$ and $\pi_{c}$ are
scaled binomials as they represent sample proportions. Assuming equal person-
time exposure in the two groups, the efficacy is often summarised in terms of
the numbers of cases in the vaccinated and unvaccinated groups, $t_{v}$ and
$t_{c}$ respectively:

$\alpha=1-\frac{\pi_{v}}{\pi_{c}}\simeq 1-\frac{t_{v}}{t_{c}}.$ (1)

It appears in the literature that only approximate methods are available for
the variance of the ratio of two binomial parameters [2, 3]. The consensus
method that is commonly used to assign confidence intervals to the risk ratio
is credited to Katz et al [3]. The method is based on asymptotic normality of
logarithm of the ratio of two binomial variables. Assuming independence of the
incidence rates, it follows that $\textrm{var}(\log(\pi_{v}/\pi_{c}))$ =
$\textrm{var}(\log(\pi_{v}))+\textrm{var}(\log(\pi_{c}))$. Using a Taylor
series, the variances are approximated as $\textrm{var}(\log(\pi))\approx$
$\textrm{var}(\pi)/\pi^{2}$ where Wald method is often used to set
$\textrm{var}(\pi)$. Then two-sided 95% confidence intervals on the efficacy
(e.g. see [4, 5, 6, 7]) can be written as

$95\%\textrm{CL}:1-\exp\bigg{(}\ln(\textrm{RR})\pm
1.96\sqrt{\frac{1-\pi_{v}}{t_{v}}+\frac{1-\pi_{c}}{t_{c}}}\bigg{)}.$ (2)

Hereafter we refer to equation 2 as pooled Wald approximation. We will show
that the method underestimates the variance espcially when the incidence rate
is low.

Equation 2 sets out the large sample asymptotic variance of the risk ratio.
However, Wald method used to define $\textrm{var}(\pi)$ is known to be
unreliable when $\pi$ is small. One may use alternative binomial proportion
confidence intervals, however, log normality of the ratio might not hold and
the variance of (the logarithm of) the ratio may be irreducible. Hightower et
al [5] raised question about the credibility of the confidence limits when the
efficacy is high and the disease risk is low. Also, O’Neill [8] noted that,
when $t\ll n$, the variance of $\ln({\textrm{RR}})$ in equation 2 remains
fairly stable and quickly converges to $1/t_{v}+1/t_{c}$.

Ratio distributions are known to have heavy tails and often no finite
variance. If one were to model the likelihood function for the efficacy
defined in equation 1 in terms of independent incidence rates, the choice of
the prior probabilities for $\pi_{v}$ and $\pi_{c}$ would be critical. One can
readily verify that the variance of the ratio of two binomial distributions
increases as binomial probabilities decrease. Uninformative priors could
simply cancel out by the division and the dependence of the posterior on the
prevalence would not become obvious. Analytical solutions using independent
incidence rates may also be hard to obtain.

For an analytical solution, we model the efficacy in terms of conditional
probabilities of the disease risks. Independence of the probabilities of the
incidence rates is neither necessary nor ideal when calculating the efficacy,
as equation 1 imposes a constraint on the two variables. Under a binomial
model with overall prevalence of $\pi=t/n$ in both groups and total population
size of $n$, overall number of cases $t=t_{c}+t_{v}$ follows
$t\sim\mathrm{Bin}(n,\pi)$, then, from equation 1 assuming
$t_{c}\sim\mathrm{Bin}(t,1/(2-\alpha))$, we expect
$t_{c}\sim\mathrm{Bin}(n,\pi/(2-\alpha))$. Were we to use Poisson
distributions for $t$ and $t_{c}$, $t_{c}$ conditional on $t$ would still
follow a binomial distribution. Modeling the efficacy in terms of conditional
probabilities has previously been suggested [4]. This notation enables to
explicitly parametrise the likelihood function in terms of the prevalence,
irrespective of the priors for $\pi_{v}$ and $\pi_{c}$.

For a general solution accounting for classification bias we assume an
imperfect diagnostic procedure with sensitivity Se and specificity Sp. Then
fraction of individuals who test positive for the disease is sum of true
positive rate and false positive rate:

$\displaystyle T$
$\displaystyle=\textrm{Se}\times\pi+(1-\textrm{Sp})\times(1-\pi)$
$\displaystyle=c_{1}+c_{2}\pi,$ (3)

where $c_{1}$=1-Sp is the false positive rate and $c_{2}$=Se+Sp-1. The
posterior distribution of $\alpha$ given that $t_{c}$ is binomial follows as

$\displaystyle p(\alpha|t_{c},\pi,c_{1},c_{2})$
$\displaystyle=\frac{p(t_{c}|\alpha,\pi,c_{1},c_{2})p(\pi)p(\alpha)}{g(\alpha)}$
$\displaystyle\propto\frac{1}{g(\alpha)}\binom{n}{t_{c}}\left(\frac{c_{1}+c_{2}\pi}{2-\alpha}\right)^{t_{c}}\left(1-\frac{c_{1}+c_{2}\pi}{2-\alpha}\right)^{n-t_{c}}f(\pi),$
(4)

where $f(\pi)$ is the prior on $\pi$ and we have assumed uniform prior on the
efficacy $\alpha\sim\mathrm{unif}\\{0,1\\}$. For a complete solution, the
marginal likelihood $g(\alpha)$ can be written in terms of the incomplete beta
function (see e.g. [9]):

$\displaystyle g(\alpha)$
$\displaystyle=f(\pi)\binom{n}{t_{c}}(c_{1}+c_{2}\pi)\big{\\{}B(c_{1}+c_{2}\pi;t_{c}-1,n-t_{c}+1)$
$\displaystyle-B((c_{1}+c_{2}\pi)/2;t_{c}-1,n-t_{c}+1)\big{\\}}.$

As we do not intend to impose a prior on the prevalence, $f(\pi)$ in equation
4 cancels out and our analysis, in essence, is likelihood based. One needs to
remember that, the posterior in equation 4, as it was derived from the second
equality in equation 1, is valid only when the individuals are equally divided
between the two groups.

The mode of the posterior of $\alpha$ is obtained by setting the derivative of
the log likelihood to zero i.e.
$\partial\ell/\partial\alpha=\partial\mathrm{ln}(p(\alpha|t_{c},\pi))/\partial\alpha=0$.
This leads to

$\alpha_{mode}=2-\frac{n(c_{1}+c_{2}\pi)}{t_{c}},$ (5)

which corresponds to the maximum likelihood estimator (MLE). Cramér–Rao bound
expresses a lower bound on the variance of any unbiased estimator of $\alpha$
in terms of the inverse of the Fisher information

$\textrm{Var}(\alpha_{mode})\geq\frac{1}{\mathcal{I}(\alpha)},$ (6)

where the Fisher information $\mathcal{I(\alpha)}$ is obtained as

$\displaystyle\mathcal{I}(\alpha)=\mathbb{E}\Big{[}\Big{(}\frac{\partial}{\partial\alpha}\ell(\alpha|t_{c},\pi)\Big{)}^{2}\Big{]}$
$\displaystyle=n\times\mathbb{E}\Big{[}\Big{(}-\frac{1-t_{c}}{2-\alpha-\pi}-\frac{-1}{2-\alpha}\Big{)}^{2}\Big{]}$
$\displaystyle=\frac{n(c_{1}+c_{2}\pi)}{(2-\alpha)^{2}(2-\alpha-(c_{1}+c_{2}\pi))}.$
(7)

Here $\mathbb{E}$ denotes ‘expected’ over $t_{c}$, where we have substituted
$\mathbb{E}[t_{c}]=\mathbb{E}[t_{c}^{2}]=(c_{1}+c_{2}\pi)/(2-\alpha)$. We will
show that the conditional binomial model has a more subtle dependence on $\pi$
compared to the pooled Wald method.

Under certain regularity conditions and assuming asymptotic normality near
MLE, 95% confidence intervals on $\alpha_{mode}$ can be estimated as

$\alpha_{mode}\pm\frac{1.96}{\sqrt{\mathcal{I}(\alpha_{mode})}}.$ (8)

However, as the posterior distribution is asymmetric, especially when the
efficacy is high, and the intervals could lie outside [0,1], we will estimate
the credible intervals computationally.

## Results

### Effect of incidence rate on vaccine efficacy

The posterior probability of vaccine efficacy given in its simplest form in
equation 4 is ready for inspection. Using binomial notation is particularly
useful in enabling us to directly plug in the numbers $n$, $t_{c}$ in the
estimation of $\alpha$. In this section we evaluate the impact of the
incidence rate on the efficacy and assign new confidence bounds to the
efficacy of COVID-19 vaccines.

Firstly, we assume a diagnostic test with perfect sensitivity and specificity
i.e. Se=Sp=1. In the absence of misclassification, mode of the posterior in
equation 5 corresponds to the expectation $\hat{\alpha}=1-t_{v}/t_{c}$. The
larger $n$ the smaller the variance of the posterior, however, for a fixed
$n$, the variance depends on $\pi$. Figure 1 shows the posterior probability
of $\alpha$ plotted over a range of $\pi$, for a fixed $n$ on the left hand,
and for a fixed $t$ on the right hand, assuming true vaccine efficacy of 70%
and 90% respectively. Also plotted in vertical lines are the independent 95%
confidence intervals from equation 2. As the event rate falls, the posterior
distributions and the confidence intervals become wider, however, for a fixed
$t$ (right plot) Wald intervals are stable over a wide range of $\pi$, and
more so when the efficacy is high. The proposed conditional binomial model
better represents the variability at low prevalence.

Figure 1:  Posterior distribution of vaccine efficacy. Blue lines represent
the normalised posterior probabilities, while vertical lines show the
independent pooled Wald confidence intervals. Left hand plot assumes a fixed
$n$=50,000 while right hand plot is for a fixed $t$=2,000. The general trend
holds for different values of the parameters. Wald method overstates the
confidence in the efficacy when $t\ll n$.

Three clinical trials of the vaccines designed to prevent COVID-19 recently
published their interim phase 3 analysis results [10, 11, 12] with two of them
reporting incredibly narrow 95% confidence bounds on the efficacy. The
reported case numbers and the efficacy rates for the primary end points are
provided in Table 1. Firstly, we note that, although the trials used different
models and priors on the efficacy, the reported confidence intervals almost
perfectly correspond with those obtained from equation 2. At large $n$ the
posterior is clearly dominated by the data and the Bayesian and the
frequentist are equivalent. Furthermore, especially where the efficacy is
high, pooled Wald confidence intervals hardly vary by the choice of $n$. If
one were to use different values for $n$ in Table 1, over a large range of the
values equation 2 would still give the same confidence intervals. Therefore,
the uncertainty caused by $n_{v}$ and $n_{c}$ is not accounted for.

Table 1: Estimated efficacy of COVID-19 vaccine trials Trial | Case numbers and reported efficacy rates | Estimated efficacy rate
---|---|---
case rate in vaccinated |  case rate in control |  reported efficacy and 95% CI |  estimated mode and 95% credible interval
AZ-Oxford (combined) | 30/5,807 | 101/5,829 | 70·4% [54·8, 80·6] | 70.3% [39.1, 90.9]
Pfizer-BioNTech | 8/18,198 | 162/18,325 | 95.0% [90.3, 97.6] | 95.1% [74.9, 99.6]
Moderna-NIH | 11/14,134 | 185/14,073 | 94.1% [89.3, 96.8] | 94.1% [75.4, 99.5]

We re-estimate the confidence intervals using the conditional binomial model
presented in the Methods. Using the case numbers reported, the likelihood of
the data in equation 4 is obtained by setting the prevalence to $\pi=T=t/n$.
Then maximum a posteriori (MAP) and 95% credible intervals for the efficacy
rates are calculated computationally. The results shown in Table 1 are
contrasted with those reported. Although estimated modes are the same, our
credible intervals are wider. Incorporating the incidence rates has removed
the overwhelming confidence originally assigned to the point estimates. Note
that, our approach requires the trial participants to be equally divided
between the vaccinated and unvaccinated groups which is roughly the case here.

Figure 2, in red, shows the posterior probabilities and the credible intervals
for COVID-19 vaccines. Of note is that, if we were to hypothetically assume
$\pi=t/n=1$, the posterior in equation 4 would produce the same intervals as
those reported by the vaccine trials and Wald approximation. Moreover, an
independent binomial model with uninformative (e.g. uniform) priors for
$\pi_{v}$ and $\pi_{c}$ would produce the pooled Wald intervals.

Figure 2:  Estimated efficacy of COVID-19 vaccines. Posterior probabilities
for the conditional binomial model are plotted in red, with shaded areas
representing the 95% credible intervals. Blue curves are for when $n$ is set
to $t_{v}+t_{c}$ and correspond with pooled Wald approximation.

### Bias in case classification

So far we have assumed no bias in classification of the cases, however,
imperfect diagnostic procedure could lead to misclassification of the infected
and uninfected individuals. In this section we examine the effect of
classification bias on estimation of the efficacy.

It is worth noting that equation 3 requires the observed infection rate $T$ to
be greater than the false positive rate $c_{1}=1-Sp$. This relates to the
‘false positive paradox’ which implies that the accuracy of a diagnostic test
is compromised if the test is used in a population where the incidence of the
disease is lower than the false positive rate of the test itself. Furthermore,
false negatives could dominate at low incidence rates. When the disease risk
is low, as the majority of the tests are negative, a small false negative rate
could lead to a situation where false negatives outnumber the positive cases.
These concepts are further explained in Note 1.

Figure 3:  Effect of imperfect diagnostic procedure. Misclassification error
biases the vaccine efficacy rate. Left plot shows the distributions for Se=1
and Sp=0.999, while the right plot is for Se=0.95 and Sp=1, with $n$=50,000 in
both. True efficacy rate is assumed at 70%. Imperfect specificity, however
small, could have disastrous effects when incidence rate is low, whereas lack
of sensitivity consistently inflates the efficacy rate.

Figure 3 illustrates the effect of classification bias on the posterior
probability of the vaccine efficacy. The left plot shows the impact of a very
small reduction in specificity to 0.999 (or increase in false positive rate),
while the right hand plot shows the effect of reduction in sensitivity to 0.95
(or increase in false negative rate). A small loss of specificity could lead
to serious underestimation of the effect size as noted by [6, 7], but it could
further lead to complete loss of precision when the incidence rate is low.
Loss of sensitivity results in overestimation of the efficacy irrespective of
the disease rate. In these plots, we have considered a larger reduction in
sensitivity, not only because reduction in specificity has a more dramatic
effect, but also as diagnostic assays typically have relatively higher
specificity than sensitivity, not least due to specimen collection,
insufficient viral load, stage of the disease, etc. [13] However, the effect
of loss of sensitivity is consistently toward shifting the mode in equation 5,
or MAP, to higher values of $\alpha$, even at low incidence rates where
negative predictive value is high.

[htb]  Note 1 J. Balayla [14] noted that there exists a prevalence threshold
below which the positive predictive value (PPV) of a diagnostic test drops
precipitously relative to the prevalence. This means that at too low a
prevalence a positive test result could more likely be a false positive than a
true positive. More underappreciated is the impact of the negative predictive
value (NPV). Though, at low incidence rates, the negative predictive value is
nearly 100%, a small loss in sensitivity could still have a marked effect as
the negative tests vastly outnumber the positive tests. We could even have a
situation where the false negatives are more than the true and false
positives. To avoid these pitfalls, the participants are pre-selected for
their symptoms before confirmation with the assay. Though this raises the pre-
test probability, it could cause collider bias [15]. Figure 4: Positive (red)
and negative (blue) predictive values are plotted in terms of population
prevalence. Solid lines are for to a diagnostic test with Se=Sp=0.99; dashed
lines are for Se=Sp=0.95. Vertical lines show the prevalence thresholds.

## Discussion

Base rate fallacy happens in situations where base rate information is ignored
in favour of individuating information. In probability terms, it often occurs
when $P(A|B)$ is confused or interchangeably used with $P(B|A)$ ignoring the
prior probability $P(A)$, e.g. probability of having a rare disease given a
positive test is wrongly equated to probability of a positive test given the
disease (or diagnostic sensitivity) ignoring the low prior probability of the
disease itself. We showed, in estimation of the vaccine efficacy when the
disease rate is low, not only diagnostic error could have deleterious effects,
but also failure to appropriately integrate the information about the base
rate or incidence rate of the disease in the calculation could lead to
underestimation of the uncertainty.

Vaccine efficacy is defined in terms of the risk ratio $\pi_{v}/\pi_{c}$, that
is the ratio of two binomial proportions. Ratio distributions are known to
have undefined variances, conversely, pooled Wald method has been
traditionally used to approximate the variance of the risk ratio. In this
article, we used a parametrisation that makes the dependence of the efficacy
on the disease prevalence explicit, without recourse to priors for $\pi_{v}$
and $\pi_{c}$. Particularly, improper priors for $\pi_{v}$ and $\pi_{c}$ could
lead to underestimation of the variance. We conditioned $t_{c}$ on
$t=t_{c}+t_{v}$ and treated $t$ as another random variable. The resulting
compound probability $t_{c}\sim\mathrm{Bin}(n,\pi/(2-\alpha))$ is over-
dispersed and better captures the variability of the variance with $\pi$,
whereas pooled Wald confidence intervals are largely insensitive to $\pi$ when
$\pi$ is small.

Wald method is intended as large sample approximation, however, the bulk of
the life sciences deals with small sample sizes. Therefore, it is likely that
the confidence intervals reported in the literature for the risk ratio (and
odds ratio) are overly optimistic. By analogy of equations 6 and 8, one could
define new confidence intervals for the risk ratio by substituting
RR=$(n_{c}/n_{v})(1-\alpha)$ for unequal sized groups in the Fisher
information. The results can be written as

$95\%\textrm{CL}:\textrm{RR}\pm
1.96\frac{n_{c}}{n_{v}}\big{(}1+\frac{t_{v}}{t_{c}}\big{)}\sqrt{\frac{1+t_{v}/t_{c}-\pi}{t_{v}+t_{c}}},$
(9)

where $\pi$=$(t_{v}+t_{c})/(n_{v}+n_{c})$. The above intervals on the risk
ratio are generally wider than but converge to the pooled Wald method when the
sample size is large. They may be preferred to those obtained from equation 2
when the sample size is small or the relative risk is low. Particularly, for a
fixed sample size as RR nears zero, the upper bound in equation 9 remains
conservative and the lower bound takes negative values and becomes
undetermined. On the contrary, as RR nears zero, the pooled Wald intervals
remain positive and shrink rapidly, giving the counterintuitive impression of
increased precision when the incidence rate is low (similar to figure 2).
However, as with Wald method, the confidence intervals in equation 9 were
derived using normal approximation which may not hold when RR significantly
deviates from 1.

Our findings have implications for pre-planning the sample sizes for phase 3
efficacy trials. Sample size calculation in case-control design is often
stated as “How many samples are needed to be randomised in order to conclude
with 100$(1-\beta)$% power that a treatment difference of size $\Delta$ exists
between the two groups at the level of significance of $\alpha$?”. Therefore
calculation of sample size requires specification of the null hypothesis
(expected treatment effect) and the alternative hypothesis defined in terms of
the difference in treatment outcomes. Here, $\alpha$ or type I error is the
probability of rejecting the null hypothesis where we should not, and $\beta$
or type II error is the probability of failing to reject the null hypothesis
where we should reject it. Under the assumption of normality of the treatment
outcome, a generic formula for per-group sample size is derived in terms of
the two-sample t-test: [1]

$n=\frac{2\sigma^{2}}{\Delta^{2}}(z_{1-\alpha/2}+z_{1-\beta})^{2},$ (10)

where $z$-scores determine the critical values for the standard normal
distribution. Therefore one needs to specify the variance of the measured
variable, the desired rates of error and the magnitude of the treatment
difference. Where the measured variable is binary (infected or uninfected),
the test statistic reduces to the test for the difference between two
proportions. Where the efficacy is of interest, the log normal approximation
of the risk ratio from equation 2 may be used to define the test statistic.
O’Neill [8] calculated the required sample sizes for a two-sided test given
the pooled Wald variance in equation 2. We re-write the total sample size in
this form:

$n=2\frac{(z_{1-\alpha/2}+z_{1-\beta})^{2}}{d^{2}}\Big{(}\frac{(2-\textrm{VE})^{2}}{\pi(1-\textrm{VE})}-2\Big{)},$
(11)

where

$d=\ln\Big{(}\Delta/(2(1-\textrm{VE}))+\sqrt{\big{(}\Delta/(2(1-\textrm{VE}))\big{)}^{2}+1}\Big{)}.$

Here VE is the anticipated efficacy and $\Delta$ is the expected difference in
VE in absolute terms. We showed, however, that at low prevalence rate,
equation 2 significantly underestimates the variance. Using an inadequately
small variance could lead to underestimation of the type I and type II errors,
potentially resulting in winner’s curse in underpowered studies [16, 17]. If
instead we were to use the proposed compound binomial model, one could simply
substitute the variance in equation 6. As in [8], under the assumption of
normality and assuming $\Delta$ is the difference between the upper and lower
limits of the confidence interval, substituting the margin of error as
$\Delta/2=z\sigma$ in equation 6 gives

$n\geq
4\frac{(z_{1-\alpha/2}+z_{1-\beta})^{2}}{\pi\Delta^{2}}(2-\textrm{VE})^{2}(2-\textrm{VE}-\pi).$
(12)

This equation sets out the total required sample size for a perfect diagnostic
test, to be equally divided between the two groups.

The proposed Cramér–Rao bound based formula 12 assumes normality of
distributions of the null and the alternative hypotheses, however, the
binomial likelihood function is asymmetric, as is pooled Wald intervals (see
[8]), and becomes more so as the efficacy increases. Notwithstanding the
limitations, we plug in the critical values for $\alpha=0.05$ and power of
$100(1-\beta)=80$ per cent ($z_{1-\alpha/2}=1.96$ and $z_{1-\beta}=0.84$) in
equations 11 and 12. The resulting sample sizes are plotted in Figure 5 for
$\Delta=10\%$ and different prevalence and efficacy rates.

Figure 5:  Sample size relative to disease prevalence. Total number of samples required to detect with 80% power and level of significance of $\alpha=0.05$ a difference in the efficacy of size $\Delta=10\%$. Solid lines represent Cramér–Rao bound and dashed lines represent pooled Wald approximation. On the left, x-axis is on logarithmic scale. y-axis is logarithmic in both plots. Table 2: Total sample sizes needed to conclude with 80% power and $\alpha$=0.05 a significant effect size effect size | event rate
---|---
VE |  $\Delta$ | 0.5 | 0.1 | 0.05 | 0.01 | 0.005 | 0.001 | 0.0005
0% | 10% | 37,632 | 238,336 | 489,216 | 2,496,256 | 5,005,056 | 25,075,456 | 50,163,456
0% | 20% | 9,408 | 59,584 | 122,304 | 624,064 | 1,251,264 | 6,268,864 | 12,540,864
0% | 30% | 4,181 | 26,482 | 54,357 | 277,362 | 556,117 | 2,786,162 | 5,573,717
0% | 40% | 2,352 | 14,896 | 30,576 | 156,016 | 312,816 | 1,567,216 | 3,135,216
30% | 10% | 21,751 | 145,009 | 299,080 | 1,531,654 | 3,072,371 | 15,398,105 | 30,805,273
30% | 20% | 5,438 | 36,252 | 74,770 | 382,913 | 768,093 | 3,849,526 | 7,701,318
30% | 30% | 2,417 | 16,112 | 33,231 | 170,184 | 341,375 | 1,710,901 | 3,422,808
30% | 40% | 1,359 | 9,063 | 18,693 | 95,728 | 192,023 | 962,382 | 1,925,330
60% | 10% | 11,064 | 79,905 | 165,957 | 854,372 | 1,714,890 | 8,599,037 | 17,204,221
60% | 20% | 2,766 | 19,976 | 41,489 | 213,593 | 428,723 | 2,149,759 | 4,301,055
60% | 30% | 1,229 | 8,878 | 18,440 | 94,930 | 190,543 | 955,449 | 1,911,580
60% | 40% | 691 | 4,994 | 10,372 | 53,398 | 107,181 | 537,440 | 1,075,264
90% | 10% | 4,553 | 37,946 | 79,686 | 413,607 | 831,009 | 4,170,221 | 8,344,237
90% | 20% | 1,138 | 9,486 | 19,921 | 103,402 | 207,752 | 1,042,555 | 2,086,059
90% | 30% | 506 | 4,216 | 8,854 | 45,956 | 92,334 | 463,358 | 927,137
90% | 40% | 285 | 2,372 | 4,980 | 25,850 | 51,938 | 260,639 | 521,515

In Figure 5 the relationship between the sample size and the incidence rate
looks linear on log-log scale as they have a power law relationship. However,
while the two methods coincide at high incidence rates, pooled Wald method
significantly underestimates the sample sizes at low incidence rates
especially when the efficacy is high (note that y-axis is on logarithmic
scale). Contrasting Figure 5 with the case rates in Table 1, it is clear that,
to achieve the narrow confidence bounds that Pfizer and Moderna have reported,
they would have needed several times more samples under pooled Wald method,
and an order of magnitude more under Cramér–Rao bound. If the event rate were
to differ from that in the general population or if possibility of
misclassification was non negligible, such a discrepancy in incidence rates
could cause such large variations in the variance that the trial population
could be unrepresentative of the larger population. Table 2 provides the total
sample sizes from Cramér–Rao bound formula 12 for different levels of efficacy
and effect size. It is clear that the sample size is also very sensitive to
the choice of $\Delta$, therefore an investigator must be wary of
misspecification of the anticipated treatment difference [1].

Throughout the Methods, we incorporated the misclassification error in the
calculations in order to emphasise the importance of accounting for
classification bias when the disease is rare. We showed that, while lack of
diagnostic sensitivity consistently inflates the estimated efficacy rates,
imperfect specificity results is serious loss of accuracy and precision at low
disease risks. Case definition for COVID-19 is particularly a major caveat.
The three vaccine trials broadly follow FDA definition of the disease. For
primary end points symptomatic cases are identified by surveillance or are
self-reported, and are subsequently confirmed with RT-PCR. Pre-selecting of
the participants for PCR assay could create the possibility for collider bias
[15]. Moreover, the highly non-specific symptoms of COVID-19, which include
symptoms as common as cough and congestion, could create the perfect
conditions for misclassification. False negatives due to e.g. selective
reporting, specimen collection, etc, and PCR false positives due to e.g.
remnant viral RNA, etc could be introduced if the test is not repeated [13,
18]. Much remains unknown about COVID-19 and its many symptoms and
presentations. Therefore, it is recommended to account for classification bias
in the calculation. The code for calculating the posterior probability of the
vaccine efficacy, which can simultaneously marginalise over the diagnostic
sensitivity and specificity is provided.

## Code

R code for the posterior probability of the efficacy was modified from code
published in [9]. It is provided in Appendix along with functions to calculate
the sample sizes from equations 11 and 12.

## Acknowledgments

The author’s position at the University of Cambridge is funded by CRUK grant
C60100/A23916. The author would like to appreciate the helpful comments
received from the Cancer Mutagenesis group at MRC Cancer Unit.

## References

  * 1. J. Wittes.  Sample Size Calculations for Randomized Controlled Trials.  Epidemiologic Reviews, 24(1):39–53, 07 2002.
  * 2. J. J. Gart and J. Nam.  Approximate interval estimation of the ratio of binomial parameters: a review and corrections for skewness.  Biometrics, 44(2):323–338, Jun 1988.
  * 3. D. Katz et al.  Obtaining confidence intervals for the risk ratio in cohort studies.  Biometrics, 34:469–474, 1978.
  * 4. M. Ewell.  Comparing methods for calculating confidence intervals for vaccine efficacy.  Stat Med, 15(21-22):2379–2392, 1996.
  * 5. A. W. Hightower et al.  Recommendations for the use of Taylor series confidence intervals for estimates of vaccine efficacy.  Bull World Health Organ, 66(1):99–105, 1988.
  * 6. P. A. Lachenbruch.  Sensitivity, specificity, and vaccine efficacy.  Controlled Clinical Trials, 19(6):569 – 574, 1998.
  * 7. A. Hahn et al.  Impact of diagnostic methods on efficacy estimation – a proof-of-principle based on historical examples.  Tropical Medicine & International Health, 25(3):357–363, 2020\.
  * 8. R. T. O’Neill.  On sample sizes to estimate the protective efficacy of a vaccine.  Stat Med, 7(12):1279–1288, Dec 1988.
  * 9. P. J. Diggle.  Estimating prevalence using an imperfect test.  Epidemiology Research International, 2011(608719), 2011.
  * 10. M. Voysey et al.  Safety and efficacy of the ChAdOx1 nCoV-19 vaccine (AZD1222) against SARS-CoV-2: an interim analysis of four randomised controlled trials in Brazil, South Africa, and the UK.  Lancet, Dec 2020.
  * 11. F. P. Polack et al.  Safety and efficacy of the bnt162b2 mrna covid-19 vaccine.  New England Journal of Medicine, 383(27):2603–2615, 2020.  PMID: 33301246.
  * 12. L. R. Baden et al.  Efficacy and safety of the mrna-1273 sars-cov-2 vaccine.  New England Journal of Medicine, 0(0):null, 0.
  * 13. I. Arevalo-Rodriguez et al.  False-negative results of initial RT-PCR assays for COVID-19: A systematic review.  PLoS One, 15(12):e0242958, 2020.
  * 14. J. Balayla.  Prevalence threshold and the geometry of screening curves.  PLoS One, 15(10):e0240215, 2020.
  * 15. G. J. Griffith et al.  Collider bias undermines our understanding of COVID-19 disease risk and severity.  Nat Commun, 11(1):5749, 11 2020.
  * 16. J. P. Ioannidis.  Why most discovered true associations are inflated.  Epidemiology, 19(5):640–648, Sep 2008.
  * 17. K. S. Button et al.  Power failure: why small sample size undermines the reliability of neuroscience.  Nat Rev Neurosci, 14(5):365–376, 05 2013.
  * 18. J. Balayla.  Bayesian updating and sequential testing: Overcoming inferential limitations of screening tests, 2020.

## Appendix

⬇

#

# R function for Bayesian estimation of the efficacy rate using an

# imperfect test, originally modified from code provided by

# P. J. Diggle. Estimating prevalence using an imperfect test. \

# Epidemiology Research International, 2011(608719), 2011

#

# Notes

#

# 1. Prior for efficacy rate is assumed uniform on (0,1)

# 2. Priors for sensitivity and specificity are independent scaled

# beta distributions

# 3. Function uses a simple quadrature algorithm with number of

# quadrature points as an optional argument ”ngrid” (see below);

# the default value ngrid=20 has been sufficient for all examples

# tried by the author, but is not guaranteed to give accurate

# results for all possible values of the other arguments.

#

library(Rmpfr)

.N <- function(.) mpfr(., precBits = 200)

efficacy.bayes<-function(alpha,Tc,n,pi,lowse=0.5,highse=1.0,

sea=1,seb=1,lowsp=0.5,highsp=1.0,

spa=1,spb=1,ngrid=20) {

#

# arguments

# alpha: vector of efficacy rates for which posterior is required

# (will be converted internally to increasing sequence of equally

# spaced values, see ”result” below)

# Tc: number of positive test results

# n: number of indiviudals tested

# lowse: lower limit of prior for sensitivity

# highse: upper limit of prior for sensitivity

# sea,seb: parameters of scaled beta prior for sensitivity

# lowsp: lower limit of prior for specificity

# highsp: upper limit of prior for specificity

# spa,spb: parameters of scaled beta prior for specificity

# ngrid: number of grid-cells in each dimension for quadrature

#

# result is a list with components

# alpha: vector of efficacy rate for which posterior density has

# been calculated

# post: vector of posterior densities

# mode: posterior mode

# interval: maximum a posteriori credible interval

#

ibeta<-function(x,a,b) {

pbeta(x,a,b)*beta(.N(a),.N(b))

}

nalpha<-length(alpha)

bin.width<-(alpha[nalpha]-alpha[1])/(nalpha-1)

alpha<-alpha[1]+bin.width*(0:(nalpha-1))

integrand<-array(0,c(nalpha,ngrid,ngrid))

h1<-(highse-lowse)/ngrid

h2<-(highsp-lowsp)/ngrid

for (i in 1:ngrid) {

se<-lowse+h1*(i-0.5)

pse<-(1/(highse-lowse))*dbeta((se-lowse)/(highse-lowse),sea,seb)

for (j in 1:ngrid) {

sp<-lowsp+h2*(j-0.5)

psp<-(1/(highsp-lowsp))*dbeta((sp-lowsp)/(highsp-lowsp),spa,spb)

if(ngrid==1){

se=highse; sp=highsp; pse=1; psp=1; h1=1; h2=1

}

c1<-1-sp

c2<-se+sp-1

g<-(c1+c2*pi)*chooseMpfr(n,Tc)*(ibeta(c1+c2*pi,Tc-1,n-Tc+1)-ibeta((c1+c2*pi)/2,Tc-1,n-Tc+1))

p<-(c1+c2*pi)/(2-alpha)

density<-rep(0,nalpha)

for (k in 1:nalpha) {

density[k]<-asNumeric(dbinom(Tc,n,.N(p[k]))/g)

}

integrand[,i,j]<-pse*psp*density

}

}

post<-rep(0,nalpha)

for (i in 1:nalpha) {

post[i]<-h1*h2*sum(integrand[i,,])

}

ord<-order(post,decreasing=Tc)

mode<-alpha[ord[1]]

cumpost=cumsum(post/sum(post))

interval=c(alpha[which.min(abs(cumpost-0.025))],

alpha[which.min(abs(cumpost-0.975))])

list(alpha=alpha,post=post,mode=mode,interval=interval)

}

#

# example, to reproduce Table 1 results set ngrid=1, highse=1, highsp=1

# oxford 30 of 5807 vs 101 of 5829

# pfizer 8 of 18198 vs 162 of 18325

# moderna 11 of 14134 vs 185 of 14073

N<-14134+14073 # total number of participants

Tc<- 185 # number of cases in control

pi<-(185+11)/(14134+14073) # total number of cases / N

ngrid<-1

lowse<-1

highse<- 1 #0.95

lowsp<-1

highsp<-1 #0.999

sea<-2 #5

seb<-2

spa<-2 #5

spb<-2

alpha<-seq(0,1,by=0.0005) #0.001*(1:400)

result<-efficacy.bayes(alpha,Tc,N,pi,lowse,highse,

sea,seb,lowsp,highsp,spa,spb,ngrid)

result$mode

result$interval

plot(result$alpha,result$post/sum(result$post),type=”l”,xlab=”alpha”,ylab=”p(alpha)”)

# Function to reproduce per-group sample sizes presented in O’Neill,

# 1988. Note that he used relative width (RW) and attack rate in

# unvaccinated group (ARU)

sample_size_95wald_paper <- function(ARU,RW){

VE=0.4

y<-RW*VE/(2*(1-VE))

d<-log(y+sqrt(y^2+1))

return((1.96)^2/d^2*((1+1/(1-VE))/ARU-2))

}

print(t(outer(c(0.01,0.005,0.001,0.0005),rev(seq(0.1,1,0.1)),sample_size_95wald_paper)),quote
= FALSE)

# Total sample sizes for pooled Wald variance (power 80%

# and alpha=0.05)

sample_size_wald <- function(VE,Pi,delta){

RW=delta/VE

ARU=Pi/(2-VE)

y<-RW*VE/(2*(1-VE))

d<-log(y+sqrt(y^2+1))

return(2*(1.96+0.84)^2/d^2*((1+1/(1-VE))/ARU-2))

}

# Total sample sizes for Cramer-Rao bound (power 80% and alpha=0.05)

sample_size_cramer <- function(VE,Pi,delta){

return(4*(1.96+0.84)^2*(2-VE)^2*(2-VE-Pi)/Pi/delta^2)

}