# Estimating the linear relation between variables that are never jointly
observed: an application in in vivo experiments

Polina Arsenteva1 Mohamed Amine Benadjaoud2 and Hervé Cardot1
1) Institut de Mathématiques de Bourgogne  UMR CNRS 5584
Université de Bourgogne  Dijon  France
2) IRSN PSE-SANTE/SERAMED  Fontenay aux roses  France

###### Abstract

This work is motivated by in vivo experiments in which measurement are
destructive so that the variables of interest can never be observed
simultaneously when the aim is to estimate the regression coefficients of a
linear regression. Assuming that the global experiment can be decomposed into
sub experiments (corresponding for example to different doses) with distinct
first moments, we propose different estimators of the linear regression which
take account of that additional information. We consider estimators based on
moments as well as estimators based optimal transport theory. These estimators
are proved to be consistent as well as asymptotically Gaussian under weak
hypotheses. The asymptotic variance has no explicit expression, except in some
particular cases, and specific bootstrap approaches are developed to build
confidence intervals for the estimated parameter. A Monte Carlo study is
conducted to assess and compare the finite sample performances of the
different approaches.

##  1 Motivation and context

In order to study the effects of a certain treatment on a living organism, in
vivo experiments are often conducted. In the context of a complex organism
response, scientists may be interested in studying multiple variables
describing the effect from different perspectives. In particular, such
variables of interest often include a macroscopic biomarker only available
through in vivo data on the one hand, and a microscopic biomarker that can
also be observed on a cellular level. The interest in this case lies in
predicting the former with the latter. For instance, in the context of
studying the adverse effects induced by radiotherapy on healthy tissues, the
potential outcomes of interest may manifest in the form of lesions, that are
quantified in a form a certain macroscopic biomarker, and the levels of a some
predictor such as gene expression. These measures often require sacrificing
the animal. As a result, the quantities of interest cannot be observed on the
same animals. Since the goal is to establish relationships between these
variables, a problem of statistical data fusion arises.

In this work we propose an approach to estimate the relationship between the
variables that are not simultaneously observed under the experimental setting
described above, based on a conditional model assuming linear relationship
within every component (experimental condition). An estimator derived with the
method of moments as well as optimal transport solution using Wasserstein
distance are considered for the real data problem in question.

###  1.1 Example of an in vivo experiment

Figure 1: Schematic representation of the design of an in vivo experiment
studying the effect of irradiated volume.

An illustrative example of an in vivo experiment where the variables of
interest are never observed simultaneously is presented in Figure 1. In this
experiment, mice are irradiated on the lungs with different volume, with a
goal of studying the role of irradiated volume in the appearing of radio-
induced adverse effects. The latter are assessed by measuring septal
thickening, a macroscopic biomarker of radio-induced adverse effects in the
lungs. The other variable that is measured with the purpose of predicting the
adverse effect related variable is the expression of multiple pro-inflammatory
genes. As shown in Figure 1, there are two independent cohorts in the study,
one is used to measure gene expression, whereas the other for measuring septal
thickness. This is a results of both measures being of destructive nature,
which does not allow them to be taken on the same animals.

Comparing distributions of the measurements arising from the two cohorts, such
as those presented in Figure 2, one may suspect a correlation or even a linear
relationship between the variables. In order to assert whether such a
relationship exists, one has to connect two variables that are never observed
simultaneously, which translates into solving a data fusion problem. In this
example, there are four groups indicating time points (1, 3, 6 and 12 months
after irradiation) when the corresponding animals are sacrificed and the
measurements are taken. Thus, the categorical variable indication the time
point, which is known for every observation, can be used as an additional
variable in order to link the predictor and the predicted variables between
each other.

Figure 2: Distribution of the data, collected from the irradiated patch under
SBRT with 3 mm beam size: the expression of the gene IL6 on the left, and
septal thickness on the right. The measurements were made 1, 3, 6 and 12
months after irradiation.

###  1.2 Existing research

The task of linking variables that are never jointly observed cannot be
approached as a typical missing values problem, since most methods for
inference on incomplete data require a sufficient overlap, which is completely
absent in the case treated here. As a result, any approaches applying related
frameworks such as multiple imputation in the context of data fusion are
unsuitable for our application. For example, Carrig et al., (2015) use
multiple imputation to integrate disparate datasets allowing for the absence
of the overlap, but requiring a calibration dataset, where all the variables
of interest must be jointly observed.

Other approaches to data fusion available in literature include factor
analysis (Cudeck,, 2000), statistical matching (Mitsuhiro and Hoshino,, 2020),
Bayesian network inference (Triantafillou et al.,, 2010; Tsamardinos et al.,,
2012) and Gaussian Markov combinations (Massa and Riccomagno,, 2017). These
methods are designed for linking variables that are not observed
simultaneously through covariates, present for both variables of interest.
This corresponds to the properties of the in vivo data described in Section
1.1. However, the covariates in these approaches are random variables,
typically continuous, and often assumed to be Gaussian, which is the case in
Cudeck, (2000) and Massa and Riccomagno, (2017). The groups variable available
through in vivo experiments may present in a continuous form, but the presence
of such categories as control and sham makes it impossible to assume
continuity and normality. The Bayesian network approaches introduced by
Triantafillou et al., (2010) and Tsamardinos et al., (2012), aimed at
inferring binary causal relationships between variables, are more suitable for
large datasets with a high number of covariates. Finally, currently available
research in statistical matching addresses such aspects as not-at-random
missingness (Mitsuhiro and Hoshino,, 2021) and high dimensionality (Mitsuhiro
and Hoshino,, 2020). This approach is based on the idea of comparing distances
between the covariates from the datasets of interest, which cannot be done by
taking the group variable as the covariate. It can be noted that the goal of
the aforementioned examples in statistical matching is to group individuals
before imputation, which is not necessary in our case since the groups are
already known.

##  2 Estimating the linear relation between variables that are never jointly
observed

###  2.1 Problem and presentation of the different identification approaches

We consider a real random variable $Y$ and a vector of $d$ real valued random
regressors $\mathbf{X}=(X_{1},\ldots,X_{d})$ and suppose that the following
linear regression hold

$\displaystyle Y$
$\displaystyle=\beta_{0}+\sum_{j=1}^{d}\beta_{j}X_{j}+\epsilon.$ (1)

The residuals $\epsilon$ are supposed to be independent of the random
covariates $X_{1},\ldots,X_{d}$, with zero mean and variance
$\sigma^{2}_{\epsilon}$. With destructive in vivo experiments, we can never
observe simultaneously $\mathbf{X}$ and $Y$, meaning that we never have at
hand the pair $(\mathbf{X},Y)$ but only $(\mathbf{X},.)$ and $(.,Y)$. This
means that only the marginal moments of $\mathbf{X}$ and $Y$ can be estimated
in presence of sampled data.

In the absence of additional information and without any strong additional
hypothesis, the parameters $(\beta_{0},\beta_{1},\ldots,\beta_{d})$ and the
variance of the noise $\sigma^{2}_{\epsilon}$ cannot be identified. Indeed, if
for example $X_{1}$ is centered with symmetric distribution, the coefficient
$\beta_{1}$ can only be determined up to sign change since $\beta_{1}X_{1}$
and $\beta_{1}(-X_{1})$ have the same distribution.

To deal with this identification issue, we consider that we can perform
different experiments in which the mean value of $X$ is allowed to vary. For
that, we suppose that there are $K$ groups (corresponding to $K$ different
experiments), defined by a discrete variable $G$ taking values in
$\\{1,\cdots,K\\}$ observed simultaneously with $Y$ and with $X$. This means
that we have now access to $(\mathbf{X},G)$ and $(Y,G)$ but not to
$(\mathbf{X},Y,G)$. We also suppose that $\epsilon$ is independent of $G$.

Given $G=k$, for $k=1,\ldots,K$, we denote by
$\mu_{Y}^{k}={\mathbb{E}}(Y|G=k)$ and
$\mu_{X_{j}}^{k}={\mathbb{E}}(X_{j}|G=k)$, $j=1,\ldots,d$, the expected values
within each group.

We now present two different approaches developed to identify the vector
$\boldsymbol{\beta}=(\beta_{0},\ldots,\beta_{d})$ of unknown regression
coefficients and the noise variance $\sigma_{\epsilon}^{2}$, taking account of
the additional information given by the discrete variable $G$.

####  2.1.1 Moment approach

A first and simple approach is based on first moments identification. Taking
the conditional expectation, given $G=k$, in (1), we have for $k=1,\ldots,K$,

$\displaystyle\mu_{Y}^{k}$
$\displaystyle=\beta_{0}+\sum_{j=1}^{d}\beta_{j}\mu_{X_{j}}^{k},$ (2)

since the residual term $\epsilon$ is supposed to satisfy
$\mathbb{E}(\epsilon|G=k)=0$ for $k=1,\ldots,K$.

We denote by $\boldsymbol{\mu}_{1,X}$ the $K\times(d+1)$ design matrix, whose
$k$th row is equal to $(1,\boldsymbol{\mu}_{X}^{k\top})$ with
$\boldsymbol{\mu}_{X}^{k}=(\mu_{X_{1}}^{k},\cdots,\mu_{X_{d}}^{k})^{\top}$, by
$\boldsymbol{\mu}_{Y}$ the $K$ dimensional vector with elements
$(\mu_{Y}^{1},\ldots,\mu_{Y}^{K})$. The $K$ linear equations in (2) can we
written in a matrix way,
$\boldsymbol{\mu}_{Y}=\boldsymbol{\mu}_{1,X}\boldsymbol{\beta}$.

We introduce the following assumption, guaranteeing the identifiability of the
model parameters:

$\displaystyle\mathbf{H}_{1}$
$\displaystyle\quad\mbox{rank}(\boldsymbol{\mu}_{1,X})=d+1,$

meaning that there are at least $K\geq d+1$ groups and that the $d+1$ column
vectors of $\boldsymbol{\mu}_{1,X}$ span a vector space of dimension $d+1$ in
$\mathbb{R}^{K}$.

######  Lemma 2.1.

If model (1) holds and assumption $\mathbf{H}_{1}$ is fulfilled,
$\boldsymbol{\beta}$ is uniquely identified in terms of the conditional first
order moments of $\mathbf{X}$ and $Y$ given $G$,

$\displaystyle\boldsymbol{\beta}$
$\displaystyle=\left(\boldsymbol{\mu}_{1,X}^{\top}\boldsymbol{\mu}_{1,X}\right)^{-1}\boldsymbol{\mu}_{1,X}^{\top}\boldsymbol{\mu}_{Y}.$

Additionally, the noise variance $\sigma_{\epsilon}^{2}$ satisfies

$\displaystyle\sigma_{\epsilon}^{2}$
$\displaystyle=\sigma_{Y}^{2}-\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}\boldsymbol{\beta}_{-0},$

where $\sigma_{Y}^{2}$ is the variance of $Y$, $\boldsymbol{\Gamma}_{X}$ is
the covariance matrix of $\mathbf{X}$, and
$\boldsymbol{\beta}_{-0}=(\beta_{1},\ldots,\beta_{d})$.

The proof of Lemma 2.1 is direct and thus omitted.

####  2.1.2 Optimal transport and minimum Wasserstein distance approach

The second approach is based on optimal transport, in particular on the idea
of estimating the linear transformation of the distribution of $\mathbf{X}$
that is the closest to that of $Y$ with respect to the Wasserstein distance
(see Panaretos and Zemel, (2019) for a general introduction for
statisticians).

The optimal transport map T between Gaussian measures on $\mathbb{R}^{d}$ is
linear, and the Wasserstein distance of order 2 between two Gaussian
distributions $D_{1}$ and $D_{2}$, with
$D_{1}\sim\mathcal{N}(\boldsymbol{\mu}_{1},\boldsymbol{\Gamma}_{1})$ and
$D_{2}\sim\mathcal{N}(\boldsymbol{\mu}_{2},\boldsymbol{\Gamma}_{2})$, is equal
to

$\displaystyle W_{2}^{2}(D_{1},D_{2})$
$\displaystyle=\|\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{2}\|^{2}+\mbox{tr}\left(\boldsymbol{\Gamma}_{1}+\boldsymbol{\Gamma}_{2}-2\left(\boldsymbol{\Gamma}_{2}^{1/2}\boldsymbol{\Gamma}_{1}\boldsymbol{\Gamma}_{2}^{1/2}\right)^{1/2}\right),$

where $\|.\|$ denotes the Euclidean norm and $\mbox{tr}(\mathbf{A})$ the trace
of matrix $\mathbf{A}$.

