# Interpreting DESI’s evidence for evolving dark energy

Marina Cortês  and Andrew R. Liddle

###### Abstract

The latest results on baryon acoustic oscillations from DESI (Dark Energy
Spectroscopic Instrument), when combined with cosmic microwave background and
supernova data, show indications of a deviation from a cosmological constant
in favour of evolving dark energy. Use of a pivot scale for the equation of
state $w$ shows that this evidence is concentrated in the derivative of $w$
rather than its mean offset from $-1$, indicating a new cosmic coincidence
where the mean equation of state matches that of the $\Lambda$CDM model
precisely in the region probed by the observations. An equivalent way to
express this is to say that the dark energy hits the maximum value that it
will ever achieve within the observed window. We argue that conclusions on
dark energy evolution are strongly driven by the assumed parameter priors and
that this coincidence, which we are naming the PhantomX coincidence (where X
stands for crossing), may be a signature of this.

##  1 Introduction

After over a decade of relentless verification of the standard $\Lambda$CDM
model through experiment, cosmology is a field eagerly seeking any signs of
deviation from this simple paradigm. Hence the excitement that measurements of
baryon acoustic oscillations (BAO) by the Dark Energy Spectroscopic Instrument
(DESI) [1, 2], when combined with cosmic microwave background (CMB) and
supernova observations, show tantalising hints of dark energy evolution [3].
The deviations within a two-parameter dark energy model are estimated in Ref.
[3] variously as 2.5 sigma, 3.5 sigma, and 3.9 sigma depending on the
supernova dataset included in the compilation (respectively, PantheonPlus [4],
Union3 [5], and DESY5 [6]).

In this article we examine the robustness of this claim, focussing almost
entirely on the case of the flat Universe two-parameter equation of state
model labelled w0waCDM.

##  2 Dark energy evolution in the w0waCDM model

###  2.1 Models, priors, and pivots

The w0waCDM model is a phenomenological model of dark energy which features a
two-parameter dark energy equation of state [7, 8]

$w(a)=w_{0}+w_{a}(1-a)\,,$ (2.1)

where $a$ is the scale factor normalized to unity at present, and $w_{0}$ and
$w_{a}$ are constants. Crucial for our discussion is the choice of prior
ranges for these parameters. In Ref. [3] they are taken to be uniform in the
ranges $[-3,1]$ and $[-3,2]$ respectively, with the additional condition
$w_{0}+w_{a}<0$ to allow early matter domination.111In practice the data go
nowhere near this last condition, which would be entirely off the top of the
constraint plot shown later. The cosmological constant, or $\Lambda$CDM, model
corresponds to $w_{0}=-1$ and $w_{a}=0$, and a regime in which $w(a)<-1$ is
called a ‘phantom’ regime.

Theoretical modellers typically find the non-phantom regime much simpler to
model than a phantom regime [9]. For instance, scalar fields with canonical
kinetic terms, normally called quintessence models, are necessarily non-
phantom, and have already been discussed in the context of the DESI results
[10, 11] (other physically-motivated scenarios being discussed in Ref. [12]).
Ref. [13] provides a sobering assessment of how useful the $w_{0}$-$w_{a}$
parametrization is in developing understanding of the fundamental nature of
dark energy. Models that cross between phantom and non-phantom pose additional
modelling difficulties due to singularities in the perturbation equations
[14].

The raw $w_{0}$–$w_{a}$ parametrization has the disadvantage that constraints
on the two parameters are typically strongly correlated, because there is no
useful data constraining the value of $w$ at the present epoch, which is what
$w_{0}$ represents. Fortunately, the model has a reparametrization invariance
under a change of pivot scale, i.e. the choice of scale factor $a_{\rm p}$ at
which to specify $w$; we replace $w_{0}$ with $w_{\rm p}$ according to

$w_{\rm p}=w_{0}+w_{a}(1-a_{\rm p})\,.$ (2.2)

where $a_{\rm p}$ is chosen to decorrelate $w_{\rm p}$ and $w_{a}$. The pivot
scale depends on the data being fit, the main disadvantage of the method being
that datasets with different pivots should not be overlaid on the same plot.
Marginalized uncertainties on $w_{a}$ are unchanged by the pivot
transformation, but those on $w_{\rm p}$ will be diminished compared to those
on $w_{0}$. The constraint area in the 2D parameter plane is preserved by the
transformation, while in a Gaussian approximation the uncertainty on $w_{\rm
p}$ will match that of a constant equation of state model ($w_{a}$ assumed
zero).

Pivot scales were originally utilized in CMB studies, capitalising on the
freedom to choose a wavenumber $k$ on which to specify the amplitude, slope,
etc. of perturbations ${\cal P}(k)$. They were first exploited in Ref. [15],
to decorrelate measurements of the spectral index and its running. Later they
featured notably in the WMAP Collaboration analyses, which set fixed scales
for defining the scalar and tensor amplitudes from its three-year science
release onwards [16], while optimizing the pivot choice for best power
spectrum constraints was demonstrated in Ref. [17].

The use of pivots for dark energy, where the pivot is at a time rather than a
length, was initiated by Huterer and Turner [18] who used a pivot within a
parametrization of $w$ with redshift. In its modern form using the w0waCDM
parametrization, it was popularized by the Dark Energy Task Force Report of
Albrecht et al. [19] who used it in defining a much-used Figure of Merit of
observational configurations (see also Refs. [20, 21]). The pivot scale
corresponds to the scale factor at which a given data compilation best
constrains the normalization of $w(a)$, which can be considered its average
value over the observational window.

###  2.2 Evolving dark energy constraints

We first note that the DESI analysis for models with constant $w$ ($w_{a}$
assumed zero) finds in all cases results consistent with $\Lambda$CDM [3].
There is no evidence for a constant offset of $w$ from the $\Lambda$CDM value
[3].

Figure 1:  Observational constraints in the $w_{0}$–$w_{a}$ plane from Ref.
[3], combining DESI BAO and CMB constraints with three different choices of
supernova sample. The magenta and red lines partition models into phantom and
non-phantom behaviour at early times and today, respectively. In combination
they cut the plane into four zones. The blue and orange lines mark parameter
values where $w$ crosses $-1$ at redshifts 0.26 and 0.33 respectively. These
correspond to the pivot redshifts for the PantheonPlus and DESY5 supernova
samples (blue) and Union3 (orange). This shows that all three choices have $w$
close to $-1$ at the pivot scale. [Adapted from Figure 6 of Ref. [3], under
Creative Commons BY 4.0 License.]

We thus move onto the two-parameter w0waCDM model. In Figure 1 we show the
results from the DESI paper (Ref. [3], their Figure 6 right panel), upon which
we have superimposed a number of lines to guide our interpretation. The DESI
data are BAO measurements from a variety of object types [1, 2]. The CMB data
is mostly from the Planck Satellite [22] with additional information from the
Atacama Cosmology Telescope [23] in measures of CMB lensing. As described
above, three options for catalogues of supernova luminosity–redshift relations
are deployed for comparison, PantheonPlus [4], Union3 [5], and DESY5 [6].

Before considering the data, note the red and purple lines which divide the
parameter plane into four regions. To the lower left are models which are in
the phantom regime throughout cosmic history. On the upper right are models
that are never phantom, where in particular all quintessence models lie, while
the small area in the top left contains models which started non-phantom but
transitioned to the phantom regime by the present. All the observational
contours lie in the fourth region (the largest region as shown in the figure,
though the plot shows only a small part of the entire prior domain); these
models transition from an early-time phantom regime to a present-day non-
phantom regime.

The most striking feature of the constraints in Figure 1 is that the elongated
ellipses point closely in the direction of the $\Lambda$CDM point. This
orientation implies that the line defining the pivot scale (the major axis of
the ellipses) nearly coincides with the line defining models which have
$w_{\rm p}=-1$. That line is given by $-1-w_{0}=w_{a}(1-a_{\rm p})$ for a
chosen $a_{\rm p}$. The pivot redshifts found in Ref. [3] of $z=0.26$ (for
PantheonPlus and DESY5) and $z=0.33$ (for Union3) correspond to $a_{\rm p}$ of
$0.79$ and $0.75$ respectively, and the lines in model space where $w_{p}=-1$
at those redshifts are shown by the blue and orange lines. The observational
constraints indeed lie on top of those lines and are oriented along them.

In reality the observations are not probing the whole redshift range, but
rather a window centred on the pivot scale. That’s true even for the CMB
probes, which are principally from the dark energy’s late-time effect on the
angular-diameter distance to last-scattering. At high redshift the dark energy
density is too low to have any observable effect, while at low redshift there
is too little volume to take constraining data.

The best-fitting models all have the following characteristic. They start deep
in the phantom regime, with $w$ increasing rapidly. Just as we reach the
redshifts which are most strongly constrained by the observations, $w$ reaches
the $\Lambda$CDM value $w=-1$, to a precision around $\pm 0.02$. At later
redshifts, passing outside the observational window, $w$ continues its assumed
linear ascent to a value significantly above $-1$. The best-fit $w(a)$ models
are shown in Figure 2, with the pivot values indicated by blobs.

Figure 2:  The best-fit $w(a)$ evolutions for the three choices of supernova
dataset, colour coded as in Figure 1. The pivot scale factors are $a_{\rm
p}=0.79$ for the green and blue lines and $a_{\rm p}=0.75$ for the orange,
with the pivot values of $w$ indicated by the blobs. The PhantomX Coincidence
is that the blobs are so close to $w=-1$ when most of the evolution is not.

The outcome is that the preferred regions establish a new cosmic coincidence.
Over cosmic history the equation of state of the preferred models exhibits
order one variation, but at the epoch of observation it is within the special
value of $w=-1$ by a few hundredths. Within the actual range constrained by
observations the likely variation is a few tenths in the more extreme models,
though there are currently no good estimates of the redshift range that is
well constrained by these observations. Gaussian process modelling of $w(a)$
[24, 25] would be one way to achieve this. Since the coincidence is that
phantom crossing occurs at or near the centre of the observational window, we
call this the PhantomX Coincidence (X stands for crossing). It is separate
from the coincidence that dark energy domination ($\Omega_{\rm de}=0.5$) also
takes place during this epoch, the onset of acceleration being somewhat
earlier.

Since the phantom transition marks the point at which the dark energy density
$\rho_{\rm de}$ stops increasing and starts to decrease, an equivalent
statement of the coincidence is that the maximum value that the dark energy
density will ever reach happens to lie in the observed window.

