# Running in circles: is practical application feasible for data fission and
data thinning in post-clustering differential analysis?

Benjamin Hivert
Univ. Bordeaux, INSERM, INRIA, SISTM team, BPH,
Vaccine Research Institute, VRI, Hôpital Henri Mondor,
Denis Agniel
RAND Corporation,
Rodolphe Thiébaut
Univ. Bordeaux, INSERM, INRIA, SISTM team, BPH,
CHU Bordeaux, Service d’information médicale,
Vaccine Research Institute, VRI, Hôpital Henri Mondor,
and
Boris P. Hejblum
Univ. Bordeaux, INSERM, INRIA, SISTM team, BPH,
Vaccine Research Institute, VRI, Hôpital Henri Mondor

###### Abstract

The standard pipeline to analyse single-cell RNA sequencing (scRNA-seq) often
involves two steps : clustering and Differential Expression Analysis (DEA) to
annotate cell populations based on gene expression. However, using clustering
results for data-driven hypothesis formulation compromises statistical
properties, especially Type I error control. Data fission was introduced to
split the information contained in each observation into two independent parts
that can be used for clustering and testing. However, data fission was
originally designed for non-mixture distributions, and adapting it for
mixtures requires knowledge of the unknown clustering structure to estimate
component-specific scale parameters. As components are typically unavailable
in practice, scale parameter estimators often exhibit bias. We explicitly
quantify how this bias affects subsequent post-clustering differential
analysis Type I error rate despite employing data fission. In response, we
propose a novel approach that involves modeling each observation as a
realization of its distribution, with scale parameters estimated non-
parametrically. Simulations study showcase the efficacy of our method when
component are clearly separated. However, the level of separability required
to reach good performance presents complexities in its application to real
scRNA-seq data.

Keywords: Unsupervised learning, Mixture Model, Post-clustering inference,
Type I error, Non-parametric estimation, local variance

##  1 Introduction

Clustering encompasses all unsupervised statistical methods that group
observations into homogeneous and separated clusters. Widely used in various
application fields such as biology or genomics (Jaeger & Banks, 2023),
clustering plays a significant role to uncover or summarise signals contained
in different kind of multivariate data. In the context of single-cell RNA-seq
(scRNA-seq) technologies which are high-throughput genomics techniques that
measure gene expression at the single-cell level, providing insights into
cellular heterogeneity and functional diversity within complex biological
tissues, cluster analysis is an integral component of the traditional pipeline
for data analysis (Amezquita et al., 2020). Combined with differential
expression analysis, it allows to identify and annotate cellular
subpopulations by simply grouping cells based on their gene expression.
Differential expression analysis tests each individual gene for differential
expression between two clusters, leading to the identification of marker genes
that could potentially serve as cell-type marker genes (Pullin & McCarthy,
2024).

However, such a two-step pipeline for post-clustering differential analysis
requires using the same data twice: first to estimate the clusters, and then
again to estimate the differences between them and perform significance
testing – a procedure sometimes refered to as “double-dipping” (Kriegeskorte
et al., 2009). In post-clustering differential analysis, it has been shown
that the primary risk of double dipping is to compromise the control of the
Type I error rate of otherwise well-calibrated testing procedure for
traditional differential analysis, leading to false positives (Zhang et al.,
2019). Fundamentally, uncertainties stemming from cluster analysis outcomes,
particularly related with determining the optimal number of clusters, could
create artificial differences between homogeneous groups of observations.
Failure to account for the double use of the data during the analysis may lead
differential analysis tools to erroneously identify those artificial
differences (Hivert, Agniel, Thiébaut & Hejblum, 2024). Although challenges
related to double dipping are increasingly studied, emphasizing their critical
importance, they remain unresolved in the context of scRNA-seq data analysis
(Lähnemann et al., 2020).

Leveraging the selective inference framework (Fithian et al., 2014), which
involves conditioning on the clustering event during the derivation of the
test statistic and the associated p-value, various methodological efforts have
been carried to address this double dipping issue in post-clustering inference
(Gao et al., 2022). Zhang et al. (2019) introduced the Truncated-Normal test
(TN-test), a four-step post-clustering differential analysis procedure that
involves splitting the data into two parts, followed by cluster analysis on
the first part of the data. Subsequently, a support vector machine (SVM)
classifier is applied to the clustered data to learn the clustering structure
and predict the cluster labels on the remaining data. Finally, a differential
analysis is conducted between the predicted clusters using a truncated-normal
test. This truncation, on either side of the hyperplane fitted by the SVM,
allows to correct for the double dipping issue. However, information loss may
occur as an inherent consequence of the data-splitting process. Moreover, the
method involves multiple steps, introducing complexity into the overall
analytical framework, and Song et al. (2023) highlighted the poor performances
of the TN-test in their numerical simulations. More recently, Hivert, Agniel,
Thiébaut & Hejblum (2024) and Chen & Gao (2023) have introduced selective
tests tailored to detect mean differences between two clusters on individual
variables. These methods explicitly condition on the clustering event to
derive their p-value, relying on the set of all perturbed data sets that would
yield the same partition when subjected to the same clustering algorithm.
Hivert, Agniel, Thiébaut & Hejblum (2024) rely on a a Monte Carlo approach to
suit any clustering algorithm, which comes at the expense of extensive
computational times. On the other hand, Chen & Gao (2023) explicitly describe
this set specifically for K-means and hierarchical clustering applied to the
squared distance matrix. Bachoc et al. (2023) introduced a selective tests
designed for convex clustering. All these proposed selective methods designed
for post-clustering inference either require specific clustering algorithms or
introduce new specific test statistics. This consequently increases their
complexity and makes their application less straightforward in the context of
scRNA-seq data analysis.

Leiner et al. (2023) have drawn inspiration from data splitting (which
effectively addresses overfitting issues in machine learning) to break free
from the selective inference framework. They propose a method called “data
fission”, wherein the information within each individual observation is split
into two independent parts. The first part could be used for cluster analysis,
and labeling observations on the second part simply by matching them to their
first counterpart. Differential analysis could then be conducted on the
remaining information, i.e. the second part. However, in the context of post-
clustering differential analysis, it is imperative that the two parts be
independent, as each analysis (cluster analysis and differential analysis)
must be performed independently to effectively prevent double dipping and the
associated inflation of Type I error. Data fission requires strong parametric
distributional assumptions, with only the Gaussian and Poisson (Neufeld et
al., 2024) distributions ensuring independence between the two fissioned
parts. Expanding upon this concept, Neufeld, Dharamshi, Gao & Witten (2023)
have generalized the fission process by introducing “data thinning”. Building
on the same foundational idea, they not only succeed in developing a process
capable of decomposing data into more than two parts, but also broaden the
spectrum of distributions where independence between each part is provided.
This includes the negative binomial distribution, widely used when modeling
RNA-seq data.

Data fission or data thinning, although initially appealing, present
significant limitations. First, these methods lack results and justification
when applied to mixture distributions, which are commonly used to model data
with a clustered structure (MacQueen et al., 1967). Consequently, the absence
of such results inherently assumes a global null hypothesis of no clusters in
the data. Additionally, these methods assume prior knowledge of parameters
(e.g., variance for the Gaussian distribution or overdispersion for the
negative binomial distribution). Although robust estimators for these
parameters could theoretically ensure the validity of the method, this further
adds complexity for clustered data where each cluster has a different
parameter value. And in the absence of knowledge about the data structure,
specifically the clusters, the only justifiable estimator is the full-sample
one (i.e. computed using all the observations regardless of their mixture
component), even though it fails to correctly estimate the intra-component
parameter value.

Data fission and data thinning have emerged as promising alternatives to
selective inference for post-clustering inference, owing to their
compatibility with any clustering methods and differential analysis tests. In
this article, we propose a strategy to address their inherent limitations by
modeling each observation as a realization from distinct random variables with
their own parameter value. This approach facilitates the application of
fission or thinning processes within the context of mixture distributions
although the accurate estimation of scale parameters still remains
challenging. We establish a link between the bias in the estimation of the
variance parameter in the Gaussian distribution and the expected Type I error
rate of the two-sample $t$-test (Welch, 1947), underscoring the critical
importance of a robust estimator and the practical complexities associated
with implementing these methods. To tackle this issue, we introduce a non-
parametric method for estimating the local individual variance to be used in
the Gaussian setting.

##  2 Methods

In the following, capital letters represent random variables with a
probability distribution function denoted as $p$, $x_{i}$ denotes a set of $n$
realizations of $X$, and matrices and vectors are bolden.

###  2.1 Data fission and data thinning

Let $\boldsymbol{X}$ be a random variable with a known distribution. Both data
fission (Leiner et al., 2023) and its generalization, data thinning (Neufeld,
Dharamshi, Gao & Witten, 2023), aim to decompose the random variable
$\boldsymbol{X}$ into two (or more in the case of data thinning) new random
variables $\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$. These new
variables are designed to i) retain information from the original variable
$\boldsymbol{X}$, and ii) be independent. Of note, data fission can generate
pairs of $\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$ that are not
independent, but her we will focus on the independent cases only. The balance
between the amount of information from $\boldsymbol{X}$ kept into either
$\boldsymbol{X}^{(1)}$ or $\boldsymbol{X}^{(2)}$ is tuned by a hyperparameter
$\tau$. Such a decomposition can be performed for various probability
distributions of $\boldsymbol{X}$. For data fission, Leiner et al. (2023)
identified only two distributions, the Gaussian and the Poisson, that satisfy
the independence requirement. In contrast, Neufeld, Dharamshi, Gao & Witten
(2023) established a comprehensive decomposition that ensures independence for
all convolution-closed distributions. This encompasses Gaussian, Poisson
(Neufeld et al., 2024) and Negative Binomial (Neufeld, Popp, Gao, Battle &
Witten, 2023) distributions.