If we assume that, given $G=k$, $\mathbf{X}$ is a Gaussian random vector and
$\epsilon$ is Gaussian, we have, if the linear model (1) holds, that $Y$ is
also Gaussian given $G=k$, with expectation
$\mu_{Y}^{k}=\beta_{0}+\sum_{j=1}^{d}\beta_{j}\mu_{X_{j}}^{k}$ and variance
$\sigma_{\epsilon}^{2}+\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}$,
where $\boldsymbol{\Gamma}_{X}^{k}$ is the variance matrix of $\mathbf{X}$
given $G=k$. Thus, the Wasserstein distance between $D_{\gamma}$, the
distribution of
$\gamma_{0}+\boldsymbol{\gamma}_{-0}^{\top}\mathbf{X}+\epsilon$, and $D_{Y}$,
the distribution of $Y$, is equal to

$\displaystyle W_{2}^{2}(D_{\gamma},D_{Y})$
$\displaystyle=\mathbb{E}\left[W_{2}^{2}(D_{\gamma},D_{Y})|G\right]$
$\displaystyle=\sum_{k=1}^{K}\pi_{k}\left[(\mu_{Y}^{k}-\alpha_{0}-\boldsymbol{\gamma}_{-0}^{\top}\boldsymbol{\mu}_{X}^{k})^{2}+\left(\sigma_{Y,k}-\sqrt{\boldsymbol{\gamma}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\gamma}_{-0}+\sigma_{\epsilon}^{2}}\right)^{2}\right],$

where $\pi_{k}=\mathbb{P}[G=k]$ and $\sigma_{Y,k}^{2}=\mbox{Var}(Y|G=k)$.

We introduce the following loss criterion

$\displaystyle\varphi(\boldsymbol{\gamma},\sigma^{2})$
$\displaystyle=\sum_{k=1}^{K}\pi_{k}\left[(\mu_{Y}^{k}-\gamma_{0}-\boldsymbol{\gamma}_{-0}^{\top}\boldsymbol{\mu}_{X}^{k})^{2}+\left(\sigma_{Y,k}-\sqrt{\boldsymbol{\gamma}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\gamma}_{-0}+\sigma^{2}}\right)^{2}\right].$
(3)

which evaluates the weighted Wasserstein distance between $Y$ and a linear
combination of the $\mathbf{X}$ variables, contaminated by a Gaussian noise,
with variance $\sigma^{2}$. We state, without proof, the following Lemma
ensuring that the parameters $\boldsymbol{\beta}$ and $\sigma_{\epsilon}^{2}$
can be identified under general conditions.

######  Lemma 2.2.

If model (1) holds and if assumption $\mathbf{H}_{1}$ is fulfilled,
$\varphi(\boldsymbol{\gamma},\sigma^{2})$ has its unique minimum at
$\boldsymbol{\gamma}=\boldsymbol{\beta}$ and
$\sigma^{2}=\sigma_{\epsilon}^{2}$.

###  2.2 Sampled data and estimators

We suppose that experiments are made for $K\geq 2$ different groups and that
for each group $k$, for $k=1,\ldots,K$, we have two independent samples
$(Y^{k}_{1},\ldots,Y^{k}_{n_{y}^{k}})$ and
$(X^{k}_{j,1},\ldots,X^{k}_{j,n_{x}^{k}})_{j=1,\ldots,d}$, with sizes
$n_{y}^{k}$ and $n_{x}^{k}$. For each unit $i=1,\ldots,n_{x}^{k}$ from group
$k$, the vector of covariates is denoted by
$\mathbf{X}_{i}^{k}=(X_{1,i},\ldots,X_{d,i})$. We also define
$N_{x}=\sum_{k=1}^{K}n_{x}^{k}$ and $N_{y}=\sum_{k=1}^{K}n_{y}^{k}$, the total
number of observations of the response $Y$ and the covariates
$X_{1},\ldots,X_{d}$.

As shown in Lemma 2.1 and Lemma 2.2, identification of parameter
$\boldsymbol{\beta}$ relies on the knowledge of the first two conditional
moments, given $G$, of $Y$ and $\mathbf{X}$. For $k=1,\ldots,K$, we denote by
$\widehat{\mu}_{Y}^{k}=\frac{1}{n_{y}^{k}}\sum_{i=1}^{n_{y}^{k}}Y^{k}_{i}$ and
$\widehat{\mu}_{X_{j}}^{k}=\frac{1}{n_{x}^{k}}\sum_{i=1}^{n_{x}^{k}}X^{k}_{j,i}$
the empirical means within each group $k$ and by
$\widehat{\sigma}^{2}_{Y,k}=\frac{1}{n_{y}^{k}}\sum_{i=1}^{n_{y}^{k}}\left(Y_{i}^{k}\right)^{2}-\left(\widehat{\mu}_{Y}^{k}\right)^{2}$,
and
$\widehat{\boldsymbol{\Gamma}}^{k}_{X}=\frac{1}{n_{x}^{k}}\sum_{i=1}^{n_{x}^{k}}\mathbf{X}_{i}^{k}(\mathbf{X}_{i}^{k})^{\top}-\widehat{\boldsymbol{\mu}}_{X}^{k}(\widehat{\boldsymbol{\mu}}_{X}^{k})^{\top}$,
the empirical variances, with
$\widehat{\boldsymbol{\mu}}_{X}^{k}=(\widehat{\mu}_{X_{1}}^{k},\ldots,\widehat{\mu}_{X_{d}}^{k})$.
We also define the overall empirical mean and variance
$\widehat{\mu}_{Y}=\frac{1}{N_{y}}\sum_{k=1}^{K}{n_{y}^{k}}\widehat{\mu}_{Y}^{k}$,
$\widehat{\boldsymbol{\mu}}_{X}=\frac{1}{N_{x}}\sum_{k=1}^{K}{n_{x}^{k}}\widehat{\boldsymbol{\mu}}_{X}^{k}$,
$\widehat{\sigma}^{2}_{Y}=\frac{1}{N_{y}}\sum_{k=1}^{K}\sum_{i=1}^{n_{y}^{k}}\left(Y_{i}^{k}\right)^{2}-\left(\widehat{\mu}_{Y}\right)^{2}$,
$\widehat{\boldsymbol{\Gamma}}_{X}=\frac{1}{N_{x}}\sum_{k=1}^{K}\sum_{i=1}^{n_{x}^{k}}\mathbf{X}_{i}^{k}(\mathbf{X}_{i}^{k})^{\top}-\widehat{\boldsymbol{\mu}}_{X}\widehat{\boldsymbol{\mu}}_{X}^{\top}$.

We denote by $\widehat{\boldsymbol{\mu}}_{1,X}$ the $K\times(d+1)$ matrix,
with the first column consisting of ones, and the rest equal to
$\widehat{\boldsymbol{\mu}}_{X}$.

####  2.2.1 Moment estimators

If $\widehat{\boldsymbol{\mu}}_{1,X}$ is full rank, moment estimators of
$\boldsymbol{\beta}=(\beta_{0},\beta_{1},\ldots,\beta_{d})$ can be built by
considering the empirical counterpart of identification equations given in
Lemma 2.1,

$\displaystyle\widehat{\boldsymbol{\beta}}^{M}$
$\displaystyle=\left(\widehat{\boldsymbol{\mu}}_{1,X}^{\top}\widehat{\boldsymbol{\mu}}_{1,X}\right)^{-1}\widehat{\boldsymbol{\mu}}_{1,X}^{\top}\widehat{\boldsymbol{\mu}}_{Y}$

where
$\widehat{\boldsymbol{\mu}}_{Y}=(\widehat{\mu}_{Y}^{1},\ldots,\widehat{\mu}_{Y}^{K})$.
We consider in the following a slightly more general moment estimator of
$\boldsymbol{\beta}$, by introducing a weight $w_{k}$ given to each group $k$
of observations, with $w_{k}>0$ and $\sum_{k=1}^{K}w_{k}=1$. The weighted
moment estimator $\boldsymbol{\beta}$ is defined as the minimizer of

$\displaystyle\psi(\boldsymbol{\gamma})$
$\displaystyle=\sum_{k=1}^{K}w_{k}\left[\widehat{\mu}_{Y}^{k}-\left(\gamma_{0}+\sum_{j=1}^{d}\gamma_{j}\widehat{\mu}_{X_{j}}^{k}\right)\right]^{2},$

which is unique if $\widehat{\boldsymbol{\mu}}_{1,X}$ is full rank and defined
by

$\displaystyle\widehat{\boldsymbol{\beta}}^{M}$
$\displaystyle=\left(\widehat{\boldsymbol{\mu}}_{1,X}^{\top}\mathbf{w}\widehat{\boldsymbol{\mu}}_{1,X}\right)^{-1}\widehat{\boldsymbol{\mu}}_{1,X}^{\top}\mathbf{w}\widehat{\boldsymbol{\mu}}_{Y}$
(4)

with $\mathbf{w}=(w_{1},\ldots,w_{K})$. The first moment estimator considered
just before corresponds to equal weights $w_{k}=K^{-1}$.

We can then define the following estimator of the noise variance,

$\displaystyle\widehat{\sigma}_{\epsilon}^{2,M}$
$\displaystyle=\widehat{\sigma}_{Y}^{2}-(\widehat{\boldsymbol{\beta}}_{-0}^{M})^{\top}\widehat{\boldsymbol{\Gamma}}_{X}\widehat{\boldsymbol{\beta}}_{-0}^{M}.$
(5)

####  2.2.2 Minimum Wasserstein distance estimators

Estimators of $\boldsymbol{\beta}$ and $\sigma^{2}_{\epsilon}$ based on an
optimal transport criterion are derived by minimizing the empirical version
$\varphi_{n}(\boldsymbol{\gamma},\sigma^{2})$ of functional
$\varphi(\boldsymbol{\gamma},\sigma^{2})$ defined by

$\displaystyle\varphi_{n}(\boldsymbol{\gamma},\sigma^{2})$
$\displaystyle=\sum_{k=1}^{K}\pi_{k}\left[(\widehat{\mu}_{Y}^{k}-\gamma_{0}-\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\mu}}_{X}^{k})^{2}+\left(\widehat{\sigma}_{Y,k}-\sqrt{\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\gamma}_{-0}+\sigma^{2}}\right)^{2}\right].$
(6)

Note that in absence of an a priori information on the probability $\pi_{k}$
of observing group $k$, we can set $\pi_{k}=K^{-1}$. We denote by
$(\widehat{\boldsymbol{\beta}}^{W},\widehat{\sigma}^{2,W})$ minimizers of
$\varphi_{n}(\boldsymbol{\gamma},\sigma^{2})$ which are obtained with
iterative optimization algorithms based on gradient descent. The algorithm can
be initialized with
$(\widehat{\boldsymbol{\beta}}^{M},\widehat{\sigma}^{2,M})$ .

###  2.3 Consistency and asymptotic distribution

To study the asymptotic behavior of the estimators of $\boldsymbol{\beta}$
defined in previous Section, we suppose that the number $K$ of groups is kept
fixed and that for all groups and all variables $\mathbf{X}$ and $Y$, the
number of observations tends to infinity. This means that
$n_{\min}=\min(n_{y}^{1},\cdots,n_{y}^{K},n_{x}^{1},\cdots,n_{x}^{K})$ the
smallest sample size among all experiments is also supposed to tend to
infinity.

Under general moment conditions, without any particular assumption on the
probability law of $Y$ and $\mathbf{X}$ we get in the following two Lemmas
that the moment estimators and the Wasserstein distance based estimators are
consistent. The proofs are postponed to the Appendix.

######  Lemma 2.3.