Shifting to the pivot scale corresponds to sliding (not rotating) the points
on the blue or orange lines of Figure 1 horizontally onto the $w_{0}=-1$ axis,
dragging the corresponding contours with them. This shifts the contours to
have an accurately vertical major axis (corresponding to decorrelation of
$w_{\rm p}$ and $w_{a}$), with all three 95% contours overlapping the $w_{\rm
p}=-1$ axis line.222The chains used to plot Figure 1 have not yet been made
public so we are unable to make the transformation ourselves to replot the
figure, though the outcome is easy to estimate by eye. This mapping preserves
contour areas, i.e. the induced prior on $w_{\rm p}$ is still flat.

In a couple of places, Ref. [3] refers to a need for extra data to break the
degeneracy between $w_{0}$ and $w_{a}$. While indeed those parameters are
degenerate, this is not really the correct picture. The freedom to make a
pivot transformation alone is sufficient to completely remove this degeneracy
for any dataset, decorrelating measurements of the amplitude and slope of
$w(a)$. What they really mean is that extra data are needed to better
constrain $w_{a}$ in the pivot parametrization.

On a separate point, Figure 1 of Ref. [3] shows only a single DESI BAO point,
the $z=0.51$ LRG $D_{\rm M}/D_{\rm H}$ ratio of transverse and line-of-sight
distances, to be significantly displaced from the best-fit $\Lambda$CDM model.
From this one might infer that the significance of the evolution is largely
driven by this datapoint. This has seen much discussion since the results were
released, especially since the SDSS estimates at similar redshifts [26], seen
in Figure 15 of Ref. [1], do not show this discrepancy from $\Lambda$CDM.
However, the analysis in Appendix A of Ref. [3], which swaps out the low-
redshift DESI BAO points for their SDSS equivalents, shows only a modest
reduction in the significance of around half a sigma in each of the data
combinations. The result is therefore due to a more complex intersection of
the different datasets in the full parameter space. Nevertheless, given the
significant offset between SDSS and DESI estimates at these redshifts, further
data in this regime is crucial.

###  2.3 Interpretation

Let’s look at those initial significances again. Depending on which supernova
dataset is combined with DESI BAO and the CMB, there is a claimed discrepancy
from $\Lambda$CDM of 2.5 sigma, 3.5 sigma, or 3.9 sigma [3], based on a
likelihood-ratio comparison. These numbers look so precise, lending credence
to the discrepancy. But the interpretation carries a substantial unstated
dependence on the choice of prior.

What motivates that prior? Why is the prior uniform, and with such a wide
range? Who thought, before the data came along, that the models with rapidly
escalating $w(a)$ at the bottom were as plausible as the quintessence models
at the top? That phantom models were a priori more likely than non-phantom
models? Who rates models that have to be patched across the phantom divide to
be as reasonable as those which do not? Clearly, the chosen prior is
substantially subjective. Certainly it bears little resemblance to the
distribution of models considered plausible when the searches for dark energy
evolution began in earnest around the time of the Dark Energy Task Force
report [19].

We could change the prior, for instance by tapering its extension into the
deep phantom regime, and in doing so we would lessen this new modest cosmic
coincidence that we have described. But if we alter the prior, and downweight
the models that lie furthest from the original theoretical expectations, we
will also reduce the discrepancy from $\Lambda$CDM. Since the location that
the observations actually constrain lies so close to $w=-1$, the evidence
against $\Lambda$CDM would diminish or perhaps disappear entirely.

While we have concentrated on the flat case, the results allowing varying
spatial curvature are extremely similar [3]. Indeed, curiously the results are
all not only consistent with zero curvature, but practically centred upon zero
far within the uncertainties. Presuming this to be just another coincidence,
our interpretations regarding dark energy apply equal well in the curved case.

To recap, our principal point is to highlight a so-far-unquantified dependence
on model priors. It is clearly discomforting to find the preferred models have
a substantially varying $w(a)$ that just happens to be equal (on average, but
accurately so) to $-1$ in exactly the range probed by the observations (or,
equivalently, that the maximum dark energy density over cosmic history
occurred within the observed window). This new cosmic coincidence of phantom
crossing in the observed window, which we have called the PhantomX
Coincidence, hints that the chosen prior is not a good representation of the
likely underlying physics.

It is not that a particular choice of prior is right or wrong. Bayesian priors
are always a matter of judgement, and sufficiently good data will ultimately
overturn poor assumptions. Indeed, setting of priors is the one place where
the Bayesian method allows the input of personal knowledge and expertise. It’s
a theorist’s chance to show that their reading of Nature is more intuitive
than that of their colleagues. But what is imperative is that any conclusions
drawn from this kind of analysis must assess the robustness against the large
changes that reasonable changes of prior may bring. In this case, we contend
that the robustness has yet to be demonstrated.

##  3 Conclusions

Despite the specific analysis methodology we are highlighting in this
manuscript, the first reaction, by any of us in the community, to the release
of cosmology results from the first year of DESI’s data, can be nothing other
than outstanding and unabashed acclaim. The announcement of this set of
articles, from the first ever Stage-IV survey,333As defined by the Dark Energy
Task Force in 2006 [19]. is the result of consistent efforts by a large team
of researchers, who have persevered over more than a decade of the survey’s
planning, funding acquisition, operational design, and analysis pipeline. We
also wish to highlight that one of the authors in this work was a member of
the DESI survey.444At the time named BigBOSS. Specifically we computed and
produced the cosmology forecasts obtained through the estimation of both BAO
and broadband power spectrum measurements featured in Ref. [27].

We know well that, at the end of the day, every scientist’s worst dilemma is
the ruthless impartiality with which one must approach scientific procedure.
The motivation and revision of one’s full set of initial assumptions and
convictions is amongst the first priorities in the process of data analysis.
It is this scientific rigour that distinguishes our community of enquiry from
other disciplines. Therefore, a more comprehensive understanding of the role
of current assumptions in the $w_{0}$–$w_{a}$ plane will be vital in fully
meeting the high standards of excellence upheld by the US Department of
Energy-led DESI survey collaboration.

## Acknowledgments

We thank Cary Sneider for encouragement, Joe Lovano for hospitality during the
initial phase of this work, Dragan Huterer and Seshadri Nadathur for
clarifications. This work was supported by the Fundação para a Ciência e a
Tecnologia (FCT) through the research grants UIDB/04434/2020 and
UIDP/04434/2020. M.C. acknowledges support from the FCT through the
Investigador FCT Contract No. CEECIND/02581/2018 and POPH/FSE (EC), and A.R.L.
from the Investigador FCT Contract No. CEECIND/02854/2017 and POPH/FSE (EC).
M.C. and A.R.L. are supported by the FCT through the research project
PTDC/FIS-AST/0054/2021. This article/publication is based upon work from COST
Action CA21136 – “Addressing observational tensions in cosmology with
systematics and fundamental physics (CosmoVerse)”, supported by COST (European
Cooperation in Science and Technology).

## References

  * [1] A. G. Adame et al. (DESI Collaboration), _DESI 2024 III: Baryon Acoustic Oscillations from Galaxies and Quasars_ , arXiv:2404.03000 [astro-ph.CO] (2024).
  * [2] A. G. Adame et al. (DESI Collaboration), DESI _2024 IV: Baryon Acoustic Oscillations from the Lyman Alpha Forest_ , arXiv:2404.03001 [astro-ph.CO] (2024).
  * [3] A. G. Adame et al. (DESI Collaboration), _DESI 2024 VI: Cosmological Constraints from the Measurements of Baryon Acoustic Oscillations_ , arXiv:2404.03002 [astro-ph.CO] (2024).
  * [4] D. Scolnic et al., _The Pantheon+ Analysis: The Full Data Set and Light-curve Release_ , Astrophys. J. 938, 113 (2022) [arXiv:2112.03863 [astro-ph.CO]].
  * [5] D. Rubin et al., _Union Through UNITY: Cosmology with 2,000 SNe Using a Unified Bayesian Framework_ , arXiv:2311.12098 [astro-ph.CO] (2023).
  * [6] T. M. C. Abbott et al. (DES collaboration), _The Dark Energy Survey: Cosmology Results With $\sim$1500 New High-redshift Type Ia Supernovae Using The Full 5-year Dataset_, arXiv:2401.02929 [astro-ph.CO] (2024).
  * [7] M. Chevallier and D. Polarski, _Accelerating Universes with Scaling Dark Matter_ , Int. J. Mod. Phys. D10, 213 (2001) [arXiv:gr-qc/0009008].
  * [8] E. V. Linder, _Exploring the Expansion History of the Universe_ Phys. Rev. Lett. 90, 091301 (2003) [arXiv:astro-ph/0208512].
  * [9] E. J. Copeland, M. Sami, and S. Tsujikawa, _Dynamics of dark energy_ , Int. J. Mod. Phys. D15, 1753 (2006) [arXiv:hep-th/0603057].
  * [10] Y. Tada and T. Terada, _Quintessential interpretation of the evolving dark energy in light of DESI_ , arXiv:2404.05722 [astro-ph.CO] (2024).
  * [11] W. Yin, _Cosmic Clues: DESI, Dark Energy, and the Cosmological Constant Problem_ , arXiv:2404.06444 [astro-ph.CO] (2024).
  * [12] D. Wang, _Constraining Cosmological Physics with DESI BAO Observations_ , arXiv:2404.06796 [astro-ph.CO] (2024).
  * [13] W. J. Wolf and P. G. Ferreira, _Underdetermination of dark energy_ , Phys. Rev. D108, 103519 (2023) [arXiv:2310.07482 [astro-ph.CO]].
  * [14] W. Fang, W. Hu, and A. Lewis, _Crossing the phantom divide with parametrized post-Friedmann dark energy_ , Phys. Rev. D78, 087303 (2008) [arXiv:0808.3125].
  * [15] E. J. Copeland, I. J. Grivell, and A. R. Liddle, _Cosmological parameter estimation and the spectral index from inflation_ , Mon. Not. Roy. Astron. Soc. 298, 1233 (1998) [arXiv:astro-ph/9712028].
  * [16] D. N. Spergel et al. (WMAP Collaboration), _Wilkinson Microwave Anisotropy Probe (WMAP) Three Year Results: Implications for Cosmology_ , Astrophys. J. Suppl. 170, 377 (2007) [arXiv:astro-ph/0603449].
  * [17] M. Cortês, A. R. Liddle, and P. Mukherjee, _On what scale should inflationary observables be constrained?_ , Phys. Rev. D75, 083520 (2007) [arXiv:astro-ph/0702170].
  * [18] D. Huterer and M. S. Turner, _Probing the dark energy: methods and strategies_ , Phys. Rev. D64, 123527 (2001) [arXiv:astro-ph/0012510].
  * [19] A. Albrecht et al., _Report of the Dark Energy Task Force_ , arXiv:astro-ph/0609591 (2006).
  * [20] E. V. Linder, _Biased Cosmology: Pivots, Parameters, and Figures of Merit_ , Astropart.Phys. 26, 102 (2006) [arXiv:astro-ph/0604280].
  * [21] D. Martin and A. Albrecht, _Talk about Pivots_ , [arXiv:astro-ph/0604401] (2006).
  * [22] N. Aghanim et al. (Planck Collaboration), _Planck 2018 results. VI. Cosmological parameters_ , A&A 641, A6 (2020) [arXiv:1807.06209 [astro-ph.CO]].
  * [23] S. Aiola et al., _The Atacama Cosmology Telescope: DR4 maps and cosmological parameters_ , JCAP 2020, 647 (2020) [arXiv:2007.07288 [astro-ph.CO]].
  * [24] T. Holsclaw, U. Alam, B. Sanso, H. Lee, K. Heitmann, S. Habib, and D. Higdon, _Nonparametric Dark Energy Reconstruction from Supernova Data_ , Phys. Rev. Lett. 105, 241302 (2010) [arXiv:1011.3079 [astro-ph.CO]].
  * [25] M. Seikel, C. Clarkson, and M. Smith, _Reconstruction of dark energy and expansion dynamics using Gaussian processes_ , JCAP 06, 036 (2012) [arXiv:1204.2832 [astro-ph.CO]].
  * [26] S. Alam et al. (eBOSS Collaboration), _Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory_ , Phys. Rev. D103, 083533 (2021), [arXiv:2007.08991 [astro-ph.CO]].
  * [27] D. Schlegel et al., _The BigBOSS Experiment_ , arXiv:1106.1706 [astro-ph.IM] (2011).