Distribution of $\boldsymbol{X}$ | $\boldsymbol{\tau}$ | Data fission | Data thinning
---|---|---|---
$\mathcal{P}(\lambda)$ | $\tau\in[0,1]$ |  | $Z\sim\operatorname{Binom}(X,\tau)$
---
$X^{(1)}=Z$
$X^{(2)}=X-Z$
| $X^{(1)}\,|\,X=x\sim\operatorname{Binom}(x,\tau)$
---
$X^{(2)}=X-X^{(1)}$
$\mathcal{N}(\boldsymbol{\mu}_{p},\boldsymbol{\Sigma}_{p\times p})$ | $\tau\in]0,+\infty)$ |  | $\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0}_{p},\boldsymbol{\Sigma}_{p\times p})$
---
$\boldsymbol{X}^{(1)}=\boldsymbol{X}+\tau\boldsymbol{Z}$
$\boldsymbol{X}^{(2)}=\boldsymbol{X}-\frac{1}{\tau}\boldsymbol{Z}$
|
$\boldsymbol{X}^{(1)}\,|\,\boldsymbol{X}=\boldsymbol{x}\sim\mathcal{N}(\tau\boldsymbol{x},\tau(1-\tau)\boldsymbol{\Sigma}_{p\times
p})$
---
$\boldsymbol{X}^{(2)}=\boldsymbol{X}-\boldsymbol{X}^{(1)}$
$\mathcal{NB}(\mu,\theta)$ | $\tau\in[0,1]$ | No fission |  | $X^{(1)}\,|\,X=x\sim\text{BetaBin}(x,\tau\theta,(1-\tau)\theta)$
---
$X^{(2)}=X-X^{(1)}$

Table 1: Data fission and data thinning decompositions for three usual
distributions: Poisson, Gaussian and Negative Binomial.

Decompositions for the three distributions into two independent random
variables $\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$, are detailed in
Table 1. For proofs of independence, please refer to the Appendix A for the
Gaussian case, and to Neufeld et al. (2024) and Neufeld, Popp, Gao, Battle &
Witten (2023) for the Poisson and negative binomial distribution,
respectively. Both Gaussian and negative binomial data fissions/thinnings
require knowledge of scale parameters (namely $\boldsymbol{\Sigma}$ or
$\theta$) for practical feasibility, as highlighted in Table 1. Theoretical
guarantees, and especially the independence between $\boldsymbol{X}^{(1)}$ and
$\boldsymbol{X}^{(2)}$, are based on using the true values of these
parameters. Yet, in real-life applications these are unknown and need to be
estimated, most likely from the data as described in Neufeld, Popp, Gao,
Battle & Witten (2023).

###  2.2 Limits in practical application of data fission an data thinning

####  2.2.1 Mixture distributions

Neufeld et al. (2024) and Neufeld, Popp, Gao, Battle & Witten (2023) have
proposed the application of data thinning for post-clustering inference.
Existing decompositions are currently limited to non-mixture distributions.
Unfortunately, these distributions can only describe a global null hypothesis,
assuming the absence of separated clusters within the data. In scenarios with
true clusters, a more appropriate modeling approach involves the use of
mixture models (Bouveyron et al., 2019) – each component of the mixture
representing a distinct cluster. Given $n$ i.i.d. realizations of random
variables following a mixture of distributions, the density function for an
observation $\mathbf{x}_{i}$ is:
$p(\mathbf{x}_{i})=\sum\limits_{k=1}^{K}\pi_{k}f(\mathbf{x}_{i}|\boldsymbol{\theta}_{k})$,
where $\pi_{k}$ is the probability that an observation was generated by the
$k\textsuperscript{th}$ component and $f(.|\boldsymbol{\theta}_{k})$ is the
density of the $k\textsuperscript{th}$ component with its specific parameters
$\boldsymbol{\theta}_{k}$ (for simplicity, we only consider the case where all
components belong to the same parametric distribution with density $f$). Data
fission and data thinning assume a homogeneous distribution for all
observations, so they are only applicable to one component. In the case of a
mixture, they can thus only be applied to each component separately.

Figure 1: Schematic view illustrating the circularity induced by data fission
for post-clustering differential analysis

Data fission and data thinning therefore create a circular situation, as
illustrated in the schematic Figure 1. To be applied for post-clustering
inference, these methods require knowledge of the intra-component parameters
$\boldsymbol{\theta}_{k}$ which, in turn, depend themselves on the components
that we would estimate. An effective approach could be to estimate a global
parameter $\hat{\boldsymbol{\theta}}$ based on all observations. However, this
would make the critical assumption that the parameter value is the same across
all components, i.e. $\boldsymbol{\theta}=\boldsymbol{\theta}_{k}$; this means
that the data are distributed according to the same distribution globally,
regardless of their component:

$p(\mathbf{x}_{i})=f(\mathbf{x}_{i}|\boldsymbol{\theta}_{k})=f(\mathbf{x}_{i}|\boldsymbol{\theta}).$
(1)

This can thus only be verified for $K=1$, so under the global null of no
cluster. Yun & Foygel Barber (2023) highlighted similar challenges for
existing post-clustering selective tests.

####  2.2.2 Scale parameter prior knowledge and estimation

We will restrict our analysis to the Gaussian setting, that is considering:

$f\left(\mathbf{x}_{i}|\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k}\right)=\frac{1}{(2\pi)^{p/2}|\boldsymbol{\Sigma}_{g}|^{1/2}}\exp\left(-\frac{1}{2}(\mathbf{x}_{i}-\boldsymbol{\mu}_{k})^{T}\boldsymbol{\Sigma}_{k}^{-1}(\mathbf{x}_{i}-\boldsymbol{\mu}_{k})\right)$

with $k=1,\dots,K$, $K\geq 1$ the number of components in the mixture, and
$\boldsymbol{\mu}_{k}\in\mathbb{R}^{p}$ and
$\boldsymbol{\Sigma}_{k}\in\mathbb{R}^{p\times p}$ respectively the mean
vector and the covariance matrix of the $k^{\textsuperscript{th}}$ components.
Much of the results apply to the negative binomial case, with the
overdispersion parameter being analogous to the variance parameter in the
Gaussian case. In this setting, the challenge lies in estimating
$\boldsymbol{\Sigma}_{k}$ for data fission or data thinning.

Figure 2: Toy example illustrating the impact of variance estimation on data
fission p-values Panel A: A two-dimensional Gaussian distribution incorrectly
clustered into $2$ clusters. Panel B: QQ-plot of the $t$-test p-values for the
comparison between the two estimated clusters across $1,\\!000$ simulations
when data fission is performed with $3$ different variance estimators. Panel
C: Extension of the problem to a two-component, two-dimensional Gaussian
mixture incorrectly clustered into $3$ clusters. $C_{1}$ and $C_{2}$ originate
from the same component, which is erroneously split into two. Panel D:
$t$-test p-values for the comparison between $C_{1}$ and $C_{2}$ over
$1,\\!000$ data simulations using the same $3$ variance estimators for the
data fission.

Figure 2 provides a comprehensive overview of the challenges associated with
variance estimation in data fission. Panel A presents an illustrative example
involving $300$ realizations of a multivariate Gaussian distribution ($K=1$)
with a mean vector $\boldsymbol{\mu}=\left(\begin{smallmatrix}0\\\
0\end{smallmatrix}\right)$ and a covariance matrix
$\boldsymbol{\Sigma}=\left(\begin{smallmatrix}4&0\\\
0&1\end{smallmatrix}\right)$. The k-means algorithm was applied to this
dataset, resulting in the estimation of two clusters, $C_{1}$ and $C_{2}$.
Data fission performance was evaluated when using different estimations of
$\boldsymbol{\Sigma}$. First, we considered the true intra-components
covariance matrix $\boldsymbol{\Sigma}_{k}$ (see section 2.3.1 for a proposed
approach to data fission assuming several mixture component). We also
considered the overall sample covariance matrix
$\widehat{\boldsymbol{\Sigma}}=\frac{1}{n-1}\sum_{i=1}^{n}\left(\boldsymbol{x}_{i}-\overline{\boldsymbol{X}}\right)\left(\boldsymbol{x}_{i}-\overline{\boldsymbol{X}}\right)^{t}$
where $\overline{\boldsymbol{X}}$ is the sample mean vector. We finally use
the $k-$means results to compute an intra-cluster covariance matrix defined
as: $\widehat{\boldsymbol{\Sigma}}_{\hat{k}}=\frac{1}{|C_{k}|-1}\sum_{i=1\in
C_{k}}\left(\boldsymbol{x}_{i}-\overline{\boldsymbol{X}}_{C_{k}}\right)\left(\boldsymbol{x}_{i}-\overline{\boldsymbol{X}}_{C_{k}}\right)^{t}$
where $\overline{\boldsymbol{X}}_{C_{k}}$ is the sample mean vector of the
cluster $C_{k}$. Since both clusters originate from the same component, there
is no inherent differences between them, implying that the $t$-test p-values
should exhibit a uniform distribution. Panel B presents a QQ-plot illustrating
the resulting p-values against the Uniform distribution for the test on
$X_{1}$ in $1,\\!000$ replications of the experiment. In this scenario where
the mixture has only one component ($K=1$), the sample covariance matrix
$\widehat{\boldsymbol{\Sigma}}$ provides an unbiased estimate of the true
$\boldsymbol{\Sigma}$, resulting in uniformly distributed p-values. However,
when considering the intra-cluster covariance matrices
$\widehat{\boldsymbol{\Sigma}}_{\hat{k}}$ (for $k=1,2$) derived from the
k-means results, the p-values no longer exhibit a uniform distribution. In
this case, the estimated matrices $\widehat{\boldsymbol{\Sigma}}_{\hat{k}}$
for each cluster drastically underestimate the true covariance matrix
$\boldsymbol{\Sigma}$, compromising the independence between
$\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$. This deviation from
independence indicates that it becomes straightforward to replicate the
clusters identified in $\boldsymbol{X}^{(1)}$ onto $\boldsymbol{X}^{(2)}$.