If ${\mathbb{E}}(Y^{2})<+\infty$ and
${\mathbb{E}}(\|\mathbf{X}\|^{2})<+\infty$, and assumption $\mathbf{H}_{1}$ is
fulfilled, the sequence of estimators
$(\widehat{\boldsymbol{\beta}}^{M},\widehat{\sigma}_{\epsilon}^{2,M})$ defined
by (4) and (5) converges in probability to
$(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$ when $n_{\min}$ tends to
infinity.

For the Wasserstein minimum distance estimator, since there is no explicit
expression for the estimators, a compactness assumption is also made to obtain
the consistency.

######  Lemma 2.4.

If ${\mathbb{E}}(Y^{2})<+\infty$ and
${\mathbb{E}}(\|\mathbf{X}\|^{2})<+\infty$,
$(\boldsymbol{\beta},\sigma^{2}_{\epsilon})\in\Theta$ and $\Theta$ is a
compact set that does not contain 0, suppose that model (1) holds and
hypothesis $\mathbf{H}_{1}$ is fulfilled, then the sequence of estimators
$(\widehat{\boldsymbol{\beta}}^{W},\widehat{\sigma}^{2,W}_{\epsilon})$ that
minimize (6) converges in probability to $(\beta,\sigma^{2}_{\epsilon})$ when
$n_{\min}$ tends to infinity.

As far as the asymptotic distribution of the estimators is concerned, and for
sake of simplicity and lighter notations, we suppose from now on that the
number of experiments is the same for all groups and all variables, that is to
say $n=n_{y}^{1}=\cdots n_{y}^{K}=n_{x}^{1}\cdots=n_{x}^{K}$.

######  Proposition 2.1.

Assume that the assumptions of Lemma 2.3 are fulfilled. Then as $n$ tends to
infinity,

$\sqrt{n}\left(\widehat{\boldsymbol{\beta}}^{M}-\boldsymbol{\beta}\right)\rightsquigarrow\mathcal{N}\left(0,\boldsymbol{\Gamma}_{\beta_{M}}\right)$

where the expression of the asymptotic covariance matrix
$\boldsymbol{\Gamma}_{\beta_{M}}$ is given in the proof.

This result is based on the CLT for empirical means and the application of the
delta method, which involves computing the Jacobian of the inverse of
matrices, making obtaining the explicit expression for
$\boldsymbol{\Gamma}_{\beta_{M}}$ a difficult task when $d>1$.

######  Remark 2.1.

The weak convergence toward a Gaussian distribution presented in Proposition
2.1 would remain true, at the expense of heavier notations and a different
asymptotic covariance matrix $\boldsymbol{\Gamma}_{\beta_{M}}$, provided that
there exist two constants, $0<c\leq C$ such that

$0<c\leq\frac{\max(n_{y}^{1},\cdots,n_{y}^{K},n_{x}^{1},\cdots,n_{x}^{K})}{\min(n_{y}^{1},\cdots,n_{y}^{K},n_{x}^{1},\cdots,n_{x}^{K})}\leq
C<+\infty,$ (7)

and $\min(n_{y}^{1},\cdots,n_{y}^{K},n_{x}^{1},\cdots,n_{x}^{K})\to\infty$.

The asymptotic normality of $\widehat{\boldsymbol{\beta}}^{W}$ relies on
classical results for M-estimators recalled in the Appendix (see Theorem A.4).
Note that since we estimate simultaneously $\beta$ and
$\sigma_{\epsilon}^{2}$, an additional condition on the existence of the
moments of order four for the covariates is required.

######  Proposition 2.2.

If model (1) holds and hypothesis $\mathbf{H}_{1}$ is fulfilled,
${\mathbb{E}}(Y^{2})<+\infty$ and ${\mathbb{E}}(\|\mathbf{X}\|^{4})<+\infty$,
$(\boldsymbol{\beta},\sigma^{2}_{\epsilon})\in\Theta$ and $\Theta$ is a
compact set that does not contain $(0,0)$, then, as $n$ tends to infinity,

$\displaystyle\sqrt{n}\left(\begin{pmatrix}\widehat{\boldsymbol{\beta}}^{W}\\\
\widehat{\sigma}^{2,W}_{\epsilon}\end{pmatrix}-\begin{pmatrix}\boldsymbol{\beta}\\\
\sigma^{2}_{\epsilon}\end{pmatrix}\right)$
$\displaystyle\rightsquigarrow\mathcal{N}\left(0,\boldsymbol{\Gamma}_{W}\right),$

for some covariance matrix $\boldsymbol{\Gamma}_{W}$.

Similarly to $\boldsymbol{\Gamma}_{\beta_{M}}$, the expression of the
asymptotic covariance matrix $\boldsymbol{\Gamma}_{W}$ is almost impossible to
explicitly derive manually, with the exception of some particularly simple
cases.

##  3 Bootstrapping for confidence intervals

Since, as noted in the previous section, it is complicated to compute
explicitly the asymptotic variance matrix of
$\widehat{\boldsymbol{\beta}}^{M}$ and $\widehat{\boldsymbol{\beta}}^{W}$, we
consider stratified bootstrap approaches in order to build confidence sets for
$\boldsymbol{\beta}$.

Our bootstrap procedure takes account of the independence between the
different groups $k=1,\ldots,K$, as well as the independence of inputs
$(X_{1}^{k},\ldots,X_{d}^{k})$ and output $Y^{k}$ within each group, meaning
more formally that, given $G=k$, the joint probability measure
$\mathbb{P}^{k}$ of $Y$ and $\mathbf{X}$ is a product measure of the marginal
measures
$\mathbb{P}^{k}=\mathbb{P}_{Y}^{k}\otimes\mathbb{P}_{\mathbf{X}}^{k}$.

Within each group $k$, we draw, with equal probability and with replacement,
$n_{y}^{k}$ observations among $Y_{1}^{k},\ldots,Y_{n_{y}^{k}}^{k}$, We also
draw independently, with equal probability and with replacement, $n_{x}^{k}$
observations among $\mathbf{X}_{1}^{k},\ldots,\mathbf{X}_{n_{x}^{k}}^{k}$. We
denote by $\mu_{Y}^{k*}$ and by $\boldsymbol{\mu}_{X}^{k*}$ the empirical
means and by $\sigma_{Y,k}^{2,*}$ and by $\boldsymbol{\Gamma}_{X}^{k*}$ the
empirical variances of $Y$ and $\mathbf{X}$ in these these bootstrap samples.
Bootstrapped estimators $\boldsymbol{\beta}^{M,*}$ and
$\boldsymbol{\beta}^{W,*}$ of $\boldsymbol{\beta}$ can now be computed by
replacing the empirical moments by the bootstrap moments in (4) and (6).

To build confidence sets for the components of $\boldsymbol{\beta}$ based on
previous bootstrap procedure, the bootstrap percentile technique, described in
Chapter 4 of Shao and Tu, (1995), is simple to use.

It can be noted that our estimators are smooth functions of sample means so
that classical bootstrap theory applies (see for example Shao and Tu, (1995),
Chapter 3). For simplicity, as in Proposition 2.1, we suppose that
$n=n_{y}^{1}=\ldots=n_{y}^{K}=n_{x}^{1}=\ldots=n_{x}^{K}$. Because of the
considered experimental design, our global "empirical distribution" is made of
products of marginal empirical distributions, the bootstrap for means is
almost surely consistent for the Kolmogorov metric, and with Theorem 3.1 in
Shao and Tu, (1995) the same result holds for the estimators of
$\boldsymbol{\beta}$ considered in this work. Then, the application of Theorem
4.1 in Shao and Tu, (1995) allows to conclude that bootstrap percentile method
gives consistent confidence sets for each component of $\boldsymbol{\beta}$.

######  Proposition 3.1.

Suppose that $\mathbb{E}(Y^{2})<\infty$ and
$\mathbb{E}\|\mathbf{X}\|^{2}<\infty$ and hypothesis $\mathbf{H}_{1}$ is
fulfilled. Then as $n\to+\infty$, the bootstrap estimator
$\boldsymbol{\beta}^{M,*}$ is strongly consistent for $\boldsymbol{\beta}$ in
the Kolmogorov metric. Consider a risk $\alpha\in(0,1)$, for each component of
$\boldsymbol{\beta}$, the bootstrap percentile approach provides, for a given
nominal level $1-\alpha$, a consistent confidence set.

We can also state a similar result for the estimators minimizing the
Wasserstein distance, under slightly more restrictive moment conditions and
compactness assumption.

######  Proposition 3.2.

Suppose that the assumptions of Proposition 2.2 as fulfilled. As
$n\to+\infty$, the bootstrap estimator
$(\boldsymbol{\beta}^{W,*},\sigma_{\epsilon}^{2,W,*})$ is strongly consistent
for $(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$ in the Kolmogorov metric.
Consider a risk $\alpha\in(0,1)$, for each component of $\boldsymbol{\beta}$,
the bootstrap percentile approach provides, for a given nominal level
$1-\alpha$, a consistent confidence set.

Propositions 3.1 and 3.2 are very important for practical applications since
they ensure that even if we are not able to compute the explicit expression of
the asymptotic variance of our estimators, we can still construct asymptotic
confidence intervals thanks to simple bootstrap procedures.

##  4 A simulation study

Some simulations are performed in order to evaluate the finite sample
performances of the proposed approaches on data that resemble in vivo data
originating from real experiments on mice.

###  4.1 Simulation design

####  4.1.1 General setting

Let $K\in\mathbb{N}$ be the number of subpopulations,
$g_{k}=\\{0,1,\dots,k-1\\}$ be the level of the grouping variable for
subpopulation $k\in\\{1,\dots,K\\}$, and the number of animals observed per
group was choose for simplicity
$n=n_{y}^{1}=\cdots=n_{y}^{K}=n_{x}^{1}\cdots=n_{x}^{K}$. For every animal
$i\in\\{1,\dots,n\\}$ and every subpopulation $k\in\\{1,\dots,K\\}$ we choose
the predictor variable $X^{k}_{i}$ in the Gaussian univariate setting:

$X^{k}_{i}\sim\mathcal{N}(\mu_{X}^{k},\sigma_{X,k}^{2})\text{, where
}\mu_{X}^{k}=\Delta_{\mu}g_{k}+C_{\mu},$

where $C_{\mu}\in\mathbb{R}$ is the global mean scale parameter, and
$\Delta_{\mu}\in\mathbb{R}$ is the parameter that determines the difference
between groups in terms of the mean. We simulate the predicted variable
independently as follows, with regression parameters $\beta_{0}$ and
$\beta_{1}$:

$Y^{k}_{i}=\beta_{0}+\beta_{1}X^{\prime k}_{i}+\epsilon^{k}_{i}\text{, where
}X^{\prime k}_{i}\sim X^{k}_{i}\text{ and
}\epsilon^{k}_{i}\sim\mathcal{N}(0,\sigma_{\epsilon,k}^{2}).$

It should be noted that $X^{k}_{i}$ and $X^{\prime k}_{i}$ are independent, we
thus recreate the situation of the predictor and predicted variables not being
simultaneously observed.

The variable sets $X=(X^{k}_{i})_{1\leq i\leq n,1\leq k\leq K}$ and
$Y=(Y^{k}_{i})_{1\leq i\leq n,1\leq k\leq K}$ are simulated
$N_{sim}\in\mathbb{N}$ times. For each simulation, $N_{boot}\in\mathbb{N}$
bootstrap samples of the size $n$ are generated from $X^{k}=(X^{k}_{i})_{1\leq
k\leq K}$ and $Y^{k}=(Y^{k}_{i})_{1\leq k\leq K}$ for each subpopulation
$k\in\\{1,\dots,K\\}$ independently, then moment estimators $\mu_{X}^{k*}$ and
$\mu_{Y}^{k*}$ are calculated. Finally, the bootstrap sample-based estimators
$\boldsymbol{\beta}^{M,*}=(\beta_{0}^{M,*},\beta_{1}^{M,*})$ and
$\boldsymbol{\beta}^{W,*}=(\beta_{0}^{W,*},\beta_{1}^{W,*})$ are calculated.
Based on $N_{boot}$ estimates, we calculate 95% confidence intervals using the
function quantile of the Python library NumPy. As a result, we obtain
$N_{sim}$ confidence intervals for each regression parameter, that we use to
calculate the following quantities of interest: coverage rate of the
intervals, their average amplitude, and power, i.e. the proportion of
intervals not including 0. Along with the bootstrap procedure, we estimate the
confidence intervals for the parameter estimators of the regression on the
means by group with a naive method, assuming that the deviation of the
parameter estimator from the true value divided by the estimator’s standard
error follows a Student’s t-distribution:

$\frac{\widehat{\beta}_{j}-\beta_{j}}{SE(\widehat{\beta_{j}})}\rightsquigarrow
t_{K-(d+1)}\text{ for }j\in\\{0,\dots,d\\}.$

####  4.1.2 Parameter levels

We fix the following parameters throughout all simulations: the regression
parameters $\beta_{0}=1$ and $\beta_{1}=2$, the number of simulations
$N_{sim}=500$, the number of bootstrap samples $N_{boot}=500$, the location
parameters of the predictor variable $\Delta_{\mu}=1$ and $C_{\mu}=10$.
Multiple parameters are varied to study their effect. We take the number of
animals $n\in\\{10,30\\}$, in particular to test whether the inference is
significantly impaired in case of a small number of animals, which is often
the case for real experimental data. We consider the number of groups
$K\in\\{4,10\\}$, 4 being the number of groups that is often observed in real
data, and 10 being a higher number that may produce sufficiently good results
with the naive approach to approximating confidence intervals with Student
distribution. Additionally, the parameter $\sigma^{2}_{X,k}$ can be adjusted
to control the extent to which the observations per group can be easily
distinguish one from another. We set
$\sigma_{X,1}^{2}=\dots=\sigma_{X,K}^{2}$, and
$\sigma^{2}_{X,k}\in\\{0.75,2\\}$, the first value corresponding to lower
overlap between groups, and the second to higher overlap. Finally, we
introduce an additional parameter $\rho\in\mathbb{R}^{+}$ controlling the
variance of the response to the variance of noise ratio, i.e.
$\rho=\cfrac{\sigma_{Y,k}}{\sigma_{\epsilon,k}}$. The choice of adjusting the
noise to signal ratio instead of the quantity of noise itself through
$\sigma_{\epsilon,k}^{2}$ is motivated by the fact that $\sigma_{Y,k}$ depends
on $\sigma_{X,k}$, hence the same level of $\sigma_{\epsilon,k}$ cannot be
interpreted the same way for different values of $\sigma_{X,k}$. The variance
of the noise can be expressed as follows:
$\sigma_{\epsilon,k}^{2}=\cfrac{\beta_{1}^{2}\sigma_{X,k}^{2}}{\rho^{2}-1}$.
The values of $\rho$ are chosen to correspond to the realistic situation,
namely the very noisy case and a slightly less noisy one:
$\rho\in\\{1.01,1.1\\}$. The effect of different values of $\rho$ on the
simulated response variable is illustrated on Figure 3.

(a)

(b)

(c)

Figure 3: The effect of different values of $\rho$ on the data, with $K=4$ and
$\sigma_{X,1}^{2}=\dots=\sigma_{X,4}^{2}=0.75$. a) Boxplots constructed from
the simulated values of $X^{k}_{i}$. b) Boxplots constructed from the
simulated values of $Y^{k}_{i}$ with lower relative noise level, i.e.
$\rho=1.1$. c) Boxplots constructed from the simulated values of $Y^{k}_{i}$
with higher relative noise level, i.e. $\rho=1.01$.