Safe Sequential Testing and Effect Estimation in Stratified Count Data

Rosanne J. Turner Peter D. Grünwald

CWI, UMC Utrecht CWI, Leiden University

###### Abstract

Sequential decision making significantly speeds up research and is more cost-
effective compared to fixed-$n$ methods. We present a method for sequential
decision making for stratified count data that retains Type-I error guarantee
or false discovery rate under optional stopping, using _e-variables_. We
invert the method to construct stratified anytime-valid confidence sequences,
where cross-talk between subpopulations in the data can be allowed during data
collection to improve power. Finally, we combine information collected in
separate subpopulations through pseudo-Bayesian averaging and switching to
create effective estimates for the minimal, mean and maximal treatment effects
in the subpopulations.

##  1 INTRODUCTION

Fixed-n hypothesis tests and confidence intervals limit research opportunities
and quick decision making, as they rely on static research designs where data
are only evaluated at one time point. We aim to develop hypothesis tests for
conditional independence and anytime-valid confidence sequences for stratified
treatment effects in subpopulations that retain a guarantee on the probability
of falsely rejecting the null hypothesis and coverage of the true effect under
continuous monitoring of data. To this end we use e-values, tools for
constructing tests that keep the _type-I error rate_ (or _false positive
rate_) controlled under sequential testing with optional stopping. Over the
last four years, e-values have become the standard tools (essentially, the
appropriate alternative for $p$-values) for dealing with such settings. Below
we summarize the essentials; for much more background on the budding field of
e-processes (also known as ‘testing by betting’ and ‘safe testing’) see the
recent overview (Ramdas et al.,, 2022) and specifically for details on
e-values refer to Grünwald et al., (2022); Vovk and Wang, (2021). In this
paper, we develop e-processes for stratified $2\times 2$ tables, enabling, in
Section 2, anytime-valid (i.e. valid under optional stopping) conditional
independence (CI) tests for Bernoulli streams for two groups $a$ and $b$ (e.g.
$a$ is control, $b$ is treatment), where the test is conditional on a third
variable, the stratum. Based on these CI tests, we then, in Section 3, develop
anytime-valid confidence sequences (henceforth just called ‘confidence
sequences’) for a notion of effect size representing divergence from CI. The
importance of our tests is ubiquitous in e.g. medical statistics — we can
think of the CI test in Section 2 as an an anytime-valid sequential version of
the Cochran-Mantel-Haenzel test, a work-horse in the field of epidemiology.
Our e-processes are generalizations of those designed for $2\times 2$ tables
(same setting as ours, but with just a single stratum) by Turner et al.,
(2021); Turner and Grünwald, (2022). To achieve the generalization, we employ
tools from the theoretical machine learning literature, most notably the
literature on prediction with expert advice (Cesa-Bianchi and Lugosi,, 2006),
which extends Bayesian learning techniques with ideas such as ‘sleeping’,
‘switching’ and the like. Moreover, inspired by these ideas, we develop the
novel notion of cross-talk between strata, which allows us to make confidence
intervals narrower if outcomes in various strata are interrelated, while
nevertheless remaining valid even if they are not. While for many statistical
models, anytime-valid tests need more data to reach a desired conclusion than
fixed $n$ methods and anytime-valid confidence intervals are somewhat wider
than standard ones (Ramdas et al.,, 2022; Grünwald et al.,, 2022), we find in
this paper that we can partially counteract this difference by employing the
cross-talk strategy (which is not available for fixed-$n$ methods), as is
illustrated by comparing our confidence sequences to fixed-$n$ confidence
intervals for Mantel-Haenszel risk differences in Section 3.

#### E-Processes

Consider a random process $Y_{1},Y_{2},\ldots$ and let $\mathcal{H}_{0}$, the
null hypothesis, be a set of distributions for this process. An e-variable for
$Y_{j},Y_{j+1},\ldots,Y_{m}$ conditional on $Y^{(j-1)}=(Y_{1},\ldots,Y_{j-1})$
for testing $\mathcal{H}_{0}$ is any nonnegative random variable $S$ that can
be written as function of $Y^{(m)}=(Y_{1},\ldots,Y_{m})$ such that

$\forall P\in\mathcal{H}_{0}:\mathbb{E}_{P}[S\mid S^{(j-1)}]\leq 1;$ (1)

for $j=1$ we set $\mathbb{E}_{P}[S\mid S^{(0)}]:=\mathbb{E}_{P}[S]$ and call
$S$ an unconditional e-variable; an e-value is the value an e-variable takes
on a realized sample. It is easily shown that for any sequence
$S_{1},S_{2},\ldots$ where $S_{j}$ is an e-variable for $Y_{(j)}$ conditional
on $Y^{(j-1)}$, the product $E^{(m)}:=\prod_{j=1}^{m}S_{j}$ is an
unconditional e-variable for $Y^{(m)}$. $E^{(1)},E^{(2)},\ldots$ is called a
test martingale or e-process (see Ramdas et al., (2022) on how e-processes
strictly generalize test martingales). Via Ville’s inequality, it is shown
that e-processes have the remarkable property that, for any $0<\alpha<1$, the
probability that there exists an $m$ such that $E^{(m)}\geq 1/\alpha$ is
bounded by $\alpha$. As a consequence, if we look at the data at some time $m$
and reject if $E^{(m)}\geq 1/\alpha$, the probability under the null of
falsely rejecting the null is at most $\alpha$ no matter how we chose $m$; it
may be determined by external circumstances (do we have money to experiment
further?) or by aggressive stopping rules such as ‘keep sampling until you can
reject the null’, or even by peeking into the future. Tests with this property
are called safe under optional stopping and Ramdas et al., (2020) show that,
in essence, all reasonable such tests should be based on e-processes. Just
like p-values can be converted into confidence intervals, e-process can be
converted into anytime-valid confidence tests, also known as confidence
sequences — we will explore these in Section 3.

#### Setting

We consider the _stratified contingency table_ setting/model. Under the global
null hypothesis (we consider more complicated nulls later), outcomes
$Y\in\\{0,1\\}$ are independent of groups $X\in\\{a,b\\}$ (e.g. representing
interventions) given their stratum $k\in[K]:=\\{1,...,K\\}$. We formalize this
by measuring time in terms of blocks: we assume that at each time
$j=1,2,\ldots$, we are given a stratum indicator $k_{j}\in[K]$ and we observe
a block of $n=n_{a}+n_{b}$ outcomes, with $n_{a}$ outcomes in group $a$ and
$n_{b}$ in group $b$, all in the same stratum $k_{j}$. We write
$Y^{(m)}=(Y_{1},\ldots,Y_{m})$ with $Y_{j}$ the data vector corresponding to
the $j$-th block arriving. Hence
$Y_{j}=(Y_{j,a,1},\ldots,Y_{j,a,n_{a}},Y_{j,b,1},\ldots,Y_{j,b,n_{b}})$ is a
vector in $\\{0,1\\}^{n}$ denoting $n=n_{a}+n_{b}$ outcomes in $k_{j}$. Under
both null and alternative, all blocks are assumed independent, with each
outcome in group $x$ in stratum $k$ independently
$\sim\text{Bernoulli}(\theta_{x,k})$. Formally, the null hypothesis then
expresses that

$\mathcal{H}_{0}:\theta_{a,k}=\theta_{b,k}\text{ for all }k.$ (2)

We will assume $n_{a}=n_{b}=1$ for all strata in simulation examples in this
paper, but these can be chosen freely in practice and can even be adapted in
between data blocks — as long as they are set at or before the beginning of a
data block, they are allowed to depend on the past. Of course, in practice, we
often deal with $2K$ i.i.d. streams of data, one for each group-stratum
combination, with data not necessarily coming in at the same rate for
different strata/groups. While superficially different, we can still recast
this setting in terms of blocks: for example, participant may sequentially
enter a study and are each independently randomized with probability $1/2$ to
receive ‘treatment’ (group $b$) or ‘control/placebo’ ($a$). We then wait until
the first time $t_{1}$ that we have seen $n_{a}$ outcomes in group $a$ and
$n_{b}$ outcomes in group $b$ in the same stratum; we call this stratum
$k_{1}$, denote these $n$ outcomes $Y_{1}$, and proceed observing outcomes in
the various streams until the first time $t_{2}$ that there is another stratum
$k_{2}$ (potentially $k_{2}=k_{1}$) so that we have seen $n_{a}$ outcomes in
group $a$, $n_{b}$ in group $b$ in stratum $k_{2}$; we denote these $n$
outcomes $Y_{2}$, and so on. If we want to stop at any time $t$, we take as
data all blocks that have been completed so far, and ignore all started-yet-
unfinished blocks.

#### Related Work