Panel C introduces a scenario with two true clusters generated using a mixture
of two Gaussian distributions ($K=2$). The K-means identifies $3$ clusters
(for illustration purposes), incorrectly splitting one mixture component into
clusters $C_{1}$ and $C_{2}$. We compared data fission performance was again
compared when the decomposition is performed using the same $3$ covariance
estimators as in the previous scenario. Panel D presents the resulting
p-values for the $t$-test on $C_{1}$ and $C_{2}$ on $X_{1}$ over $1,\\!000$
replications of the analysis. Efficient Type I error control is achievable
only by considering the true intra-components covariance matrix. Both global
and estimated intra-cluster covariances are biased estimators, leading to
correlations between $\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$,
compromising control of the Type I error rate. This underscores the intricate
challenges associated with covariance estimation for data fission in practical
scenarios, mainly due to the misspecification of the generative model for data
fission arising when genuine clusters exist within the dataset.

###  2.3 Practical solutions

####  2.3.1 Individual fission (or thining) for mixture models

In Gaussian mixture models, parameters are typically component-specific,
meaning that they are assumed to be shared across all individuals within each
component. As explained above and in Figure 1, this assumption poses a
circular challenge as it requires knowledge of the true data structure to be
able to accurately estimate the component-specific covariance matrices that
are then needed to perform data fission or data thinning. To address this
limitation, we propose an alternative approach in which each observation is
modeled as a realization of its own Gaussian distribution, i.e.
$p_{i}(\boldsymbol{x}_{i})=f\left(\boldsymbol{x_{i}}|\boldsymbol{\mu}_{i},\boldsymbol{\Sigma}_{i}\right)$.
Thus, observations are no longer assumed to be identically distributed).
Consequently, the covariance matrix $\boldsymbol{\Sigma}_{i}$ is no longer
specific to the components but instead to the individual observations. Despite
this individual-level definition, two individuals drawn from the same
component are expected to have very close variance parameters. By bypassing
the components in the definition of $\boldsymbol{\Sigma}_{i}$, this modeling
strategy theoretically encompasses both the global null ($K=1$) and the
mixture ($K\geq 1$) settings, and opens up individual-level data fission and
data thinning. As highlighted above, variance estimation remains crucial for
the practical feasibility of these methods. This new modeling assumes
individual variances that still need to be known (ideally), or precisely
estimated in real-life settings.

####  2.3.2 Non-parametric local variance estimation

To estimate $\boldsymbol{\Sigma}_{i}$, we propose to use weighted variances
where the weights are determined through non-parametric estimations using
kernel smoothing. The underlying principle behind this approach is that,
despite the individual-specific nature of the variance, two observations
within the same component of the mixture (i.e., two neighbors) should display
similar variance patterns. The non-parametric weights assigned to each
individual reflect their contribution to the estimation of the variance for
the $i\textsuperscript{th}$ observation, effectively capturing the proximity
between individuals. First, let’s assume that we are under the univariate
setting, that is: $X_{i}\sim\mathcal{N}\left(\mu_{i},\sigma^{2}_{i}\right)$
for $i=1,\dots,n$. We define $\widehat{\sigma^{2}}_{i}$, the resulting
estimate of $\sigma^{2}_{i}$, as:

$\widehat{\sigma^{2}}_{i}=\frac{\sum\limits_{j=1}^{n}w_{ij}\left(x_{i}-\widehat{m}_{i}\right)^{2}}{\left(\sum\limits_{j=1}^{n}w_{ij}\right)-1}$
(2)

where $w_{ij}$ are individual-weights and
$\widehat{m}_{i}=\sum\limits_{j=1}w_{ij}x_{i}\bigg{/}\sum\limits_{j=1}^{n}w_{ij}$
are individual-specific weighted means. Ideally, $w_{ij}$ should be zero (or
very small) for all the observations $x_{j}$ that are not in the same
component as $x_{i}$.

The catch of the variance computation in (2) lies in determining each
individual weight $w_{ij}$. Since it is essential for $w_{ij}$ to
appropriately capture the proximity between observations, we opted for a
kernel-based definition: $w_{ij}=K\left(x_{i}-x_{j}\right)$, where
$K(u)=\frac{1}{\sqrt{2\pi}h}e^{-\frac{u^{2}}{2h^{2}}}$ represents the uniform
kernel providing the proximity measurement between $x_{i}$ and $x_{j}$. This
ensures that individual weights reflect the local relationships within the
data. This kernel choice focuses on nearby points, emphasizing their influence
on the weighted variance estimate. The parameter $h$ in the definition of $K$
serves as the bandwidth parameter, controlling the width of the kernel and,
consequently, the neighborhood around each observation $x_{i}$ that
contributes to its weighted variance estimate. A smaller $h$ results in more
localized estimations, emphasizing nearby points, but can lead to
underestimate the variances. Conversely, a larger $h$ includes a broader range
of observations, and an excessively large value might consider almost all
observations in the variance estimation, yielding an estimator close to the
full sample one. Therefore, an optimal choice for $h$ would involve
considering only the observations in the same components of the underlying
mixture. However, achieving this ideal scenario is impractical as it is
equivalent to knowing the true components of the mixture.

Bandwidth calibration is a crucial step in any kernel method (Heidenreich et
al., 2013). We propose to use individual-specific bandwidth $h_{i}=h(x_{i})$
to reduce the bias in kernel density estimation. To accomplish this, we first
estimate the changepoint, $i^{*}$ in the spread of distances from observation
$x_{i}$ to all other observations. This changepoint delineates the change in
components of the mixture: observations with distances preceding this
changepoint are deemed part of the same component as $x_{i}$, whereas those
following it are considered part of different components. Subsequently, we
define $h_{i}=|x_{i}-x_{i^{*}}|$ as the distance between $x_{i}$ and the
observation $x_{i^{*}}$ that is furthest from $x_{i}$ before the break in the
mixture. Thus, $h_{i}$ is determined in such a way that the individual
bandwidth accommodates only the observations in the same component as $x_{i}$,
ensuring its size is sufficient to encompass all observations of the component
(Chacón & Duong, 2020).

##  3 Results

###  3.1 Quantification of Type I error rate as a function of the bias in the
variance estimation

Independence between $\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$ is
guaranteed solely when the true variance is used for decomposition.
Substituting the true variance $\boldsymbol{\Sigma}_{i}$ with an estimate
$\widehat{\boldsymbol{\Sigma}}_{i}$ introduces correlations between these new
random variables. Consider a Gaussian random variable
$\boldsymbol{X}\sim\mathcal{N}\left(\boldsymbol{\mu}_{p},\boldsymbol{\Sigma}_{p\times
p}\right)$ that we aim to decompose using the data fission or data thinning
processes outlined in Table 1, but with an estimated covariance matrix,
$\widehat{\boldsymbol{\Sigma}}$, in place of $\boldsymbol{\Sigma}$. Table 2
provides the covariance values between $\boldsymbol{X}^{(1)}$ and
$\boldsymbol{X}^{(2)}$ as a function of the bias of
$\widehat{\boldsymbol{\Sigma}}$ (for proofs, refer to Appendix B for data
fission and to Neufeld, Dharamshi, Gao & Witten (2023) for data thinning). If
$\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$ are not independent (as with
a biased estimator of $\boldsymbol{\Sigma}$ that induces significant
covariance), a split of a single component into two estimated clusters in
$\boldsymbol{X}^{(1)}$ is easily transferred to $\boldsymbol{X}^{(2)}$ and
results in false positives during the inference step (even if the latter is
carried out on $\boldsymbol{X}^{(2)}$).

| Data Fission | Data Thinning
---|---|---
$\operatorname{Cov}\left(\boldsymbol{X}^{(1)},\boldsymbol{X}^{(2)}\right)$ | $\boldsymbol{\Sigma}-\widehat{\boldsymbol{\Sigma}}$ | $\tau(1-\tau)\left(\boldsymbol{\Sigma}-\widehat{\boldsymbol{\Sigma}}\right)$
Table 2: Bias in the estimation of $\boldsymbol{\Sigma}$ induces correlations
between $\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$