###  4.2 Results

The results of the simulation study are presented in Figure 4, namely the
obtained coverage rates for the linear regression slope confidence intervals
on Figure 4(a), the amplitudes on Figure 4(b), and the associated test powers
on Figure 4(c). Every figure contains four tables, each corresponding to a
combination of parameters with respect to the number of animals $N_{anim}$ and
the number of groups $K$. The columns of the tables represent combinations of
parameters with respect to the group dispersion parameter $\sigma^{2}_{X,k}$
and the noise to signal ratio $\rho$. Finally, the lines of each table
indicate the method used to estimate the confidence intervals: the proposed
method of moments and optimal transport based bootstrap estimators, and the
classical linear regression procedure based on estimated means within each
group.

(a)

(b)

(c)

Figure 4: a) Coverage rates, b) amplitudes, and c) powers of the confidence
intervals for the estimators of $\beta_{1}$ obtained from $N_{sim}$
simulations for each parameter combination. The method labels on the left: "mm
(boot)" stands for the method of moments with bootstrap, "ot (boot)" for the
optimal transport method with bootstrap, and "mm (student)" for the naive
linear regression on means approach based on Student’s distribution.

In general, it can be observed that the bootstrap estimators produce
confidence intervals with smaller amplitudes and with higher power, but with
lower coverage rate. Whereas the coverage rates in all cases remain within
91-96%, the extent to which the amplitudes are smaller and the powers are
bigger is significantly more important for the bootstrap estimators in almost
all cases. This implies that the naive approach based on the Student’s
distribution is more likely to produce false negatives in terms of
significance.This trend of further amplified by the number of groups
parameter: whereas the results are overall worsened with the decrease in
either the number of animals or groups, it is the case with the small number
of groups that demonstrates the biggest difference in the approaches. Indeed,
in almost all cases within the tables with $K=4$ we observe the amplitudes
approximately twice as important for the naive approach, and a similar trend
in terms of lower powers. The latter result is important since lower power
implies bigger probability of not detecting a significant relationship between
the predictor and the predicted variables, when it is actually present.
Overall, these results mean that the proposed bootstrap estimators are more
effective if the experimental design entails a small number of groups.

Concerning the remaining two parameters, as expected, in general the best
results are obtained with lower $\sigma^{2}_{X,k}$ and higher $\rho$. In most
case, the results for the naive and the bootstrap estimators are either both
good or both bad in terms of power, with the latter being slightly better. A
particularly complicated case can be distinguished, with high dispersion and
high noise level, with few groups and few animals, where all estimators
drastically fail: we observe almost equally bad powers (0.1 for both bootstrap
estimators and 0.09 for the naive estimator), despite the significant
difference in amplitudes. On the other hand, we can also distinguish two case
where the powers of the bootstrap estimators are more than 90%, whereas those
of the naive estimator are under 50%: in both cases there are 4 groups and
high noise, in the first case there are only 10 animals but lower dispersion,
in the second case high dispersion level is compensated but a higher number of
animals. This implies that if the underlying distributions per group are
characterized by a reasonable amount of overlap, or a significant overlap is
compensated by having more observations, the bootstrap estimators manage to
detect the significant relationship in most cases, unlike the naive estimator.

Lastly, it can be observed that the estimator based on optimal transport
produces confidence intervals with slightly smaller amplitudes compared to the
method of moments estimator. The difference appears to be relatively more
important in the cases with higher group overlap $\sigma^{2}_{X,k}=2$.
However, the powers are not affected by this difference. This may be explained
by a more important bias associated with the optimal transport estimator. The
estimator will likely produce better results in terms of the power than the
method of moments estimator if the bias is corrected.

##  5 Application to real data

In order to illustrate the proposed estimators on real data, we studied the
data from experiments conducted on mice in order to assess the adverse effects
induced in the context of different irradiated volume. In the course of these
experiments, mice were exposed to either stereotactic body radiation therapy
(SBRT) with different beam sizes at 90 Gy on the left lung, or whole-thorax
irradiation (WTI) at 19 Gy. For one cohort, the expressions of three pro-
inflammatory genes (IL1 $\alpha$ and $\beta$, IL6 and TNF) were measured, for
the other cohort the measurements of the thickness of the alveolar septas were
made due to its role as a macroscopic biomarker of radio-induced pulmonary
lesions. In the case of the SBRT, the measurements were taken in multiple
locations: the irradiated patch (inside the irradiation field), the remaining
part of the left lung referred to as ipsilateral lung, and the right lung
(contralateral lung). The goal of this statistical analysis is to determine
whether there is a statistical association between the gene expression as
predictors and the septal thickening as as outcome. Our approach is applied
since the variables are measured on different animals, but within each
irradiation condition there are shared groups in terms of the measurement time
points.

Table 1: Results of estimation of the linear regression slope predicting septal thickening with the pro-inflammatory genes expression, with three methods, for WTI and SBRT with different beam sizes, with measurements taken in different parts of lungs. | Method of moments | Optimal transport | Lin. Reg. on means
---|---|---|---
| (bootstrap) | (bootstrap) | (Student)
Loc. | Vol. | Gene | ${\widehat{\beta}_{1}}$ | 95% C.I. | Signif. | ${\widehat{\beta}_{1}}$ | 95% C.I. | Signif. | ${\widehat{\beta}_{1}}$ | 95% C.I. | Signif.
|  | IL1b | -0.35 | (-1.03, 0.4) | ✗ | -0.19 | (-0.44, 0.23) | ✗ | -0.35 | (-1.99, 1.28) | ✗
| 1 mm | IL6 | 0.43 | (-0.23, 1.2) | ✗ | 0.2 | (-0.2, 0.84) | ✗ | 0.43 | (-1.32, 2.17) | ✗
Ipsilateral |  | TNF | 0.2 | (-0.23, 0.84) | ✗ | 0.16 | (-0.18, 0.84) | ✗ | 0.2 | (-1.26, 1.66) | ✗
lung |  | IL1b | 0.9 | (0.0, 2.11) | ✔ | 1.05 | (0.01, 1.44) | ✔ | 0.9 | (-0.43, 2.24) | ✗
| 3 mm | IL6 | 0.05 | (-0.34, 0.46) | ✗ | 0.06 | (-0.33, 0.49) | ✗ | 0.05 | (-1.65, 1.76) | ✗
|  | TNF | 0.65 | (-0.12, 1.6) | ✗ | 0.63 | (-0.12, 1.46) | ✗ | 0.65 | (-1.57, 2.87) | ✗
|  | IL1b | 2.31 | (-1.73, 3.87) | ✗ | 0.82 | (-0.48, 0.94) | ✗ | 2.31 | (-0.73, 5.35) | ✗
| 1 mm | IL6 | 0.85 | (-0.7, 2.38) | ✗ | 0.61 | (-0.56, 1.6) | ✗ | 0.85 | (-3.75, 5.45) | ✗
Irradiated |  | TNF | 0.85 | (-0.6, 2.3) | ✗ | 0.69 | (-0.54, 1.53) | ✗ | 0.85 | (-3.96, 5.66) | ✗
patch |  | IL1b | 2.85 | (0.84, 5.14) | ✔ | 2.59 | (0.73, 3.58) | ✔ | 2.85 | (-1.07, 6.78) | ✗
| 3 mm | IL6 | 1.35 | (0.22, 2.47) | ✔ | 1.3 | (0.21, 2.26) | ✔ | 1.35 | (-1.8, 4.5) | ✗
|  | TNF | 3.81 | (1.01, 6.37) | ✔ | 3.37 | (0.99, 5.33) | ✔ | 3.81 | (-1.86, 9.48) | ✗
|  | IL1b | -0.92 | (-2.2, 1.06) | ✗ | -0.42 | (-0.89, 0.51) | ✗ | -0.92 | (-3.97, 2.14) | ✗
| 0 mm | IL6 | 2.23 | (-1.99, 2.28) | ✗ | 0.83 | (-0.83, 0.97) | ✗ | 2.23 | (-2.65, 7.12) | ✗
|  | TNF | 2 | (-2.0, 2.82) | ✗ | 0.88 | (-1.0, 1.11) | ✗ | 2 | (-1.72, 5.72) | ✗
|  | IL1b | 0.63 | (-0.82, 2.32) | ✗ | 0.23 | (-0.33, 0.56) | ✗ | 0.63 | (-2.95, 4.21) | ✗
| 1 mm | IL6 | 1.03 | (-0.57, 2.19) | ✗ | 0.46 | (-0.3, 1.0) | ✗ | 1.03 | (-0.88, 2.94) | ✗
Right |  | TNF | 1.05 | (-0.47, 2.43) | ✗ | 0.66 | (-0.31, 1.26) | ✗ | 1.05 | (-1.01, 3.11) | ✗
lung |  | IL1b | 1.07 | (0.17, 1.97) | ✔ | 0.92 | (0.17, 1.13) | ✔ | 1.07 | (-1.96, 4.1) | ✗
| 3 mm | IL6 | 0.3 | (-1.45, 1.12) | ✗ | 0.37 | (-0.74, 0.86) | ✗ | 0.3 | (-4.94, 5.53) | ✗
|  | TNF | 2.02 | (0.27, 4.1) | ✔ | 1.05 | (0.04, 1.44) | ✔ | 2.02 | (0.03, 4.02) | ✔
|  | IL1b | 6.41 | (-11.6, 15.4) | ✗ | 1.62 | (-2.28, 3.1) | ✗ | 6.41 | (-16.7, 29.52) | ✗
| 7 mm | IL6 | -0.58 | (-1.26, 0.09) | ✗ | -0.6 | (-1.23, 0.1) | ✗ | -0.58 | (-6.23, 5.08) | ✗
|  | TNF | 0.51 | (-1.7, 3.86) | ✗ | 0.43 | (-1.45, 1.89) | ✗ | 0.51 | (-23.21, 24.22) | ✗
Whole | IL1a | 2.93 | (-1.18, 6.63) | ✗ | 2.47 | (-1.49, 3.02) | ✗ | 2.93 | (-11.08, 16.94) | ✗
thorax | IL6 | 3.7 | (1.53, 5.99) | ✔ | 2.51 | (1.39, 3.43) | ✔ | 3.7 | (1.11, 6.29) | ✔
irradiation | TNF | 2.35 | (0.57, 4.73) | ✔ | 1.67 | (0.53, 1.88) | ✔ | 2.35 | (-3.86, 8.57) | ✗

The linear regression parameters were estimated with three estimators the same
way it was done in the simulation study in Section 4.2. The results are
presented in Table LABEL:tab:vivo_results. The focus is placed on estimating
the slope parameter $\beta_{1}$ in particular. The table contains the
estimations of $\beta_{1}$ as well as the estimated confidence intervals for
the slope estimator, and the corresponding test result on the significance of
the estimated relationship.

Figure 5: Distribution of the data, collected from the irradiated patch under
SBRT with 7 mm beam size: the expression of the gene TNF on the left, and
septal thickness on the right. The measurements were made 30, 45 and 60 days
after irradiation.

Major differences between the results for the naive and the bootstrap
estimators can be observed in terms of the detected significance. On the one
hand, the relationship between the pro-inflammatory genes and septal
thickening has been identified by all methods in the case of whole-thorax
irradiation (two out of three genes for the bootstrap estimators, and only one
for the naive estimator), which is an expected result. On the other hand, we
also expect to identify a strong correlation in the case of the measurements
taken directly from the irradiated patch. This is only the case for the
bootstrap estimators, but not for the naive one. This results is in accordance
with the results we obtain with simulated data: the confidence intervals are
often overestimated with the naive approach, which may result in false
negatives in terms of significance.