The first paper to use e-processes for conditional independence testing is
(Lindon and Malek,, 2020), but their tests are very different from ours and
involve a simple null hypothesis, allowing them to use Bayes factors for their
e-processes. Further, Turner et al., (2021); Turner and Grünwald, (2022)
develop independence tests and confidence sequence for $2\times 2$ tables; our
paper is a direct extension of theirs, extending their techniques to the
stratum-conditional case. Very recently four other related papers have
appeared: (Pandeva et al.,, 2022; Grünwald et al.,, 2022; Shaer et al.,, 2022;
Duan et al.,, 2022): these papers all differ from ours in that they assume
data are jointly i.i.d. (i.e. one observes a single i.i.d. stream
$(X_{1},Y_{1},K_{1}),(X_{2},Y_{2},K_{2}),\ldots$). The latter three also make
the so-called Model-X assumption (the distribution of $X_{i}\mid K_{i}$ is
assumed known). Our paper is complementary: we do not need the i.i.d. or
Model-X assumption and as explained above, our setting does not just capture
data in blocks (such as paired data) but also data in the form of $2K$ i.i.d.
streams, one for each group in each stratum, with no stochastic assumptions
about what group or what stratum arrives at what time. The price we have to
pay is that we can only deal with a small number of strata and with finite
sets of outcomes and number of groups (in this paper we focus on $2$ but
extension to the finite case is straightforward); aforemetioned references can
deal with arbitrary covariate and outcome random variables $K_{i}$ and
$Y_{i}$. Nevertheless, small-strata-count-studies are highly common in the
medical statistics world, and we show here how to construct efficient
sequential tests for them.

The code used for experiments in this paper will initially be placed on the
repository linked to this publication (Turner,, 2023), and will later be
integrated in the safestats R package (Ly et al.,, 2022).

##  2 E-VARIABLES FOR TESTING THE GLOBAL NULL

We first consider the case where there is only one stratum, $k_{j}=k^{*}$ for
each each $j$. The problem is then reduced to testing whether two Bernoulli
data streams come from the same source. Turner et al., (2021) showed that in
this case, for arbitrary estimators
$\breve{\theta}_{a}|Y^{(j-1)},\breve{\theta}_{b}|Y^{(j-1)}$, the following is
an e-variable for $Y_{j}$ conditional on $Y^{(j-1)}$, i.e. (1) holds with
$S:=S_{j}$ given by

$\displaystyle
S_{j}=\prod_{i=1}^{n_{a}}\frac{p_{\breve{\theta}_{a}|Y^{(j-1)}}(Y_{j,a,i})}{p_{\breve{\theta}_{0}|Y^{(j-1)}}(Y_{j,a,i})}\prod_{i=1}^{n_{b}}\frac{p_{\breve{\theta}_{b}|Y^{(j-1)}}(Y_{j,b,i})}{p_{\breve{\theta}_{0}|Y^{(j-1)}}(Y_{j,b,i})},$
(3)

where $p_{\theta}(Y)=\theta^{Y}(1-\theta)^{1-Y}$ denotes the
Bernoulli$(\theta)$ probability of $Y\in\\{0,1\\}$), as long as we pick
$\breve{\theta}_{0}\in\Theta_{0}=[0,1]$ as follows:

$\displaystyle\breve{\theta}_{0}=\breve{\theta}_{0}|Y^{(j-1)}$
$\displaystyle:=\arg\min_{\theta\in[0,1]}D(P_{\breve{\theta}_{a},\breve{\theta}_{b}}\|P_{\theta,\theta})$
$\displaystyle\overset{(a)}{=}\frac{n_{a}}{n}\breve{\theta}_{a}|Y^{(j-1)}+\frac{n_{b}}{n}\breve{\theta}_{b}|Y^{(j-1)}.$
(4)

Here and in the sequel, $P_{\theta_{a},\theta_{b}}$ represents the
distribution on $n_{a}+n_{b}$ independent binary outcomes with the first
$n_{a}$ outcomes $\sim$ Bernoulli$(\theta_{a})$ and the subsequent $n_{b}$
outcomes $\sim$ Bernoulli$(\theta_{b})$, i.e. the distribution of outcomes in
a single block according to $(\theta_{a},\theta_{b})$, and
$D(P_{\theta_{a},\theta_{b}}\|P_{\theta^{\prime}_{a},\theta^{\prime}_{b}})$
abbreviates the KL divergence between two such distributions. Equality (a)
follows by simple calculus.

Importantly, in (3),
_$(\breve{\theta}_{a},\breve{\theta}_{b})\in\Theta_{1}=[0,1]^{2}$ can be
chosen as a function of past data anyway we like, not affecting the Type-I
error guarantee_. Nevertheless, if we were given the true probabilities
$\theta^{*}_{a}$ and $\theta^{*}_{b}$ of the two groups in block $j$, then we
could set $\breve{\theta}_{a}=\theta^{*}_{a}$ and
$\breve{\theta}_{b}=\theta^{*}_{b}$ and this choice is special: the e-variable
(3) then has, among all e-variables, the largest expected logarithm under the
true alternative $P_{\theta^{*}_{a},\theta^{*}_{b}}$. We then say it is
_growth-rate optimal_ (GRO) for collecting evidence against the null
hypothesis (Grünwald et al.,, 2022). Formally, we define

$\textsc{gro}(\theta^{*}_{a},\theta^{*}_{b}):=\sup_{S}{\bf E}_{Y_{j}\sim
P_{\theta^{*}_{a},\theta^{*}_{b}}}[\log S]$ (5)

where the supremum is over all random variables $S$ that are e-variables for
$Y_{j}$ under ${\cal H}_{0}$. It directly follows from (Grünwald et al.,,
2022, Theorem 1) that, if we plug in $\breve{\theta}_{a}=\theta^{*}_{a}$ and
$\breve{\theta}_{b}=\theta^{*}_{b}$ into (3), then the resulting $S_{j}$ is
GRO and its growth rate is equal to the KL divergence, i.e.

${\bf E}_{Y_{j}\sim P_{\theta^{*}_{a},\theta^{*}_{b}}}[\log
S_{j}]=\textsc{gro}(\theta^{*}_{a},\theta^{*}_{b})=D(P_{\theta^{*}_{a},\theta^{*}_{b}}\|P_{\tilde{\theta},\tilde{\theta}}),$
(6)

where $\tilde{\theta}=(n_{a}/n)\theta^{*}_{a}+(n_{b}/n)\theta^{*}_{b}$.
Growth-rate optimality is the analogue of statistical power in the sequential
setting: if we plug in these ‘true’
$\breve{\theta}_{a}=\theta^{*}_{a},\breve{\theta}_{b}=\theta^{*}_{b}$, we
expect the product $E^{(m)}$ to increase as fast as possible in $m$, enabling
us to reach $1/\alpha$ and reject the null hypothesis as fast as possible,
compared with all other possible e-processes. In practice though,
$\theta^{*}_{a}$ and $\theta^{*}_{b}$ are unknown, but to get near-grow-
optimal e-variables, we can estimate $\breve{\theta}_{a}$ and
$\breve{\theta}_{b}$ based on all data seen before data block $j$ — then
$\breve{\theta}_{a}$ and $\breve{\theta}_{b}$ converge to
$\theta^{*}_{a},\theta^{*}_{b}$ and our e-variables $S_{j}$ get better and
better in the GRO sense. We follow Turner et al., (2021) who successfully
chose to place a beta prior on the parameter space and took the Bayesian
posterior mean as an estimate.

In treatment/ control test settings, there often exists prior knowledge of a
minimal clinically relevant or expected odds ratio
$\text{OR}(\theta_{a},\theta_{b}):=(\theta_{b}/(1-\theta_{b}))((1-\theta_{a})/\theta_{a})$,
i.e. it is known that $\text{OR}(\theta_{a},\theta_{b})=\phi$ for some given
$\phi$. In that case, one can restrict estimating $\breve{\theta}_{a}$ and
$\breve{\theta}_{b}$ to
$\Theta_{1}(\phi)=\\{(\theta_{a},\theta_{b});\text{OR}(\theta_{a},\theta_{b})=\phi\\}$,
possibly improving power and growth-rate of the test (Turner et al.,, 2021).
Both search spaces are illustrated in Figure 1.

Figure 1: Parameter space $\breve{\theta}_{a}|Y^{(j-1)}$ and
$\breve{\theta}_{b}|Y^{(j-1)}$ are estimated in, in $2\times 2$ table without
strata; either through placing a beta prior on the entire unit square (in
light orange) and calculating the posterior mean with all data up to and
including time $j-1$ or through restricting the posterior estimation to a
particular odds ratio value $\phi$ and placing a beta prior on all pairs
$(\theta_{a},\theta_{b})$ corresponding to this odds ratio value (for example
the red curve, for $\phi=2$).

#### Combining e-variables from individual strata

We can use the e-variable in (3) to calculate e-process values $E^{(m),k}$ for
each stratum $k$ separately. To be precise, we set $S^{k}_{j}$ to the
equivalent of (3) if $k=k_{j}$,

$S^{k}_{j}=\prod_{i=1}^{n_{a}}\frac{p_{\breve{\theta}_{a,k}|Y^{(j-1)}}(Y_{j,a,i})}{p_{\breve{\theta}_{0,k}|Y^{(j-1)}}(Y_{j,a,i})}\prod_{i=1}^{n_{b}}\frac{p_{\breve{\theta}_{b,k}|Y^{(j-1)}}(Y_{j,b,i})}{p_{\breve{\theta}_{0,k}|Y^{(j-1)}}(Y_{j,b,i})},$
(7)

and $S^{k}_{j}=1$ otherwise, i.e. if $k_{j}\neq k$, and
$E^{(m),k}:=\prod_{j=1}^{m}S^{k}_{j}$ — note that at each ‘time $j$’, the
product e-variable only changes for the $k$ such that $j$-th block was a block
of outcomes in stratum $k$.

We now need to combine the e-processes-per-stratum into a single e-process for
(2) to measure evidence against ${\cal H}_{0}$ and allowing tests with type-I
error probability guarantee on (2), the global null hypothesis that the odds
ratio of the success probabilities equals $1$ in each stratum. There are
several ways to do this. The first and most straightforward option is to
_multiply_ the individual e-values across the strata:

$E^{(m)}=\prod_{j=1}^{m}S_{j}^{k_{j}}=\prod_{j=1}^{m}\prod_{k=1}^{K}S_{j}^{k}.$
(8)

To see that $E^{(1)},E^{(2)},\ldots$ is an e-process, simply note that each
$\smash{S_{j}^{k_{j}}}$ is a conditional e-variable (i.e. it satisfies (1)
with $S=\smash{S_{j}^{k_{j}}}$) since, given that $S_{j}$ in (3) is a
conditional e-variable, $\smash{S_{j}^{k_{j}}}$ must be an e-variable as well.
When $\theta_{a,k}\approx\theta_{b,k}$ in a few of the strata, this might be a
data-inefficient approach, as one would need to collect a lot of extra
evidence in the strata where the success probabilities are substantially
different to counteract the expected small e-values in the other strata. A
second option that possibly better handles these cases is to create a _convex
combination_ , i.e. a mixture, of e-values at each time point $j$ (any convex
combination of e-variables is also an e-variable (Vovk and Wang,, 2021)): A
simple first option is to pick some prior distribution on the strata $\pi(k)$,
and to use that distribution for calculating the mixture after each batch
comes in:

$\displaystyle S_{j}:=\sum_{k=1}^{K}\pi(k)S^{k}_{j}\ ;\
E^{(m)}=\prod_{j=1}^{m}S_{j}\text{\ so that also}$
$\displaystyle\text{$E^{(m)}=\prod_{k=1}^{K}E^{(m),k}$ with
$E^{(m),k}=\prod_{j=1}^{m}S_{j}^{k}$}.$ (9)

Extending the simple averaging above, we could replace the prior $\pi(k)$ in
(2) with a distribution $\pi(k|y^{(j-1)})$ that depends on previous data
$y^{(j-1)}$, since, since we assume the data itself in each block are
independent, dependency of $\pi$ on past data will not affect guarantee (1).
Such an approach is called the method of mixtures in the anytime-valid testing
literature (Ramdas et al.,, 2022). Thus, any distribution on $[K]$ that
depends on the past is allowed here, but an intuitive choice is a pseudo-
Bayesian posterior

$\pi(k|y^{(j-1)}):=\frac{\pi(k)(E^{(j-1),k})^{\eta}}{\sum_{k^{\prime}}\pi(k^{\prime})(E^{(j-1),k^{\prime}})^{\eta}},$
(10)

where by definition, $E^{(0)}=1$ and we pick $\eta$ beforehand as a _learning
rate_ : if we set it to a higher value, we will focus on strata with higher
e-values more quickly; with $\eta=1$, (10) becomes similar to a Bayesian
posterior. Just as the beta-posterior used to determine $\breve{\theta}_{x,k}$
in (7) allows us to learn $\theta^{*}_{x,k}$, this new posterior allows us to
learn which strata can help us most to reject the null. However, even for
$\eta=1$ the analogy to Bayes only goes so far — for example, at each $j$,
only the e-variable $S^{k_{j}}$ for stratum $k_{j}$ changes; the other $S^{k}$
‘sleep’ (Koolen and van Erven,, 2010) and thus $E^{(j-1),k}$ behaves
differently from a likelihood. This more general past-determined updating
originates in the area of machine learning called prediction with expert
advice where many other such ‘posterior’-updates have been considered
(Herbster and Warmuth,, 1998; van Erven et al.,, 2007; Koolen and de Rooij,,
2013). These include the more extreme approach called _switching_. With this
approach, we calculate (2) with $\pi(k)$ replaced by any distribution we like
(the choice is again allowed to depend on $Y^{(j-1)}$) up to and including a
particular batch $j^{*}$. Thereafter, for $j\geq j^{*}$, we set

$\pi^{*}(k|y^{(j)})=\begin{cases}1&\text{if $k=k^{*}$ with
$k^{*}=\arg\max_{k}E^{(j^{*}),k}$}\\\ 0&\text{otherwise}\end{cases}$ (11)

creating a new E-process $E^{(1)}_{[j^{*}]},E^{(2)}_{[j^{*}]},\ldots$ such
that, for $m\leq j^{*}$, $E^{(m)}_{[j^{*}]}=E^{(m)}$ and, for $m>j^{*}$,

$E^{(m)}_{[j^{*}]}=E^{(j^{*})}\cdot\prod_{j=j^{*}+1}^{m}E^{(j),k^{*}}$ (12)

$j^{*}$ could arbitrarily be picked prior to the study, or we could also place
a prior on the moment of switching and take a weighted average over (12) for
various values of $j^{*}$ for each $\delta$, thereby obtaining yet another
e-process with $j^{*}$ ‘integrated out’ (see Figure S2 in the supplementary
material for a more elaborate comparison of switch priors in a simulation
experiment for confidence sequences).

In Figure 2, the three different methods for combining e-variables for testing
$\mathcal{H}_{0}$ are compared with respect to _power_ : the expected
probability of rejecting $\mathcal{H}_{0}$ under some fixed data generating
distribution. For Figure 2, data were sampled from a distribution where risk
differences and control group rates all differed between strata. It can be
observed that all methods that took the stratification into account
outperformed the unstratified approach, where just one sequential e-variable
was calculated for all strata combined. The three different methods will be
re-compared for confidence sequences in Figure 6.

Figure 2: Power for rejecting the null at level $\alpha=0.05$ that the odds
ratio in all strata equals 1 estimated with 1000 repeated experiments for
various e-variable combination methods. $40$ batches were collected in each of
three strata (so maximum sample size was $m=120$) and sampling was stopped as
soon as $E^{(m)}\geq\frac{1}{\alpha}$. Real control group success rates were
$0.1,0.2,0.8$ and real risk differences were $0.05,0.4,-0.6$. Pseudo-Bayesian
approaches were implemented with learning rates (LR) $1$ and $2$. Switch
approaches were implemented for switching at point $j^{*}=10$, or with a
uniform prior on switch times $j=5$ until $m-5$.

#### Cross-talk between strata

To further improve power of the hypothesis test, we will allow for cross-talk
between strata while estimating $\breve{\theta}_{a,k}$ and
$\breve{\theta}_{b,k}$ based on data seen so far. In the current simple
setting of testing the global null, ‘cross-talk’ simply amounts to design
$S_{j}$ that grow faster (allowing for faster rejecting of the null) if the
alternative satisfies certain constraints. For example, if one expects
treatment effects (say, measured as odds ratios) to be stable (identical)
throughout different strata, but control group recovery rates to vary, one
would like cross-talk about the odds ratios between strata. Practically, this
means that to arrive at the estimates $\breve{\theta}_{x,k}\mid Y^{(j-1)}$, we
first limit the parameter space to $\Theta_{1}(\smash{\hat{\phi}^{(j-1)}})$,
i.e. all vectors $\theta_{x,k}$ with odds ratio $\hat{\phi}^{(j-1)}$, set to
be the maximum likelihood odds ratio based on all previous data in all strata,
i.e. calculated by ignoring strata. We then calculate
$\breve{\theta}_{x,k}\mid Y^{(j-1)}$ as posterior means using beta priors
conditioned on the parameters being in $\Theta_{1}(\hat{\phi}^{(j-1)})$.
Similarly, when one expects control group recovery rates to be stable, but the
treatment effects to vary because of a possible interaction with stratum
characteristics, allowing cross-talk about control group recovery rates might
improve power. In practice, we achieve this by using as beta prior parameters
for the control group rate $\breve{\theta}_{a,k}|Y^{(j-1)}$ the total counts
of failures and successes aggregated over all strata (summed with some initial
prior parameters to ensure stable estimates at time point $j=1$; we set
initial prior values $0.18$ for both the fail and success rate based on a
suggestion by (Turner et al.,, 2021)). In the odds-ratio cross-talk scenario,
we effectively constrain the parameters of the alternative
$\breve{\theta}_{x,k}\mid Y^{(j-1)}$ at each $j$ to share the same odds-ratio;
in the control-group cross-talk scenario, we constrain these parameters to
share the same $\theta_{a}$, i.e.
$\breve{\theta}_{a,k}|Y^{(j-1)}=\breve{\theta}_{a,k^{\prime}}|Y^{(j-1)}$ for
each $k,k^{\prime}$. Would one be unsure whether cross-talk would improve
power at all, and if so, whether one should cross-talk on the odds ratios or
the cross ratios, one could put prior mass $1/3$ on each of the corresponding
three e-values, say $E^{(m)}_{\rho}$ for
$\rho\in\\{\textsc{none},\textsc{odds},\textsc{control rate}\\}$, where none
stands for the standard e-variable without cross-talk. One could then, for
each block $j$, use a mixture e-variable, where the three e-values are mixed
as in (10) with $\eta=1$, $k$ replaced by $\rho$ and ‘$E^{(j-1),k}$’ replaced
by ’$E^{(j-1)}_{\rho}$’ giving a new ‘mix’ e-process. All four cross-talk
scenarios are explored in simulations in Figure 3, where data were generated
from strata with similar control group success rates, but different risk
differences, and different control group success rates, but similar odds
ratios showing that allowing for cross-talk on control rate or odds ratio
improves power in the respective scenarios. The cross-talk mixture performs
comparably to the optimal cross-talk options in both cases. Cross-talk can be
expected to improve power even if, in truth, under the alternative, the odds-
ratio resp. control-group rate is just similar, but not exactly the same under
all groups; and the confidence sequences of the next section remain valid (but
will get wider) even if the odds-ratios resp. control-group rates happen to be
completely different. Thus, the method described here cannot really be viewed
as a constraint on the model, and we chose to call it cross-talk instead: data
in one stratum informs, ‘talks to’ estimates for other strata.

Figure 3: Power for rejecting the null hypothesis at level $\alpha=0.05$ that
the odds ratio in all strata equals 1 estimated with 100 repeated experiments
for various types of cross-talk. $40$ batches were collected in each stratum
and sampling was stopped as soon as $E^{(m)}\geq\frac{1}{\alpha}$. On the
left, real control group success rates were $0.49,0.5$ and $0.51$ in each
stratum; risk differences were $-0.09,-0.49,0.39$. On the right, real odds
ratios were $4,4.01,2.95$.

(a) all different

(b) same control group rate

(c) same OR

(d) all different

(e) same control group rate

(f) same OR

Figure 4: Examples of $95$% stratified confidence intervals ((a), (b) and (c))
and mean confidence interval widths estimated over $100$ runs ((d), (e) and
(f)) with different types of cross-talk. In (a), (b) and (c) the true risk
difference of the data generating distribution in each stratum is indicated by
a dashed line. For (a) and (d), the data were generated by distributions with
different control group success rates ($0.1$, $0.2$ and $0.8$) and risk
differences ($0.05$, $0.4$ and $-0.6$) in each stratum. For (b) and (e),
strata sizes were unbalanced: as can be seen for stratum 1, the red points,
data collection stopped after $10$ batches. Control group success rates were
all $0.5$ and risk differences were different ($-0.49$, $-0.25$ and $0.1$).
For (c) and (f), strata sizes were unbalanced as well, and now odds ratios
were the same in each stratum ($2$), but control group rates differed again
($0.2$, $0.25$ and $0.85$).

#### A GRO-Sanity Check

While the simulations above and below show encouraging empirical results
regarding the power of our methods, it is still useful to have some
theoretical assurance that, no matter the ‘true’ alternative generating the
data, all methods we consider produce e-values that grow fast (i.e. achieve
good power) under this alternative. We now provide a simple theorem to this
end. As usual in the e-value and safe-testing literature, and for reasons
explained by Grünwald et al., (2022), we concentrate on GRO (5) rather than
power.

######  Theorem 2.1.