To further describe the repercussions of variance estimation in the context of
data fission, we derived an analytical expression for the Type I error rate of
the $t$-test as a function of the bias in estimating this parameter. Let
$X\sim\mathcal{N}(\mu,\sigma^{2})$ be a random variable. Given that $X$ is not
generated from a mixture of distributions, it becomes evident that distinct
clusters are not inherent in the variable. Our aim is to perform data fission
in $X$, using $X^{(1)}$ for $k-$means clustering (with $K=2$) and $X^{(2)}$
for differential testing between the two inferred clusters. Notably, we employ
an estimate $\widehat{\sigma^{2}}$ of $\sigma^{2}$ in the fission process.
Given the generation process of $X$, it is established that regardless of the
clustering on $X^{(1)}$, there should be no mean difference between the
estimated clusters on $X^{(2)}$ as long as independence is achieved. Let
$\widehat{C}_{1}$ and $\widehat{C}_{2}$ be the two estimated clusters on
$X^{(1)}$ with the same intra-cluster variance, which is a reasonable
hypothesis with $k-$means clustering as explained in Appendix C. Since we are
under the null hypothesis of no mean difference between the clusters, the $T$
statistic for the $t$-test between $\widehat{C}_{1}$ and $\widehat{C}_{2}$ on
$X^{(2)}$ is given by:

$T_{\widehat{C}_{1},\widehat{C}_{2}}=\frac{\widehat{\delta}_{X^{(2)}}}{\sqrt{\frac{4s^{2}(X^{(2)})}{n}}}$
(3)

Here, $\widehat{\delta}_{X^{(2)}}$ represents the observed mean difference
between $\widehat{C}_{1}$ and $\widehat{C}_{2}$ on $X^{(2)}$ and
$s^{2}\left(X^{(2)}\right)$ is their shared intra-cluster variance. Of note,
$\widehat{\delta}_{X^{(2)}}$ and $s^{2}\left(X^{(2)}\right)$ both depend on
the ability to reproduce the clustering results obtained on $X^{(1)}$ into
$X^{(2)}$. This is quantified by their correlation
$\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)$, which can be easily computed
from their covariance given in Table 2. Based on this correlation, it can be
demonstrated that:

$T_{\widehat{C}_{1},\widehat{C}_{2}}\stackrel{{\scriptstyle\mathcal{H}_{0}}}{{\sim}}\mathcal{N}\left(\frac{\sqrt{n}\sqrt{\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}}{\sqrt{1-\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}},1\right)$

where
$\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)=\frac{\left(\sigma^{2}-\widehat{\sigma^{2}}\right)}{\sqrt{\left(\sigma^{2}+\tau^{2}\widehat{\sigma^{2}}\right)\times\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)}}$.
Appendix C gives details on the derivation of this test statistic and its
distribution. The associated Type I error rate for this test is given by
$1-F(q_{\alpha/2})+F(-q_{\alpha/2})$, where $F$ is the cumulative distribution
function of
$\mathcal{N}\left(\frac{\sqrt{n}\sqrt{\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}}{\sqrt{1-\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}},1\right)$
(or the corresponding non-central $\mathcal{T}(n-2)$ distribution), and
$q_{\alpha/2}$ is the quantile of a standard Gaussian distribution
$\mathcal{N}(0,1)$ (or of a $\mathcal{T}(n-2)$ distribution).

We validated this theoretical result through numerical simulations, and
conducted a detailed exploration of the influence of variance values and
sample size on the resulting Type I error rate. We generated $n$ realizations
of a Gaussian random variable with a mean $\mu=0$ and variances $\sigma^{2}$,
with $1,\\!000$ Monte Carlo repetitions. Subsequently, we used data fission
with varying values of $\widehat{\sigma^{2}}$, obtaining $X^{(1)}$ for k-means
clustering with $K=2$ and $X^{(2)}$ for testing mean differences between the
two estimated clusters. Initially, we examine the impact of the original true
variance $\sigma^{2}$ by considering $\sigma^{2}$ values of $\\{0.01,1,4\\}$
for a fixed sample size of $n=100$. Figure 3A, illustrates the relationship
between the bias in estimating $\sigma^{2}$ and the Type I error rate,
demonstrating a consistent agreement with the derived theoretical error rate.
We further documented the behavior of the Type I error for a fixed
$\sigma^{2}=1$ with varying sample sizes $n\in\\{50,100,250,100,1000\\}$ in
Figure 3B, which shows the expected impact of the sample size on the Type I
error rate. These findings collectively underscore the critical importance of
accurately estimating the variance to achieve a well-calibrated Type I error
rate with data fission.

Figure 3: Impact of variance estimation on Type I error rate in data fission.
Panel A: Evolution of the estimated Type I error rate in data fission as a
function of the relative bias and the original variance of the data. Panel B:
Evolution of the estimated Type I error rate in data fission as a function of
the relative bias and the sample size.

###  3.2 Performances of the local variance estimator

We conducted simulation studies to assess the performance of the non-
parametric variance estimator defined in (2). Generating univariate data akin
to the motivating example in Figure 2 with $n=100$ realizations from a two-
component univariate Gaussian mixture:
$0.5\mathcal{N}\left(0,\sigma^{2}\right)+0.5\mathcal{N}\left(\delta,\sigma^{2}\right)$.
We explored a range of ratio $\delta/\sigma$ values from $0$ (i.e. no
separation between the components, representing the global null of no clusters
in the data) to $100$ (indicating an extreme separation). Different values of
$\sigma^{2}$ were considered: $\sigma^{2}\in\\{0.01,1,4\\}$. For each pair
$\left(\delta,\sigma^{2}\right)$ defining the mixture, variance was estimated
from the data using our proposed weighted local variance, and the result was
used for individual data fission. We then considered only the observations
coming from the first component of the mixture for $k-$means clustering on the
corresponding $X^{(1)}$ with $K=2$. This resulted in one true component
splitting into two incorrect clusters. We then tested the mean differences
between those two clusters using the $t$-test on $X^{(2)}$. This scenario was
replicated $1,\\!000$ times, and we computed the empirical Type I error rate
at the $\alpha=5\%$ level.

Figure 4: Performance evaluation of the non-parametric variance estimator in a
simulated univariate setting. All empirical results were obtained through
$1,\\!000$ simulations of the data. Panel A: Median relative bias, defined as
$\left(\widehat{\sigma^{2}}-\sigma^{2}\right)/\sigma^{2}$ and its associated
inter-quantile range, against the signal-to-noise ratio $\delta/\sigma$. Panel
B: Type I error rate at the $\alpha=5\%$ level against $\delta/\sigma$. Panel
C: Median relative bias and its associated inter-quantile range as a function
of the sample size for three degrees of separation between the two components
(informed by the ratio $\delta/\sigma$). Panel D: Type I error rate at the
$\alpha=5\%$ level as a function of sample size.

Figure 4A shows that local variances are underestimated (manifested as a
negative relative bias) until the signal vs. noise ratio $\delta/\sigma$
reaches approximately $3$ – a threshold value for separation in a Gaussian
mixture model previously reported in the literature (Siffer et al., 2018;
Hivert, Agniel, Thiébaut & Hejblum, 2024). Consequently, there are no clear
change-points in the spread of distances for individual observations and all
observations should contribute significantly to the variance calculation.
However, this is not the case here (with a bandwidth $h_{i}$ that is too
small), leading to underestimated local variances. As the ratio increases
within the range $3.5\leq\delta/\sigma\leq 10$, component separation becomes
clearer but the identification of change-points remains challenging. Change-
points detection is intricate for observations in the tails in-between the
components, leading to non-zero weight for observations in both components.
Consequently, an overestimation of the local variances (positive relative
biases) is observed together with an increase in their associated inter-
quantile range. Finally, for $\delta/\sigma>10$, sufficient component
separation ensures consistent local variance estimation, highlighting the
critical importance of clear-cut component separation for accurate variance
estimation. Thus, the methodology outlined in 2.3.2 requires well-separated
observations to be able to provide accurate variance estimations. In Figure
4B, the Type I error rate at $5\%$ remains well calibrated only for values of
$\delta/\sigma$ that ensures robust variance estimations (that is
$\delta/\sigma>10$), emphasizing again the necessity of good component
separations for reliable testing. Figure 4C demonstrates how increasing the
sample size fails to improve non-parametric local variance estimations for
three representative values of the ratio $\delta/\sigma\in\\{0.5,5,10\\}$,
indicating that performance is more dependent upon component separation rather
than sample size. The corresponding Type I error rates depicted in Figure 4B
thus align with the observed relative bias results, underscoring the critical
role of accurate estimation of local variances to ensure Type I error control.

###  3.3 Application to single-cell RNA-seq data analysis