For the irradiated patch, as well as in all other cases with identified
significance, it is only the case for the beam size of 3 mm, the results that
is consistent with literature, indicating it as the beam size starting with
which the long-term lesions start appearing (Bertho et al.,, 2020). Multiple
significant associations have been identified with the bootstrap estimators in
the ipsilateral lung and in the right lung for the beam size of 3 mm. However,
none of the genes has been identified as significantly linked to lesions in
the case of 7 mm. This may be caused by the relationship being of non-linear
nature, as appears to be the case for the gene TNF, presented in Figure 5.

Finally, among the cases where a significant relationship has been detected,
the estimated values of the slope are always positive, which indicates a
general radio-induced up-regulation trend. These values are in general bigger
in case of the whole-thorax irradiation and within the patch than for the
ipsilateral or right lung for SBRT, which is also in line with biological
knowledge, since the genes are expected to be more strongly up-regulated in
the areas of injury than further away. This effect is illustrated in Figure 6
on the example of the linear model prediction made for the gene IL6 .

Figure 6: Linear model prediction of septal thickness based on IL6 expression,
plotted for different locations and beam sizes, with the results from two
bootstrap estimators.

##  6 Discussion

This work focuses on a statistical framework designed for extracting
dependencies form in vivo experiments, specifically introducing linear
regression estimators in the context where the predictor and the predicted
variables are never jointly observed. In this work we chose the basic linear
multivariate setting, prioritizing simplicity and computational feasibility.
Particularly, the estimator based on the method of moments makes no hypotheses
on data distribution and can be calculated explicitly. The estimator based on
optimal transport includes a simple optimization problem, and is based on the
Gaussian form of the Wasserstein distance but does not technically require the
data to be Gaussian, seeking to approximate them with Gaussian variables in
whatever case.

However, these approaches are inapplicable in the cases where linear
relationship hypothesis cannot be satisfied. For instance, it is the case with
predicting survival data with some continuous biomarker, which is of
particular interest in research into radio-induced adverse effects. To be able
to consider such scenarios, our model can be extended to a more general case,
namely with generalized linear model. The optimal transport estimator appears
promising in this context given the fact that Wasserstein distance allows to
compare probability distributions of different nature (for example, continuous
and discrete).

Another further research direction lies in investigating alternative methods
based on integrated likelihood and Bayesian approaches, which are likely to
produce better results in many cases but require putting priors on
distributions. Concerning a likelihood-based approach, the joint likelihood to
be maximized can be written in a factorized form due to the independence of
$X$ and $Y$:
$L(x,y;\beta_{0},\beta_{1},\mu_{x}^{k},\sigma_{\epsilon},\sigma_{x})=L_{X}(x;\mu_{x},\sigma_{x})L_{Y}(y;\beta_{0},\beta_{1},\mu_{x},\sigma_{\epsilon},\sigma_{x})$
for a group $k$ in the univariate setting. The marginal likelihoods can be
expressed explicitly, assuming that the model (1) is true, and that
$X|(G=k)\sim\mathcal{N}(\mu_{x}^{k},\sigma_{x,k}^{2})$, and thus
$Y|(G=k)\sim\mathcal{N}(\beta_{0}+\beta_{1}\mu_{x}^{k},\sigma_{\epsilon}^{2}+\beta_{1}^{2}\sigma_{x,k}^{2})$,
with the noise $\epsilon\sim\mathcal{N}(0,\sigma_{\epsilon}^{2})$ independent
of $X$.

Finally, it would be of interest to work on improving theoretical properties
of our estimators, namely correcting the negative bias that appears for both
estimators. The latter is particularly important in the case of the optimal
transport estimator, which may be arising naturally with the Wasserstein
distance (Wasserman,, 2023). Correcting this bias would considerably improve
the estimator, making it competitive with the approaches mentioned earlier
that make numerous assumptions on the data.

## Funding

This work is supported by the European Union through the PO FEDER-FSE
Bourgogne 2014/2020 programs as part of the project ModBioCan2020, and by
Institut de Radioprotection et de Sûreté nucléaire as part of the project
ROSIRIS.

Conflict of Interest: none declared.

Acknowledgements: The authors warmly thank Dr. Patrick Tardivel for his
numerous remarks that helped to improve the presentation of the manuscript.

## References

  * Bertho et al., (2020) Bertho, A., Santos, M. D., Buard, V., Paget, V., Guipaud, O., Tarlet, G., Milliat, F., and François, A. (2020).  Preclinical model of stereotactic ablative lung irradiation using arc delivery in the mouse: Effect of beam size changes and dose effect at constant collimation.  International Journal of Radiation Oncology, Biology, Physics, 107(3):548–562.  Publisher: Elsevier.
  * Carrig et al., (2015) Carrig, M. M., Manrique-Vallier, D., Ranby, K. W., Reiter, J. P., and Hoyle, R. H. (2015).  A nonparametric, multiple imputation-based method for the retrospective integration of data sets.  Multivariate Behavioral Research, 50(4):383–397.
  * Cudeck, (2000) Cudeck, R. (2000).  An estimate of the covariance between variables which are not jointly observed.  Psychometrika, 65(4):539–546.
  * Magnus and Neudecker, (2019) Magnus, J. R. and Neudecker, H. (2019).  Matrix differential calculus with applications in statistics and econometrics.  Wiley Ser. Probab. Stat. Hoboken, NJ: John Wiley & Sons, 3rd updated edition edition.
  * Massa and Riccomagno, (2017) Massa, M. S. and Riccomagno, E. (2017).  Algebraic representations of Gaussian Markov combinations.  Bernoulli, 23(1):626–644.  Publisher: Bernoulli Society for Mathematical Statistics and Probability.
  * Mitsuhiro and Hoshino, (2020) Mitsuhiro, M. and Hoshino, T. (2020).  Kernel canonical correlation analysis for data combination of multiple-source datasets.  Japanese Journal of Statistics and Data Science, 3(2):651–668.
  * Mitsuhiro and Hoshino, (2021) Mitsuhiro, M. and Hoshino, T. (2021).  Bayesian data combination model with Gaussian process latent variable model for mixed observed variables under NMAR missingness.  arXiv:2109.00462 [stat].
  * Newey and McFadden, (1994) Newey, W. K. and McFadden, D. (1994).  Large sample estimation and hypothesis testing.  In Handbook of econometrics, Vol. IV, volume 2 of Handbooks in Econom., pages 2111–2245. North-Holland, Amsterdam.
  * Panaretos and Zemel, (2019) Panaretos, V. M. and Zemel, Y. (2019).  Statistical aspects of Wasserstein distances.  Annu. Rev. Stat. Appl., 6:405–431.
  * Shao and Tu, (1995) Shao, J. and Tu, D. (1995).  The Jackknife and Bootstrap.  Springer Ser. Stat. New York, NY: Springer-Verlag.
  * Triantafillou et al., (2010) Triantafillou, S., Tsamardinos, I., and Tollis, I. (2010).  Learning causal structure from overlapping variable sets.  In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 860–867. JMLR Workshop and Conference Proceedings.  ISSN: 1938-7228.
  * Tsamardinos et al., (2012) Tsamardinos, I., Triantafillou, S., and Lagani, V. (2012).  Towards Integrative Causal Analysis of Heterogeneous Data Sets and Studies.  Journal of Machine Learning Research, 13(39):1097–1157.
  * van der Vaart, (1998) van der Vaart, A. W. (1998).  Asymptotic statistics.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge.
  * Wasserman, (2023) Wasserman, L. (2023).  Optimal transport in statistics.  54th Journées de Statistique of the SFdS, Université libre de Bruxelles.

##  Appendix A Some classical theorems in asymptotic statistics

A proof of the classical continuous mapping theorem can be found in van der
Vaart, (1998) (Theorem 2.3).

######  Theorem A.1.

(Continuous mapping theorem).
Let $g:\mathbb{R}^{d}\to\mathbb{R}^{m}$ be continuous at every point of $C$
such that $\mathbb{P}[X\in C]=1$.
If the sequence of random variables $(X_{n})_{n\geq 1}$ converges in
distribution (resp. probability, resp. almost surely) to $X$ then
$(g(X_{n}))_{n\geq 1}$ converges in distribution (resp. probability, resp.
almost surely) to $g(X)$.

We also recall some well known results that are useful to show the consistency
of estimators $\widehat{\theta}_{n}$ defined as the minimizers of functionals
$Q_{n}(\theta)$ which have some regularity properties at the limit.

######  Theorem A.2.

(Lemma 2.9 in Newey and McFadden, (1994))
Suppose that $\theta\in\Theta$ and $\Theta$ is compact, $Q_{0}(\theta)$ is
continuous and $\forall\theta\in\Theta$, $Q_{n}(\theta)\to Q_{0}(\theta)$ in
probability as $n$ tends to infinity. If there is $\alpha>0$ and
$B_{n}=O_{p}(1)$ such that

$\forall(\widetilde{\theta},\theta)\in\Theta\times\Theta,\
|Q_{n}(\widetilde{\theta})-Q_{n}(\theta)|\leq
B_{n}\|\widetilde{\theta}-\theta\|^{\alpha}$

then

$\sup_{\theta\in\Theta}\left|Q_{n}(\theta)-Q_{0}(\theta)\right|\to 0\ \mbox{
in probability.}$

######  Theorem A.3.

(Theorem 2.1 in Newey and McFadden, (1994))
Suppose that $\theta\in\Theta$ and $\Theta$ is compact, $Q_{0}(\theta)$ is
continuous $\forall\theta\in\Theta$. If $Q_{0}(\theta)$ is uniquely maximized
at $\theta_{0}$ and, as $n$ tends to infinity,
$\sup_{\theta\in\Theta}\left|Q_{n}(\theta)-Q_{0}(\theta)\right|\to 0$ in
probability, then $\widehat{\theta}_{n}\to\theta_{0}$ in probability.

Under additional hypotheses, we also get the asymptotic normality of the
sequence of estimators $\widehat{\theta}_{n}$ of $\theta_{0}$. we denote by
$\nabla_{00}Q_{n}(\theta)$ the Hessian matrix of functional $Q_{n}$ evaluated
at $\theta$.

######  Theorem A.4.

(Theorem 3.1 in Newey and McFadden, (1994))
Suppose that $\widehat{\theta}_{n}\to\theta_{0}$ in probability, (i)
$\theta_{0}$ is an interior point of $\Theta$, (ii) $Q_{n}(\theta)$ is twice
differentiable in a neighborhood $\mathcal{N}$ of $\theta_{0}$, (iii)
$\sqrt{n}\nabla_{0}Q_{n}(\theta_{0})\rightsquigarrow\mathcal{N}\left(0,\boldsymbol{\Sigma}\right)$,
(iv) there is $\mathbf{H}(\theta)$ continuous at $\theta_{0}$ and
$\sup_{\theta\in\mathcal{N}}\|\nabla_{00}Q_{n}(\theta)-\mathbf{H}(\theta)\|\to
0$ in probability (v) $\mathbf{H}=\mathbf{H}(\theta_{0})$ is non singular.
Then

$\sqrt{n}\left(\widehat{\theta}_{n}-\theta_{0}\right)\rightsquigarrow\mathcal{N}\left(0,\mathbf{H}^{-1}\boldsymbol{\Sigma}\mathbf{H}^{-1}\right)$

We also recall the central limit theorem for bootstrap means (see Theorem 23.4
in van der Vaart, (1998) for a proof).

######  Theorem A.5.

(CLT for bootstrap means)
Let $X_{1},X_{2},\ldots$ be i.i.d. random vectors with mean $\boldsymbol{\mu}$
and covariance matrix $\boldsymbol{\Gamma}$. Then conditionally on
$X_{1},X_{2},\ldots$, for almost every sequence $X_{1},X_{2},\ldots$

$\sqrt{n}\left(\overline{X}^{*}_{n}-\overline{X}_{n}\right)\rightsquigarrow\mathcal{N}\left(0,\boldsymbol{\Gamma}\right)$

where $\overline{X}_{n}$ is the empirical mean and $\overline{X}^{*}_{n}$ is
the empirical mean of $n$ independent observations drawn from the empirical
distribution.

##  Appendix B Proofs

###### Proof.

of Lemma 2.3
First note that the assumptions ${\mathbb{E}}(Y^{2})<+\infty$ and
${\mathbb{E}}(\|\mathbf{X}\|^{2})<+\infty$ ensure the existence of
$\sigma_{Y}^{2}$ and $\boldsymbol{\Gamma}_{X}$. From the law of large numbers,
we have that for all $k\in\\{1,\ldots,K\\},\
\widehat{\boldsymbol{\mu}}_{1,X}^{k}\to\boldsymbol{\mu}_{X}^{k}$ and
$\widehat{\mu}_{Y}^{k}\to\mu_{Y}^{k}$ in probability when $n_{\min}$ tends to
infinity.