Suppose that we observe $m=m_{1}+\ldots+m_{K}$ blocks, with $m_{k}$ blocks
lying in stratum $k$, each such block sampled independently from
$P_{\theta^{*}_{a,k},\theta^{*}_{b,k}}$. Then, with ${\bf E}$ denoting
expectation under this distribution, the e-process $E^{(m)}$ defined by
multiplication as in (8) and the mix e-process $E^{(m)}$ as above with
constituent e-processes defined multiplicatively as in (8) both achieve:

$\sum_{k=1}^{K}m_{k}\textsc{gro}(\theta^{*}_{a,k},\theta^{*}_{b,k})={\bf
E}\left[\log E^{(m)}\right]+O(\log m).$ (13)

To interpret the result, note that, if an oracle were to supply us with
$\theta^{*}_{a}=(\theta^{*}_{a,1},\ldots,\theta^{*}_{a,k}),\theta^{*}_{b}=(\theta^{*}_{b,1},\ldots,\theta^{*}_{b,k})$
i.e. if we were told ‘if the alternative were true, then its parameters would
be $P_{\theta^{*}_{a},\theta^{*}_{b}}$’, then we could use the GRO (growth
optimal e-variable) which, conditional on observing a block in stratum $k$,
would obtain the optimal, largest possible expected growth
$\textsc{gro}(\theta^{*}_{a,k},\theta^{*}_{b,k})$. Since we assume data to be
independent, the best growth we could obtain with such an oracle is given by
the left-hand side of (13). The theorem expresses that the price for learning
(via Bayes predictive distributions $\breve{\theta}_{x,k}$ based on beta-
priors) rather than knowing $\theta^{*}_{a},\theta^{*}_{b}$ is modest, namely
logarithmic in $m$ whereas the growth itself is linear in $m$; this is the
standard situation for parametric settings, described in detail by Grünwald et
al., (2022). We may expect the constant hidden in the $O(\log m)$ to become
substantially smaller if the preconditions for effective cross-talk hold as
described above, e.g. odds ratios or group recovery rates are identical or
similar across strata; but determining this constant precisely across cases,
as well as extending the analysis to pseudo-Bayesian and switch e-processes,
is complicated and will be left for future research. The proof of this theorem
can be found in the appendix.

##  3 EXTENSION TO CONFIDENCE SEQUENCES

Turner and Grünwald, (2022) showed that (3) in the $2\times 2$-table (single
stratum) can be generalized, to test null hypotheses
$\mathcal{H}_{0}:=\\{P_{(\theta_{a},\theta_{b})};(\theta_{a},\theta_{b})\in\Theta_{0}\\}$
beyond ‘$\theta_{a}=\theta_{b}$’:

$\displaystyle
S_{j,[\Theta_{0}]}=\prod_{i=1}^{n_{a}}\frac{p_{\breve{\theta}_{a}|Y^{(j-1)}}(Y_{j,a,i})}{p_{\breve{\theta}_{a}^{\circ}|Y^{(j-1)}}(Y_{j,a,i})}\prod_{i=1}^{n_{b}}\frac{p_{\breve{\theta}_{b}|Y^{(j-1)}}(Y_{j,b,i})}{p_{\breve{\theta}^{\circ}_{b}|Y^{(j-1)}}(Y_{j,b,i})}$
(14)

is an e-variable, as long as $\Theta_{0}\subset[0,1]^{2}$ is convex and
closed. Here $(\breve{\theta}^{\circ}_{a}\mid
Y^{(j-1)},\breve{\theta}^{\circ}_{b}\mid Y^{(j-1)})$ is defined to minimize KL
divergence, i.e. is the pair $(\theta_{a},\theta_{b})\in\Theta_{0}$ that
minimizes, over $\Theta_{0}$,
$D(P_{\breve{\theta}_{a}\mid Y^{(j-1)},\breve{\theta}_{b}\mid
Y^{(j-1)}}(Y^{n_{a}}_{a},Y^{n_{b}}_{b})\|P_{\theta_{a},\theta_{b}}(Y^{n_{a}}_{a},Y^{n_{b}}_{b}))$.
(3) is a special case since with
$\Theta_{0}=\\{(\theta,\theta):\theta\in[0,1]\\}$, this KL divergence is
minimized by $(\breve{\theta}^{\circ}_{0},\breve{\theta}^{\circ}_{0})$ with
$\breve{\theta}^{\circ}$ as defined underneath (3). Again,
$\breve{\theta}_{a}$ and $\breve{\theta}_{b}$ are estimated based on past data
$Y^{(j-1)}$ as in (3). Based on (14) one can construct an exact
(nonasymptotic) confidence sequence (CS)

$\text{\sc
CS}_{\alpha,(m)}=\left\\{\delta:E^{(m)}_{[\Theta_{0}(\delta)]}\leq\frac{1}{\alpha}\right\\},$
(15)

with $\Theta_{0}(\delta)\subset[0,1]^{2}$ a null hypothesis determined by a
divergence measure. By construction, such a confidence sequence is always-
valid (Ramdas et al.,, 2022) in the sense that for any $\delta$, any
$\theta\in\Theta_{0}(\delta)$, the $P_{\theta}$-probability that there will
ever be an $m$ such that $\delta\not\in\text{\sc CS}_{\alpha,(m)}$ is at most
$\alpha$. This means that we can take the _running intersection_ of the
confidence sequence while retaining coverage, which will be used throughout
the simulation experiments in this paper. In this paper, we are going to
construct confidence sequences for risk differences as examples, where we are
going to test hypotheses of the form
$\Theta_{0}(\delta):=\left\\{(\theta_{a},\theta_{b})\in[0,1]^{2}:\theta_{b}-\theta_{a}=\delta\right\\}$
— below we extend this to the case that differentiates in terms of the strata.
Still, everything could also easily be adapted to construct confidence
intervals for other divergence measures, such as odds and risk ratios (Turner
and Grünwald,, 2022).

###  3.1 One CS per stratum

If we expect the effect size values to differ between the strata, one could
decide to report a separate confidence sequence for each stratum using (15)
above. To reach a better estimate sooner, we could however still allow cross-
talk on control group success rates or odds ratios between subpopulations, as
described in section 2 above. In this setup, we would end up with a
_collection_ of $k$ confidence sequences:

$\text{\sc
CS}^{k}_{\alpha,(m)}=\left\\{\delta:E^{(m),k}_{[\Theta_{0}(\delta)]}\leq\frac{1}{\alpha}\right\\},$
(16)

with $\breve{\theta}_{a}$ and $\breve{\theta}_{b}$ in $E^{(m),k}$ estimated
based on data seen up to time $m$ and $E^{(m),k}$ defined as in (2) with
$S^{k}_{j}$ replaced by $S^{k}_{j,[\Theta_{0}]}$ as in (14), calculated for
stratum $k$. Illustrations of confidence intervals over time with the three
options for cross-talk are depicted in Figure 4. As can be observed there, not
allowing cross-talk gives the best results when the true data generating
distributions in the strata have different control group success rates and
odds ratios (see the circle-shaped points in Figure 4d, especially in the
third stratum, where the effect size has a different sign). However, when
control group rates or odds ratios are similar across strata, allowing cross-
talk improves results. See for example Figure 4e, where interval width
decreases much faster in the smaller stratum 1 while allowing cross-talk about
the control group rate. Similar experiments for comparing confidence sequences
with and without the mixture of cross-talk methods can be found in the
supplementary material, Figure S1.

###  3.2 CS for the minimum or maximum

In some scenarios, for example when we do not have the means to collect a
large data sample, or when data is very unbalanced in one or more strata, it
could be more informative to create one CS for the minimum or maximum effect
size value over all strata. To achieve this, we introduce two new forms of
null hypotheses and corresponding e-variables that will subsequently be
inverted to create two one-sided confidence sequences, for lower and upper
bounds on the minimum or maximum.

#### One-sided CS: upper bound

We will first illustrate how to estimate an upper bound on some minimal effect
size value over strata111Analogously, with this method a lower bound on some
maximal effect size value can be estimated by reversing all signs.. To this
end, we consider a null hypothesis of the form $\mathcal{H}_{0,\delta}:\forall
k:\theta_{k}\in\Theta_{0}(\geq\delta)$ (i.e. for risk difference effect size,
$\Theta_{0}(\geq\delta)=\\{(\theta_{a},\theta_{b})\in[0,1]^{2}:\theta_{b}-\theta_{a}\geq\delta\\}$)
and aim to design e-variables to test it. E.g. in the example depicted in
Figure 5(a), we aim to design an e-variable that will reject
$\mathcal{H}_{0,\delta^{\prime\prime}}$ at any batch $j$ with probability less
than $\alpha$ (i.e., that offers type-I error guarantee), when the data in the
strata are in reality generated by $(\theta_{a,1},\theta_{b,1})$ and
$(\theta_{a,2},\theta_{b,2})$. We do eventually want to reject
$\mathcal{H}_{0,\delta^{\prime}}$ as
$\delta((\theta_{a,2},\theta_{b,2}))<\delta^{\prime}$. As we collect more and
more data, we can reject null hypotheses corresponding to values of
$\delta^{\prime}$ for which
$\delta^{\prime}-\delta((\theta_{a,2},\theta_{b,2}))$ gets closer and closer
to $0$.

(a) Examples of parameter spaces for $\mathcal{H}_{0,\delta}:\forall
k:\theta_{k}\in\Theta_{0}(\geq\delta)$.

(b) Examples of parameter spaces for $\mathcal{H}_{0,\delta}:\exists
k:\theta_{k}\in\Theta_{0}(\leq\delta)$.

Figure 5: Parameter space examples for hypotheses tested to construct upper
and lower bounds on minima and maxima of effect size values

Let us denote the e-process consisting of the e-variables for testing
$\theta_{k}\in\Theta_{0}(\geq\delta)$ in each stratum combined, using any of
the methods described above in Section 2, as $E^{*(m)}_{\delta}$. The one-
sided confidence interval for the minimum effect can be defined as:

$\text{CS}^{+}_{\alpha,(m)}:=\left[-1,\min\left\\{\delta:E^{*(m)}_{\delta}\geq\frac{1}{\alpha}\right\\}\right].$
(17)

All possible approaches for combining e-variables from separate strata, as
described in Section 2 above, to find an upper bound for the minimal effect
size value are compared in the confidence intervals in the paragraph below.

Figure 6: Example of confidence sequences for the lower- (LB) and upper (UB)
bounds of the minimum effect. $30$ observations were made in each stratum, and
the real differences were $0.5$, $0.4$ and $0.05$. With the switch method, a
uniform prior ranging from $m_{\text{switch}}=5$ until $30$ was applied. With
the pseudo-Bayesian approach, the learning rate $\eta$ was set to $1$ and $2$.
$\alpha$ was set to $0.05$.

#### One-sided CS: lower bound