Single-cell RNA-seq (scRNA-seq) data analysis pipelines often involve an
initial cluster analysis followed by differential analysis to identify marker
genes and annotate cell populations based on gene expression. Given their
overdispersed count nature, the negative binomial distribution is favored over
the Gaussian distribution to model scRNA-seq data. Unfortunately, the negative
binomial distribution also raises challenges for data thinning. The
overdispersion parameter $\theta$ plays a role similar to the variance
parameter in the Gaussian distribution, as it is considered known. Therefore,
the quality of its estimation is directly impacts
$\operatorname{Cov}\left(X^{(1)},X^{(2)}\right)=\tau(1-\tau)\frac{\mu^{2}}{\theta}\left(1-\frac{\theta+1}{\widehat{\theta}+1}\right)$,
where $\widehat{\theta}$ is an estimated of $\theta$ (Neufeld, Dharamshi, Gao
& Witten, 2023). Also, for negative binomial mixtures, data thinning
decomposition is once again feasible only at the component level. As the
overdispersion parameter is component-specific (Li et al., 2018), providing an
estimator that ensures independence is here also a harduous and circular task,
given that the components themselves are again unknown and thus require
estimation through data thinning.

We further illustrate the necessity of applying intra-component data thinning
(with the associated intra-components overdispersion) to ensure Type I error
control in the post-clustering inference setting with numerical simulations.
We generated $n=100$ observations from a two-component negative binomial
mixture:
$0.5\text{NegBin}(\mu_{1},\theta_{1})+0.5\text{NegBin}(\mu_{2},\theta_{2})$
with component parameters $(\mu_{1},\theta_{1})=(5,5)$ and
$(\mu_{2},\theta_{2})=(60,40)$. We conducted similar post-clustering inference
as in Figure 2. In Figure 5A, applying the $k-$means with $K=3$ clusters
reveals that the first mixture component is erroneously split into 2 clusters
(labeled $C_{1}$ and $C_{3}$). We then assessed the Type I error rate
associated with the Wilcoxon test between these incorrect clusters when
applying data thinning with various overdispersion estimators. First, we
applied intra-component data thinning using oracle estimates
$\widetilde{\theta}_{k}$, $k=1,2$, representing true intra-component
overdispersion parameters (unfeasible in real-life application where the true
cluster structure is unknown). We compared those results with two alternatives
that are feasible in practice: i) applying intra-cluster data thinning based
on the $k-$means results from Figure 5A with their associated
$\widehat{\theta}_{\hat{k}}$, $\hat{k}=1,2,3$, and global data thinning with
its associated $\widehat{\theta}$. Of note, all the overdispersion estimations
were performed using Maximum Likelihood. Figure 5B presents the QQ-plot
against the Uniform distribution of the associated Wilcoxon p-values across
$1,\\!000$ simulations. Similarly to the Gaussian setting, achieving uniformly
distributed p-values is only possible through intra-component data thinning
using their oracle overdispersion estimates $\widetilde{\theta}_{k}$. All
other data thinning approaches are performed with biased estimators of the
overdispersion compromising the independence between $X^{(1)}$ and $X^{(2)}$,
and therefore leading to failure in the control of the post-clustering Type I
error rate.

Figure 5: Challenges in estimating overdispersion, a gene-specific parameter,
for negative binomial data thinning. Panel A showcases erroneous clustering
results on a simulated dataset, while Panel B presents a QQ-plot against the
Uniform distribution, displaying Wilcoxon p-values from various data thinning
approaches over $1,\\!000$ simulations when testing between the two erroneous
clusters. In Panel C, the estimated overdispersions of genes within different
cell populations (macrophages, monocytes, granulocytes, and CD4) are plotted
against those in neutrophils, along with their associated Root Mean Square
Error (RMSE) for comparison.

Extending our investigation from simulated scenarios to real-life
applications, we used a single-cell RNA-seq dataset from the Tabula Sapiens
Consortium (Consortium* et al., 2022) to delve into the practical challenges
of estimating overdispersion. Our analysis focused on five distinct cell
populations: $2,\\!560$ neutrophils, $105$ macrophages, $386$ monocytes, $454$
granulocytes, and $833$ CD4 T cells – all collected from a single donor. In
this controlled setting, where cell types were known, we succeeded in
estimating the overdispersion of $8,\\!333$ genes for each cell type using the
variance stabilizing transformation implemented in the sctransform package
(Choudhary & Satija, 2022). Figure 5C illustrates the comparison of gene
overdispersion when estimated solely in neutrophils versus the overdispersion
estimated for the same genes in the other four cellular populations. Root Mean
Squared Error (RMSE) values were computed to quantify the agreement between
estimations. Our findings reveal that overdispersion is specific to each cell
population, as evidenced by a notable deviation from the diagonal and
relatively high RMSE values. This underscores the challenge of accurately
estimating this parameter without prior knowledge of the true mixture
underlying the data. Combined with our simulation studies, these results
demonstrate how Type I error can easily be inflated in real-life applications
of data thining for scRNA-seq data analysis, particularly due to the
difficulty in providing an unbiased estimate of gene overdispersion.

##  4 Discussion

We highlight here the practical limitations inherent in data fission and its
extension, data thinning, for post-clustering inference challenges. A crucial
issue is the assumption of a homogeneous data distribution, which implies an
absence of true clusters in the data. To address this limitation and adapt to
scenarios with true classes, a shift towards mixture models becomes
imperative. However, these models lack a predefined decomposition through data
fission or thinning.

We have proposed an intra-component decomposition for data fission and data
thining, and demonstrated its theoretical validity. It relies on a priori
knowledge of the mixture components scale parameters, such as variances in the
Gaussian distribution or the overdispersion in the negative binomial
distribution. However in real-life applications, these parameters are unknown.
Adequately estimating these parameters becomes intricate in the presence of
true clusters, given their component-specific nature ; meanwhile the quality
of the estimation of those parameters is directly linked with the covariance
between the new random variables, $X^{(1)}$ and $X^{(2)}$, decomposing the
original data. Only unbiased estimation of these parameters ensures the
independence between $X^{(1)}$ and $X^{(2)}$. That independence is paramount
for post-clustering inference to adequately control the Type I error rate.

In the Gaussian framework, we theoretically quantified the relationship
between the relative bias in variance estimation and the associated Type I
error in post-clustering $t$-tests. In practice, our simulation results
suggest that a small relative bias can be acceptable while still achieving
effective Type I error control. These first results pave the way for defining
a principled approach to tuning the hyperparmater $tau$ in data fission and
data thining for post-clustering inference to optimize statistical power.

As a solution to avoid the need for prior knowledge of component-specific
variance parameters in the Gaussian mixture setting, we propose a
heteroscedastic model with individual variances, that can be estimated through
a non-parametric weighted estimator. This approach aligns more closely with
the distributional assumptions made by data fission and thinning. However, the
performance of this non-parametric approach relies heavily on the choice of
its bandwidth. The best bandwidth would be the one capturing only the
observations originating from a same components, but it would again require
knowledge of the true mixture components whereas their estimation is part of
the clustering first step of the method.

Finally, we demonstrated that the results derived the Gaussian distribution
context are readily applicable to the negative binomial distribution commonly
employed in modeling RNA-seq data. Specifically, we illustrate on real data
that overdispersion is also component-specific. Therefore, without knowing the
true cluster structure data, data thinning can not ensure the needed
independence between clustering and differential testing.

In practice, the application of data fission or data thinning for post-
clustering inference appears to be akin to a self-referential loop generating
circular reasoning. While introduced as a solution for addressing challenges
in post-clustering inference, all strategies that could theoretically ensure
independence between the two stages of the analysis ultimately rely on
knowledge of the true, but unknown, clusters. Despite its conceptual appeal,
the practical utility of these methods for post-clustering inference is
remains limited to extreme cases with extremely high signal vs. noise ratios,
emphasizing the need for alternative methodologies that can navigate the
complexities of unknown class structures more effectively.

All codes and data needed to reproduce the results presented here are openly
accessible from Zenodo with DOI 10.5281/zenodo.11207777 (Hivert, Agniel,
Thiébaut & Hejblum, 2024).

##  5 Acknowledgements