We deduce from the continuous mapping theorem that
$\widehat{\boldsymbol{\mu}}_{1,X}^{\top}\widehat{\boldsymbol{\mu}}_{1,X}\to\boldsymbol{\mu}_{1,X}^{\top}\mathbf{w}\boldsymbol{\mu}_{1,X}$
and
$\widehat{\boldsymbol{\mu}}_{1,X}^{\top}\mathbf{w}\widehat{\boldsymbol{\mu}}_{Y}\to\boldsymbol{\mu}_{1,X}^{\top}\mathbf{w}\boldsymbol{\mu}_{Y}$
in probability. Under hypothesis $\mathbf{H}_{1}$, the inverse being
continuous in a neighborhood of
$\boldsymbol{\mu}_{1,X}^{\top}\mathbf{w}\boldsymbol{\mu}_{1,X}$ another
application of the continuous mapping theorem gives that
$\left(\widehat{\boldsymbol{\mu}}_{1,X}^{\top}\mathbf{w}\widehat{\boldsymbol{\mu}}_{1,X}\right)^{-1}\to\left(\boldsymbol{\mu}_{1,X}^{\top}\mathbf{w}\boldsymbol{\mu}_{1,X}\right)^{-1}$
and

$\widehat{\boldsymbol{\beta}}^{M}=\left(\widehat{\boldsymbol{\mu}}_{1,X}^{\top}\mathbf{w}\widehat{\boldsymbol{\mu}}_{1,X}\right)^{-1}\widehat{\boldsymbol{\mu}}_{1,X}^{\top}\mathbf{w}\widehat{\boldsymbol{\mu}}_{Y}\to\left(\boldsymbol{\mu}_{1,X}^{\top}\mathbf{w}\boldsymbol{\mu}_{1,X}\right)^{-1}\boldsymbol{\mu}_{1,X}^{\top}\mathbf{w}\boldsymbol{\mu}_{Y}=\boldsymbol{\beta}$

in probability as $n_{\min}$ tends to infinity.

The law of large numbers gives that
$\widehat{\boldsymbol{\Gamma}}_{X}^{2}\to\boldsymbol{\Gamma}_{X}^{2}$ and
$\widehat{\sigma}_{Y}^{2}\to\sigma_{Y}^{2}$ in probability and we deduce, with
another application of the continuous mapping theorem, that
$\widehat{\sigma}_{Y}^{2}-\widehat{\boldsymbol{\beta}}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{2}\widehat{\boldsymbol{\beta}}\to\sigma_{Y}^{2}-\boldsymbol{\beta}^{\top}\boldsymbol{\Gamma}_{X}^{2}\boldsymbol{\beta}=\sigma_{\epsilon}^{2}$
in probability as $n_{\min}$ tends to infinity. ∎

###### Proof.

of Lemma 2.4
The proof is based on Lemma 2.9 and Theorem 2.1 in Newey and McFadden, (1994),
which are recalled in Appendix A.

The law of large numbers and the continuous mapping theorem give us that for
all $(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\in\Theta$,
$\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\to\varphi(\boldsymbol{\gamma},\sigma^{2}_{\gamma})$
in probability, when $n_{\min}$ tends to infinity.

Consider now $(\boldsymbol{\alpha},\sigma_{\alpha}^{2})\in\Theta$. We have,

$\displaystyle\left|(\widehat{\mu}_{Y}^{k}-\gamma_{0}-\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\mu}}_{X}^{k})^{2}-(\widehat{\mu}_{Y}^{k}-\alpha_{0}-\boldsymbol{\alpha}_{-0}^{\top}\widehat{\boldsymbol{\mu}}_{X}^{k})^{2}\right|$
$\displaystyle=\left|(\boldsymbol{\alpha}-\boldsymbol{\gamma})^{T}\begin{pmatrix}1\\\
\widehat{\boldsymbol{\mu}}_{X}^{k}\end{pmatrix}\left(2\widehat{\mu}_{Y}^{k}-(\boldsymbol{\alpha}+\boldsymbol{\gamma})^{T}\begin{pmatrix}1\\\
\widehat{\boldsymbol{\mu}}_{X}^{k}\end{pmatrix}\right)\right|$
$\displaystyle\leq\left\lVert\boldsymbol{\alpha}-\boldsymbol{\gamma}\right\rVert
A_{n}^{k},$

with Cauchy-Schwarz inequality and $A_{n,k}=\mathcal{O}_{p}(1)$ because
$\|\widehat{\boldsymbol{\mu}}_{X}^{k}\|=O_{p}(1)$,
$\widehat{\mu}_{Y}^{k}=O_{p}(1)$ and for some constant $C_{1}$ that does not
depend on $\boldsymbol{\alpha}$ and $\boldsymbol{\gamma}$,
$\|\boldsymbol{\alpha}+\boldsymbol{\gamma}\|\leq C_{1}<\infty$ because
$\Theta$ is supposed to be compact.

On the other hand, we have

$\displaystyle\left|\left(\widehat{\sigma}_{Y,k}-\sqrt{\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\gamma}_{-0}+\sigma_{\gamma}^{2}}\right)^{2}-\left(\widehat{\sigma}_{Y,k}-\sqrt{\boldsymbol{\alpha}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\alpha}_{-0}+\sigma_{\alpha}^{2}}\right)^{2}\right|$
$\displaystyle=\left|\sqrt{\boldsymbol{\alpha}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\alpha}_{-0}+\sigma_{\alpha}^{2}}-\sqrt{\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\gamma}_{-0}+\sigma_{\gamma}^{2}}\right|\left(2\widehat{\sigma}_{Y,k}+\sqrt{\boldsymbol{\alpha}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\alpha}_{-0}+\sigma_{\alpha}^{2}}+\sqrt{\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\gamma}_{-0}+\sigma_{\gamma}^{2}}\right)$
$\displaystyle=\left|\sqrt{\boldsymbol{\alpha}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\alpha}_{-0}+\sigma_{\alpha}^{2}}-\sqrt{\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\gamma}_{-0}+\sigma_{\gamma}^{2}}\right|O_{p}(1)$

since $\Theta$ is compact and
$\|\widehat{\boldsymbol{\Gamma}}_{X}^{k}\|_{sp}=O_{p}(1)$, where $\|.\|_{sp}$
denotes the spectral norm. Because
$\boldsymbol{\alpha}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\alpha}_{-0}-\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\gamma}_{-0}=\boldsymbol{\alpha}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\left(\boldsymbol{\alpha}_{-0}-\boldsymbol{\gamma}_{-0}\right)+\left(\boldsymbol{\alpha}_{-0}-\boldsymbol{\gamma}_{-0}\right)^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\gamma}_{-0}$
we have, for some constant $C_{2,k}>0$,

$\displaystyle\left|\boldsymbol{\alpha}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\alpha}_{-0}-\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\gamma}_{-0}\right|$
$\displaystyle\leq
C_{2,k}\left\|\widehat{\boldsymbol{\Gamma}}_{X}^{k}\right\|_{sp}\|\boldsymbol{\alpha}-\boldsymbol{\gamma}\|.$
(8)

Using now the fact that function $x\mapsto\sqrt{x}$ is concave and
differentiable, we have for $x>0$ and $y>0$ that
$\sqrt{y}\leq\sqrt{x}+\frac{y-x}{2\sqrt{x}}$. Thus, if $y>x>0$ then
$0<\sqrt{y}-\sqrt{x}\leq\frac{y-x}{2\sqrt{x}}$ and if $x>y>0$, then
$0<\sqrt{x}-\sqrt{y}\leq\frac{x-y}{2\sqrt{y}}$. Consequently, we have
$|\sqrt{y}-\sqrt{x}|\leq\frac{|x-y|}{2\min(\sqrt{x},\sqrt{y})}$ and we deduce
that,

$\displaystyle\left|\sqrt{\boldsymbol{\alpha}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\alpha}_{-0}+\sigma_{\alpha}^{2}}-\sqrt{\boldsymbol{\gamma}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\gamma}_{-0}+\sigma_{\gamma}^{2}}\right|$
$\displaystyle\leq
B_{n}^{k}\left(\|\boldsymbol{\alpha}-\boldsymbol{\gamma}\|+|\sigma_{\alpha}^{2}-\sigma_{\gamma}^{2}|\right)$
(9)

where $B_{n}^{k}=O_{p}(1)$.

Combining previous inequalities, we get

$\displaystyle\left|\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})-\varphi_{n}(\boldsymbol{\alpha},\sigma_{\alpha}^{2})\right|$
$\displaystyle\leq\left(\|\boldsymbol{\alpha}-\boldsymbol{\gamma}\|+|\sigma_{\alpha}^{2}-\sigma_{\gamma}^{2}|\right)\sum_{k=1}^{K}\pi_{k}\left(B_{n}^{k}+A_{n}^{k}\right),$
(10)

with $\sum_{k=1}^{K}\pi_{k}\left(B_{n}^{k}+A_{n}^{k}\right)=O_{p}(1)$. As a
result, it can be deduced from Lemma 2.9 in Newey and McFadden, (1994) that

$\sup_{(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\in\Theta}\left|\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})-\varphi(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\right|\to
0\quad\mbox{ in probability.}$

We conclude the proof by recalling that
$\varphi(\boldsymbol{\gamma},\sigma^{2}_{\gamma})$ attains its unique minimum
at $(\boldsymbol{\beta},\sigma^{2}_{\epsilon})\in\Theta$ if assumption
$\mathbf{H}_{1}$ is fulfilled, so that
$(\widehat{\boldsymbol{\beta}}^{W},\widehat{\sigma}^{2,W})\to(\boldsymbol{\beta},\sigma^{2}_{\epsilon})$
in probability in view of Theorem 2.1 in Newey and McFadden, (1994). ∎

###### Proof.

of Proposition 2.1

The central limit theorem applies directly to the independent sequences of
independent random variables
$(\mathbf{X}_{1}^{1},\cdots,\mathbf{X}_{n}^{1}),\ldots,(\mathbf{X}_{1}^{K},\cdots,\mathbf{X}_{n}^{K})$
and $(Y_{1}^{1},\cdots,Y_{n}^{1}),,\ldots,(Y_{1}^{K},\cdots,Y_{n}^{K})$ so
that, as $n$ tends to infinity

$\displaystyle\sqrt{n}\begin{pmatrix}\widehat{\boldsymbol{\mu}}_{X}^{1}-\boldsymbol{\mu}_{X}^{1}\\\
\vdots\\\ \widehat{\boldsymbol{\mu}}_{X}^{K}-\boldsymbol{\mu}_{X}^{K}\\\
\widehat{\mu}_{Y}^{1}-\mu_{Y}^{1}\\\ \vdots\\\
\widehat{\mu}_{Y}^{K}-\mu_{Y}^{K}\end{pmatrix}$
$\displaystyle\rightsquigarrow\mathcal{N}\left(0,\boldsymbol{\Gamma}_{\mu}\right)$
(11)

where $\boldsymbol{\Gamma}_{\mu}$ is a block diagonal matrix, with diagonal
elements
$(\boldsymbol{\Gamma}_{X}^{1},\ldots,\boldsymbol{\Gamma}_{X}^{K},\sigma_{Y,1}^{2},\ldots,\sigma_{Y,K}^{2})$,
with
$\boldsymbol{\Gamma}_{X}^{k}=\mbox{Var}(\mathbf{X}|G=k)={\mathbb{E}}\left(\mathbf{X}^{k}(\mathbf{X}^{k})^{\top}\right)-\boldsymbol{\mu}_{X}^{k}(\boldsymbol{\mu}_{X}^{k})^{\top}$
and $\sigma_{Y,k}^{2}=\mbox{Var}(Y|G=k)$. Consider the application
$g:\mathbb{R}^{dK+K}\to\mathbb{R}^{d+1}$ defined by

$g(\boldsymbol{\mu}_{X}^{1},\ldots,\boldsymbol{\mu}_{X}^{K},\mu_{Y}^{1},\ldots,\mu_{Y}^{K})=\left(\boldsymbol{\mu}_{1,X}^{\top}\mathbf{w}\boldsymbol{\mu}_{1,X}\right)^{-1}\boldsymbol{\mu}_{1,X}^{\top}\mathbf{w}\boldsymbol{\mu}_{Y}.$

Application $g$ is differentiable at
$\boldsymbol{\theta}=(\boldsymbol{\mu}_{X}^{1},\ldots,\boldsymbol{\mu}_{X}^{K},\mu_{Y}^{1},\ldots,\mu_{Y}^{K})$,
with non null Jacobian matrix denoted by $\mathbf{J}_{\theta}$ (see Chapter 8
and more particularly Theorem 8.3 in Magnus and Neudecker, (2019)). The
application of the Delta method (see Theorem 3.1 in van der Vaart, (1998))
permits to get the asymptotic normality convergence result,

$\displaystyle\sqrt{n}\left(\widehat{\boldsymbol{\beta}}^{M}-\boldsymbol{\beta}\right)$
$\displaystyle\rightsquigarrow\mathcal{N}\left(0,\boldsymbol{\Gamma}_{\beta_{M}}\right),$

where
$\boldsymbol{\Gamma}_{\beta_{M}}=\mathbf{J}_{\theta}\boldsymbol{\Gamma}_{\mu}\mathbf{J}_{\theta}^{\top}$.
∎

###### Proof.

of Proposition 2.2

The proof consists in checking the different points of Theorem A.4. Point (i)
is satisfied by the hypotheses, and the point (ii) follows directly from the
fact $\varphi_{n}(\boldsymbol{\gamma},\sigma^{2})$ is twice-differentiable in
a neighborhood of $(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$. To show that
(iii) is fulfilled, we consider the following expansion, based on the
empirical version of the gradient of $\varphi$,

$\displaystyle\nabla\varphi_{n}$
$\displaystyle=\begin{pmatrix}-2\sum_{k=1}^{K}\pi_{k}\left(\widehat{\mu}_{Y}^{k}-\beta_{0}-\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\mu}}_{X}^{k}\right)\\\
-2\sum_{k=1}^{K}\pi_{k}\left[\left(\widehat{\mu}_{Y}^{k}-\beta_{0}-\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\mu}}_{X}^{k}\right)\widehat{\boldsymbol{\mu}}_{X}^{k}+\left(\frac{\widehat{\sigma}_{Y,k}}{\sqrt{\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}}}-1\right)\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}\right]\\\
\sum_{k=1}^{K}\pi_{k}\left(1-\frac{\widehat{\sigma}_{Y,k}}{\sqrt{\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}}}\right)\end{pmatrix}$
(12)