We now also aim to estimate a lower bound for the minimal effect size value
(or, analogously, an upper bound for the maximal effect size value). To
achieve this, we now consider a null hypothesis of the form
$\mathcal{H}_{0,\delta}:\exists k:\theta_{k}\in\Theta_{0}(\leq\delta)$.
Looking at Figure 5(b) as an example, where data are generated by
$(\theta_{a,1},\theta_{b,1})$ and $(\theta_{a,2},\theta_{b,2}),$ we aim to
design an e-variable that will reject $\mathcal{H}_{0,\delta^{\prime}}$ at any
batch $j$ with probability less than $\alpha$ (i.e., we again want type-I
error guarantee if $\mathcal{H}_{0,\delta^{\prime}}$ is true), as
$\delta((\theta_{a,2},\theta_{b,2}))<\delta^{\prime}$. We do want to reject as
quickly as possible $\mathcal{H}_{0,\delta^{\prime\prime}}$, as $\forall
k,\delta(\theta^{(k)})>\delta^{\prime\prime}$. As we collect more data, we can
reject null hypotheses with values of $\delta^{\prime\prime}$ for which
$\delta((\theta_{a,2},\theta_{b,2}))-\delta^{\prime\prime}$ gets closer and
closer to 0.

To build our one-sided confidence interval $\text{CS}^{-}_{\alpha,(m)}$, we
again want to construct a compound e-variable $E^{*(m)}_{\delta}$ testing the
null hypothesis corresponding to each value of $\delta$, but now take
$\max\\{\delta:E^{*(m)}_{\delta}\geq 1/\alpha\\}$ as our lower bound. To test
$\mathcal{H}_{0,\delta}$ we will use the _minimum_ of
$E^{(j),k}_{\Theta_{0}(\leq\delta)}$ over all $k$, which provides an
e-variable for ${\cal H}_{0,\delta}$. To see this, let us assume
$\mathcal{H}_{0,\delta}$ is true an that for some $k^{*}$,
$\theta_{k^{*}}\in\Theta_{0}(\leq\delta)$; the other data generating
distributions might or might not come from $\Theta_{0}(\leq\delta)$. Then:
$\mathbb{E}(\min_{k}S^{k})\leq\min_{k}\mathbb{E}(S^{k})\leq\mathbb{E}(S^{k^{*}})\leq
1.$

#### Combining into confidence interval

We now combine the lower bound and upper bound estimation methods established
above to build confidence intervals for the minimal effect size value. This
can be achieved through taking the intersection of the one-sided confidence
sequences introduced above:
$\text{CS}_{\alpha,(m)}:=\text{CS}^{-}_{\alpha,(m)}\cap\text{CS}^{+}_{\alpha,(m)}.$
Results from an experiment where in one of the strata the treatment effect was
substantially smaller than in the others are depicted in Figure 6 (with
average interval widths in the supplementary material, S3). In early phases of
data collection, multiplication gives the quickest convergence, but as more
data is collected, the “sequential learning” methods converge quicker. When
risk differences where about the same across all strata, multiplication
converged the quickest (see Figure S4 in the Supplementary material).

###  3.3 CS for the mean effect

In addition to estimating the minimum or maximum effect in one of the strata,
one might be interested in estimating the mean effect an intervention will
have on an entire population, given the existence of subpopulations. For
example, one might want to estimate the effect a vaccination will have on the
probability of people being contaminated with a disease, taking into account
that a certain proportion of the population concerns elderly or
immunocompromised citizens.

Assuming we have a trustworthy estimate of the proportion of subjects
belonging to each stratum $k$ in the population of interest, $\pi_{k}$, we aim
to estimate the mean risk difference (mean expected effect of the
intervention) $\delta^{*}:=\sum_{k}\pi_{k}\delta_{k}$. We can build a
confidence sequence for $\delta^{*}$ by constructing an e-variable for the set
of all possible success probability distributions satisfying this
$\delta^{*}$,
$\mathcal{H}_{0,\delta^{*}}:\\{P_{\vec{\theta}};d(\vec{\theta})=\sum_{k}\pi_{k}d((\theta_{a,k},\theta_{b,k}))=\delta^{*}\\}$.
It is not directly clear what an optimal e-variable could look like; one
option that offers both the type-I error guarantee with potentially good power
is to combine the growth-rate optimal e-variable (3) for a specific
$\delta_{k}$ in each stratum with the _universal inference_ (Wasserman et
al.,, 2020) method for designing e-processes. Based on this strategy, we look
at the set of all vectors $\vec{\delta}:=(\delta_{1},...,\delta_{K})$ that
satisfy $\sum_{k}\pi_{k}\delta_{k}=\delta^{*}$. For one member of the set, we
can calculate the e-variable _based on all batches of data seen up to and
including time $m$_ according to (3):

$E^{(m)}_{[\vec{\delta}]}=\prod_{k}E^{(m),k}_{[\Theta_{0}(\delta_{k})]},$

where $\smash{E^{(m),k}_{[\Theta_{0}(\delta_{k})]}}$ can be calculated using
estimates for $\breve{\theta}_{a,k}$ and $\breve{\theta}_{b,k}$ as before,
only including data seen up to _and not including_ batch $m$. The e-variable
for $\mathcal{H}_{0,\delta^{*}}$ can then be calculated as (Wasserman et al.,,
2020):
$E^{*(m)}_{\delta^{*}}=\min_{\vec{\delta}}\smash{E^{(m)}_{[\vec{\delta}]}}$,
and the corresponding confidence sequence can be constructed as before,
analogously to (17).

Figure 7: Simulated example of $95$% confidence sequences for the mean effect
across subpopulations. $25$ observations were made in each stratum, and the
real risk difference of $0.4$ was homogeneous across subpopulations. The
confidence sequence for the mean effect is plotted alongside the Miettinen-
Nuninen confidence interval, a fixed-n confidence interval method, at batch
number $50$ (the purple triangles). In the supplementary materials, figure S5,
the mean effect CS is further illustrated for heterogeneous risk differences
in strata.

#### Comparison to fixed-n CI for Mantel-Haenszel risk difference

Much of the research into estimating stratified risk differences with coverage
guarantee has considered Mantel-Haenszel risk differences, where risk
differences or odds ratios are _homogeneous_ across strata but control group
rates can vary (see for example (Qiu et al.,, 2019)), with fixed-n designs.
This is a strong assumption, and we do not make it ourselves; but we can use
cross-talk on the risk difference to tailor our confidence sequences so that
they adapt (get narrow) if the risk difference is indeed homogeneous. One
recent fixed-n approach for this setting was described and implemented by
Klingenberg, (2014). In Figure 7, our confidence sequence for the mean effect
is compared to the Miettinen-Nuninen (MN) confidence interval from
Klingenberg, (2014) at fixed time $50$ in a setting where risk differences
were homogeneous. The MN-interval is slightly narrower, but because we are
allowed to continuously monitor the confidence interval while retaining
coverage with the confidence sequence, we can exclude 0 from the CS
considerably earlier than with the fixed-n method — which is remarkable
because unlike the MN fixed-n confidence interval, our anytime-valid
confidence sequences are also valid if in fact risk differences are not
homogeneous.

##  4 CONCLUSION AND FUTURE WORK

We have introduced a new method for global null hypothesis testing and
constructing exact anytime-valid confidence sequences in stratified count
data. Our method is complementary to previously proposed methods for similar
settings as we need no stochastic assumptions about the arrival times of the
subgroups or strata, and no Model-X assumptions. We have shown that our tests
and estimates are efficient in terms of power, and that precise effect size
estimations can be reached with less strong model assumptions compared to pre-
existing fixed-n methods, while retaining coverage guarantees and allowing
sequential decision making. We have also shown that we can improve the
traditional model of global null testing in the CMH-setting through
incorporating ideas from machine-learning: allowing for cross-talk between
strata, and incorporating pseudo-Bayesian learning and switching between
strata for learning compound effect measures.

Our work extends that of Turner et al., (2021) and Turner and Grünwald, (2022)
to incorporate strata for count data. Their methods, however, are generally
implementable for any convex null hypothesis, and future work should explore
if they also can feasibly be extended to stratified sequential effect
estimation for continuous outcome variables.

###  4.1 Acknowledgements

This work is part of the Enabling Personalized Interventions (EPI) project,
which is supported by the Dutch Research Council (NWO) in the Commit2Data
–Data2Person program, contract 628.011.028.

## References

  * Cesa-Bianchi and Lugosi, (2006) Cesa-Bianchi, N. and Lugosi, G. (2006).  Prediction, Learning and Games.  Cambridge University Press, Cambridge, UK.
  * Duan et al., (2022) Duan, B., Ramdas, A., and Wasserman, L. (2022).  Interactive rank testing by betting.  In Proceedings of the First Conference on Causal Learning and Reasoning, volume 177 of Proceedings of Machine Learning Research, pages 201–235.
  * Grünwald, (2007) Grünwald, P. (2007).  The Minimum Description Length Principle.  MIT Press, Cambridge, MA.
  * Grünwald et al., (2022) Grünwald, P., de Heide, R., and Koolen, W. (2022).  Safe testing.  accepted, pending minor revision, for publication in Journal of the Royal Statistical Society: Series B.
  * Grünwald et al., (2022) Grünwald, P., Henzi, A., and Lardy, T. (2022).  Anytime valid tests of conditional independence under model-x.  arXiv preprint arXiv:2209.12637.
  * Herbster and Warmuth, (1998) Herbster, M. and Warmuth, M. K. (1998).  Tracking the best expert.  Machine learning, 32(2):151–178.
  * Klingenberg, (2014) Klingenberg, B. (2014).  A new and improved confidence interval for the mantel–haenszel risk difference.  Statistics in Medicine, 33(17):2968–2983.
  * Koolen and de Rooij, (2013) Koolen, W. M. and de Rooij, S. (2013).  Universal codes from switching strategies.  IEEE Transactions on Information Theory, 59(11):7168–7185.
  * Koolen and van Erven, (2010) Koolen, W. M. and van Erven, T. (2010).  Freezing and sleeping: Tracking experts that learn by evolving past posteriors.  arXiv preprint arXiv:1008.4654.
  * Lindon and Malek, (2020) Lindon, M. and Malek, A. (2020).  Anytime-valid inference for multinomial count data.  arXiv preprint arxiv:2011.03567.
  * Ly et al., (2022) Ly, A., Turner, R., and Ter Schure, J. (2022).  R-package safestats.  CRAN.
  * Pandeva et al., (2022) Pandeva, T., Bakker, T., Naesseth, C. A., and Forré, P. (2022).  E-valuating classifier two-sample tests.  arXiv preprint arxiv:2210.13027.
  * Qiu et al., (2019) Qiu, S.-F., Poon, W.-Y., Tang, M.-L., and Tao, J.-R. (2019).  Construction of confidence intervals for the risk differences in stratified design with correlated bilateral data.  Journal of Biopharmaceutical Statistics, 29(3):446–467.
  * Ramdas et al., (2022) Ramdas, A., Grünwald, P., Vovk, V., and Shafer, G. (2022).  Game-theoretic statistics and safe anytime-valid inference.  arXiv preprint arXiv:2210.01948.
  * Ramdas et al., (2020) Ramdas, A., Ruf, J., Larsson, M., and Koolen, W. (2020).  Admissible anytime-valid sequential inference must rely on nonnegative martingales.  arXiv preprint arXiv:2009.03167.
  * Shaer et al., (2022) Shaer, S., Maman, G., and Romano, Y. (2022).  Model-free sequential testing for conditional independence via testing by betting.  arXiv preprint arXiv:2210.00354.
  * Turner and Grünwald, (2022) Turner, R. and Grünwald, P. (2022).  Exact anytime-valid confidence intervals for contingency tables and beyond.  arXiv preprint arxiv:2203.09785.
  * Turner et al., (2021) Turner, R., Ly, A., and Grünwald, P. (2021).  Generic e-variables for exact sequential k-sample tests that allow for optional stopping.  arXiv preprint arxiv:2106.02693.
  * Turner, (2023) Turner, R. J. (2023).  safeSequentialTestingAISTATS2023.  Code corresponding to this AISTATS Paper, accessible at https://github.com/rosanneturner/safeSequentialTestingAISTATS2023.
  * van Erven et al., (2007) van Erven, T., Grünwald, P., and de Rooij, S. (2007).  Catching up faster in bayesian model selection and model averaging.  In Advances in Neural Information Processing Systems, volume 20.
  * Vovk and Wang, (2021) Vovk, V. and Wang, R. (2021).  E-values: Calibration, combination, and applications.  Annals of Statistics.
  * Wasserman et al., (2020) Wasserman, L., Ramdas, A., and Balakrishnan, S. (2020).  Universal inference.  Proceedings of the National Academy of Sciences, 117(29):16880–16890.