BH is supported partly by the Digital Public Health Graduate’s school, funded
by the PIA $3$ (Investments for the Future – Project reference: 17-EURE-0019).
The work was supported through the DESTRIER Inria Associate-Team from the
Inria@SiliconValley program (analytical code: DRI-012215), and by the CARE
project funded from the Innovative Medicines Initiative 2 Joint Undertaking
(JU) under grant agreement no IMI2-101005077. The JU receives support from the
European Union’s Horizon $2020$ research and innovation programme and EFPIA
and Bill $\&$ Melinda Gates Foundation, Global Health Drug Discovery
Institute, University of Dundee. This work has received funding from the
European Union’s Horizon $2020$ research and innovation programme under EHVA
grant agreement no H2020-681032. This study was carried out in the framework
of the University of Bordeaux’s France 2030 program / RRI PHDS. Computer time
for this study was provided by the computing facilities MCIA (Mésocentre de
Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université
de Pau et des Pays de l’Adour.

## References

  * (1)
  * Amezquita et al. (2020) Amezquita, R. A., Lun, A. T., Becht, E., Carey, V. J., Carpp, L. N., Geistlinger, L., Marini, F., Rue-Albrecht, K., Risso, D., Soneson, C. et al. (2020), ‘Orchestrating single-cell analysis with bioconductor’, Nature methods 17(2), 137–145.
  * Bachoc et al. (2023) Bachoc, F., Maugis-Rabusseau, C. & Neuvial, P. (2023), ‘Selective inference after convex clustering with $\ell_{1}$ penalization’, arXiv preprint arXiv:2309.01492 .
  * Bouveyron et al. (2019) Bouveyron, C., Celeux, G., Murphy, T. B. & Raftery, A. E. (2019), Model-based clustering and classification for data science: with applications in R, Vol. 50, Cambridge University Press.
  * Chacón & Duong (2020) Chacón, J. E. & Duong, T. (2020), Multivariate Kernel Smoothing and Its Applications SMOOTHING AND ITS APPLICATIONS, CRC PRESS.
  * Chen & Gao (2023) Chen, Y. T. & Gao, L. L. (2023), ‘Testing for a difference in means of a single feature after clustering’, arXiv preprint arXiv:2311.16375 .
  * Choudhary & Satija (2022) Choudhary, S. & Satija, R. (2022), ‘Comparison and evaluation of statistical error models for scrna-seq’, Genome Biology 23, 20.
https://doi.org/10.1186/s13059-021-02584-9

  * Consortium* et al. (2022) Consortium*, T. T. S., Jones, R. C., Karkanias, J., Krasnow, M. A., Pisco, A. O., Quake, S. R., Salzman, J., Yosef, N., Bulthaup, B., Brown, P. et al. (2022), ‘The tabula sapiens: A multiple-organ, single-cell transcriptomic atlas of humans’, Science 376(6594), eabl4896.
  * Fithian et al. (2014) Fithian, W., Sun, D. & Taylor, J. (2014), ‘Optimal inference after model selection’, arXiv preprint arXiv:1410.2597 .
  * Gao et al. (2022) Gao, L. L., Bien, J. & Witten, D. (2022), ‘Selective inference for hierarchical clustering’, Journal of the American Statistical Association pp. 1–11.
  * Heidenreich et al. (2013) Heidenreich, N.-B., Schindler, A. & Sperlich, S. (2013), ‘Bandwidth selection for kernel density estimation: a review of fully automatic selectors’, AStA Advances in Statistical Analysis 97, 403–433.
  * Hivert, Agniel, Thiébaut & Hejblum (2024) Hivert, B., Agniel, D., Thiébaut, R. & Hejblum, B. P. (2024), ‘Post-clustering difference testing: Valid inference and practical considerations with applications to ecological and biological data’, Computational Statistics & Data Analysis 193, 107916.
  * Hivert, Agniel, Thiébaut & Hejblum (2024) Hivert, B., Agniel, D., Thiébaut, R. & Hejblum, B. P. (2024), ‘Reproducible codes and results for Running in circles: is practical application feasible for data fission and data thinning in post-clustering differential analysis? (Version v1) [Data set]’, Zenodo .  DOI: 10.5281/zenodo.11207777.
https://doi.org/10.5281/zenodo.11207777

  * Jaeger & Banks (2023) Jaeger, A. & Banks, D. (2023), ‘Cluster analysis: A modern statistical review’, Wiley Interdisciplinary Reviews: Computational Statistics 15(3), e1597.
  * Kriegeskorte et al. (2009) Kriegeskorte, N., Simmons, W. K., Bellgowan, P. S. & Baker, C. I. (2009), ‘Circular analysis in systems neuroscience: the dangers of double dipping’, Nature neuroscience 12(5), 535–540.
  * Lähnemann et al. (2020) Lähnemann, D., Köster, J., Szczurek, E., McCarthy, D. J., Hicks, S. C., Robinson, M. D., Vallejos, C. A., Campbell, K. R., Beerenwinkel, N., Mahfouz, A. et al. (2020), ‘Eleven grand challenges in single-cell data science’, Genome biology 21(1), 1–35.
  * Leiner et al. (2023) Leiner, J., Duan, B., Wasserman, L. & Ramdas, A. (2023), ‘Data fission: splitting a single data point’, Journal of the American Statistical Association pp. 1–12.
  * Li et al. (2018) Li, Q., Noel-MacDonnell, J. R., Koestler, D. C., Goode, E. L. & Fridley, B. L. (2018), ‘Subject level clustering using a negative binomial model for small transcriptomic studies’, BMC bioinformatics 19(1), 1–10.
  * MacQueen et al. (1967) MacQueen, J. et al. (1967), ‘Some methods for classification and analysis of multivariate observations’, 1(14), 281–297.
  * Neufeld, Dharamshi, Gao & Witten (2023) Neufeld, A., Dharamshi, A., Gao, L. L. & Witten, D. (2023), ‘Data thinning for convolution-closed distributions’, arXiv preprint arXiv:2301.07276 .
  * Neufeld et al. (2024) Neufeld, A., Gao, L. L., Popp, J., Battle, A. & Witten, D. (2024), ‘Inference after latent variable estimation for single-cell rna sequencing data’, Biostatistics 25(1), 270–287.
  * Neufeld, Popp, Gao, Battle & Witten (2023) Neufeld, A., Popp, J., Gao, L. L., Battle, A. & Witten, D. (2023), ‘Negative binomial count splitting for single-cell rna sequencing data’, arXiv preprint arXiv:2307.12985 .
  * Pullin & McCarthy (2024) Pullin, J. M. & McCarthy, D. J. (2024), ‘A comparison of marker gene selection methods for single-cell rna sequencing data’, Genome Biology 25(1), 56.
  * Siffer et al. (2018) Siffer, A., Fouque, P.-A., Termier, A. & Largouët, C. (2018), Are your data gathered?, in ‘Proceedings of the 24th acm sigkdd international conference on knowledge discovery & data mining’, pp. 2210–2218.
  * Song et al. (2023) Song, D., Li, K., Ge, X. & Li, J. J. (2023), ‘Clusterde: a post-clustering differential expression (de) method robust to false-positive inflation caused by double dipping’, Research Square .
  * Welch (1947) Welch, B. L. (1947), ‘The generalization of ‘student’s’problem when several different population varlances are involved’, Biometrika 34(1-2), 28–35.
  * Yun & Foygel Barber (2023) Yun, Y.-J. & Foygel Barber, R. (2023), ‘Selective inference for clustering with unknown variance’, Electronic Journal of Statistics 17(2), 1923–1946.
  * Zhang et al. (2019) Zhang, J. M., Kamath, G. M. & David, N. T. (2019), ‘Valid post-clustering differential analysis for single-cell rna-seq’, Cell systems 9(4), 383–392.

##  Appendix A Independence proof of the Gaussian process

Let
$\boldsymbol{X}\sim\mathcal{N}\left(\boldsymbol{\mu},\boldsymbol{\Sigma}_{p\times
p}\right)$. Considering the fission process for Gaussian data described in
Table 1, we can decompose $\boldsymbol{X}$ into two new random variables
$\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$ using a new random variable
$\boldsymbol{Z}\sim\mathcal{N}\left(\boldsymbol{0},\boldsymbol{\Sigma}\right)$.
It follows form this decomposition that:

$\boldsymbol{X}^{(1)}\sim\mathcal{N}\left(\boldsymbol{\mu},\left(1+\tau^{2}\right)\boldsymbol{\Sigma}\right)\quad\text{and}\quad\boldsymbol{X}^{(2)}\sim\mathcal{N}\left(\boldsymbol{\mu},\left(1+\frac{1}{\tau^{2}}\right)\boldsymbol{\Sigma}\right)$

Moreover, we have
$\boldsymbol{X}^{(1)}\perp\\!\\!\\!\perp\boldsymbol{X}^{(2)}$. Indeed, we
have:

$\displaystyle\text{Cov}(\boldsymbol{X}^{(1)},\boldsymbol{X}^{(2)})$
$\displaystyle=\mathbb{E}\left[\left(\boldsymbol{X}^{(1)}-\mathbb{E}[\boldsymbol{X}^{(1)}]\right)\left(\boldsymbol{X}^{(2)}-\mathbb{E}[\boldsymbol{X}^{(2)}]\right)^{t}\right]$
$\displaystyle=\mathbb{E}\left[\left(\boldsymbol{X}^{(1)}-\boldsymbol{\mu}\right)\left(\boldsymbol{X}^{(2)}-\boldsymbol{\mu}\right)^{t}\right]\qquad\mbox{since}\qquad\mathbb{E}[\boldsymbol{X}^{(1)}]=\mathbb{E}[\boldsymbol{X}^{(2)}]=\boldsymbol{\mu}$
$\displaystyle=\mathbb{E}\left[\boldsymbol{X}^{(1)}{\boldsymbol{X}^{(2)}}^{t}-\boldsymbol{X}^{(1)}\boldsymbol{\mu}^{t}-\boldsymbol{\mu}{\boldsymbol{X}^{(2)}}^{t}+\boldsymbol{\mu}\boldsymbol{\mu}^{t}\right]$
$\displaystyle=\mathbb{E}[\boldsymbol{X}^{(1)}{\boldsymbol{X}^{(2)}}^{t}]-\mathbb{E}[{\boldsymbol{X}^{(1)}}]\boldsymbol{\mu}^{t}-\boldsymbol{\mu}\mathbb{E}[{\boldsymbol{X}^{(2)}}^{t}]+\boldsymbol{\mu}\boldsymbol{\mu}^{t}$
$\displaystyle=\mathbb{E}[\boldsymbol{X}^{(1)}{\boldsymbol{X}^{(2)}}^{t}]-\boldsymbol{\mu}\boldsymbol{\mu}^{t}-\boldsymbol{\mu}\boldsymbol{\mu}^{t}+\boldsymbol{\mu}\boldsymbol{\mu}^{t}$
$\displaystyle=\mathbb{E}[\boldsymbol{X}^{(1)}{\boldsymbol{X}^{(2)}}^{t}]-\boldsymbol{\mu}\boldsymbol{\mu}^{t}$

Moreover, we have:

$\displaystyle\mathbb{E}[\boldsymbol{X}^{(1)}{\boldsymbol{X}^{(2)}}^{t}]$
$\displaystyle=\mathbb{E}\left[\left(\boldsymbol{X}+\tau\boldsymbol{Z}\right)\left(\boldsymbol{X}-\frac{1}{\tau}\boldsymbol{Z}\right)^{t}\right]$
$\displaystyle=\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}^{t}-\frac{1}{\tau}\boldsymbol{X}\boldsymbol{Z}^{t}+\tau\boldsymbol{Z}\boldsymbol{X}^{t}-\boldsymbol{Z}\boldsymbol{Z}^{t}\right]$
$\displaystyle=\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^{t}]-\frac{1}{\tau}\mathbb{E}[\boldsymbol{X}\boldsymbol{Z}^{t}]+\tau\mathbb{E}[\boldsymbol{Z}\boldsymbol{X}^{t}]-\mathbb{E}[\boldsymbol{Z}\boldsymbol{Z}^{t}]$
$\displaystyle=\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^{t}]-\frac{1}{\tau}\mathbb{E}[\boldsymbol{X}]\mathbb{E}[\boldsymbol{Z}^{t}]+\tau\mathbb{E}[\boldsymbol{Z}]\mathbb{E}[\boldsymbol{X}^{t}]-\mathbb{E}[\boldsymbol{Z}\boldsymbol{Z}^{t}]\qquad\mbox{since}\qquad\boldsymbol{X}\perp\\!\\!\\!\perp\boldsymbol{Z}$
$\displaystyle=\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^{t}]-\mathbb{E}[\boldsymbol{Z}\boldsymbol{Z}^{t}]\qquad\mbox{since}\qquad\mathbb{E}[\boldsymbol{Z}]=0$