Since model (1) holds, $\nabla\varphi=0$ and
$\widehat{\mu}_{Y}^{k}-\beta_{0}-\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\mu}}_{X}^{k}=(\widehat{\mu}_{Y}^{k}-{\mu}_{Y}^{k})-\boldsymbol{\beta}_{-0}^{\top}\left(\widehat{\boldsymbol{\mu}}_{X}^{k}-\boldsymbol{\mu}_{X}^{k}\right)$,
we thus deduce with (11) the asymptotic normality of the first component of
the gradient $\nabla\varphi_{n}$, that is to say
$\sqrt{n}\left(-2\sum_{k=1}^{K}\pi_{k}\left(\widehat{\mu}_{Y}^{k}-\beta_{0}-\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\mu}}_{X}^{k}\right)\right)$
converges in distribution to a centered Gaussian distribution. As far as the
second component is concerned, it can be noted that
$\widehat{\boldsymbol{\Gamma}}_{X}^{k}$ converges in probability to
$\boldsymbol{\Gamma}_{X}^{k}$ and by the continuous mapping theorem,
$\sqrt{\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}}\to\sigma_{Y,k}$
in probability. It can also be noted that, under the moment condition
${\mathbb{E}}\left[\|\mathbf{X}\|^{4}|G=k\right]<\infty$, the central limit
theorem gives that
$\sqrt{n}\left(\widehat{\boldsymbol{\Gamma}}_{X}^{k}-\boldsymbol{\Gamma}_{X}^{k}\right)$
converges in distribution to a centered Gaussian multivariate distribution,
and we deduce with the Cramer-Wold device, the continuous mapping theorem and
Slutsky’s theorem that the second component of $\nabla\varphi_{n}$ multiplied
by $\sqrt{n}$ also in distribution to a centered Gaussian random vector. It is
immediate to deduce that the same convergence result holds for the third
component, which is to say that
$\sqrt{n}\left(\sum_{k=1}^{K}\pi_{k}\left(1-\frac{\widehat{\sigma}_{Y,k}}{\sqrt{\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}}}\right)\right)$
converges in distribution to a centered Gaussian random variable. We finally
deduce, with the Cramer-Wold device, that (iii) is fulfilled.