##  Appendix S1 PROOFS

###### Proof.

(of theorem 2.1). First consider the basic case with $E^{(m)}$ as in (8). As
we show below, we have, with ${\bf E}\equiv{\bf E}_{P_{\theta^{*}}}$,

$\displaystyle{\bf E}\left[\log E^{(m)}\right]={\bf E}\left[\sum_{j=1}^{m}\log
S_{j}\right]={\bf
E}\left[\sum_{j=1..m}\sum_{x\in\\{a,b\\}}\sum_{i=1..n_{x}}\log\frac{p_{\breve{\theta}_{x,k_{j}}|Y^{(j-1)}}(Y_{j,x,i})}{p_{\breve{\theta}_{0,k_{j}}|Y^{(j-1)}}(Y_{j,x,i})}\right]\geq$
$\displaystyle{\bf
E}\left[\sum_{j=1..m}\sum_{x\in\\{a,b\\}}\sum_{i=1..n_{x}}\log\frac{p_{\breve{\theta}_{x,k_{j}}|Y^{(j-1)}}(Y_{j,x,i})}{p_{\tilde{\theta}_{0,k_{j}}}(Y_{j,x,i})}\right]\geq{\bf
E}\left[\sum_{\begin{subarray}{c}j=1..m\\\ x\in\\{a,b\\}\\\
i=1..n_{x}\end{subarray}}\log\frac{p_{\theta^{*}_{x,k_{j}}}(Y_{j,x,i})}{p_{\tilde{\theta}_{0,k_{j}}}(Y_{j,x,i})}-\sum_{\begin{subarray}{c}k=1..K\\\
x\in\\{a,b\\}\end{subarray}}\log\left(n_{x}m_{k}\right)\right]+O(1)=$
$\displaystyle\sum_{k=1..K}m_{k}\cdot
D(P_{\theta^{*}_{a,k},\theta^{*}_{b,k}}\|P_{\tilde{\theta}_{0,k},\tilde{\theta}_{0,k}}))+O(\log
m)$ (18)

where we use notation
$D(P_{\theta^{*}_{a},\theta^{*}_{b}}\|P_{\theta_{0},\theta_{0}})$ as in (2);
and $\tilde{\theta}_{0,k}$ is defined as
$\arg\min_{\theta\in[0,1]}D(P_{\theta^{*}_{a,k},\theta^{*}_{b,k}}\|P_{\theta,\theta})$
which by the same calculation as the one leading up to (2, is given by
$\tilde{\theta}_{0,k}=(n_{a}/n)\theta^{*}_{a,k}+(n_{b}/n)\theta^{*}_{b,k}$,
and $m_{k}$ denotes the number of times that an instance of block $k$ was
observed in the first $m$ blocks, and we remind the reader that $+O(\log m)$
may also indicate a negative difference of order $\log m$. (S1) immediately
implies the result, using (6).

The first two equalities in (S1) are immediate. The first inequality follows
because $P_{\tilde{\theta}_{0,k_{j}},\tilde{\theta}_{0,k_{j}}}$ minimizes KL
divergence to $P_{\theta^{*}_{a,k_{j}},\theta^{*}_{b,k_{j}}}$ among all
$\theta\in[0,1]$, within each block $j$. The final equality follows by
independence and basic calculus. It remains to show the second inequality.
This one follows because we use a prior $W(\theta_{a,k},\theta_{b,k}$ under
which $\theta_{a}$ and $\theta_{b}$ are independently beta distributed with
strictly positive densities on $(0,1)$. We can then use a standard Laplace
approximation of the Bayesian marginal likelihood to obtain, for each fixed
$k\in\\{1,\ldots,K\\}$, where the expectation ${\bf E}$ is over
$Y^{\prime}_{(1)},\ldots,Y^{\prime}_{(m^{\prime})}\sim
P_{\theta^{*}_{a,k},\theta^{*}_{b,k}}$:

$\displaystyle{\bf
E}\left[-\log\prod_{j=1}^{m^{\prime}}\prod_{x\in\\{a,b\\}}\prod_{i=1}^{n_{x}}{p_{\breve{\theta}_{x,k}|Y^{(j-1)}}(Y_{j,x,i})}\right]={\bf
E}\left[-\log\left(\int\prod_{j=1}^{m^{\prime}}\prod_{x\in\\{a,b\\}}\prod_{i=1}^{n_{x}}{p_{\theta_{x,k}}(Y_{j,x,i})}\right)dW(\theta_{a,k},\theta_{b,k})\right]$
$\displaystyle\leq{\bf E}\left[\sum_{j=1}^{m^{\prime}}-\log
p_{\theta^{*}_{a,k},\theta^{*}_{b,k}}(Y_{(j)})\right]+\log(n_{a}+n_{b})m^{\prime}+O(1).$

Here the equality is standard telescoping of the Bayesian marginal likelihood,
and the inequality is the Laplace approximation, i.e. the same calculation as
the one leading up to the $(d/2)\log n$ BIC approximation of Bayesian marginal
likelihood for a $d$-parameter exponential family; here $d=2$ since we have
two free parameters, $\theta^{*}_{a,k}$ and $\theta^{*}_{b,k}$; see
(Grünwald,, 2007, Chapter 8) for proof and detailed explanation).

This shows the result for the basic case that $E^{(m)}$ is arrived at by
multiplication, (8). The case for $E^{(m)}_{\textsc{mix}}$ follows similarly
by noting that, by construction, $E^{(m)}_{\textsc{mix}}\geq
E^{(m)}_{\textsc{none}}/3$, where $E^{(m)}_{\textsc{none}}$ denotes the
standard e-process with multiplication and without cross-talk, for which we
have already (just) shown the result. ∎

##  Appendix S2 ADDITIONAL EXPERIMENTS

(a) all different

(b) same control group rate

(c) same OR

(d) all different

(e) same control group rate

(f) same OR

Figure S1: Examples of $95$% stratified confidence intervals ((a), (b) and
(c)) and mean confidence interval widths estimated over $100$ runs ((d), (e)
and (f)) with different types of cross-talk, including mixing different types
of cross-talk. In (a), (b) and (c) the true risk difference of the data
generating distribution in each stratum is indicated by a dashed line. For (a)
and (d), the data were generated by distributions with different control group
success rates ($0.1$, $0.2$ and $0.8$) and risk differences ($0.05$, $0.4$ and
$-0.6$) in each stratum. For (b) and (e), strata sizes were unbalanced: as can
be seen for stratum 1, the red points, data collection stopped after $10$
batches. Control group success rates were all $0.5$ and risk differences were
different ($-0.49$, $-0.25$ and $0.1$). For (c) and (f), strata sizes were
unbalanced as well, and now odds ratios were the same in each stratum ($2$),
but control group rates differed again ($0.2$, $0.25$ and $0.85$).

(a) Upper bound sequence example

(b) Average difference with true minimum

Figure S2: Example of a confidence sequence and average difference from upper
bound to true minimal effect size value through 100 simulations, for different
switch priors on $j^{*}$. $30$ observations were made in each stratum, and the
real differences were $0.5$, $0.4$ and $0.05$. For the priors on early switch
times, all prior mass was distributed between batch numbers $5$ up to
$10$.$\alpha$ was set to $0.05$. Figure S3: Average interval width (upper
bound for the respective methods minus lower bound estimated with the minimum
method) of confidence sequences for the lower- (LB) and upper (UB) bounds of
the minimum effect and estimated through 100 simulations. $30$ observations
were made in each stratum, and the real differences were $0.5$, $0.4$ and
$0.05$. With the switch method, a uniform prior ranging from $j^{*}=5$ until
$30$ was applied. With the pseudo-Bayesian approach, the learning rate $\eta$
was set to $1$ and $2$. $\alpha$ was set to $0.05$.

(a) Confidence sequence example

(b) Average width

Figure S4: Example of confidence sequences for the lower- (LB) and upper (UB)
bounds of the minimum effect, and average interval width (upper bound for the
respective methods minus lower bound estimated with the minimum method). $30$
observations were made in each stratum, and the real differences were $0.4$,
$0.4$ and $0.5$. With the switch method, a uniform prior ranging from
$m_{\text{switch}}=5$ until $30$ was applied. With the pseudo-Bayesian
approach, the learning rate $\eta$ was set to $1$ and $2$. $\alpha$ was set to
$0.05$. Figure S5: Simulated example of a confidence sequence for the mean
effect across subpopulations. $25$ observations were made in each stratum, and
the real risk differences were $0.2$ and $0.5$. The confidence sequence for
the mean difference is plotted alongside the confidence sequence for the
minimum of the differences, estimated with pseudo-Bayesian averaging and a
uniform switch prior. $\alpha$ was set to $0.05$.