But we also have:

$\displaystyle\mathbb{V}\text{ar}(\boldsymbol{X})=\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^{t}]-\mathbb{E}[\boldsymbol{X}]\mathbb{E}[\boldsymbol{X}^{t}]=\boldsymbol{\Sigma}$
$\displaystyle\iff\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^{t}]=\boldsymbol{\Sigma}+\mathbb{E}[\boldsymbol{X}]\mathbb{E}[\boldsymbol{X}^{t}]$
$\displaystyle\iff\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^{t}]=\boldsymbol{\Sigma}+\boldsymbol{\mu}\boldsymbol{\mu}^{t}\qquad\mbox{since}\qquad\mathbb{E}[\boldsymbol{X}]=\boldsymbol{\mu}$

and

$\displaystyle\mathbb{V}\text{ar}(\boldsymbol{Z})=\mathbb{E}[\boldsymbol{Z}\boldsymbol{Z}^{t}]-\mathbb{E}[\boldsymbol{Z}]\mathbb{E}[\boldsymbol{Z}^{t}]=\boldsymbol{\Sigma}$
$\displaystyle\iff\mathbb{E}[\boldsymbol{Z}\boldsymbol{Z}^{t}]=\boldsymbol{\Sigma}+\mathbb{E}[\boldsymbol{Z}]\mathbb{E}[\boldsymbol{Z}^{t}]$
$\displaystyle\iff\mathbb{E}[\boldsymbol{Z}\boldsymbol{Z}^{t}]=\boldsymbol{\Sigma}\qquad\mbox{since}\qquad\mathbb{E}[\boldsymbol{Z}]=0$

So finally,

$\text{Cov}(\boldsymbol{X}^{(1)},\boldsymbol{X}^{(2)})=\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^{t}]-\mathbb{E}[\boldsymbol{Z}\boldsymbol{Z}^{t}]-\boldsymbol{\mu}\boldsymbol{\mu}^{t}=\boldsymbol{\Sigma}+\boldsymbol{\mu}\boldsymbol{\mu}^{t}-\boldsymbol{\Sigma}-\boldsymbol{\mu}\boldsymbol{\mu}^{t}=0$
(4)

##  Appendix B Impact of covariance estimation under the independence between
$\boldsymbol{X}^{(1)}$ and $\boldsymbol{X}^{(2)}$

Now, let suppose that we use
$\boldsymbol{Z}\sim\mathcal{N}\left(\boldsymbol{0},\hat{\boldsymbol{\Sigma}}\right)$
to perform data fission. It follows from equation (4) that, since
$\operatorname{\mathbb{V}ar}\left(\boldsymbol{Z}\right)=\hat{\boldsymbol{\Sigma}}$,

$\text{Cov}(\boldsymbol{X}^{(1)},\boldsymbol{X}^{(2)})=\boldsymbol{\Sigma}-\hat{\boldsymbol{\Sigma}}$

##  Appendix C Derivation of $t$-test statistic under the null of no-cluster
in the univariate data fission post-clustering setting

The fission process, as outlined in Table 1, gives:

$X^{(1)}\sim\mathcal{N}(\mu,\sigma^{2}+\tau^{2}\widehat{\sigma^{2}})\quad\text{and}\quad
X^{(2)}\sim\mathcal{N}\left(\mu,\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)$

After applying k-means clustering to the $n$ realizations of $X^{(1)}$, two
clusters $\widehat{C}_{1}$ and $\widehat{C}_{2}$ are obtained. As $X^{(1)}$
follows a Gaussian distribution, the k-means clustering splits it into two
parts $X^{(1)}|C_{k}$, $k=1,2$, with the same size $n/2$ and the same
variance, forming two half-Gaussian distributions. It follows from the half-
Gaussian distribution that:

$\mathbb{E}[X^{(1)}|C_{k}]=\pm\sqrt{\frac{2\left(\sigma^{2}+\tau^{2}\widehat{\sigma^{2}}\right)}{\pi}}$
(5)

where the sign depends on whether $\hat{C}_{k}$ represents the truncation on
the right or left side of $\mu$. Using (5), the mean differences between
$\widehat{C}_{1}$ and $\widehat{C}_{2}$ is then given by
$\hat{\delta}_{X^{(1)}}=\left|\mathbb{E}\left[X^{(1)}|\widehat{C}_{2}\right]-\mathbb{E}\left[X^{(1)}|\widehat{C}_{1}\right]\right|=2\sqrt{\frac{2}{\pi}\left(\sigma^{2}+\tau^{2}\widehat{\sigma^{2}}\right)}$.
Due to the k-means clustering, the two estimated clusters share the same
intra-cluster variance. Following the same half-Gaussian arguments, this
shared intra-cluster variance in $X^{(1)}$ is given by:

$s_{X^{(1)}}^{2}=(\sigma^{2}+\tau^{2}\widehat{\sigma^{2}})\left(1-\frac{2}{\pi}\right)$

In the post-clustering data fission setting, hypothesis testing is applied on
the realizations of $X^{(2)}$ using $\widehat{C}_{1}$ and $\widehat{C}_{2}$.
Our focus is on the mean differences $\hat{\delta}_{X^{(2)}}$ between these
two clusters and their shared intra-cluster variances $s^{2}(X^{(2)})$ on
$X^{(2)}$. These values depend on the ability to reproduce the clustering
results on $X^{(2)}$ which is quantified by
$\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)$. Assuming the same clustering
for $X^{(2)}$ as for $X^{(1)}$ implies the same truncation on $X^{(2)}$. In
this case, because
$X^{(2)}\sim\mathcal{N}\left(\mu,\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)$,
we have:

$\hat{\delta}_{X^{(2)}}=2\sqrt{\frac{2}{\pi}\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)}$

However, considering that the cluster analysis is applied first to $X^{(1)}$
and then transferred to $X^{(2)}$, we need to account for the correlation
between $X^{(1)}$ and $X^{(2)}$, leading to:

$\displaystyle\hat{\delta}_{X^{(2)}}$
$\displaystyle=2\sqrt{\frac{2}{\pi}\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)}\times\left|\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)\right|$
$\displaystyle=2\sqrt{\frac{2}{\pi}\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)}\times\left|\frac{(\sigma^{2}-\widehat{\sigma^{2}})}{\sqrt{(\sigma^{2}+\tau^{2}\widehat{\sigma^{2}})\times(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}})}}\right|$

This mean difference $\hat{\delta}_{X^{(2)}}$ can be interpreted as a trade-
off between how easily $\widehat{C}_{1}$ and $\widehat{C}_{2}$ can be
reproduced on $X^{(2)}$ and the mean differences that would be observed if
they were perfectly reproduced on $X^{(2)}$. For the shared intra-cluster
variance, $s^{2}(X^{(2)})$ should approximate
$(\sigma^{2}+\frac{1}{\tau}^{2}\widehat{\sigma^{2}})(1-\frac{2}{\pi})$ when
$\left|\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)\right|$ is close to $1$
(the scenario where we accurately replicate $\widehat{C}_{1}$ and
$\widehat{C}_{2}$ on $X^{(2)}$). In this situation, as the two estimated
clusters are fully preserved, their variances precisely follow the half-
Gaussian variance on $X^{(2)}$. Conversely, when
$\left|\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)\right|$ is close to $0$
(indicating the inability to reproduce the clustering results on $X^{(2)}$ due
to its independence with $X^{(1)}$), $s^{2}(X^{(2)})$ should approximate
$\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}$, representing the entire
variance of $X^{(2)}$. This is a consequence of the independence between
$X^{(1)}$ and $X^{(2)}$ that leads to random clusters on $X^{(2)}$. Thus,
these two considerations give:

$s^{2}(X^{(2)})=\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)\left(1-\frac{2}{\pi}\right)\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}+\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)\left(1-\left(\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)\right)^{2}\right)$
(6)

In the post-clustering inference setting, we are interested in testing a mean
difference between the two estimated clusters $\widehat{C}_{1}$ and
$\widehat{C}_{2}$. So, the null hypothesis we are interested in is:

$\mathcal{H}_{0}:\mu_{C_{1}}=\mu_{C_{2}}\qquad\mbox{vs}\qquad\mathcal{H}_{1}:\mu_{C_{1}}\neq\mu_{C_{2}}$

Under the null,
$T_{\widehat{C}_{1},\widehat{C}_{2}}=\frac{\overline{X^{(2)}_{\widehat{C}_{1}}}-\overline{X^{(2)}_{\widehat{C}_{2}}}}{\sqrt{\frac{\hat{s}^{2}_{\widehat{C}_{1}}}{n_{1}}+\frac{\hat{s}^{2}_{\widehat{C}_{2}}}{n_{2}}}}$
where $\overline{X^{(2)}_{C_{k}}}=\frac{1}{n_{k}}\sum\limits_{i\in
C_{k}}x_{i}^{(2)}$.

We assumed that the $2-$means clustering results in clusters such that
$n_{1}=n_{2}=\frac{n}{2}$. Moreover, using the half Gaussian argument, we
supposed that $\widehat{C}_{1}$ and $\widehat{C}_{2}$ shares the same intra-
cluster variances given by $s^{2}(X^{(2)})$. So,

$\displaystyle T_{\widehat{C}_{1},\widehat{C}_{2}}$
$\displaystyle=\frac{\overline{X^{(2)}_{\widehat{C}_{1}}}-\overline{X^{(2)}_{\widehat{C}_{2}}}}{\sqrt{\frac{s^{2}(X^{(2)})}{\frac{n}{2}}+\frac{s^{2}(X^{(2)})}{\frac{n}{2}}}}$
$\displaystyle=\frac{\overline{X^{(2)}_{\widehat{C}_{1}}}-\overline{X^{(2)}_{\widehat{C}_{2}}}}{\sqrt{\frac{4s^{2}(X^{(2)})}{n}}}$

We want to compute the null distribution of
$T_{\widehat{C}_{1},\widehat{C}_{2}}$. Using the Central Limit theorem, we
know that:

$\overline{X^{(2)}}_{\widehat{C}_{1}}-\overline{X^{(2)}}_{\widehat{C}_{2}}\sim\mathcal{N}\left(\underbrace{\mathbb{E}[X^{(2)}|\widehat{C}_{1}]-\mathbb{E}[X^{(2)}|\widehat{C}_{2}]}_{\hat{\delta}_{X^{(2)}}},\operatorname{\mathbb{V}ar}\left(\overline{X^{(2)}_{\widehat{C}_{1}}}-\overline{X^{(2)}_{\widehat{C}_{2}}}\right)\right)$

But,

$\displaystyle\operatorname{\mathbb{V}ar}\left(\overline{X^{(2)}_{\widehat{C}_{1}}}-\overline{X^{(2)}_{\widehat{C}_{2}}}\right)$
$\displaystyle=\operatorname{\mathbb{V}ar}\left(\overline{X^{(2)}_{\widehat{C}_{1}}}\right)+\operatorname{\mathbb{V}ar}\left(\overline{X^{(2)}_{\widehat{C}_{2}}}\right)\qquad\mbox{by
independence}$
$\displaystyle=2\operatorname{\mathbb{V}ar}\left(\overline{X^{(2)}_{\widehat{C}_{1}}}\right)\qquad\mbox{because
we assumed same intra-cluster variance}$
$\displaystyle=\frac{2}{\left(\frac{n}{2}\right)^{2}}\sum\limits_{i\in\widehat{C}_{1}}\operatorname{\mathbb{V}ar}\left(X^{(2)}_{i}|C_{1}\right)$
$\displaystyle=\frac{8}{n^{2}}\frac{n}{2}\operatorname{\mathbb{V}ar}\left(X^{(2)}_{1}|C_{1}\right)$
$\displaystyle=\frac{4}{n}\operatorname{\mathbb{V}ar}\left(X^{(2)}_{1}|C_{1}\right)$
$\displaystyle=\frac{4}{n}s^{2}(X^{(2)})$

So,

$\overline{X^{(2)}}_{\widehat{C}_{1}}-\overline{X^{(2)}}_{\widehat{C}_{2}}\sim\mathcal{N}\left(\widehat{\delta}_{X^{(2)}},\frac{4}{n}s^{2}(X^{(2)})\right)$

Finally, we have:

$\displaystyle\mathbb{E}\left[T_{\widehat{C}_{1},\widehat{C}_{2}}\right]$
$\displaystyle=\mathbb{E}\left[\frac{\overline{X^{(2)}_{\widehat{C}_{1}}}-\overline{X^{(2)}_{\widehat{C}_{2}}}}{\sqrt{\frac{4s^{2}(X^{(2)})}{n}}}\right]$
$\displaystyle=\frac{1}{\sqrt{\frac{4s^{2}(X^{(2)})}{n}}}\mathbb{E}\left[\overline{X^{(2)}_{\widehat{C}_{1}}}-\overline{X^{(2)}_{\widehat{C}_{2}}}\right]$
$\displaystyle=\frac{\widehat{\delta}_{X^{(2)}}}{\sqrt{\frac{4s^{2}(X^{(2)})}{n}}}$

and:

$\displaystyle\operatorname{\mathbb{V}ar}\left(T_{\widehat{C}_{1},\widehat{C}_{2}}\right)$
$\displaystyle=\operatorname{\mathbb{V}ar}\left(\frac{\overline{X^{(2)}_{\widehat{C}_{1}}}-\overline{X^{(2)}_{\widehat{C}_{2}}}}{\sqrt{\frac{4s^{2}(X^{(2)})}{n}}}\right)$
$\displaystyle=\frac{1}{\frac{4s^{2}(X^{(2)})}{n}}\operatorname{\mathbb{V}ar}\left(\overline{X^{(2)}_{\widehat{C}_{1}}}-\overline{X^{(2)}_{\widehat{C}_{2}}}\right)$
$\displaystyle=\frac{1}{\frac{4s^{2}(X^{(2)})}{n}}\frac{4s^{2}(X^{(2)})}{n}$
$\displaystyle=1$

So,

$T_{\widehat{C}_{1},\widehat{C}_{2}}\mathop{\sim}\limits^{\mathcal{H}_{0}}\mathcal{N}\left(\frac{\widehat{\delta}_{X^{(2)}}}{\sqrt{\frac{4s^{2}(X^{(2)})}{n}}},1\right)$

Since we have,
$\widehat{\delta}_{X^{(2)}}=2\sqrt{\frac{2}{\pi}\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)}\times\left|\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)\right|$
and
$s^{2}(X^{(2)})=(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}})(1-\frac{2}{\pi})\left(\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)\right)^{2}+(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}})\left(1-\left(\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)\right)^{2}\right)$,
which can be simplified as
$\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)\left(1-\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}\right)$,
then:

$\displaystyle\frac{\widehat{\delta}_{X^{(2)}}}{\sqrt{\frac{4s^{2}(X^{(2)})}{n}}}$
$\displaystyle=\frac{2\sqrt{\frac{2}{\pi}\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)}\times\left|\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)\right|}{\sqrt{\frac{4\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)\left(1-\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}\right)}{n}}}$
$\displaystyle=\frac{2\sqrt{\frac{2}{\pi}}\sqrt{\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)}\sqrt{\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}}{\frac{2\sqrt{\left(\sigma^{2}+\frac{1}{\tau^{2}}\widehat{\sigma^{2}}\right)}\sqrt{1-\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}}{\sqrt{n}}}$
$\displaystyle=\frac{\sqrt{n}\sqrt{\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}}{\sqrt{1-\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}}$

So, under $\mathcal{H}_{0}$,

$T_{\widehat{C}_{1},\widehat{C}_{2}}\mathop{\sim}\limits^{\mathcal{H}_{0}}\mathcal{N}\left(\frac{\sqrt{n}\sqrt{\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}}{\sqrt{1-\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}},1\right)$

The associated Type I error rate for this test is given by
$1-F(q_{\alpha/2})+F(-q_{\alpha/2})$, where $F$ is the cumulative distribution
function of
$\mathcal{N}\left(\frac{\sqrt{n}\sqrt{\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}}{\sqrt{1-\frac{2}{\pi}\operatorname{Cor}\left(X^{(1)},X^{(2)}\right)^{2}}},1\right)$
(or the corresponding non-central $\mathcal{T}(n-2)$ distribution), and
$q_{\alpha/2}$ is the quantile of a standard Gaussian distribution
$\mathcal{N}(0,1)$ (or of a $\mathcal{T}(n-2)$ distribution).