To prove that (iv) also holds, consider the Hessian matrix of functional
$\varphi_{n}$, evaluated at $(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$:

$\displaystyle\nabla_{00}\varphi_{n}$
$\displaystyle=\begin{pmatrix}2&2\left(\sum_{k=1}^{K}\pi_{k}\widehat{\boldsymbol{\mu}}_{X}^{k}\right)^{\top}&0\\\
2\sum_{k=1}^{K}\pi_{k}\widehat{\boldsymbol{\mu}}_{X}^{k}&\widehat{\mathbf{H}}(\boldsymbol{\beta}_{-0})&\sum_{k=1}^{K}\pi_{k}\widehat{\sigma}_{Y,k}\left(\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{-3/2}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}\\\
0&\sum_{k=1}^{K}\pi_{k}\widehat{\sigma}_{Y,k}\left(\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{-3/2}\left(\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{\top}&\frac{1}{2}\sum_{k=1}^{K}\pi_{k}\widehat{\sigma}_{Y,k}\left(\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{-3/2}\end{pmatrix},$

where

$\displaystyle\widehat{\mathbf{H}}(\boldsymbol{\beta}_{-0})=$ $\displaystyle
2\sum_{k=1}^{K}\pi_{k}\bigg{[}\widehat{\sigma}_{Y,k}\left(\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{-3/2}\bigg{[}\left(\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}\right)\left(\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{\top}-\left(\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)\widehat{\boldsymbol{\Gamma}}_{X}^{k}\bigg{]}$
$\displaystyle+\widehat{\boldsymbol{\mu}}_{X}^{k}\left(\widehat{\boldsymbol{\mu}}_{X}^{k}\right)^{\top}+\widehat{\boldsymbol{\Gamma}}_{X}^{k}\bigg{]}.$

By similar arguments as those used to show that
$\varphi_{n}(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$ converges in
probability to $\varphi(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$, we deduce
that $\nabla_{00}\varphi_{n}$ converges in probability to some matrix
$\mathbf{H}(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$, defined as follows

$\displaystyle\mathbf{H}(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$
$\displaystyle=\begin{pmatrix}2&2\left(\sum_{k=1}^{K}\pi_{k}\boldsymbol{\mu}_{X}^{k}\right)^{\top}&0\\\
2\sum_{k=1}^{K}\pi_{k}\boldsymbol{\mu}_{X}^{k}&\mathbf{H}(\boldsymbol{\beta}_{-0})&\sum_{k=1}^{K}\pi_{k}\sigma_{Y,k}\left(\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{-3/2}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\\\
0&\sum_{k=1}^{K}\pi_{k}\sigma_{Y,k}\left(\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{-3/2}\left(\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{\top}&\frac{1}{2}\sum_{k=1}^{K}\pi_{k}\sigma_{Y,k}\left(\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{-3/2}\end{pmatrix}$

where

$\displaystyle\mathbf{H}(\boldsymbol{\beta}_{-0})=$ $\displaystyle
2\sum_{k=1}^{K}\pi_{k}\bigg{(}\sigma_{Y,k}\left(\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{-3/2}\bigg{[}\left(\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)\left(\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{\top}-\left(\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)\boldsymbol{\Gamma}_{X}^{k}\bigg{]}$
$\displaystyle+\boldsymbol{\mu}_{X}^{k}\left(\boldsymbol{\mu}_{X}^{k}\right)^{\top}+\boldsymbol{\Gamma}_{X}^{k}\bigg{)}.$

We now must check that $\mathbf{H}(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$
is a positive definite matrix. For that we show that at the minimizer value
$(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$ its determinant is strictly
positive. We first note that
$\sigma_{Y,k}=\left(\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{1/2}$
so that
$\sigma_{Y,k}\left(\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}\right)^{-3/2}=\frac{1}{\sigma_{Y,k}^{2}}$
and $\mathbf{H}(\boldsymbol{\beta}_{-0})$ can be written in a simpler form,

$\displaystyle\mathbf{H}(\boldsymbol{\beta}_{-0})=$ $\displaystyle
2\sum_{k=1}^{K}\pi_{k}\left[\boldsymbol{\mu}_{X}^{k}\left(\boldsymbol{\mu}_{X}^{k}\right)^{\top}+\frac{1}{\sigma_{Y,k}^{2}}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\left(\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{\top}\right],$
(13)

which is a definite positive matrix under hypothesis $\mathbf{H}_{1}$. Using a
block matrix determinant formula, we have

$\displaystyle\left|\mathbf{H}(\boldsymbol{\beta},\sigma_{\epsilon}^{2})\right|$
$\displaystyle=\begin{vmatrix}2&0\\\
0&\frac{1}{2}\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\end{vmatrix}\
\begin{vmatrix}\mathbf{H}(\boldsymbol{\beta}_{-0})-\mathbf{C}\begin{pmatrix}\frac{1}{2}&0\\\
0&\frac{2}{\sum_{k}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}}\end{pmatrix}\mathbf{C}^{\top}\end{vmatrix}$
(14)

where
$\mathbf{C}=\begin{pmatrix}2\sum_{k=1}^{K}\pi_{k}\boldsymbol{\mu}_{X}^{k}&\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\end{pmatrix}$
and it only has to be verified that the second determinant at the righthand
side of (14) is strictly positive. We now have to show that

$\displaystyle\mathbf{H}(\boldsymbol{\beta}_{-0})-\mathbf{C}\begin{pmatrix}\frac{1}{2}&0\\\
0&\frac{2}{\sum_{k}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}}\end{pmatrix}\mathbf{C}^{\top}=2\sum_{k=1}^{K}\pi_{k}\boldsymbol{\mu}_{X}^{k}\left(\boldsymbol{\mu}_{X}^{k}\right)^{\top}-2\left(\sum_{k=1}^{K}\pi_{k}\boldsymbol{\mu}_{X}^{k}\right)\left(\sum_{k=1}^{K}\pi_{k}\boldsymbol{\mu}_{X}^{k}\right)^{\top}$
$\displaystyle+2\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\left(\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{\top}-\frac{2}{\sum_{k}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}}\left(\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)\left(\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{\top}$
(15)

is a positive matrix. We can remark that by Cauchy Schwarz inequality, for
$\mathbf{u}\in\mathbb{R}^{d}$,

$\displaystyle\mathbf{u}^{\top}\left(\sum_{k=1}^{K}\pi_{k}\boldsymbol{\mu}_{X}^{k}\right)\left(\sum_{k=1}^{K}\pi_{k}\boldsymbol{\mu}_{X}^{k}\right)^{\top}\mathbf{u}=\left(\sum_{k=1}^{K}\pi_{k}\mathbf{u}^{\top}\boldsymbol{\mu}_{X}^{k}\right)^{2}$
$\displaystyle\leq\sum_{k=1}^{K}\pi_{k}\left(\mathbf{u}^{\top}\boldsymbol{\mu}_{X}^{k}\right)^{2}$
$\displaystyle=\mathbf{u}^{\top}\left(\sum_{k=1}^{K}\pi_{k}\boldsymbol{\mu}_{X}^{k}\left(\boldsymbol{\mu}_{X}^{k}\right)^{\top}\right)\mathbf{u}$

using the fact that $\sum_{k}(\sqrt{\pi_{k}})^{2}=1$. It can be noted that if
$\mathbf{u}\neq 0$, previous inequality is strict unless
$\mathbf{u}^{\top}\boldsymbol{\mu}_{X}^{1}=\cdots=\mathbf{u}^{\top}\boldsymbol{\mu}_{X}^{K}$,
which can not happen under hypothesis $\mathbf{H}_{1}$. The second part at the
righthand side of (15) is handled the same way. We have

$\displaystyle\mathbf{u}^{\top}\left(\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)\left(\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{\top}\mathbf{u}$
$\displaystyle=\left(\mathbf{u}^{\top}\left(\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)\right)^{2}$
$\displaystyle\leq\sum_{k=1}^{K}\sqrt{\frac{\pi_{k}}{\sigma_{Y,k}^{2}}}^{2}\sum_{k=1}^{K}\left(\sqrt{\frac{\pi_{k}}{\sigma_{Y,k}^{2}}}^{2}\mathbf{u}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{2}$
$\displaystyle=\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\sum_{k=1}^{K}\frac{\pi_{k}}{\sigma_{Y,k}^{2}}\mathbf{u}^{\top}\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\left(\boldsymbol{\Gamma}_{X}^{k}\boldsymbol{\beta}_{-0}\right)^{\top}\mathbf{u},$

and consequently the determinant of
$\mathbf{H}(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$ is strictly positive.

To finish the proof, it remains to check that in a neighborhood $\mathcal{N}$
of $(\boldsymbol{\beta},\sigma_{\epsilon}^{2})$, we have

$\sup_{(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\in\mathcal{N}}\|\nabla_{00}\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})-\mathbf{H}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\|\to
0\ \mbox{ in probability.}$

This is a direct consequence of the continuous mapping theorem, which gives us
that for all $(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\in\mathcal{N}$,
$\|\nabla_{00}\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})-\mathbf{H}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\|\to
0$ in probability, and the fact that third order partial derivatives of
$\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})$ are bounded in
probability for $(\boldsymbol{\gamma},\sigma^{2}_{\gamma})$ so that Theorem
A.2 can apply. ∎

###### Proof.

of Proposition 3.1

The fact that the bootstrap estimator $\boldsymbol{\beta}^{M,*}$ is strongly
consistent for $\boldsymbol{\beta}$ is a direct consequence of Theorem 3.1 in
Shao and Tu, (1995), noting that

$\widehat{\boldsymbol{\beta}}^{M}=g(\widehat{\boldsymbol{\mu}}_{X}^{1},\ldots,\widehat{\boldsymbol{\mu}}_{X}^{K},\widehat{\mu}_{Y}^{1},\ldots,\widehat{\mu}_{Y}^{K})$

is a continuously differentiable function of means at
$(\boldsymbol{\mu}_{X}^{1},\ldots,\boldsymbol{\mu}_{X}^{K},\mu_{Y}^{1},\ldots,\mu_{Y}^{K})$.
The fact that confidence sets based on the percentile approach are consistent
is proved by checking the assumptions in Theorem 4.1 (iii) Shao and Tu,
(1995), namely the bootstrap estimator $\boldsymbol{\beta}^{M,*}$ is
consistent, $\widehat{\boldsymbol{\beta}}^{M}$ is consistent (Lemma 2.3), with
asymptotic Gaussian distribution (Proposition 2.1). ∎

###### Proof.

of Proposition 3.2

We denote by
$\boldsymbol{\theta}_{0}=(\boldsymbol{\beta},\sigma_{\epsilon}^{2,W})$ the
vector of true parameters, by
$\widehat{\boldsymbol{\theta}}=(\boldsymbol{\beta}^{W},\sigma_{\epsilon}^{2})$
the sequence of minimum Wasserstein distance estimators and by
$\boldsymbol{\theta}^{*}=(\boldsymbol{\beta}^{W,*},\sigma_{\epsilon}^{2,W,*})$
bootstrap estimators of $\boldsymbol{\theta}_{0}$. The vector of parameters
$\boldsymbol{\theta}^{*}$ is the minimizer of functional $\varphi_{n}^{*}$
defined as follows,

$\displaystyle\varphi_{n}^{*}(\boldsymbol{\gamma},\sigma^{2})$
$\displaystyle=\sum_{k=1}^{K}\pi_{k}\left[(\mu_{Y}^{k,*}-\gamma_{0}-\boldsymbol{\gamma}_{-0}^{\top}\boldsymbol{\mu}_{X}^{k,*})^{2}+\left(\sigma_{Y,k}^{*}-\sqrt{\boldsymbol{\gamma}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k,*}\boldsymbol{\gamma}_{-0}+\sigma^{2}}\right)^{2}\right].$
(16)

We first show with arguments similar to those employed in the proof of Lemma
2.4, that $\boldsymbol{\theta}^{*}$ is a consistent estimator for
$\boldsymbol{\theta}_{0}$, based on the fact that $\varphi_{n}^{*}$ is a
smooth function converging to $\varphi$ and the sample mean theorem for
bootstrap (see for example Theorem 23.4 in van der Vaart, (1998)). Indeed, we
first recall that for all
$(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\in\Theta$,
$\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\to\varphi(\boldsymbol{\gamma},\sigma^{2}_{\gamma})$
in probability, when $n_{\min}$ tends to infinity and

$\displaystyle\left|\varphi_{n}^{*}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})-\varphi(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\right|$
$\displaystyle\leq\left|\varphi_{n}^{*}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})-\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\right|+\left|\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})-\varphi(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\right|.$
(17)

Since the bootstrap means converge to the empirical ones we deduce with the
continuous mapping theorem that
$\varphi_{n}^{*}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\to\varphi_{n}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})$
in probability, when $n_{\min}$ tends to infinity, so that
$\varphi_{n}^{*}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\to\varphi(\boldsymbol{\gamma},\sigma^{2}_{\gamma})$.
We also have, as in (10), where empirical means are replaced by bootstrap
means,

$\displaystyle\left|\varphi_{n}^{*}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})-\varphi_{n}^{*}(\boldsymbol{\alpha},\sigma_{\alpha}^{2})\right|$
$\displaystyle\leq\left(\|\boldsymbol{\alpha}-\boldsymbol{\gamma}\|+|\sigma_{\alpha}^{2}-\sigma_{\gamma}^{2}|\right)\sum_{k=1}^{K}\pi_{k}\left(B_{n}^{k,*}+A_{n}^{k,*}\right),$
(18)

for any $(\boldsymbol{\alpha},\sigma_{\alpha}^{2})\in\Theta$, with
$\sum_{k=1}^{K}\pi_{k}\left(B_{n}^{k,*}+A_{n}^{k,*}\right)=O_{p}(1)$. As a
result, we deduce from Lemma 2.4, inequality (17) and Lemma 2.9 in Newey and
McFadden, (1994) that

$\sup_{(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\in\Theta}\left|\varphi_{n}^{*}(\boldsymbol{\gamma},\sigma^{2}_{\gamma})-\varphi(\boldsymbol{\gamma},\sigma^{2}_{\gamma})\right|\to
0\quad\mbox{ in probability.}$

We conclude that $\boldsymbol{\theta}^{*}\to\boldsymbol{\theta}_{0}$ in
probability in view of Theorem 2.1 in Newey and McFadden, (1994).

We now prove that
$\sqrt{n}\left(\boldsymbol{\theta}^{*}-\widehat{\boldsymbol{\theta}}\right)$
and
$\sqrt{n}\left(\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\right)$
have the same asymptotic distribution. By definition of
$\widehat{\boldsymbol{\theta}}$ and Taylor expansion we have

$\nabla\varphi_{n}(\widehat{\boldsymbol{\theta}})=\nabla\varphi_{n}(\boldsymbol{\theta}_{0})+\nabla_{00}\varphi_{n}(\overline{\boldsymbol{\theta}})\left(\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\right)=0$
(19)

where $\overline{\boldsymbol{\theta}}$ belongs, componentwise, to the segment
between $\boldsymbol{\theta}_{0}$ and $\widehat{\boldsymbol{\theta}}$. We have
a similar expansion for boostrap estimators, as well as

$\nabla\varphi_{n}^{*}(\boldsymbol{\theta}^{*})=\nabla\varphi_{n}^{*}(\boldsymbol{\theta}_{0})+\nabla_{00}^{*}\varphi_{n}(\overline{\boldsymbol{\theta}}^{*})\left(\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_{0}\right)=0$
(20)

where $\overline{\boldsymbol{\theta}}^{*}$ belongs, componentwise, to the
segment between $\boldsymbol{\theta}_{0}$ and $\boldsymbol{\theta}^{*}$.
Combining (19) and (20), we deduce

$\displaystyle\boldsymbol{\theta}^{*}-\widehat{\boldsymbol{\theta}}$
$\displaystyle=\left(\nabla_{00}^{*}\varphi_{n}(\overline{\boldsymbol{\theta}}^{*})\right)^{-1}\nabla\varphi_{n}^{*}(\boldsymbol{\theta}_{0})-\left(\nabla_{00}\varphi_{n}(\overline{\boldsymbol{\theta}})\right)^{-1}\nabla\varphi_{n}(\boldsymbol{\theta}_{0})$
$\displaystyle=\left(\left(\nabla_{00}^{*}\varphi_{n}(\overline{\boldsymbol{\theta}}^{*})\right)^{-1}-\left(\nabla_{00}\varphi_{n}(\overline{\boldsymbol{\theta}})\right)^{-1}\right)\nabla\varphi_{n}^{*}(\boldsymbol{\theta}_{0})+\left(\nabla_{00}\varphi_{n}(\overline{\boldsymbol{\theta}})\right)^{-1}\left(\nabla\varphi_{n}^{*}(\boldsymbol{\theta}_{0})-\nabla\varphi_{n}(\boldsymbol{\theta}_{0})\right)$
(21)

Noticing that $\nabla_{00}^{*}\varphi_{n}(\overline{\boldsymbol{\theta}}^{*})$
and $\nabla_{00}\varphi_{n}(\overline{\boldsymbol{\theta}})$ both tend in
probability to the same limit
$\mathbf{H}(\boldsymbol{\beta},\sigma^{2}_{\epsilon})$ and we have, with
similar arguments as those used in the proof of Proposition 2.2, that
$\nabla\varphi_{n}^{*}(\boldsymbol{\theta}_{0})$ is $O_{p}(n^{-1/2})$, it can
be deduce that

$\displaystyle\boldsymbol{\theta}^{*}-\widehat{\boldsymbol{\theta}}$
$\displaystyle=\left(\nabla_{00}\varphi_{n}(\overline{\boldsymbol{\theta}})\right)^{-1}\left(\nabla\varphi_{n}^{*}(\boldsymbol{\theta}_{0})-\nabla\varphi_{n}(\boldsymbol{\theta}_{0})\right)+o_{P}(n^{-1/2}).$
(22)

Using arguments similar to those employed in the expansion of
$\nabla\varphi_{n}$ in the proof of Proposition 2.2, we make appear the
difference between bootstrap means and empirical means or a differentiable
functional of these quantities:

$\displaystyle\nabla\varphi_{n}^{*}(\boldsymbol{\theta}_{0})-\nabla\varphi_{n}(\boldsymbol{\theta}_{0})=\footnotesize\begin{pmatrix}2\sum_{k=1}^{K}\pi_{k}\left((\widehat{\mu}_{Y}^{k}-{\mu}_{Y}^{k,*})-\beta_{0}-\boldsymbol{\beta}_{-0}^{\top}\left(\widehat{\boldsymbol{\mu}}_{X}-\boldsymbol{\mu}_{X}^{k,*}\right)\right)\\\
\begin{pmatrix}2\sum_{k=1}^{K}\pi_{k}\left[\left(\widehat{\mu}_{Y}^{k}-\mu_{Y}^{k}-\beta_{0}-\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\mu}}_{X}^{k}\right)\widehat{\boldsymbol{\mu}}_{X}^{k}+\left(\frac{\widehat{\sigma}_{Y,k}}{\sqrt{\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}}}-1\right)\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}\right]\\\
-2\sum_{k=1}^{K}\pi_{k}\left[\left(\mu_{Y}^{k,*}-\mu_{Y}^{k}-\beta_{0}-\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\mu}_{X}^{k,*}\right)\boldsymbol{\mu}_{X}^{k,*}+\left(\frac{\sigma_{Y,k}^{*}}{\sqrt{\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k,*}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}}}-1\right)\boldsymbol{\Gamma}_{X}^{k,*}\boldsymbol{\beta}_{-0}\right]\end{pmatrix}\\\
\sum_{k=1}^{K}\left(\frac{\widehat{\sigma}_{Y,k}}{\sqrt{\boldsymbol{\beta}_{-0}^{\top}\widehat{\boldsymbol{\Gamma}}_{X}^{k}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}}}-\frac{\sigma_{Y,k}^{*}}{\sqrt{\boldsymbol{\beta}_{-0}^{\top}\boldsymbol{\Gamma}_{X}^{k,*}\boldsymbol{\beta}_{-0}+\sigma_{\epsilon}^{2}}}\right)\end{pmatrix},$
(23)

which satisfies the central limit theorem for bootstrap means, or the Delta
method for bootstrap estimators (see the Appendix as well as Theorem 23.4 and
Theorem 23.5 in van der Vaart, (1998)). Consequently,
$\nabla\varphi_{n}^{*}(\boldsymbol{\theta}_{0})-\nabla\varphi_{n}(\boldsymbol{\theta}_{0})$
and
$\nabla\varphi_{n}(\boldsymbol{\theta}_{0})-\nabla\varphi(\boldsymbol{\theta}_{0})$
have the same asymptotic distribution. By Slustky’s theorem, the asymptotic
distribution of
$\sqrt{n}\left(\boldsymbol{\theta}^{*}-\widehat{\boldsymbol{\theta}}\right)$
is thus the same as the asymptotic distribution of
$\mathbf{H}(\boldsymbol{\beta},\sigma^{2}_{\epsilon})\sqrt{n}\nabla\varphi_{n}(\boldsymbol{\theta}_{0})$,
and we can conclude that
$\sqrt{n}\left(\boldsymbol{\theta}^{*}-\widehat{\boldsymbol{\theta}}\right)$
and
$\sqrt{n}\left(\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\right)$
have also the same asymptotic Gaussian distribution.

∎