# Asymptotic properties of cross-classified sampling designs

Jean Rubin, Guillaume Chauvet

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jean Rubin, Guillaume Chauvet. Asymptotic properties of cross-classified sampling designs. 2024. [hal-04579261](https://hal.science/hal-04579261)

HAL Id: hal-04579261

https://hal.science/hal-04579261

Preprint submitted on 17 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Asymptotic properties of cross-classified sampling designs

Jean Rubin¹,² and Guillaume Chauvet¹

¹Ensai (Irmar), Campus de Ker Lann, Bruz - France

²Insee

May 17, 2024

## Abstract

We investigate the family of cross-classified sampling designs across an arbitrary number of dimensions. We introduce a variance decomposition that enables the derivation of general asymptotic properties for these designs and the development of straightforward and asymptotically unbiased variance estimators. Additionally, we demonstrate the suitability of weighted bootstrap techniques for CCS, given the availability of a weighted bootstrap technique in each dimension. Our conclusions are supported by an extensive simulation study. Finally, we apply the proposed methods to a French longitudinal survey conducted among children.

## 1 Introduction

Multistage sampling designs are commonly used in household and health surveys. In case of two-stage sampling, the population units are grouped into large blocks (e.g. municipalities or counties), called Primary Sampling Units (PSUs), which are sampled at the first stage. At the second stage, a list of population units is obtained inside the selected PSUs, and a sample of these units is drawn. A detailed treatment of multistage sampling may be found in Cochran (1977), Särndal et al. (1992) and Fuller (2011). In some situations, a population unit $k$ is more easily represented as a $D$-uple vector ($k_1, \dots, k_D$). For example, in the ELFE maternity survey (Juillard et al., 2017), a sample is obtained by selecting a sample of maternity units ($k_1$), a sample of days (unit $k_2$), and by crossing the two samples. In this case, a population unit is therefore given by a day-maternity couple ($k_1, k_2$).
# Asymptotic properties of cross-classified sampling designs

Jean Rubin¹,² and Guillaume Chauvet¹

¹Ensai (Irmar), Campus de Ker Lann, Bruz - France

²Insee

May 17, 2024

## Abstract

We investigate the family of cross-classified sampling designs across an arbitrary number of dimensions. We introduce a variance decomposition that enables the derivation of general asymptotic properties for these designs and the development of straightforward and asymptotically unbiased variance estimators. Additionally, we demonstrate the suitability of weighted bootstrap techniques for CCS, given the availability of a weighted bootstrap technique in each dimension. Our conclusions are supported by an extensive simulation study. Finally, we apply the proposed methods to a French longitudinal survey conducted among children.

## 1 Introduction

Multistage sampling designs are commonly used in household and health surveys. In case of two-stage sampling, the population units are grouped into large blocks (e.g. municipalities or counties), called Primary Sampling Units (PSUs), which are sampled at the first stage. At the second stage, a list of population units is obtained inside the selected PSUs, and a sample of these units is drawn. A detailed treatment of multistage sampling may be found in Cochran (1977), Särndal et al. (1992) and Fuller (2011). In some situations, a population unit $k$ is more easily represented as a $D$-uple vector ($k_1, \dots, k_D$). For example, in the ELFE maternity survey (Juillard et al., 2017), a sample is obtained by selecting a sample of maternity units ($k_1$), a sample of days (unit $k_2$), and by crossing the two samples. In this case, a population unit is therefore given by a day-maternity couple ($k_1, k_2$).
In such situations, Ohlsson (1996) introduced the two-dimensional cross-classified sampling (CCS)
design, under which independent samples $S_1$ and $S_2$ are selected in each dimension. By taking the
Cartesian product of $S_1$ and $S_2$, the final sample $S = S_1 \times S_2$ is naturally obtained. Other examples
of cross-classified sampling designs include consumer price index surveys (Dalén and Ohlsson, 1995)
and business surveys (Skinner, 2015). To produce reliable estimators with associated confidence
intervals, some basic statistical properties are needed for cross-classified sampling designs, including
the consistency and the asymptotic normality of Horvitz-Thompson estimators. It is also desirable
to provide (at least approximately) unbiased variance estimators. Thereupon, it would also be of
both theoretical and practical interest to derive appropriate bootstrap methods for cross-classified
sampling. With the notable exception of Skinner (2015), who proposed a bootstrap algorithm for
with-replacement sampling designs in each dimension, this last topic has not been studied in the
literature.

In this work, we extend cross-classified sampling to an arbitrary number of dimensions. Under some
mild conditions, we prove the consistency and the asymptotic normality of the Horvitz-Thompson
(HT)-estimator. Using the Hoeffding-Sobol decomposition (Hoeffding, 1948), we generalize the
variance formula given in Ohlsson (1996) and prove that the variance of the HT-estimator can
be decomposed into a sum of multiple terms with different orders of magnitude. By identifying
the leading terms in this variance decomposition, we obtain simple, consistent variance estimators.
This decomposition is also used to derive bootstrap methods suitable for cross-classified sampling.

The article is structured as follows. In Section 2, we introduce our notations and outline our main
assumptions. Section 3 utilizes the Hoeffding-Sobol decomposition to derive a comprehensive vari-
ance decomposition and establish the consistency and asymptotic normality of the HT-estimator.
We also present straightforward and consistent variance estimators, which are exemplified in Sec-
tion 4. Section 5 introduces and examines a weighted bootstrap method, while its application is
illustrated in Section 6. Section 7 presents the results of a simulation study. In Section 8, we apply
the proposed methods to the ELFE survey. The proofs are provided in the Appendix.
In such situations, Ohlsson (1996) introduced the two-dimensional cross-classified sampling (CCS)
design, under which independent samples $S_1$ and $S_2$ are selected in each dimension. By taking the
Cartesian product of $S_1$ and $S_2$, the final sample $S = S_1 \times S_2$ is naturally obtained. Other examples
of cross-classified sampling designs include consumer price index surveys (Dalén and Ohlsson, 1995)
and business surveys (Skinner, 2015). To produce reliable estimators with associated confidence
intervals, some basic statistical properties are needed for cross-classified sampling designs, including
the consistency and the asymptotic normality of Horvitz-Thompson estimators. It is also desirable
to provide (at least approximately) unbiased variance estimators. Thereupon, it would also be of
both theoretical and practical interest to derive appropriate bootstrap methods for cross-classified
sampling. With the notable exception of Skinner (2015), who proposed a bootstrap algorithm for
with-replacement sampling designs in each dimension, this last topic has not been studied in the
literature.

In this work, we extend cross-classified sampling to an arbitrary number of dimensions. Under some
mild conditions, we prove the consistency and the asymptotic normality of the Horvitz-Thompson
(HT)-estimator. Using the Hoeffding-Sobol decomposition (Hoeffding, 1948), we generalize the
variance formula given in Ohlsson (1996) and prove that the variance of the HT-estimator can
be decomposed into a sum of multiple terms with different orders of magnitude. By identifying
the leading terms in this variance decomposition, we obtain simple, consistent variance estimators.
This decomposition is also used to derive bootstrap methods suitable for cross-classified sampling.

The article is structured as follows. In Section 2, we introduce our notations and outline our main
assumptions. Section 3 utilizes the Hoeffding-Sobol decomposition to derive a comprehensive vari-
ance decomposition and establish the consistency and asymptotic normality of the HT-estimator.
We also present straightforward and consistent variance estimators, which are exemplified in Sec-
tion 4. Section 5 introduces and examines a weighted bootstrap method, while its application is
illustrated in Section 6. Section 7 presents the results of a simulation study. In Section 8, we apply
the proposed methods to the ELFE survey. The proofs are provided in the Appendix.
# 2 Multi-dimensional cross-classified sampling

## 2.1 Notation

Suppose that we are interested in a finite population $\mathcal{U} = \prod_{d=1}^{D} \mathcal{U}_d$ which can be seen as the cartesian product of $D$ finite populations of respective sizes $N_d$, $d = 1, \dots, D$. The size of the product population $\mathcal{U}$ is therefore $N = \prod_{d=1}^{D} N_d$. For each unit $\mathbf{k} = (k_1, \dots, k_D) \in \mathcal{U}$, a variable of interest $y$ takes the value $y_{\mathbf{k}}$. We are interested in estimating the population total

$$ Y = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}}. \qquad (2.1) $$

For computations considered in the following sections, it is convenient to introduce some notations for sub-totals. Let $I \subseteq \{1, \dots, D\}$ denote a subset of dimensions, and let $\mathcal{U}_I = \prod_{d \in I} \mathcal{U}_d$ be the product population associated to these dimensions. For any $\mathbf{k}' \in \mathcal{U}_I$, we let

$$ Y_{\mathbf{k}'} = \sum_{\substack{1 \in \mathcal{U} \\ \forall d \in I, l_d = k'_d}} y_1 \qquad (2.2) $$

denote the sub-total of $y$ when the set of coordinates in $\mathcal{U}_I$ remains fixed and equal to $\mathbf{k}'$.

For $d \in \{1, \dots, D\}$, we let $p_d(\cdot)$ denote a sampling design used in the population $\mathcal{U}_d$. Under a $D$-dimensional cross-classified sampling design, $D$ independent samples $S_d$ are selected in the populations $\mathcal{U}_d$, $d = 1, \dots, D$, and their cartesian product $S := \prod_{d=1}^{D} S_d$ is the overall sample. Therefore, the resulting sampling design $p(\cdot)$ is such that

$$ \forall d \in \{1, \dots, D\}, \forall s_d \in \mathcal{P}(\mathcal{U}_d), p\left(\prod_{d=1}^{D} s_d\right) = \prod_{d=1}^{D} p_d(s_d). \qquad (2.3) $$

In the particular case when $D = 2$, we obtain the usual two-dimensional cross-classified design introduced by Ohlsson (1996). It is somewhat similar to a two-stage sampling design where each unit $k_1$ in $\mathcal{U}_1$ would be a PSU, while $\mathcal{U}_2$ would be the sub-population of SSUs for any $k_1$. The main difference is that, due to the independence in the selection of $S_1$ and $S_2$, the same subsample of SSUs is used inside any PSU: in other words, the cross-classified design does not verify the so-called *independence property* for two-stage sampling (Särndal et al., 1992, pages 134-135). However, the *invariance property* is satisfied, namely the sampling design used in a given dimension does not depend on the outcome of the sampling performed in another dimension, see Särndal et al. (1992, Section 4.3.1) and Juillard et al. (2017).
# 2 Multi-dimensional cross-classified sampling

## 2.1 Notation

Suppose that we are interested in a finite population $\mathcal{U} = \prod_{d=1}^{D} \mathcal{U}_d$ which can be seen as the cartesian product of $D$ finite populations of respective sizes $N_d$, $d = 1, \dots, D$. The size of the product population $\mathcal{U}$ is therefore $N = \prod_{d=1}^{D} N_d$. For each unit $\mathbf{k} = (k_1, \dots, k_D) \in \mathcal{U}$, a variable of interest $y$ takes the value $y_{\mathbf{k}}$. We are interested in estimating the population total

$$Y = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}}. \qquad (2.1)$$

For computations considered in the following sections, it is convenient to introduce some notations for sub-totals. Let $I \subseteq \{1, \dots, D\}$ denote a subset of dimensions, and let $\mathcal{U}_I = \prod_{d \in I} \mathcal{U}_d$ be the product population associated to these dimensions. For any $\mathbf{k}' \in \mathcal{U}_I$, we let

$$Y_{\mathbf{k}'} = \sum_{\substack{1 \in \mathcal{U} \\ \forall d \in I, l_d = k'_d}} y_1 \qquad (2.2)$$

denote the sub-total of $y$ when the set of coordinates in $\mathcal{U}_I$ remains fixed and equal to $\mathbf{k}'$.

For $d \in \{1, \dots, D\}$, we let $p_d(\cdot)$ denote a sampling design used in the population $\mathcal{U}_d$. Under a $D$-dimensional cross-classified sampling design, $D$ independent samples $S_d$ are selected in the populations $\mathcal{U}_d$, $d = 1, \dots, D$, and their cartesian product $S := \prod_{d=1}^{D} S_d$ is the overall sample. Therefore, the resulting sampling design $p(\cdot)$ is such that

$$\forall d \in \{1, \dots, D\}, \forall s_d \in \mathcal{P}(\mathcal{U}_d), p\left(\prod_{d=1}^{D} s_d\right) = \prod_{d=1}^{D} p_d(s_d). \qquad (2.3)$$

In the particular case when $D = 2$, we obtain the usual two-dimensional cross-classified design introduced by Ohlsson (1996). It is somewhat similar to a two-stage sampling design where each unit $k_1$ in $\mathcal{U}_1$ would be a PSU, while $\mathcal{U}_2$ would be the sub-population of SSUs for any $k_1$. The main difference is that, due to the independence in the selection of $S_1$ and $S_2$, the same subsample of SSUs is used inside any PSU: in other words, the cross-classified design does not verify the so-called *independence property* for two-stage sampling (Särndal et al., 1992, pages 134-135). However, the *invariance property* is satisfied, namely the sampling design used in a given dimension does not depend on the outcome of the sampling performed in another dimension, see Särndal et al. (1992, Section 4.3.1) and Juillard et al. (2017).
In the population $\mathcal{U}_d$, we let $\delta_k^d$ denote the sample membership indicator of the unit $k$ in the sample $S_d$. We let $\pi_k^d$ denote the probability that $k$ is selected in $S_d$, then $n_d = \sum_{k \in \mathcal{U}_d} \pi_k^d$ denote the expected size of the sample $S_d$, and $\pi_{kl}^d$ denote the probability that units $k$ and $l$ are jointly selected in $S_d$. We will assume in the following that these probabilities are all strictly positive. Finally we also use the notation

$$ \Delta_{kl}^d = \mathrm{Cov}(\delta_k^d, \delta_l^d) = \pi_{kl}^d - \pi_k^d \pi_l^d \quad \text{for any } k, l \in \mathcal{U}_d. \qquad (2.4) $$

## 2.2 Estimation

We consider weighted total estimators as:

$$ \hat{Y} = \sum_{\mathbf{k} \in S} y_{\mathbf{k}} \prod_{d=1}^{D} w_{k_d}^{d}(S_d), \qquad (2.5) $$

where $\{w_k^d(S_d)\}_{k \in \mathcal{U}_d}$ is a set of estimation weights available for each dimension $d \in \{1, \dots, D\}$. We suppose that $\{w_k^d(S_d)\}_{k \in \mathcal{U}_d}$ depends only on the sample $S_d$, and not on the samples selected in the other dimensions. When there is no risk of confusion, we simplify the notation as $w_k^d(S_d) = w_k^d$. We also suppose that for any $d \in \{1, \dots, D\}$ and $k \in \mathcal{U}_d$, $w_k^d = 0$ if $k \notin S_d$. This implies that each unit $\mathbf{k} \in \mathcal{U}$ has an associated weight $w_k = \prod_{d=1}^D w_{k_d}^d$ that can be decomposed as a product of independent weights, leading to the following compact expression:

$$ \hat{Y} = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} w_{\mathbf{k}}. \qquad (2.6) $$

In the important specific case when $w_k^d = \{\pi_k^d\}^{-1}\delta_k^d$ for $k \in \mathcal{U}_d$, we obtain the Horvitz-Thompson (HT) estimator

$$ \hat{Y}_{\pi} = \sum_{\mathbf{k} \in S} \frac{y_{\mathbf{k}}}{\prod_{d=1}^{D} \pi_{k_d}^{d}} = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d=1}^{D} \frac{\delta_{k_d}^{d}}{\pi_{k_d}^{d}}. \qquad (2.7) $$

We may further introduce for $\mathbf{k} \in \mathcal{U}$ the quantities $\delta_{\mathbf{k}} = \prod_{d=1}^{D} \delta_{k_d}^d$ and $\pi_{\mathbf{k}} = \prod_{d=1}^{D} \pi_{k_d}^d$, leading to the compact expression

$$ \hat{Y}_{\pi} = \sum_{\mathbf{k} \in S} \frac{y_{\mathbf{k}}}{\pi_{\mathbf{k}}} = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \frac{\delta_{\mathbf{k}}}{\pi_{\mathbf{k}}}. \qquad (2.8) $$

We will also use the notation $\pi_{kl} = \prod_{d=1}^{D} \pi_{kdl_d}^d$ for any $\mathbf{k}, \mathbf{l} \in \mathcal{U}$.

## 2.3 Assumptions

In this paper, we consider the following assumptions:
In the population $\mathcal{U}_d$, we let $\delta_k^d$ denote the sample membership indicator of the unit $k$ in the sample $S_d$. We let $\pi_k^d$ denote the probability that $k$ is selected in $S_d$, then $n_d = \sum_{k \in \mathcal{U}_d} \pi_k^d$ denote the expected size of the sample $S_d$, and $\pi_{kl}^d$ denote the probability that units $k$ and $l$ are jointly selected in $S_d$. We will assume in the following that these probabilities are all strictly positive. Finally we also use the notation

$$ \Delta_{kl}^d = \mathrm{Cov}(\delta_k^d, \delta_l^d) = \pi_{kl}^d - \pi_k^d \pi_l^d \quad \text{for any } k, l \in \mathcal{U}_d. \qquad (2.4) $$

## 2.2 Estimation

We consider weighted total estimators as:

$$ \hat{Y} = \sum_{\mathbf{k} \in S} y_{\mathbf{k}} \prod_{d=1}^{D} w_{k_d}^{d}(S_d), \qquad (2.5) $$

where $\{w_k^d(S_d)\}_{k \in \mathcal{U}_d}$ is a set of estimation weights available for each dimension $d \in \{1, \dots, D\}$. We suppose that $\{w_k^d(S_d)\}_{k \in \mathcal{U}_d}$ depends only on the sample $S_d$, and not on the samples selected in the other dimensions. When there is no risk of confusion, we simplify the notation as $w_k^d(S_d) = w_k^d$. We also suppose that for any $d \in \{1, \dots, D\}$ and $k \in \mathcal{U}_d$, $w_k^d = 0$ if $k \notin S_d$. This implies that each unit $\mathbf{k} \in \mathcal{U}$ has an associated weight $w_k = \prod_{d=1}^{D} w_{k_d}^d$ that can be decomposed as a product of independent weights, leading to the following compact expression:

$$ \hat{Y} = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} w_{\mathbf{k}}. \qquad (2.6) $$

In the important specific case when $w_k^d = \{\pi_k^d\}^{-1}\delta_k^d$ for $k \in \mathcal{U}_d$, we obtain the Horvitz-Thompson (HT) estimator

$$ \hat{Y}_{\pi} = \sum_{\mathbf{k} \in S} \frac{y_{\mathbf{k}}}{\prod_{d=1}^{D} \pi_{k_d}^{d}} = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d=1}^{D} \frac{\delta_{k_d}^{d}}{\pi_{k_d}^{d}}. \qquad (2.7) $$

We may further introduce for $\mathbf{k} \in \mathcal{U}$ the quantities $\delta_{\mathbf{k}} = \prod_{d=1}^{D} \delta_{k_d}^d$ and $\pi_{\mathbf{k}} = \prod_{d=1}^{D} \pi_{k_d}^d$, leading to the compact expression

$$ \hat{Y}_{\pi} = \sum_{\mathbf{k} \in S} \frac{y_{\mathbf{k}}}{\pi_{\mathbf{k}}} = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \frac{\delta_{\mathbf{k}}}{\pi_{\mathbf{k}}}. \qquad (2.8) $$

We will also use the notation $\pi_{\mathbf{k}l} = \prod_{d=1}^{D} \pi_{kdl_d}^d$ for any $\mathbf{k}, \mathbf{l} \in \mathcal{U}$.

## 2.3 Assumptions

In this paper, we consider the following assumptions:
H1. There exists some constant $\alpha$ such that

$$ \frac{1}{N} \sum_{\mathbf{k} \in \mathcal{U}} (y_{\mathbf{k}})^2 \leq \alpha. \qquad (2.9) $$

H2. For any $d = 1, \dots, D$, we have $n_d \to \infty$ and $N_d \to \infty$, and there exists some constant $f_d \in [0, 1]$ such that

$$ \frac{n_d}{N_d} \to f_d. \qquad (2.10) $$

H3. For any $d = 1, \dots, D$, there exists some constant $\lambda_d > 0$ such that

$$ \forall k \in \mathcal{U}_d, \pi_k^d \geq \lambda_d \frac{n_d}{N_d}. \qquad (2.11) $$

H4. For every $d = 1, \dots, D$, there exists a constant $\gamma_d$ such that

$$ \forall k \neq l \in \mathcal{U}_d, |\Delta_{k,l}^d| \leq \gamma_d \frac{n_d}{N_d^2}. \qquad (2.12) $$

The assumption (H1) pertains to the variable of interest by positing that it possesses a finite moment of order 2. Assumptions (H2)-(H4) are associated with the sampling design. The asymptotic framework is delineated in (H2), wherein it is explicitly assumed that the sample size $n_d \to \infty$ across all dimensions. This assumption is imperative for achieving the consistency of the Horvitz-Thompson estimator, as described in Section 3.2. Furthermore, (H3) presupposes that within each dimension $d$, the first-order inclusion probabilities maintain a lower bound of order $n_d/N_d$. On the other hand, assumption (H4) addresses the second-order inclusion probabilities. The parameter $|\Delta_{k,l}^d|$ can be interpreted as a gauge of interdependence in the selection of units $k, l \in \mathcal{U}_d$ within the sample $S_d$. These quantities assume a value of 0 if the units are independently selected in $S_d$, a concept akin to Poisson sampling. Notably, this assumption holds true for simple random sampling and rejective sampling (Hájek, 1964). In summary, assumptions (H1)-(H4) adhere to conventional standards. Specifically, in the realm of two-dimensional cross-classified sampling, they align with assumptions (H1)-(H3) outlined in Juillard et al. (2017).

# 3 Properties of total estimators

## 3.1 Hoeffding-Sobol variance decomposition of the CCS

From equation (2.5), we can derive the following general variance formula:

$$ V_p(\hat{Y}) = \sum_{\mathbf{k},l \in \mathcal{U}} y_{\mathbf{k}} y_{l} \left[ \prod_{d=1}^{D} E_p(w_{k_d}^d w_{l_d}^d) - \prod_{d=1}^{D} E_p(w_{k_d}^d) E_p(w_{l_d}^d) \right]. \qquad (3.1) $$
H1. There exists some constant $\alpha$ such that

$$ \frac{1}{N} \sum_{\mathbf{k} \in \mathcal{U}} (y_{\mathbf{k}})^2 \leq \alpha. \qquad (2.9) $$

H2. For any $d = 1, \dots, D$, we have $n_d \to \infty$ and $N_d \to \infty$, and there exists some constant $f_d \in [0, 1]$ such that

$$ \frac{n_d}{N_d} \to f_d. \qquad (2.10) $$

H3. For any $d = 1, \dots, D$, there exists some constant $\lambda_d > 0$ such that

$$ \forall k \in \mathcal{U}_d, \pi_k^d \geq \lambda_d \frac{n_d}{N_d}. \qquad (2.11) $$

H4. For every $d = 1, \dots, D$, there exists a constant $\gamma_d$ such that

$$ \forall k \neq l \in \mathcal{U}_d, |\Delta_{k,l}^d| \leq \gamma_d \frac{n_d}{N_d^2}. \qquad (2.12) $$

The assumption (H1) pertains to the variable of interest by positing that it possesses a finite moment of order 2. Assumptions (H2)-(H4) are associated with the sampling design. The asymptotic framework is delineated in (H2), wherein it is explicitly assumed that the sample size $n_d \to \infty$ across all dimensions. This assumption is imperative for achieving the consistency of the Horvitz-Thompson estimator, as described in Section 3.2. Furthermore, (H3) presupposes that within each dimension $d$, the first-order inclusion probabilities maintain a lower bound of order $n_d/N_d$. On the other hand, assumption (H4) addresses the second-order inclusion probabilities. The parameter $|\Delta_{k,l}^d|$ can be interpreted as a gauge of interdependence in the selection of units $k, l \in \mathcal{U}_d$ within the sample $S_d$. These quantities assume a value of 0 if the units are independently selected in $S_d$, a concept akin to Poisson sampling. Notably, this assumption holds true for simple random sampling and rejective sampling (Hájek, 1964). In summary, assumptions (H1)-(H4) adhere to conventional standards. Specifically, in the realm of two-dimensional cross-classified sampling, they align with assumptions (H1)-(H3) outlined in Juillard et al. (2017).

# 3 Properties of total estimators

## 3.1 Hoeffding-Sobol variance decomposition of the CCS

From equation (2.5), we can derive the following general variance formula:

$$ V_p(\hat{Y}) = \sum_{\mathbf{k},l \in \mathcal{U}} y_{\mathbf{k}} y_l \left[ \prod_{d=1}^{D} E_p(w_{k_d}^d w_{l_d}^d) - \prod_{d=1}^{D} E_p(w_{k_d}^d) E_p(w_{l_d}^d) \right]. \qquad (3.1) $$
The proof is provided in Appendix A.1. However, the unbiased variance estimator derived from equation (3.1) may not be practical in many sampling designs, as illustrated by Ohlsson (1996). To explore the statistical properties of estimators and propose simpler variance estimators, we adopt an alternative variance decomposition approach, extending Ohlsson's approach for the two-dimensional case. We express $\hat{Y}$ as a sum of uncorrelated components, as outlined in Proposition 1, with its proof given in Appendix A.2. This decomposition, inspired by the works of Hoeffding (1948), is commonly known as the Hoeffding-Sobol or functional ANOVA decomposition.

**Proposition 1.** We express $\hat{Y}$ as:

$$ \hat{Y} = \sum_{I \subseteq \{1, \dots, D\}} \hat{Y}^I, \qquad (3.2) $$

where for any subset $I \subseteq \{1, \dots, D\}$:

$$ \hat{Y}^I = \sum_{I' \in \mathcal{P}(I)} (-1)^{|I| - |I'|} E_p \{\hat{Y}|(S_d)_{d \in I'}\} \qquad (3.3) $$

$$ = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d \in I} \{w_{k_d}^d - E_p(w_{k_d}^d)\} \prod_{d \notin I} E_p(w_{k_d}^d). \qquad (3.4) $$

Additionally, the components $\{\hat{Y}^I\}_{I \subseteq \{1, \dots, D\}}$ are uncorrelated, and

$$ V_p(\hat{Y}) = \sum_{I \subseteq \{1, \dots, D\}} V_p(\hat{Y}^I), \qquad (3.5) $$

where for any $I \neq \emptyset$:

$$ V_p(\hat{Y}^I) = \sum_{\mathbf{k}, l \in \mathcal{U}} y_{\mathbf{k}} y_{l} \prod_{d \in I} \mathrm{Cov}_p(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} E_p(w_{k_d}^d) E_p(w_{l_d}^d). \qquad (3.6) $$

From Proposition 1, we derive an unbiased variance estimator for $\hat{Y}$, as summarized in Corollary 1, with its proof provided in Appendix A.3.

**Corollary 1.** Assume that the set of weights is adapted to the sampling design, such that for any $d \in \{1, \dots, D\}$:

$$ \forall k \in \mathcal{U}_d, E_p(w_k^d) = 1. \qquad (3.7) $$

Additionally, assume that for any $d \in \{1, \dots, D\}$ and for any pair of units $k, l \in \mathcal{U}_d$, $\widehat{\mathrm{Cov}}^d(w_k^d, w_l^d)$ serves as an unbiased estimator of $\mathrm{Cov}_p(w_k^d, w_l^d)$ constructed from $S_d$, and that $\widehat{\mathrm{Cov}}^d(w_k^d, w_l^d) = 0$ when either $k$ or $l$ is not in $S_d$. Consequently, $V_p(\hat{Y})$ can be estimated unbiasedly as follows

$$ \hat{V}_p(\hat{Y}) = \sum_{I \subseteq \{1, \dots, D\}} \hat{V}_p(\hat{Y}^I), \qquad (3.8) $$
The proof is provided in Appendix A.1. However, the unbiased variance estimator derived from equation (3.1) may not be practical in many sampling designs, as illustrated by Ohlsson (1996). To explore the statistical properties of estimators and propose simpler variance estimators, we adopt an alternative variance decomposition approach, extending Ohlsson's approach for the two-dimensional case. We express $\hat{Y}$ as a sum of uncorrelated components, as outlined in Proposition 1, with its proof given in Appendix A.2. This decomposition, inspired by the works of Hoeffding (1948), is commonly known as the Hoeffding-Sobol or functional ANOVA decomposition.

**Proposition 1.** We express $\hat{Y}$ as:

$$ \hat{Y} = \sum_{I \subseteq \{1, \dots, D\}} \hat{Y}^I, \qquad (3.2) $$

where for any subset $I \subseteq \{1, \dots, D\}$:

$$ \hat{Y}^I = \sum_{I' \in \mathcal{P}(I)} (-1)^{|I| - |I'|} E_p \{ \hat{Y} | (S_d)_{d \in I'} \} \qquad (3.3) $$

$$ = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d \in I} \{w_{k_d}^d - E_p(w_{k_d}^d)\} \prod_{d \notin I} E_p(w_{k_d}^d). \qquad (3.4) $$

Additionally, the components $\{\hat{Y}^I\}_{I \subseteq \{1, \dots, D\}}$ are uncorrelated, and

$$ V_p(\hat{Y}) = \sum_{I \subseteq \{1, \dots, D\}} V_p(\hat{Y}^I), \qquad (3.5) $$

where for any $I \neq \emptyset$:

$$ V_p(\hat{Y}^I) = \sum_{\mathbf{k}, l \in \mathcal{U}} y_{\mathbf{k}} y_{l} \prod_{d \in I} \mathrm{Cov}_p(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} E_p(w_{k_d}^d) E_p(w_{l_d}^d). \qquad (3.6) $$

From Proposition 1, we derive an unbiased variance estimator for $\hat{Y}$, as summarized in Corollary 1, with its proof provided in Appendix A.3.

**Corollary 1.** Assume that the set of weights is adapted to the sampling design, such that for any $d \in \{1, \dots, D\}$:

$$ \forall k \in \mathcal{U}_d, E_p(w_k^d) = 1. \qquad (3.7) $$

Additionally, assume that for any $d \in \{1, \dots, D\}$ and for any pair of units $k, l \in \mathcal{U}_d$, $\widehat{\mathrm{Cov}}^d(w_k^d, w_l^d)$ serves as an unbiased estimator of $\mathrm{Cov}_p(w_k^d, w_l^d)$ constructed from $S_d$, and that $\widehat{\mathrm{Cov}}^d(w_k^d, w_l^d) = 0$ when either $k$ or $l$ is not in $S_d$. Consequently, $V_p(\hat{Y})$ can be estimated unbiasedly as follows

$$ \hat{V}_p(\hat{Y}) = \sum_{I \subseteq \{1, \dots, D\}} \hat{V}_p(\hat{Y}^I), \qquad (3.8) $$
where

$$ \hat{V}_p(\hat{Y}^I) = \sum_{\mathbf{k},l \in S} y_{\mathbf{k}} y_{l} \prod_{d \in I} \widehat{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} \frac{1}{\pi_{k_d, l_d}^d}. \quad (3.9) $$

## 3.2 Properties of the Horvitz-Thompson estimator

We now delve into the HT-estimator, a widely used method in practical scenarios. In this case, the variance formula (3.1) simplifies to

$$ V_p(\hat{Y}_\pi) = \sum_{\mathbf{k},l \in \mathcal{U}} \frac{y_{\mathbf{k}} y_l}{\pi_{\mathbf{k}} \pi_l} (\pi_{\mathbf{k}l} - \pi_{\mathbf{k}} \pi_l). \quad (3.10) $$

By using Proposition 1 and Corollary 1, we derive an alternative variance decomposition for $\hat{Y}_\pi$ alongside an unbiased variance estimator, succinctly presented in Corollary 2. The proof is explained in Appendix B.1.

**Corollary 2.** The variance of the HT-estimator can be expressed as:

$$ V_p(\hat{Y}_\pi) = \sum_{I \subseteq \{1, \dots, D\}} V_p(\hat{Y}_\pi^I), \quad (3.11) $$

where for any $I \subseteq \{1, \dots, D\}$

$$ \hat{Y}_{\pi}^{I} = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d \in I} \left( \frac{\delta_{k_d}^{d}}{\pi_{k_d}^{d}} - 1 \right), \quad (3.12) $$

with the convention that $\hat{Y}_\pi^\emptyset = Y$. Moreover, if $I \neq \emptyset$, we have

$$ V_p(\hat{Y}_\pi^I) = \sum_{\mathbf{k}', l' \in \mathcal{U}_I} \frac{Y_{\mathbf{k}'} Y_{l'}}{\pi_{\mathbf{k}'} \pi_{l'}} \Delta_{\mathbf{k}', l'}, \quad (3.13) $$

where $\Delta_{\mathbf{k}', l'} = \prod_{d \in I} \Delta_{k_d, l_d'}^d$ and $\pi_{\mathbf{k}'} = \prod_{d \in I} \pi_{k_d}^d$ for any $\mathbf{k}', l' \in \mathcal{U}_I$. This variance is estimated unbiasedly by

$$ \hat{V}_p(\hat{Y}_\pi^I) = \sum_{\mathbf{k},l \in S} \frac{y_{\mathbf{k}} y_l}{\pi_{\mathbf{k}l}} \prod_{d \in I} \frac{\Delta_{k_d, l_d}^d}{\pi_{k_d}^d \pi_{l_d}^d}. \quad (3.14) $$

The Hoeffding-Sobol decomposition allows us to express the HT-estimator as the sum of $2^D$ uncorrelated terms, as presented in (3.12), leading to the unbiased variance estimator detailed in (3.14). In the case of $D=2$, this decomposition simplifies to:

$$ \hat{Y}_{\pi} = \hat{Y}_{\pi}^{\emptyset} + \hat{Y}_{\pi}^{\{1\}} + \hat{Y}_{\pi}^{\{2\}} + \hat{Y}_{\pi}^{\{1,2\}}, \quad (3.15) $$
where

$$ \hat{V}_p(\hat{Y}^I) = \sum_{\mathbf{k},l \in S} y_{\mathbf{k}} y_{l} \prod_{d \in I} \widehat{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} \frac{1}{\pi_{k_d, l_d}^d} \quad (3.9) $$

## 3.2 Properties of the Horvitz-Thompson estimator

We now delve into the HT-estimator, a widely used method in practical scenarios. In this case, the variance formula (3.1) simplifies to

$$ V_p(\hat{Y}_\pi) = \sum_{\mathbf{k},l \in \mathcal{U}} \frac{y_{\mathbf{k}} y_l}{\pi_{\mathbf{k}} \pi_l} (\pi_{\mathbf{k}l} - \pi_{\mathbf{k}} \pi_l). \quad (3.10) $$

By using Proposition 1 and Corollary 1, we derive an alternative variance decomposition for $\hat{Y}_\pi$ alongside an unbiased variance estimator, succinctly presented in Corollary 2. The proof is explained in Appendix B.1.

**Corollary 2.** The variance of the HT-estimator can be expressed as:

$$ V_p(\hat{Y}_\pi) = \sum_{I \subseteq \{1, \dots, D\}} V_p(\hat{Y}_\pi^I), \quad (3.11) $$

where for any $I \subseteq \{1, \dots, D\}$

$$ \hat{Y}_{\pi}^{I} = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d \in I} \left( \frac{\delta_{k_d}^{d}}{\pi_{k_d}^{d}} - 1 \right), \quad (3.12) $$

with the convention that $\hat{Y}_{\pi}^{\emptyset} = Y$. Moreover, if $I \neq \emptyset$, we have

$$ V_p(\hat{Y}_\pi^I) = \sum_{\mathbf{k}', l' \in \mathcal{U}_I} \frac{Y_{\mathbf{k}'} Y_{l'}}{\pi_{\mathbf{k}'} \pi_{l'}} \Delta_{\mathbf{k}', l'}, \quad (3.13) $$

where $\Delta_{\mathbf{k}', l'} = \prod_{d \in I} \Delta_{k_d, l_d'}^d$ and $\pi_{\mathbf{k}'} = \prod_{d \in I} \pi_{k_d}^d$ for any $\mathbf{k}', l' \in \mathcal{U}_I$. This variance is estimated unbiasedly by

$$ \hat{V}_p(\hat{Y}_\pi^I) = \sum_{\mathbf{k},l \in S} \frac{y_{\mathbf{k}} y_l}{\pi_{\mathbf{k}l}} \prod_{d \in I} \frac{\Delta_{k_d, l_d}^d}{\pi_{k_d}^d \pi_{l_d}^d}. \quad (3.14) $$

The Hoeffding-Sobol decomposition allows us to express the HT-estimator as the sum of $2^D$ uncorrelated terms, as presented in (3.12), leading to the unbiased variance estimator detailed in (3.14). In the case of $D = 2$, this decomposition simplifies to:

$$ \hat{Y}_{\pi} = \hat{Y}_{\pi}^{\emptyset} + \hat{Y}_{\pi}^{\{1\}} + \hat{Y}_{\pi}^{\{2\}} + \hat{Y}_{\pi}^{\{1,2\}}, \quad (3.15) $$
with

$$
\begin{align}
\hat{Y}_{\pi}^{\{1\}} &= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \left( \frac{\delta_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} - 1 \right), \nonumber \\
\hat{Y}_{\pi}^{\{2\}} &= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \left( \frac{\delta_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} - 1 \right), \tag{3.16} \\
\hat{Y}_{\pi}^{\{1,2\}} &= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \left( \frac{\delta_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} - 1 \right) \left( \frac{\delta_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} - 1 \right). \nonumber
\end{align}
$$

This decomposition mirrors the one proposed in Ohlsson (1996), and Corollary 2 extends his theo-
rem 3.1 to any number of dimensions. This case is further explored in Section 4, focusing particularly
on variance estimation for specific sampling designs.

It is important to specify the orders of magnitude of the terms in the decomposition (3.11). Propo-
sition 2 addresses this, with the proof provided in Appendix B.2.

**Proposition 2.** *Under assumptions (H1), (H3), and (H4), for any non-empty subset I of 1, . . . , D :*

$$
V_p(\hat{Y}_\pi^I) = O\left( \frac{N^2}{\prod_{d \in I} n_d} \right). \qquad (3.17)
$$

Proposition 2 yields significant implications. Firstly, it indicates that the dominant terms in the Hoeffding-Sobol variance decomposition correspond to subsets *I* of size 1, while other terms are of negligible order under assumption (H2). Thus, it is feasible to obtain an approximately unbiased variance estimator by focusing solely on the *D* terms in (3.14) associated with singletons. This leads to the simplified variance estimator:

$$
\hat{V}^{SIMPL}(\hat{Y}_{\pi}) = \sum_{\mathbf{k},l \in S} \frac{y_{\mathbf{k}} y_l}{\pi_{\mathbf{k}l}} \left( \sum_{d=1}^{D} \frac{\Delta_{k_d, l_d}^d}{\pi_{k_d}^d \pi_{l_d}^d} \right). \quad (3.18)
$$

However, this simplified variance estimator neglects the (non-negative) terms of higher order in the
Hoeffding-Sobol decomposition, resulting in negative bias if sample sizes $n_d$ are moderately large.
Secondly, assumption (H2) must be retained for the HT-estimator to remain consistent. In other
words, if the sample size $n_d$ in some dimension $d$ is bounded, the variance of $N^{-1}\hat{Y}_{\pi}^d$ does not
converge to zero, rendering the Horvitz-Thompson estimator inconsistent. Weak consistency and
asymptotic normality of the HT-estimator are established in Proposition 3, with the proof provided
in Appendix B.3.
with

$$
\begin{align}
\hat{Y}_{\pi}^{\{1\}} &= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \left( \frac{\delta_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} - 1 \right), \nonumber \\
\hat{Y}_{\pi}^{\{2\}} &= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \left( \frac{\delta_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} - 1 \right), \tag{3.16} \\
\hat{Y}_{\pi}^{\{1,2\}} &= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \left( \frac{\delta_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} - 1 \right) \left( \frac{\delta_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} - 1 \right). \nonumber
\end{align}
$$

This decomposition mirrors the one proposed in Ohlsson (1996), and Corollary 2 extends his theo-
rem 3.1 to any number of dimensions. This case is further explored in Section 4, focusing particularly
on variance estimation for specific sampling designs.

It is important to specify the orders of magnitude of the terms in the decomposition (3.11). Propo-
sition 2 addresses this, with the proof provided in Appendix B.2.

**Proposition 2.** *Under assumptions (H1), (H3), and (H4), for any non-empty subset I of 1, . . . , D :*

$$
V_p(\hat{Y}_\pi^I) = O\left( \frac{N^2}{\prod_{d \in I} n_d} \right). \qquad (3.17)
$$

Proposition 2 yields significant implications. Firstly, it indicates that the dominant terms in the
Hoeffding-Sobol variance decomposition correspond to subsets *I* of size 1, while other terms are of
negligible order under assumption (H2). Thus, it is feasible to obtain an approximately unbiased
variance estimator by focusing solely on the *D* terms in (3.14) associated with singletons. This
leads to the simplified variance estimator:

$$
\hat{V}^{SIMPL}(\hat{Y}_{\pi}) = \sum_{\mathbf{k},l \in S} \frac{y_{\mathbf{k}} y_l}{\pi_{\mathbf{k}l}} \left( \sum_{d=1}^{D} \frac{\Delta_{k_d, l_d}^d}{\pi_{k_d}^d \pi_{l_d}^d} \right). \quad (3.18)
$$

However, this simplified variance estimator neglects the (non-negative) terms of higher order in the
Hoeffding-Sobol decomposition, resulting in negative bias if sample sizes $n_d$ are moderately large.
Secondly, assumption (H2) must be retained for the HT-estimator to remain consistent. In other
words, if the sample size $n_d$ in some dimension $d$ is bounded, the variance of $N^{-1}\hat{Y}_{\pi}^d$ does not
converge to zero, rendering the Horvitz-Thompson estimator inconsistent. Weak consistency and
asymptotic normality of the HT-estimator are established in Proposition 3, with the proof provided
in Appendix B.3.
**Proposition 3.** Suppose that assumptions (H1)-(H4) hold. Then

$$V_p(N^{-1}\hat{Y}_\pi) = O(n_m^{-1}), \qquad (3.19)$$

where $n_m = \min_{d=1,...,D} n_d$. Additionally, suppose that:

H5. There exists some constant $C > 0$ such that

$$V_p(\hat{Y}_\pi) \geq C N^2 n_m^{-1}. \qquad (3.20)$$

H6. For any $d = 1, ..., D$:

$$\frac{\hat{Y}_{\pi}^{\{d\}}}{\sqrt{V_p(\hat{Y}_{\pi}^{\{d\}})}} \to_{\mathcal{L}} \mathcal{N}(0, 1), \qquad (3.21)$$

where $\to_{\mathcal{L}}$ denotes convergence in distribution.

H7. For any $d = 1, ..., D$, there exists some constant $\gamma_d \ge 0$ such that

$$\frac{V_p(\hat{Y}_{\pi}^{\{d\}})}{V_p(\hat{Y}_{\pi})} \to (\gamma_d)^2. \qquad (3.22)$$

Then

$$\frac{\hat{Y}_{\pi} - Y}{\sqrt{V_p(\hat{Y}_{\pi})}} \to_{\mathcal{L}} \mathcal{N}(0, 1). \qquad (3.23)$$

Assumption (H5) posits that the variance of the HT-estimator is non-vanishing, following the typical order of magnitude $O(N^2n_m^{-1})$. Proposition 3 establishes that if the HT-estimator is asymptotically normally distributed in each dimension, as stated in (H6), then it is also asymptotically normally distributed under CCS. This generalizes Theorem 1 in Juillard et al. (2017) to an arbitrary number of dimensions.

### 3.3 Plug-in variance estimation

Another viable variance estimator arises from the decomposition in equation (3.11). For any term $V_p(\hat{Y}_\pi^I)$ of the variance decomposition, a plug-in estimator based on the expression in (3.13) is obtained:

$$\hat{V}^{\text{PLUG}}(\hat{Y}_{\pi}^{I}) = \sum_{\mathbf{k}', \mathbf{l}' \in S_I} \frac{\hat{Y}_{\mathbf{k}'} \hat{Y}_{\mathbf{l}'} \Delta_{\mathbf{k}'\mathbf{l}'}}{\pi_{\mathbf{k}'} \pi_{\mathbf{l}'} \pi_{\mathbf{k}'\mathbf{l}'}}, \qquad (3.24)$$
**Proposition 3.** Suppose that assumptions (H1)-(H4) hold. Then

$$V_p(N^{-1}\hat{Y}_\pi) = O(n_m^{-1}), \qquad (3.19)$$

where $n_m = \min_{d=1,...,D} n_d$. Additionally, suppose that:

H5. There exists some constant $C > 0$ such that

$$V_p(\hat{Y}_\pi) \geq C N^2 n_m^{-1}. \qquad (3.20)$$

H6. For any $d = 1, ..., D$:

$$\frac{\hat{Y}_{\pi}^{\{d\}}}{\sqrt{V_p(\hat{Y}_{\pi}^{\{d\}})}} \to_{\mathcal{L}} \mathcal{N}(0, 1), \qquad (3.21)$$

where $\to_{\mathcal{L}}$ denotes convergence in distribution.

H7. For any $d = 1, ..., D$, there exists some constant $\gamma_d \ge 0$ such that

$$\frac{V_p(\hat{Y}_{\pi}^{\{d\}})}{V_p(\hat{Y}_{\pi})} \to (\gamma_d)^2. \qquad (3.22)$$

Then

$$\frac{\hat{Y}_{\pi} - Y}{\sqrt{V_p(\hat{Y}_{\pi})}} \to_{\mathcal{L}} \mathcal{N}(0, 1). \qquad (3.23)$$

Assumption (H5) posits that the variance of the HT-estimator is non-vanishing, following the typical order of magnitude $O(N^2n_m^{-1})$. Proposition 3 establishes that if the HT-estimator is asymptotically normally distributed in each dimension, as stated in (H6), then it is also asymptotically normally distributed under CCS. This generalizes Theorem 1 in Juillard et al. (2017) to an arbitrary number of dimensions.

### 3.3 Plug-in variance estimation

Another viable variance estimator arises from the decomposition in equation (3.11). For any term $V_p(\hat{Y}_\pi^I)$ of the variance decomposition, a plug-in estimator based on the expression in (3.13) is obtained:

$$\hat{V}^{\text{PLUG}}(\hat{Y}_{\pi}^{I}) = \sum_{\mathbf{k}', \mathbf{l}' \in S_I} \frac{\hat{Y}_{\mathbf{k}'} \hat{Y}_{\mathbf{l}'} \Delta_{\mathbf{k}'\mathbf{l}'}}{\pi_{\mathbf{k}'} \pi_{\mathbf{l}'} \pi_{\mathbf{k}'\mathbf{l}'}}, \qquad (3.24)$$
where each $k' \in S_I$ in the partial sum $Y_{k'}$ is replaced with the unbiased estimator given by:

$$ \hat{Y}_{k'} = \sum_{\substack{l \in S \\ \forall d \in I, l_d = k'_d}} \frac{y_l}{\prod_{d \notin I} \pi_{l_d}}. \qquad (3.25) $$

This estimator employs a Horvitz-Thompson-like variance estimator, leading to the plug-in variance estimator:

$$ \hat{V}^{PLUG}(\hat{Y}_{\pi}) = \sum_{I \subseteq \{1, \dots, D\}} \hat{V}^{PLUG}(\hat{Y}_{\pi}^{I}). \qquad (3.26) $$

Following the result obtained in Proposition 2, this estimator may be further simplified by restricting the sum in (3.26) to the $D$ terms associated with singletons. This yields the second simplified variance estimator:

$$ \hat{V}^{SIMPL2}(\hat{Y}_{\pi}) = \sum_{d=1}^{D} \sum_{k_d, l_d \in S_d} \frac{\hat{Y}_{k_d} \hat{Y}_{l_d} \Delta_{k_d, l_d}^{d}}{\pi_{k_d}^{d} \pi_{l_d}^{d} \pi_{k_d, l_d}^{d}}. \qquad (3.27) $$

By construction, this approximation considers only first-order variance terms, simplifying higher-order interactions within these variance terms. Proposition 4 reveals that the bias of the plug-in estimators can be expressed in terms of the $V_p(\hat{Y}_{\pi}^I)$, with the proof provided in Appendix B.4.

**Proposition 4.** For every non-empty subset $I$ of $\{1, \dots, D\}$, we have

$$ E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^I)] = V_p(\hat{Y}_{\pi}^I) + \sum_{I' \supseteq I} V_p(\hat{Y}_{\pi}^{I'}). \qquad (3.28) $$

It follows from Proposition 2 that the second term on the right-hand side of (3.28) is asymptotically negligible under assumptions (H1)-(H4). Thus, $\hat{V}^{PLUG}(\hat{Y}_{\pi}^I)$ is asymptotically unbiased for $V_p(\hat{Y}_{\pi}^I)$. Likewise, $\hat{V}^{PLUG}(\hat{Y}_{\pi})$ is asymptotically unbiased for $V_p(\hat{Y}_{\pi})$. Equation (3.28) indicates that the plug-in variance estimator overestimates in particular the variance contribution associated with singletons, which constitutes the major part of the variance. If the sample sizes $n_d$ are moderately large, the plug-in variance estimator is therefore expected to be conservative.

# 4 Illustrations of the simplified estimations

In this section, we demonstrate the proposed variance estimators in the two-dimensional scenario. The general formulas for the case $D = 2$ are initially presented in Section 4.1. We delve into the scenario of simple random sampling in each dimension in Section 4.2, and the scenario of Poisson sampling in each dimension in Section 4.3.
where each $k' \in S_I$ in the partial sum $Y_{k'}$ is replaced with the unbiased estimator given by:

$$ \hat{Y}_{k'} = \sum_{\substack{l \in S \\ \forall d \in I, l_d = k'_d}} \frac{y_l}{\prod_{d \notin I} \pi_{l_d}}. \qquad (3.25) $$

This estimator employs a Horvitz-Thompson-like variance estimator, leading to the plug-in variance estimator:

$$ \hat{V}^{PLUG}(\hat{Y}_{\pi}) = \sum_{I \subseteq \{1, \dots, D\}} \hat{V}^{PLUG}(\hat{Y}_{\pi}^{I}). \qquad (3.26) $$

Following the result obtained in Proposition 2, this estimator may be further simplified by restricting the sum in (3.26) to the $D$ terms associated with singletons. This yields the second simplified variance estimator:

$$ \hat{V}^{SIMPL2}(\hat{Y}_{\pi}) = \sum_{d=1}^{D} \sum_{k_d, l_d \in S_d} \frac{\hat{Y}_{k_d} \hat{Y}_{l_d} \Delta_{k_d, l_d}^{d}}{\pi_{k_d}^{d} \pi_{l_d}^{d} \pi_{k_d, l_d}^{d}}. \qquad (3.27) $$

By construction, this approximation considers only first-order variance terms, simplifying higher-order interactions within these variance terms. Proposition 4 reveals that the bias of the plug-in estimators can be expressed in terms of the $V_p(\hat{Y}_{\pi}^I)$, with the proof provided in Appendix B.4.

**Proposition 4.** For every non-empty subset $I$ of $\{1, \dots, D\}$, we have

$$ E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^I)] = V_p(\hat{Y}_{\pi}^I) + \sum_{I' \supseteq I} V_p(\hat{Y}_{\pi}^{I'}). \qquad (3.28) $$

It follows from Proposition 2 that the second term on the right-hand side of (3.28) is asymptotically negligible under assumptions (H1)-(H4). Thus, $\hat{V}^{PLUG}(\hat{Y}_{\pi}^I)$ is asymptotically unbiased for $V_p(\hat{Y}_{\pi}^I)$. Likewise, $\hat{V}^{PLUG}(\hat{Y}_{\pi})$ is asymptotically unbiased for $V_p(\hat{Y}_{\pi})$. Equation (3.28) indicates that the plug-in variance estimator overestimates in particular the variance contribution associated with singletons, which constitutes the major part of the variance. If the sample sizes $n_d$ are moderately large, the plug-in variance estimator is therefore expected to be conservative.

# 4 Illustrations of the simplified estimations

In this section, we demonstrate the proposed variance estimators in the two-dimensional scenario. The general formulas for the case $D = 2$ are initially presented in Section 4.1. We delve into the scenario of simple random sampling in each dimension in Section 4.2, and the scenario of Poisson sampling in each dimension in Section 4.3.
## 4.1 General case

Expanding upon the specific Hoeffding-Sobol decomposition in the case $D = 2$, as derived in equation (3.15), we derive the variance decomposition:

$$V_p(\hat{Y}_\pi) = V_p(\hat{Y}_\pi^{\{1\}}) + V_p(\hat{Y}_\pi^{\{2\}}) + V_p(\hat{Y}_\pi^{\{1,2\}}), \quad (4.1)$$

where

$$V_p(\hat{Y}_\pi^{\{1\}}) = \sum_{k_1, l_1 \in \mathcal{U}_1} \frac{Y_{k_1 \bullet} Y_{l_1 \bullet}}{\pi_{k_1}^{(1)} \pi_{l_1}^{(1)}} \Delta_{k_1 l_1}^{(1)} \quad \text{with } Y_{k_1 \bullet} = \sum_{k_2 \in \mathcal{U}_2} y_{k_1 k_2}, \\ V_p(\hat{Y}_\pi^{\{2\}}) = \sum_{k_2, l_2 \in \mathcal{U}_2} \frac{Y_{\bullet k_2} Y_{\bullet l_2}}{\pi_{k_2}^{(2)} \pi_{l_2}^{(2)}} \Delta_{k_2 l_2}^{(2)} \quad \text{with } Y_{\bullet k_2} = \sum_{k_1 \in \mathcal{U}_1} y_{k_1 k_2}, \quad (4.2)$$

$$V_p(\hat{Y}_\pi^{\{1,2\}}) = \sum_{k_1, l_1 \in \mathcal{U}_1} \sum_{k_2, l_2 \in \mathcal{U}_2} \frac{y_{k_1 k_2}}{\pi_{k_1}^{(1)} \pi_{k_2}^{(2)}} \frac{y_{l_1 l_2}}{\pi_{l_1}^{(1)} \pi_{l_2}^{(2)}} \Delta_{k_1 l_1}^{(1)} \Delta_{k_2 l_2}^{(2)},$$

as described in Ohlsson (1996, Theorem 3.1). Utilizing equation (3.14), the corresponding unbiased variance estimator is:

$$\hat{V}_p(\hat{Y}_\pi) = \hat{V}_p(\hat{Y}_\pi^{\{1\}}) + \hat{V}_p(\hat{Y}_\pi^{\{2\}}) + \hat{V}_p(\hat{Y}_\pi^{\{1,2\}}), \quad (4.3)$$

where

$$\hat{V}_p(\hat{Y}_\pi^{\{1\}}) = \sum_{k_1, l_1 \in S_1} \frac{\Delta_{k_1 l_1}^{(1)}}{\pi_{k_1 l_1}^{(1)} \pi_{k_1 l_1}^{(1)}} \sum_{k_2, l_2 \in S_2} \frac{y_{k_1 k_2} y_{l_1 l_2}}{\pi_{k_2 l_2}^{(2)}}, \\ \hat{V}_p(\hat{Y}_\pi^{\{2\}}) = \sum_{k_2, l_2 \in S_2} \frac{\Delta_{k_2 l_2}^{(2)}}{\pi_{k_2 l_2}^{(2)} \pi_{k_2 l_2}^{(2)}} \sum_{k_1, l_1 \in S_1} \frac{y_{k_1 k_2} y_{l_1 l_2}}{\pi_{k_1 l_1}^{(1)}}, \quad (4.4)$$

$$\hat{V}_p(\hat{Y}_\pi^{\{1,2\}}) = \sum_{k_1, l_1 \in S_1} \sum_{k_2, l_2 \in S_2} \frac{y_{k_1 k_2} y_{l_1 l_2}}{\pi_{k_1 l_1}^{(1)} \pi_{k_2 l_2}^{(2)}} \frac{\Delta_{k_1 l_1}^{(1)} \Delta_{k_2 l_2}^{(2)}}{\pi_{k_1}^{(1)} \pi_{l_1}^{(1)} \pi_{k_2}^{(2)} \pi_{l_2}^{(2)}}.$$

The first simplified variance estimator provided in (3.18) can be reformulated as:

$$\hat{V}^{SIMPL}(\hat{Y}_{\pi}) = \hat{V}_{p}(\hat{Y}_{\pi}^{\{1\}}) + \hat{V}_{p}(\hat{Y}_{\pi}^{\{2\}}). \quad (4.5)$$

The second simplified variance estimator provided in (3.27) can be reformulated as:

$$\hat{V}^{SIMPL2}(\hat{Y}_{\pi}) = \hat{V}^{PLUG}(\hat{Y}_{\pi}^{\{1\}}) + \hat{V}^{PLUG}(\hat{Y}_{\pi}^{\{2\}}), \quad (4.6)$$

where

$$\hat{V}^{PLUG}(\hat{Y}_{\pi}^{\{1\}}) = \sum_{k_1, l_1 \in S_1} \frac{\hat{Y}_{k_1 \bullet} \hat{Y}_{l_1 \bullet}}{\pi_{k_1}^{(1)} \pi_{l_1}^{(1)}} \frac{\Delta_{k_1 l_1}^{(1)}}{\pi_{k_1 l_1}^{(1)}} \quad \text{with } \hat{Y}_{k_1 \bullet} = \sum_{k_2 \in S_2} \frac{y_{k_1 k_2}}{\pi_{k_2}^{(2)}}, \\ \hat{V}^{PLUG}(\hat{Y}_{\pi}^{\{2\}}) = \sum_{k_2, l_2 \in S_2} \frac{\hat{Y}_{\bullet k_2} \hat{Y}_{\bullet l_2}}{\pi_{k_2}^{(2)} \pi_{l_2}^{(2)}} \frac{\Delta_{k_2 l_2}^{(2)}}{\pi_{k_2 l_2}^{(2)}} \quad \text{with } \hat{Y}_{\bullet k_2} = \sum_{k_1 \in S_1} \frac{y_{k_1 k_2}}{\pi_{k_1}^{(1)}}. \quad (4.7)$$
## 4.1 General case

Expanding upon the specific Hoeffding-Sobol decomposition in the case $D = 2$, as derived in equation (3.15), we derive the variance decomposition:

$$V_p(\hat{Y}_\pi) = V_p(\hat{Y}_\pi^{\{1\}}) + V_p(\hat{Y}_\pi^{\{2\}}) + V_p(\hat{Y}_\pi^{\{1,2\}}), \quad (4.1)$$

where

$$V_p(\hat{Y}_\pi^{\{1\}}) = \sum_{k_1, l_1 \in \mathcal{U}_1} \frac{Y_{k_1 \bullet} Y_{l_1 \bullet}}{\pi_{k_1}^{(1)} \pi_{l_1}^{(1)}} \Delta_{k_1 l_1}^{(1)} \quad \text{with } Y_{k_1 \bullet} = \sum_{k_2 \in \mathcal{U}_2} y_{k_1 k_2}, \\ V_p(\hat{Y}_\pi^{\{2\}}) = \sum_{k_2, l_2 \in \mathcal{U}_2} \frac{Y_{\bullet k_2} Y_{\bullet l_2}}{\pi_{k_2}^{(2)} \pi_{l_2}^{(2)}} \Delta_{k_2 l_2}^{(2)} \quad \text{with } Y_{\bullet k_2} = \sum_{k_1 \in \mathcal{U}_1} y_{k_1 k_2}, \quad (4.2)$$

$$V_p(\hat{Y}_\pi^{\{1,2\}}) = \sum_{k_1, l_1 \in \mathcal{U}_1} \sum_{k_2, l_2 \in \mathcal{U}_2} \frac{y_{k_1 k_2}}{\pi_{k_1}^{(1)} \pi_{k_2}^{(2)}} \frac{y_{l_1 l_2}}{\pi_{l_1}^{(1)} \pi_{l_2}^{(2)}} \Delta_{k_1 l_1}^{(1)} \Delta_{k_2 l_2}^{(2)},$$

as described in Ohlsson (1996, Theorem 3.1). Utilizing equation (3.14), the corresponding unbiased variance estimator is:

$$\hat{V}_p(\hat{Y}_\pi) = \hat{V}_p(\hat{Y}_\pi^{\{1\}}) + \hat{V}_p(\hat{Y}_\pi^{\{2\}}) + \hat{V}_p(\hat{Y}_\pi^{\{1,2\}}), \quad (4.3)$$

where

$$\hat{V}_p(\hat{Y}_\pi^{\{1\}}) = \sum_{k_1, l_1 \in S_1} \frac{\Delta_{k_1 l_1}^{(1)}}{\pi_{k_1 l_1}^{(1)} \pi_{k_1 l_1}^{(1)}} \sum_{k_2, l_2 \in S_2} \frac{y_{k_1 k_2} y_{l_1 l_2}}{\pi_{k_2 l_2}^{(2)}}, \\ \hat{V}_p(\hat{Y}_\pi^{\{2\}}) = \sum_{k_2, l_2 \in S_2} \frac{\Delta_{k_2 l_2}^{(2)}}{\pi_{k_2 l_2}^{(2)} \pi_{k_2 l_2}^{(2)}} \sum_{k_1, l_1 \in S_1} \frac{y_{k_1 k_2} y_{l_1 l_2}}{\pi_{k_1 l_1}^{(1)}}, \quad (4.4)$$

$$\hat{V}_p(\hat{Y}_\pi^{\{1,2\}}) = \sum_{k_1, l_1 \in S_1} \sum_{k_2, l_2 \in S_2} \frac{y_{k_1 k_2} y_{l_1 l_2}}{\pi_{k_1 l_1}^{(1)} \pi_{k_2 l_2}^{(2)}} \frac{\Delta_{k_1 l_1}^{(1)} \Delta_{k_2 l_2}^{(2)}}{\pi_{k_1}^{(1)} \pi_{l_1}^{(1)} \pi_{k_2}^{(2)} \pi_{l_2}^{(2)}}.$$

The first simplified variance estimator provided in (3.18) can be reformulated as:

$$\hat{V}^{SIMPL}(\hat{Y}_{\pi}) = \hat{V}_{p}(\hat{Y}_{\pi}^{\{1\}}) + \hat{V}_{p}(\hat{Y}_{\pi}^{\{2\}}). \quad (4.5)$$

The second simplified variance estimator provided in (3.27) can be reformulated as:

$$\hat{V}^{SIMPL2}(\hat{Y}_{\pi}) = \hat{V}^{PLUG}(\hat{Y}_{\pi}^{\{1\}}) + \hat{V}^{PLUG}(\hat{Y}_{\pi}^{\{2\}}), \quad (4.6)$$

where

$$\hat{V}^{PLUG}(\hat{Y}_{\pi}^{\{1\}}) = \sum_{k_1, l_1 \in S_1} \frac{\hat{Y}_{k_1 \bullet} \hat{Y}_{l_1 \bullet}}{\pi_{k_1}^{(1)} \pi_{l_1}^{(1)}} \frac{\Delta_{k_1 l_1}^{(1)}}{\pi_{k_1 l_1}^{(1)}} \quad \text{with } \hat{Y}_{k_1 \bullet} = \sum_{k_2 \in S_2} \frac{y_{k_1 k_2}}{\pi_{k_2}^{(2)}}, \\ \hat{V}^{PLUG}(\hat{Y}_{\pi}^{\{2\}}) = \sum_{k_2, l_2 \in S_2} \frac{\hat{Y}_{\bullet k_2} \hat{Y}_{\bullet l_2}}{\pi_{k_2}^{(2)} \pi_{l_2}^{(2)}} \frac{\Delta_{k_2 l_2}^{(2)}}{\pi_{k_2 l_2}^{(2)}} \quad \text{with } \hat{Y}_{\bullet k_2} = \sum_{k_1 \in S_1} \frac{y_{k_1 k_2}}{\pi_{k_1}^{(1)}}. \quad (4.7)$$
## 4.2 Simple random sampling

We examine the scenario where simple random sampling of size $n_d$ is employed in each dimension $d = 1, 2$. The sampling fraction in dimension $d$ is denoted as $f_d = n_d/N_d$. The Horvitz-Thompson estimator can be expressed as:

$$ \hat{Y}_{\pi} = \frac{N_1 N_2}{n_1 n_2} \sum_{k_1 \in S_1} \sum_{k_2 \in S_2} y_{k_1 k_2} \quad (4.8) $$

By employing equation (4.1), the Hoeffding-Sobol variance decomposition becomes

$$ V_p(\hat{Y}_{\pi}) = N_1^2 N_2^2 \left[ (1 - f_1) \frac{S_1^2}{n_1} + (1 - f_2) \frac{S_2^2}{n_2} + (1 - f_1)(1 - f_2) \frac{S_{12}^2}{n_1 n_2} \right] \quad (4.9) $$

where

$$ \begin{aligned} S_1^2 &= \frac{1}{N_1 - 1} \sum_{k_1 \in \mathcal{U}_1} (\bar{Y}_{k_1 \bullet} - \bar{Y}_{\bullet\bullet})^2 \quad \text{with} \quad \bar{Y}_{k_1 \bullet} = \frac{1}{N_2} \sum_{k_2 \in \mathcal{U}_2} y_{k_1 k_2}, \\ S_2^2 &= \frac{1}{N_2 - 1} \sum_{k_2 \in \mathcal{U}_2} (\bar{Y}_{\bullet k_2} - \bar{Y}_{\bullet\bullet})^2 \quad \text{with} \quad \bar{Y}_{\bullet k_2} = \frac{1}{N_1} \sum_{k_1 \in \mathcal{U}_1} y_{k_1 k_2}, \end{aligned} \quad (4.10) $$

$$ S_{12}^2 = \frac{1}{N_1 - 1} \frac{1}{N_2 - 1} \sum_{k_1 \in \mathcal{U}_1} \sum_{k_2 \in \mathcal{U}_2} (y_{k_1 k_2} - \bar{Y}_{k_1 \bullet} - \bar{Y}_{\bullet k_2} + \bar{Y}_{\bullet\bullet})^2, $$

and where $\bar{Y}_{\bullet\bullet} = N^{-1}Y$ is the population mean. This expression matches that found in Ohlsson (1996, Example 1.1) and Skinner (2015, Theorem 2.1).

Applying equation (4.3), we derive a term-by-term unbiased estimator. As indicated by Skinner (2015), this process involves substituting each term in (4.10) with the following unbiased estimators:

$$ \begin{aligned} \hat{S}_1^2 &= \frac{1}{n_1 - 1} \sum_{k_1 \in S_1} (\bar{y}_{k_1 \bullet} - \bar{y}_{\bullet\bullet})^2 - (1 - f_2) \frac{\hat{S}_{12}^2}{n_2}, \\ \hat{S}_2^2 &= \frac{1}{n_2 - 1} \sum_{k_2 \in S_2} (\bar{y}_{\bullet k_2} - \bar{y}_{\bullet\bullet})^2 - (1 - f_1) \frac{\hat{S}_{12}^2}{n_1}, \end{aligned} \quad (4.11) $$

$$ \hat{S}_{12}^2 = \frac{1}{n_1 - 1} \frac{1}{n_2 - 1} \sum_{k_1 \in S_1} \sum_{k_2 \in S_2} (y_{k_1 k_2} - \bar{y}_{k_1 \bullet} - \bar{y}_{\bullet k_2} + \bar{y}_{\bullet\bullet})^2, $$

where

$$ \bar{y}_{k_1 \bullet} = \frac{1}{n_2} \sum_{k_2 \in S_2} y_{k_1 k_2} \quad \text{and} \quad \bar{y}_{\bullet k_2} = \frac{1}{n_1} \sum_{k_1 \in S_1} y_{k_1 k_2} \quad (4.12) $$

represent the unbiased estimators of the partial means $\bar{Y}_{k_1 \bullet}$ and $\bar{Y}_{\bullet k_2}$, respectively, with $\bar{y}_{\bullet\bullet}$ denoting the sample mean. It is worth noting that both $\hat{S}_1^2$ and $\hat{S}_2^2$ incorporate correction terms for unbiasedness, which become negligible with sufficiently large sample sizes.
## 4.2 Simple random sampling

We examine the scenario where simple random sampling of size $n_d$ is employed in each dimension $d = 1, 2$. The sampling fraction in dimension $d$ is denoted as $f_d = n_d/N_d$. The Horvitz-Thompson estimator can be expressed as:

$$ \hat{Y}_{\pi} = \frac{N_1 N_2}{n_1 n_2} \sum_{k_1 \in S_1} \sum_{k_2 \in S_2} y_{k_1 k_2} \quad (4.8) $$

By employing equation (4.1), the Hoeffding-Sobol variance decomposition becomes

$$ V_p(\hat{Y}_{\pi}) = N_1^2 N_2^2 \left[ (1 - f_1) \frac{S_1^2}{n_1} + (1 - f_2) \frac{S_2^2}{n_2} + (1 - f_1)(1 - f_2) \frac{S_{12}^2}{n_1 n_2} \right] \quad (4.9) $$

where

$$ \begin{aligned} S_1^2 &= \frac{1}{N_1 - 1} \sum_{k_1 \in \mathcal{U}_1} (\bar{Y}_{k_1 \bullet} - \bar{Y}_{\bullet\bullet})^2 & \text{with } \bar{Y}_{k_1 \bullet} &= \frac{1}{N_2} \sum_{k_2 \in \mathcal{U}_2} y_{k_1 k_2}, \\ S_2^2 &= \frac{1}{N_2 - 1} \sum_{k_2 \in \mathcal{U}_2} (\bar{Y}_{\bullet k_2} - \bar{Y}_{\bullet\bullet})^2 & \text{with } \bar{Y}_{\bullet k_2} &= \frac{1}{N_1} \sum_{k_1 \in \mathcal{U}_1} y_{k_1 k_2}, \end{aligned} \quad (4.10) $$

$$ S_{12}^2 = \frac{1}{N_1 - 1} \frac{1}{N_2 - 1} \sum_{k_1 \in \mathcal{U}_1} \sum_{k_2 \in \mathcal{U}_2} (y_{k_1 k_2} - \bar{Y}_{k_1 \bullet} - \bar{Y}_{\bullet k_2} + \bar{Y}_{\bullet\bullet})^2, $$

and where $\bar{Y}_{\bullet\bullet} = N^{-1}Y$ is the population mean. This expression matches that found in Ohlsson (1996, Example 1.1) and Skinner (2015, Theorem 2.1).

Applying equation (4.3), we derive a term-by-term unbiased estimator. As indicated by Skinner (2015), this process involves substituting each term in (4.10) with the following unbiased estimators:

$$ \begin{aligned} \hat{S}_1^2 &= \frac{1}{n_1 - 1} \sum_{k_1 \in S_1} (\bar{y}_{k_1 \bullet} - \bar{y}_{\bullet\bullet})^2 - (1 - f_2) \frac{\hat{S}_{12}^2}{n_2}, \\ \hat{S}_2^2 &= \frac{1}{n_2 - 1} \sum_{k_2 \in S_2} (\bar{y}_{\bullet k_2} - \bar{y}_{\bullet\bullet})^2 - (1 - f_1) \frac{\hat{S}_{12}^2}{n_1}, \end{aligned} \quad (4.11) $$

$$ \hat{S}_{12}^2 = \frac{1}{n_1 - 1} \frac{1}{n_2 - 1} \sum_{k_1 \in S_1} \sum_{k_2 \in S_2} (y_{k_1 k_2} - \bar{y}_{k_1 \bullet} - \bar{y}_{\bullet k_2} + \bar{y}_{\bullet\bullet})^2, $$

where

$$ \bar{y}_{k_1 \bullet} = \frac{1}{n_2} \sum_{k_2 \in S_2} y_{k_1 k_2} \quad \text{and} \quad \bar{y}_{\bullet k_2} = \frac{1}{n_1} \sum_{k_1 \in S_1} y_{k_1 k_2} \quad (4.12) $$

represent the unbiased estimators of the partial means $\bar{Y}_{k_1 \bullet}$ and $\bar{Y}_{\bullet k_2}$, respectively, with $\bar{y}_{\bullet\bullet}$ denoting the sample mean. It is worth noting that both $\hat{S}_1^2$ and $\hat{S}_2^2$ incorporate correction terms for unbiasedness, which become negligible with sufficiently large sample sizes.
The first proposed variance estimator $\hat{V}^{SIMPL}(\hat{Y}_{\pi})$, provided in (3.18), is obtained by eliminating the estimator of the interaction term in (4.9), resulting in:

$$ \hat{V}^{SIMPL}(\hat{Y}_{\pi}) = N_1^2 N_2^2 \left[ (1 - f_1) \frac{\hat{S}_1^2}{n_1} + (1 - f_2) \frac{\hat{S}_2^2}{n_2} \right]. \quad (4.13) $$

It is important to note that this estimator still incorporates correction terms within $\hat{S}_1^2$ and $\hat{S}_2^2$, necessitating the computation of $\hat{S}_{12}^2$. On the other hand, the second proposed variance estimator $\hat{V}^{SIMPL2}(\hat{Y}_{\pi})$ not only eliminates the interaction term $\hat{S}_{12}^2$ but also these correction terms, yielding:

$$ \hat{V}^{SIMPL2}(\hat{Y}_{\pi}) = N_1^2 N_2^2 \left[ (1 - f_1) \frac{\hat{S}_1^{2, \text{PLUG}}}{n_1} + (1 - f_2) \frac{\hat{S}_2^{2, \text{PLUG}}}{n_2} \right] \quad (4.14) $$

where

$$ \begin{aligned} \hat{S}_1^{2, \text{PLUG}} &= \frac{1}{n_1 - 1} \sum_{k_1 \in S_1} (\bar{y}_{k_1 \bullet} - \bar{y}_{\bullet \bullet})^2, \\ \hat{S}_2^{2, \text{PLUG}} &= \frac{1}{n_2 - 1} \sum_{k_2 \in S_2} (\bar{y}_{\bullet k_2} - \bar{y}_{\bullet \bullet})^2, \end{aligned} \quad (4.15) $$

represent the plug-in estimators of $S_1^2$ and $S_2^2$, respectively.

## 4.3 Poisson sampling

We now consider the scenario where Poisson sampling is applied independently in each dimension $d = 1, 2$. The HT-estimator takes the form:

$$ \hat{Y}_{\pi} = \sum_{k_1 \in S_1} \sum_{k_2 \in S_2} \frac{y_{k_1 k_2}}{\pi_{k_1}^{(1)} \pi_{k_2}^{(2)}}. \quad (4.16) $$

By employing equation (4.1), the Hoeffding-Sobol variance decomposition becomes:

$$ \begin{aligned} V_p(\hat{Y}_{\pi}) = & \sum_{k_1 \in \mathcal{U}_1} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) Y_{k_1 \bullet}^2 + \sum_{k_2 \in \mathcal{U}_2} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) Y_{\bullet k_2}^2 \\ & + \sum_{k_1 \in \mathcal{U}_1} \sum_{k_2 \in \mathcal{U}_2} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) y_{k_1 k_2}^2, \end{aligned} \quad (4.17) $$

where $Y_{k_1 \bullet}$ and $Y_{\bullet k_2}$ are defined in equation (4.2). A similar variance decomposition is presented in Example 3.1 of Ohlsson (1996).
The first proposed variance estimator $\hat{V}^{SIMPL}(\hat{Y}_{\pi})$, provided in (3.18), is obtained by eliminating the estimator of the interaction term in (4.9), resulting in:

$$ \hat{V}^{SIMPL}(\hat{Y}_{\pi}) = N_1^2 N_2^2 \left[ (1 - f_1) \frac{\hat{S}_1^2}{n_1} + (1 - f_2) \frac{\hat{S}_2^2}{n_2} \right]. \quad (4.13) $$

It is important to note that this estimator still incorporates correction terms within $\hat{S}_1^2$ and $\hat{S}_2^2$, necessitating the computation of $\hat{S}_{12}^2$. On the other hand, the second proposed variance estimator $\hat{V}^{SIMPL2}(\hat{Y}_{\pi})$ not only eliminates the interaction term $\hat{S}_{12}^2$ but also these correction terms, yielding:

$$ \hat{V}^{SIMPL2}(\hat{Y}_{\pi}) = N_1^2 N_2^2 \left[ (1 - f_1) \frac{\hat{S}_1^{2, \text{PLUG}}}{n_1} + (1 - f_2) \frac{\hat{S}_2^{2, \text{PLUG}}}{n_2} \right] \quad (4.14) $$

where

$$ \begin{aligned} \hat{S}_1^{2, \text{PLUG}} &= \frac{1}{n_1 - 1} \sum_{k_1 \in S_1} (\bar{y}_{k_1 \bullet} - \bar{y}_{\bullet \bullet})^2, \\ \hat{S}_2^{2, \text{PLUG}} &= \frac{1}{n_2 - 1} \sum_{k_2 \in S_2} (\bar{y}_{\bullet k_2} - \bar{y}_{\bullet \bullet})^2, \end{aligned} \quad (4.15) $$

represent the plug-in estimators of $S_1^2$ and $S_2^2$, respectively.

## 4.3 Poisson sampling

We now consider the scenario where Poisson sampling is applied independently in each dimension $d = 1, 2$. The HT-estimator takes the form:

$$ \hat{Y}_{\pi} = \sum_{k_1 \in S_1} \sum_{k_2 \in S_2} \frac{y_{k_1 k_2}}{\pi_{k_1}^{(1)} \pi_{k_2}^{(2)}}. \quad (4.16) $$

By employing equation (4.1), the Hoeffding-Sobol variance decomposition becomes:

$$ \begin{aligned} V_p(\hat{Y}_{\pi}) = & \sum_{k_1 \in \mathcal{U}_1} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) Y_{k_1 \bullet}^2 + \sum_{k_2 \in \mathcal{U}_2} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) Y_{\bullet k_2}^2 \\ & + \sum_{k_1 \in \mathcal{U}_1} \sum_{k_2 \in \mathcal{U}_2} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) y_{k_1 k_2}^2, \end{aligned} \quad (4.17) $$

where $Y_{k_1 \bullet}$ and $Y_{\bullet k_2}$ are defined in equation (4.2). A similar variance decomposition is presented in Example 3.1 of Ohlsson (1996).
By utilizing equation (4.3), the variance is estimated without bias by substituting each term in (4.17) with an unbiased counterpart, resulting in the variance estimator:

$$ 
\begin{aligned}
\hat{V}_p(\hat{Y}_\pi) = & \sum_{k_1 \in S_1} \frac{1}{\pi_{k_1}^{(1)}} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \left( \hat{Y}_{k_1 \bullet}^2 - \sum_{k_2 \in S_2} \frac{1}{\pi_{k_2}^{(2)}} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) y_{k_1 k_2}^2 \right) \\
& + \sum_{k_2 \in S_2} \frac{1}{\pi_{k_2}^{(2)}} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) \left( \hat{Y}_{\bullet k_2}^2 - \sum_{k_1 \in S_1} \frac{1}{\pi_{k_1}^{(1)}} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) y_{k_1 k_2}^2 \right) \\
& + \sum_{k_1 \in S_1} \sum_{k_2 \in S_2} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) \frac{y_{k_1 k_2}^2}{\pi_{k_1}^{(1)} \pi_{k_2}^{(2)}}
\end{aligned} 
\quad (4.18) $$

where $\hat{Y}_{k_1 \bullet}$ and $\hat{Y}_{\bullet k_2}$ are defined in equation (4.7). It is worth noting that the two first terms in the right-hand side of (4.18) incorporate correction terms for unbiasedness, which become negligible with sufficiently large sample sizes.

The first proposed variance estimator $\hat{V}^{SIMPL}(\hat{Y}_\pi)$, presented in (3.18), is derived by eliminating the third interaction term in (4.18), yielding:

$$ 
\begin{aligned}
\hat{V}^{SIMPL}(\hat{Y}_{\pi}) = & \sum_{k_1 \in S_1} \frac{1}{\pi_{k_1}^{(1)}} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \left( \hat{Y}_{k_1 \bullet}^2 - \sum_{k_2 \in S_2} \frac{1}{\pi_{k_2}^{(2)}} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) y_{k_1 k_2}^2 \right) \\
& + \sum_{k_2 \in S_2} \frac{1}{\pi_{k_2}^{(2)}} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) \left( \hat{Y}_{\bullet k_2}^2 - \sum_{k_1 \in S_1} \frac{1}{\pi_{k_1}^{(1)}} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) y_{k_1 k_2}^2 \right).
\end{aligned} 
\quad (4.19) $$

The second proposed variance estimator $\hat{V}^{SIMPL2}(\hat{Y}_\pi)$ further removes the correction terms, resulting in:

$$ 
\hat{V}^{SIMPL2}(\hat{Y}_{\pi}) = \sum_{k_1 \in S_1} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \frac{\hat{Y}_{k_1 \bullet}^2}{\pi_{k_1}^{(1)}} + \sum_{k_2 \in S_2} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) \frac{\hat{Y}_{\bullet k_2}^2}{\pi_{k_2}^{(2)}}. 
\quad (4.20) $$

# 5 Weighted bootstrap method for CCS

We now introduce a weighted bootstrap method for CSS, assuming that for each dimension $d \in 1, \dots, D$, a weighted bootstrap method tailored to the weights $(w_{k_d}^d)_{k_d \in \mathcal{U}_d}$ is available. Specifically, we consider a set of weights $(w_{k_d}^{d*})_{k_d \in S_d}$ for each $d \in 1, \dots, D$ that satisfies the first and second-order moment constraints outlined by Beaumont and Patak (2012):

$$ \forall k_d \in S_d, E_*(w_{k_d}^{d*}) = w_{k_d}^d, 
\quad (5.1) $$

$$ \forall k_d, l_d \in S_d, \mathrm{Cov}_*(w_{k_d}^{d*}, w_{l_d}^{d*}) = \tilde{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d), 
\quad (5.2) $$
By utilizing equation (4.3), the variance is estimated without bias by substituting each term in (4.17) with an unbiased counterpart, resulting in the variance estimator:

$$ 
\begin{aligned}
\hat{V}_p(\hat{Y}_\pi) = & \sum_{k_1 \in S_1} \frac{1}{\pi_{k_1}^{(1)}} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \left( \hat{Y}_{k_1 \bullet}^2 - \sum_{k_2 \in S_2} \frac{1}{\pi_{k_2}^{(2)}} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) y_{k_1 k_2}^2 \right) \\
& + \sum_{k_2 \in S_2} \frac{1}{\pi_{k_2}^{(2)}} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) \left( \hat{Y}_{\bullet k_2}^2 - \sum_{k_1 \in S_1} \frac{1}{\pi_{k_1}^{(1)}} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) y_{k_1 k_2}^2 \right) \\
& + \sum_{k_1 \in S_1} \sum_{k_2 \in S_2} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) \frac{y_{k_1 k_2}^2}{\pi_{k_1}^{(1)} \pi_{k_2}^{(2)}}
\end{aligned}
\quad (4.18) $$

where $\hat{Y}_{k_1 \bullet}$ and $\hat{Y}_{\bullet k_2}$ are defined in equation (4.7). It is worth noting that the two first terms in the right-hand side of (4.18) incorporate correction terms for unbiasedness, which become negligible with sufficiently large sample sizes.

The first proposed variance estimator $\hat{V}^{SIMPL}(\hat{Y}_\pi)$, presented in (3.18), is derived by eliminating the third interaction term in (4.18), yielding:

$$ 
\begin{aligned}
\hat{V}^{SIMPL}(\hat{Y}_\pi) = & \sum_{k_1 \in S_1} \frac{1}{\pi_{k_1}^{(1)}} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \left( \hat{Y}_{k_1 \bullet}^2 - \sum_{k_2 \in S_2} \frac{1}{\pi_{k_2}^{(2)}} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) y_{k_1 k_2}^2 \right) \\
& + \sum_{k_2 \in S_2} \frac{1}{\pi_{k_2}^{(2)}} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) \left( \hat{Y}_{\bullet k_2}^2 - \sum_{k_1 \in S_1} \frac{1}{\pi_{k_1}^{(1)}} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) y_{k_1 k_2}^2 \right).
\end{aligned}
\quad (4.19) $$

The second proposed variance estimator $\hat{V}^{SIMPL2}(\hat{Y}_\pi)$ further removes the correction terms, resulting in:

$$ \hat{V}^{SIMPL2}(\hat{Y}_\pi) = \sum_{k_1 \in S_1} \left( \frac{1 - \pi_{k_1}^{(1)}}{\pi_{k_1}^{(1)}} \right) \frac{\hat{Y}_{k_1 \bullet}^2}{\pi_{k_1}^{(1)}} + \sum_{k_2 \in S_2} \left( \frac{1 - \pi_{k_2}^{(2)}}{\pi_{k_2}^{(2)}} \right) \frac{\hat{Y}_{\bullet k_2}^2}{\pi_{k_2}^{(2)}}. \quad (4.20) $$

# 5 Weighted bootstrap method for CCS

We now introduce a weighted bootstrap method for CSS, assuming that for each dimension $d \in 1, \dots, D$, a weighted bootstrap method tailored to the weights $(w_{k_d}^d)_{k_d \in \mathcal{U}_d}$ is available. Specifically, we consider a set of weights $(w_{k_d}^{d*})_{k_d \in S_d}$ for each $d \in 1, \dots, D$ that satisfies the first and second-order moment constraints outlined by Beaumont and Patak (2012):

$$ \forall k_d \in S_d, E_*(w_{k_d}^{d*}) = w_{k_d}^d, \quad (5.1) $$

$$ \forall k_d, l_d \in S_d, \mathrm{Cov}_*(w_{k_d}^{d*}, w_{l_d}^{d*}) = \overline{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d), \quad (5.2) $$
where $E_*(\cdot)$, $V_*(\cdot)$ and $Cov_*(\cdot, \cdot)$ denote the (bootstrap) expectation, variance and covariance, conditioned on the original sample $S$.

A natural bootstrap version for CSS of the weighted estimator in (2.5) is:

$$ \hat{Y}^* = \sum_{\mathbf{k} \in S} y_{\mathbf{k}} \prod_{d=1}^{D} w_{k_d}^{d*} \quad (5.3) $$

where the weights are simulated independently in each dimension $d \in 1, \dots, D$. Applying Proposition 1 conditionally on the original sample $S$, we obtain the Hoeffding-Sobol decomposition:

$$ \hat{Y}^* = \sum_{I \subseteq \{1, \dots, D\}} \hat{Y}^{I*}, \quad (5.4) $$

$$ \text{with } \hat{Y}^{I*} = \sum_{\mathbf{k} \in S} y_{\mathbf{k}} \prod_{d \in I} (w_{k_d}^{d*} - w_{k_d}^d) \prod_{d \notin I} w_{k_d}^d \text{ for any nonempty } I \subseteq \{1, \dots, D\}. $$

By applying equation (3.6), we obtain:

$$ \begin{aligned} V_*(\hat{Y}^{I*}) &= \sum_{\mathbf{k}, l \in S} y_{\mathbf{k}} y_l \prod_{d \in I} Cov_*(w_{k_d}^{d*}, w_{l_d}^{d*}) \prod_{d \notin I} E_*(w_{k_d}^{d*}) E_*(w_{l_d}^{d*}) \\ &= \sum_{\mathbf{k}, l \in S} y_{\mathbf{k}} y_l \prod_{d \in I} \widehat{Cov}_p(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} w_{k_d}^d w_{l_d}^d. \end{aligned} \quad (5.5) $$

The last line in (5.5) follows from equations (5.1) and (5.3).

In the case of Horvitz-Thompson estimation, for any $k_d, l_d \in S_d$, we have:

$$ w_{k_d}^d = \frac{1}{\pi_{k_d}^d} \quad \text{and} \quad \widehat{Cov}^d(w_{k_d}^d, w_{l_d}^d) = \frac{\Delta_{k_d l_d}^d}{\pi_{k_d l_d}^d} \frac{1}{\pi_{k_d}^d \pi_{l_d}^d}. $$

For the bootstrap version $\hat{Y}_\pi^*$ of the HT-estimator, equations (5.4) and (5.5) lead to:

$$ \hat{Y}_{\pi}^* = \sum_{I \subseteq \{1, \dots, D\}} \hat{Y}_{\pi}^{I*} \quad (5.6) $$

and

$$ V_*(\hat{Y}_\pi^{I*}) = \sum_{\mathbf{k}, l \in S} y_{\mathbf{k}} y_l \prod_{d \in I} \frac{\Delta_{k_d l_d}^d}{\pi_{k_d l_d}^d} \frac{1}{\pi_{k_d}^d \pi_{l_d}^d} \prod_{d \notin I} \frac{1}{\pi_{k_d}^d \pi_{l_d}^d} \quad (5.7) $$

$$ = \sum_{\mathbf{k}', l' \in S_I} \frac{\hat{Y}_{\mathbf{k}'} \hat{Y}_{l'}}{\pi_{\mathbf{k}'} \pi_{l'}} \frac{\Delta_{\mathbf{k} l'}}{\pi_{\mathbf{k}'} \pi_{l'}}. \quad (5.8) $$

This corresponds to the plug-in variance estimator for $\hat{Y}_\pi^I$ given in (3.24). Therefore, $V_*(\hat{Y}_\pi^*)$ equals the plug-in variance estimator $\hat{V}^{PLUG}(\hat{Y}_\pi)$ as defined in (3.26), and is also asymptotically
where $E_*(\cdot)$, $V_*(\cdot)$ and $Cov_*(\cdot, \cdot)$ denote the (bootstrap) expectation, variance and covariance, conditioned on the original sample $S$.

A natural bootstrap version for CSS of the weighted estimator in (2.5) is:

$$ \hat{Y}^* = \sum_{\mathbf{k} \in S} y_{\mathbf{k}} \prod_{d=1}^{D} w_{k_d}^{d*} \qquad (5.3) $$

where the weights are simulated independently in each dimension $d \in 1, \dots, D$. Applying Proposition 1 conditionally on the original sample $S$, we obtain the Hoeffding-Sobol decomposition:

$$ \hat{Y}^* = \sum_{I \subseteq \{1, \dots, D\}} \hat{Y}^{I*}, \qquad (5.4) $$

$$ \text{with } \hat{Y}^{I*} = \sum_{\mathbf{k} \in S} y_{\mathbf{k}} \prod_{d \in I} (w_{k_d}^{d*} - w_{k_d}^d) \prod_{d \notin I} w_{k_d}^d \text{ for any nonempty } I \subseteq \{1, \dots, D\}. $$

By applying equation (3.6), we obtain:

$$ \begin{aligned} V_*(\hat{Y}^{I*}) &= \sum_{\mathbf{k}, l \in S} y_{\mathbf{k}} y_{l} \prod_{d \in I} Cov_*(w_{k_d}^{d*}, w_{l_d}^{d*}) \prod_{d \notin I} E_*(w_{k_d}^{d*}) E_*(w_{l_d}^{d*}) \\ &= \sum_{\mathbf{k}, l \in S} y_{\mathbf{k}} y_{l} \prod_{d \in I} \widehat{Cov}_p(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} w_{k_d}^d w_{l_d}^d. \end{aligned} \qquad (5.5) $$

The last line in (5.5) follows from equations (5.1) and (5.3).

In the case of Horvitz-Thompson estimation, for any $k_d, l_d \in S_d$, we have:

$$ w_{k_d}^d = \frac{1}{\pi_{k_d}^d} \quad \text{and} \quad \widehat{Cov}^d(w_{k_d}^d, w_{l_d}^d) = \frac{\Delta_{k_d l_d}^d}{\pi_{k_d}^d \pi_{l_d}^d} \frac{1}{\pi_{k_d}^d \pi_{l_d}^d}. $$

For the bootstrap version $\hat{Y}_\pi^*$ of the HT-estimator, equations (5.4) and (5.5) lead to:

$$ \hat{Y}_{\pi}^{*} = \sum_{I \subseteq \{1, \dots, D\}} \hat{Y}_{\pi}^{I*} \qquad (5.6) $$

and

$$ V_*(\hat{Y}_\pi^{I*}) = \sum_{\mathbf{k}, l \in S} y_{\mathbf{k}} y_{l} \prod_{d \in I} \frac{\Delta_{k_d l_d}^d}{\pi_{k_d}^d \pi_{l_d}^d} \prod_{d \notin I} \frac{1}{\pi_{k_d}^d \pi_{l_d}^d} \qquad (5.7) $$

$$ = \sum_{\mathbf{k}', l' \in S_I} \frac{\hat{Y}_{\mathbf{k}'} \hat{Y}_{l'}}{\pi_{\mathbf{k}'} \pi_{l'}} \frac{\Delta_{\mathbf{k} l'}}{\pi_{\mathbf{k}'} \pi_{l'}}. \qquad (5.8) $$

This corresponds to the plug-in variance estimator for $\hat{Y}_\pi^I$ given in (3.24). Therefore, $V_*(\hat{Y}_\pi^*)$ equals the plug-in variance estimator $\hat{V}^{PLUG}(\hat{Y}_\pi)$ as defined in (3.26), and is also asymptotically
unbiased by Proposition 4. This is a very important result. It indicates that if we have bootstrap procedures that satisfy conditions (5.1) and (5.3) in each dimension, then under our conditions (H1)-(H4), it suffices to independently apply these bootstrap procedures in each dimension to obtain an asymptotically unbiased bootstrap variance estimator for CCS. In the following section, we illustrate this result for two bootstrap procedures in the case of simple random sampling in each dimension.

# 6 Case study: simple random sampling in each dimension

In this section, we focus on the particularly significant case of CCS in dimension $D = 2$, where simple random sampling is employed in each dimension. We investigate the application of the pseudo-population bootstrap in Section 6.1, and the utilization of the rescaled bootstrap is examined in Section 6.2.

## 6.1 Pseudo-Population Bootstrap method

We introduce a pseudo-population bootstrap approach tailored for simple random sampling without replacement in each dimension, drawing inspiration from the method proposed by Gross (1980). For ease of exposition, let us assume that both $N_2/n_2$ and $N_1/n_1$ are integers. Once $S_1$ and $S_2$ are sampled, the methodology involves constructing a pseudo-population $U_d^*$ by replicating each unit in $S_d$ for $N_d/n_d$ times for each $d = 1, 2$. Subsequently, simple random sampling without replacement is applied in each pseudo-population to generate the bootstrap sample $S_d^*$. Consequently, the bootstrap estimator is given by:

$$ \hat{Y}_{\pi}^* = \frac{N_1 N_2}{n_1 n_2} \sum_{k_1 \in S_1^*} \sum_{k_2 \in S_2^*} y_{k_1 k_2}. \qquad (6.1) $$

We can easily check that $E_*(\hat{Y}_{\pi}^*) = \hat{Y}_{\pi}$. Furthermore, the bootstrap variance is

$$ V_*(\hat{Y}_{\pi}^*) = N_1^2 N_2^2 \left[ (1-f_1) \frac{S_1^{2*}}{n_1} + (1-f_2) \frac{S_2^{2*}}{n_2} + (1-f_1)(1-f_2) \frac{S_{12}^{2*}}{n_1 n_2} \right] \qquad (6.2) $$

with

$$ \begin{aligned} S_1^{2*} &= \frac{N_1(n_1-1)}{n_1(N_1-1)} \left[ \hat{S}_1^2 + (1-f_2) \frac{\hat{S}_{12}^2}{n_2} \right], \\ S_2^{2*} &= \frac{N_2(n_2-1)}{n_2(N_2-1)} \left[ \hat{S}_2^2 + (1-f_1) \frac{\hat{S}_{12}^2}{n_1} \right], \\ S_{12}^{2*} &= \frac{N_1(n_1-1)}{n_1(N_1-1)} \frac{N_2(n_2-1)}{n_2(N_2-1)} \hat{S}_{12}^2. \end{aligned} \qquad (6.3) $$
unbiased by Proposition 4. This is a very important result. It indicates that if we have bootstrap procedures that satisfy conditions (5.1) and (5.3) in each dimension, then under our conditions (H1)-(H4), it suffices to independently apply these bootstrap procedures in each dimension to obtain an asymptotically unbiased bootstrap variance estimator for CCS. In the following section, we illustrate this result for two bootstrap procedures in the case of simple random sampling in each dimension.

# 6 Case study: simple random sampling in each dimension

In this section, we focus on the particularly significant case of CCS in dimension $D = 2$, where simple random sampling is employed in each dimension. We investigate the application of the pseudo-population bootstrap in Section 6.1, and the utilization of the rescaled bootstrap is examined in Section 6.2.

## 6.1 Pseudo-Population Bootstrap method

We introduce a pseudo-population bootstrap approach tailored for simple random sampling without replacement in each dimension, drawing inspiration from the method proposed by Gross (1980). For ease of exposition, let us assume that both $N_2/n_2$ and $N_1/n_1$ are integers. Once $S_1$ and $S_2$ are sampled, the methodology involves constructing a pseudo-population $U_d^*$ by replicating each unit in $S_d$ for $N_d/n_d$ times for each $d = 1, 2$. Subsequently, simple random sampling without replacement is applied in each pseudo-population to generate the bootstrap sample $S_d^*$. Consequently, the bootstrap estimator is given by:

$$ \hat{Y}_{\pi}^* = \frac{N_1 N_2}{n_1 n_2} \sum_{k_1 \in S_1^*} \sum_{k_2 \in S_2^*} y_{k_1 k_2}. \qquad (6.1) $$

We can easily check that $E_*(\hat{Y}_{\pi}^*) = \hat{Y}_{\pi}$. Furthermore, the bootstrap variance is

$$ V_*(\hat{Y}_{\pi}^*) = N_1^2 N_2^2 \left[ (1-f_1) \frac{S_1^{2*}}{n_1} + (1-f_2) \frac{S_2^{2*}}{n_2} + (1-f_1)(1-f_2) \frac{S_{12}^{2*}}{n_1 n_2} \right] \qquad (6.2) $$

with

$$ \begin{aligned} S_1^{2*} &= \frac{N_1(n_1-1)}{n_1(N_1-1)} \left[ \hat{S}_1^2 + (1-f_2) \frac{\hat{S}_{12}^2}{n_2} \right], \\ S_2^{2*} &= \frac{N_2(n_2-1)}{n_2(N_2-1)} \left[ \hat{S}_2^2 + (1-f_1) \frac{\hat{S}_{12}^2}{n_1} \right], \\ S_{12}^{2*} &= \frac{N_1(n_1-1)}{n_1(N_1-1)} \frac{N_2(n_2-1)}{n_2(N_2-1)} \hat{S}_{12}^2. \end{aligned} \qquad (6.3) $$
Under our asymptotic framework where $n_1, n_2 \to \infty$, the bootstrap variance asymptotically converges to the unbiased variance estimator of $V_p(\hat{Y}_\pi)$, as illustrated in equation (4.11). However, when $n_1 \to \infty$ and $n_2$ is bounded, the component $S_1^2$ tends to be overestimated, while the components $S_2^2$ and $S_{12}^2$ tend to be underestimated. When $n_2 \to \infty$ and $n_1$ is bounded, the component $S_2^2$ tends to be overestimated whereas the components $S_1^2$ and $S_{12}^2$ tend to be underestimated.

## 6.2 Rescaled bootstrap method

We propose a rescaled bootstrap method tailored for simple random sampling without replacement in each dimension, drawing inspiration from the approach proposed in Rao et al. (1992). The Horvitz-Thompson estimator employs weights of the form $w_{k_1 k_2} = w_{k_1}^{(1)} w_{k_2}^{(2)}$ for each unit $(k_1, k_2) \in U$, where $w_{k_1}^{(1)} = 1/f_1$ and $w_{k_2}^{(2)} = 1/f_2$. The method proposed in Rao et al. (1992) involves constructing bootstrap weights $w_{k_d}^{d*}$ for $d = 1, 2$ by multiplying the original weights by an adjustment factor $a_{k_d}^{d*}$, resulting in $w_{k_d}^{d*} = a_{k_d}^{d*} w_{k_d}^d$. These adjustment factors are determined by sampling multiplicities $(m_{k_d}^*)_{k_d \in S_d}$ from a multinomial distribution $M(n_d^*; \frac{1}{n_d}, \dots, \frac{1}{n_d})$, leading to:

$$a_{k_d}^{d*} = 1 + \sqrt{\frac{n_d^*(1-f_d)}{n_d-1}} \left( \frac{n_d m_{k_d}^{d*}}{n_d^*} - 1 \right) \quad \text{for any } k_d \in S_d, \quad (6.4)$$

where $n_d^*$ represents the resampling size, commonly chosen as $n_d = n_d - 1$.

The rescaled bootstrap for CCS involves constructing bootstrap weights $w_{k_1 k_2}^* = w_{k_1}^{(1*)} w_{k_2}^{(2*)}$ by employing the rescaled bootstrap in each dimension. Performing straightforward algebra yields:

$$V_*(\hat{Y}_\pi^*) = N_1^2 N_2^2 \left[ (1-f_1) \frac{S_1^{2*}}{n_1} + (1-f_2) \frac{S_2^{2*}}{n_2} + (1-f_1)(1-f_2) \frac{S_{12}^{2*}}{n_1 n_2} \right], \quad (6.5)$$

with

$$\begin{align}
S_1^{2*} &= \hat{S}_1^2 + (1 - f_2) \frac{\hat{S}_{12}^2}{n_2}, \nonumber \\
S_2^{2*} &= \hat{S}_2^2 + (1 - f_1) \frac{\hat{S}_{12}^2}{n_1}, \tag{6.6} \\
S_{12}^{2*} &= \hat{S}_{12}^2. \nonumber
\end{align}$$

While the second-order moment condition is met in each dimension, the method by Rao and Wu does not ensure that the bootstrap variance estimator is unbiased for CCS. However, it does ensure that under our asymptotic framework where $n_1, n_2 \to \infty$, the bootstrap variance estimator converges asymptotically to the unbiased variance estimator of $V_p(\hat{Y}_\pi)$. Additionally, if $n_1$ (respectively,
Under our asymptotic framework where $n_1, n_2 \to \infty$, the bootstrap variance asymptotically converges to the unbiased variance estimator of $V_p(\hat{Y}_\pi)$, as illustrated in equation (4.11). However, when $n_1 \to \infty$ and $n_2$ is bounded, the component $S_1^2$ tends to be overestimated, while the components $S_2^2$ and $S_{12}^2$ tend to be underestimated. When $n_2 \to \infty$ and $n_1$ is bounded, the component $S_2^2$ tends to be overestimated whereas the components $S_1^2$ and $S_{12}^2$ tend to be underestimated.

## 6.2 Rescaled bootstrap method

We propose a rescaled bootstrap method tailored for simple random sampling without replacement in each dimension, drawing inspiration from the approach proposed in Rao et al. (1992). The Horvitz-Thompson estimator employs weights of the form $w_{k_1k_2} = w_{k_1}^{(1)}w_{k_2}^{(2)}$ for each unit $(k_1, k_2) \in U$, where $w_{k_1}^{(1)} = 1/f_1$ and $w_{k_2}^{(2)} = 1/f_2$. The method proposed in Rao et al. (1992) involves constructing bootstrap weights $w_{k_d}^{d*}$ for $d = 1, 2$ by multiplying the original weights by an adjustment factor $a_{k_d}^{d*}$, resulting in $w_{k_d}^{d*} = a_{k_d}^{d*}w_{k_d}^d$. These adjustment factors are determined by sampling multiplicities $(m_{k_d}^*)_{k_d \in S_d}$ from a multinomial distribution $M(n_d^*; \frac{1}{n_d}, \dots, \frac{1}{n_d})$, leading to:

$$a_{k_d}^{d*} = 1 + \sqrt{\frac{n_d^*(1-f_d)}{n_d-1}} \left( \frac{n_d m_{k_d}^{d*}}{n_d^*} - 1 \right) \quad \text{for any } k_d \in S_d, \quad (6.4)$$

where $n_d^*$ represents the resampling size, commonly chosen as $n_d = n_d - 1$.

The rescaled bootstrap for CCS involves constructing bootstrap weights $w_{k_1k_2}^* = w_{k_1}^{(1*)}w_{k_2}^{(2*)}$ by employing the rescaled bootstrap in each dimension. Performing straightforward algebra yields:

$$V_*(\hat{Y}_\pi^*) = N_1^2 N_2^2 \left[ (1-f_1) \frac{S_1^{2*}}{n_1} + (1-f_2) \frac{S_2^{2*}}{n_2} + (1-f_1)(1-f_2) \frac{S_{12}^{2*}}{n_1 n_2} \right], \quad (6.5)$$

with

$$\begin{align}
S_1^{2*} &= \hat{S}_1^2 + (1 - f_2) \frac{\hat{S}_{12}^2}{n_2}, \nonumber \\
S_2^{2*} &= \hat{S}_2^2 + (1 - f_1) \frac{\hat{S}_{12}^2}{n_1}, \tag{6.6} \\
S_{12}^{2*} &= \hat{S}_{12}^2. \nonumber
\end{align}$$

While the second-order moment condition is met in each dimension, the method by Rao and Wu does not ensure that the bootstrap variance estimator is unbiased for CCS. However, it does ensure that under our asymptotic framework where $n_1, n_2 \to \infty$, the bootstrap variance estimator converges asymptotically to the unbiased variance estimator of $V_p(\hat{Y}_\pi)$. Additionally, if $n_1$ (respectively,
$n_2$) is bounded, then the component $S_2$ (respectively, $S_1$) tends to be overestimated, rendering our proposed adaptation of the Rao-Wu-Yue bootstrap a conservative approach.

# 7 Simulation study

In this section, we conduct a simulation study in the $D = 2$ dimensional case, inspired by the setup outlined in Juillard et al. (2017). Our code to realize these simulations is accessible on github in the following link: https://github.com/Jean-Rubin/CrossClassifiedSampling. Our objective is to assess the performance of the proposed variance estimators for estimating both a total and a ratio. The underlying generation model we consider is:

$$y_{k_1 k_2} = \mu + \sigma_1 U_{k_1} + \sigma_2 V_{k_2} + \sigma_{12} W_{k_1 k_2}, \quad (7.1)$$

where the $U_{k_1}$, $V_{k_2}$, $W_{k_1 k_2}$ are independent standard normal variables.

## 7.1 Variance estimation for a total

In this simulation study, we generate data according to model (7.1) with the following parameters: $N_1 = N_2 = 1000$, $\mu = 200$, and $\sigma_1 = \sigma_2 = 5$. The parameter $\sigma_{12}$ varies in {5, 10, 50} to examine the impact of interaction variance.

We perform $T = 10,000$ iterations of sample selection using CCS with simple random sampling in each dimension, with varying sample sizes $n_1$ and $n_2$ chosen from {5, 10, 100, 500}. For each iteration, we estimate the total $Y$, and the estimator for the $t^{th}$ sample is denoted as $\hat{Y}^{(t)}$. We assess the performance of the *Unbiased* variance estimator (see equation 4.11), the simplified variance estimator $\hat{V}^{SIMPL}(\hat{Y}_{\pi})$ (Simpl) given in equation (4.13), and the simplified variance estimator $\hat{V}^{SIMPL2}(\hat{Y}_{\pi})$ (Simpl2) given in equation (4.14). We also compute the bootstrap variance estimator $\hat{V}_{Gross}^{(t)}$ associated to the pseudo-population bootstrap (Gross), $\hat{V}_{RaoWu}^{(t)}$ associated to the rescaled bootstrap (RaoWu), and $\hat{V}_{Skinner}^{(t)}$ associated to a bootstrap procedure for with-replacement sampling in each dimension proposed by Skinner (2015) (Skinner). We use $B = 1,000$ bootstrap resamples for each estimator. The bootstrap variance estimator $\hat{V}^{(t)}$ is computed using the formula:

$$\hat{V}^{(t)} = \frac{1}{B-1} \sum_{b=1}^{B} \left( \hat{Y}^{(b*)}(t) - \frac{1}{B} \sum_{b=1}^{B} \hat{Y}^{(b*)}(t) \right)^2 . \quad (7.2)$$
$n_2$) is bounded, then the component $S_2$ (respectively, $S_1$) tends to be overestimated, rendering our proposed adaptation of the Rao-Wu-Yue bootstrap a conservative approach.

# 7 Simulation study

In this section, we conduct a simulation study in the $D = 2$ dimensional case, inspired by the setup outlined in Juillard et al. (2017). Our code to realize these simulations is accessible on github in the following link: https://github.com/Jean-Rubin/CrossClassifiedSampling. Our objective is to assess the performance of the proposed variance estimators for estimating both a total and a ratio. The underlying generation model we consider is:

$$y_{k_1 k_2} = \mu + \sigma_1 U_{k_1} + \sigma_2 V_{k_2} + \sigma_{12} W_{k_1 k_2}, \quad (7.1)$$

where the $U_{k_1}$, $V_{k_2}$, $W_{k_1 k_2}$ are independent standard normal variables.

## 7.1 Variance estimation for a total

In this simulation study, we generate data according to model (7.1) with the following parameters: $N_1 = N_2 = 1000$, $\mu = 200$, and $\sigma_1 = \sigma_2 = 5$. The parameter $\sigma_{12}$ varies in $\{5, 10, 50\}$ to examine the impact of interaction variance.

We perform $T = 10,000$ iterations of sample selection using CCS with simple random sampling in each dimension, with varying sample sizes $n_1$ and $n_2$ chosen from $\{5, 10, 100, 500\}$. For each iteration, we estimate the total $Y$, and the estimator for the $t^{th}$ sample is denoted as $\hat{Y}^{(t)}$. We assess the performance of the *Unbiased* variance estimator (see equation 4.11), the simplified variance estimator $\hat{V}^{SIMPL}(\hat{Y}_{\pi})$ (Simpl) given in equation (4.13), and the simplified variance estimator $\hat{V}^{SIMPL2}(\hat{Y}_{\pi})$ (Simpl2) given in equation (4.14). We also compute the bootstrap variance estimator $\hat{V}_{Gross}^{(t)}$ associated to the pseudo-population bootstrap (Gross), $\hat{V}_{RaoWu}^{(t)}$ associated to the rescaled bootstrap (RaoWu), and $\hat{V}_{Skinner}^{(t)}$ associated to a bootstrap procedure for with-replacement sampling in each dimension proposed by Skinner (2015) (Skinner). We use $B = 1,000$ bootstrap resamples for each estimator. The bootstrap variance estimator $\hat{V}^{(t)}$ is computed using the formula:

$$\hat{V}^{(t)} = \frac{1}{B-1} \sum_{b=1}^{B} \left( \hat{Y}^{(b*)}(t) - \frac{1}{B} \sum_{b=1}^{B} \hat{Y}^{(b*)}(t) \right)^2 . \quad (7.2)$$
To evaluate these estimators, we compute the Monte-Carlo relative bias (RB) and relative stability (RS) as follows:

$$RB(\hat{V}) = 100 \times \frac{\frac{1}{T} \sum_{t=1}^{T} \hat{V}^{(t)} - V}{V},$$

$$RS(\hat{V}) = 100 \times \frac{\sqrt{\frac{1}{T} \sum_{t=1}^{T} (\hat{V}^{(t)} - V)^2}}{V},$$

where $V$ is an approximation of the true variance of $\hat{Y}$, computed from an independent run of 100,000 simulations. Additionally, confidence intervals for the bootstrap methods are computed using the reverse percentile method. More precisely, assume that the bootstrap estimates are reordered as $\hat{Y}^{(1*)}(t) \le \hat{Y}^{(2*)}(t) \le \dots \le \hat{Y}^{(B*)}(t)$. Then the reverse percentile confidence interval is obtained by considering $(\hat{Y}^* - \hat{Y})$ as an approximation of the distribution of $\hat{Y} - Y$, which leads to:

$$CI_{\alpha}^{REV,(t)}(Y) = [2\hat{Y}^{(t)} - \hat{Y}^{(U*)}(t), 2\hat{Y}^{(t)} - \hat{Y}^{(L*)}(t)]. \quad (7.3)$$

For $CI_{\alpha}^{REV,(t)}(Y)$ to have the coverage level $\alpha$, the values of $L$ and $U$ are given by $L = \lfloor \frac{1-\alpha}{2}B \rfloor$ and $U = \lfloor \frac{1+\alpha}{2}B \rfloor$. This comprehensive evaluation allows us to assess the accuracy, stability, and coverage of the proposed variance estimators under various conditions.

The simulation results are presented in Table 1. As anticipated, the relative bias of the two proposed simplified variance estimators approaches zero as both sample sizes increase. However, the first simplified variance estimator exhibits negative bias when both sample sizes are small or when $\sigma_{12}$ is large, as the dropped variance component in these scenarios is non-negligible. Similarly, in these situations, the second simplified variance estimator displays positive bias, consistent with the heuristic reasoning outlined in Section 3.3. All three analytic variance estimators demonstrate similar levels of stability.

Turning to the bootstrap methods, Skinner's approach is conservative across all cases, as expected given its design for with-replacement sampling. As anticipated, the relative bias of both the Gross bootstrap and the Rao-Wu-Yue bootstrap variance estimators diminishes as both sample sizes increase. The bias of these two bootstrap variance estimators increases with $\sigma_{12}$, aligning with the formulation of the bootstrap variance outlined in equations (6.2) and (6.5). Overall, the bootstrap variance estimator exhibits smaller bias and greater stability under Gross's method, particularly for moderate sample sizes. While the coverage rates of the bootstrap confidence intervals are well
To evaluate these estimators, we compute the Monte-Carlo relative bias (RB) and relative stability (RS) as follows:

$$RB(\hat{V}) = 100 \times \frac{\frac{1}{T} \sum_{t=1}^{T} \hat{V}^{(t)} - V}{V},$$

$$RS(\hat{V}) = 100 \times \frac{\sqrt{\frac{1}{T} \sum_{t=1}^{T} (\hat{V}^{(t)} - V)^2}}{V},$$

where $V$ is an approximation of the true variance of $\hat{Y}$, computed from an independent run of 100,000 simulations. Additionally, confidence intervals for the bootstrap methods are computed using the reverse percentile method. More precisely, assume that the bootstrap estimates are reordered as $\hat{Y}^{(1*)}(t) \le \hat{Y}^{(2*)}(t) \le \dots \le \hat{Y}^{(B*)}(t)$. Then the reverse percentile confidence interval is obtained by considering $(\hat{Y}^* - \hat{Y})$ as an approximation of the distribution of $\hat{Y} - Y$, which leads to:

$$CI_{\alpha}^{REV,(t)}(Y) = [2\hat{Y}^{(t)} - \hat{Y}^{(U*)}(t), 2\hat{Y}^{(t)} - \hat{Y}^{(L*)}(t)]. \quad (7.3)$$

For $CI_{\alpha}^{REV,(t)}(Y)$ to have the coverage level $\alpha$, the values of $L$ and $U$ are given by $L = \lfloor \frac{1-\alpha}{2}B \rfloor$ and $U = \lceil \frac{1+\alpha}{2}B \rceil$. This comprehensive evaluation allows us to assess the accuracy, stability, and coverage of the proposed variance estimators under various conditions.

The simulation results are presented in Table 1. As anticipated, the relative bias of the two proposed simplified variance estimators approaches zero as both sample sizes increase. However, the first simplified variance estimator exhibits negative bias when both sample sizes are small or when $\sigma_{12}$ is large, as the dropped variance component in these scenarios is non-negligible. Similarly, in these situations, the second simplified variance estimator displays positive bias, consistent with the heuristic reasoning outlined in Section 3.3. All three analytic variance estimators demonstrate similar levels of stability.

Turning to the bootstrap methods, Skinner's approach is conservative across all cases, as expected given its design for with-replacement sampling. As anticipated, the relative bias of both the Gross bootstrap and the Rao-Wu-Yue bootstrap variance estimators diminishes as both sample sizes increase. The bias of these two bootstrap variance estimators increases with $\sigma_{12}$, aligning with the formulation of the bootstrap variance outlined in equations (6.2) and (6.5). Overall, the bootstrap variance estimator exhibits smaller bias and greater stability under Gross's method, particularly for moderate sample sizes. While the coverage rates of the bootstrap confidence intervals are well
<table><caption>Table 1: Relative biases, relative stabilities and coverage probabilities of multiple estimators of the variance of Ŷ$_{π}$ in percentage</caption><thead><tr><th rowspan="2">σ<sub>12</sub></th><th rowspan="2">n<sub>1</sub></th><th rowspan="2">n<sub>2</sub></th><th colspan="6">RB</th><th colspan="3">RS</th><th colspan="3">CI<sub>95%</sub></th></tr><tr><th>Unbiased</th><th>Simpl</th><th>Simpl2</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Unbiased</th><th>Simpl</th><th>Simpl2</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Gross</th><th>RaoWu</th><th>Skinner</th></tr></thead><tbody><tr><td rowspan="2">5</td><td>5</td><td>5</td><td>2.2</td><td>-8.6</td><td>9.4</td><td>-7.2</td><td>17.5</td><td>18.1</td><td>55.3</td><td>56.2</td><td>55.9</td><td>44.8</td><td>57.6</td><td>58.4</td><td>91.1</td><td>94.2</td><td>94.2</td></tr><tr><td>10</td><td>10</td><td>1.2</td><td>-4.3</td><td>4.7</td><td>-1.7</td><td>8.7</td><td>10.2</td><td>34.5</td><td>35.1</td><td>35.0</td><td>31.4</td><td>36.0</td><td>37.1</td><td>93.4</td><td>94.5</td><td>94.9</td></tr><tr><td rowspan="2">10</td><td>100</td><td>0.6</td><td>0.2</td><td>2.2</td><td>-6.6</td><td>1.8</td><td>5.6</td><td>47.8</td><td>47.6</td><td>47.5</td><td>44.2</td><td>48.0</td><td>49.1</td><td>91.6</td><td>93.0</td><td>93.0</td></tr><tr><td>100</td><td>100</td><td>-0.1</td><td>-0.6</td><td>0.3</td><td>-0.2</td><td>0.7</td><td>11.9</td><td>10.1</td><td>10.1</td><td>10.2</td><td>10.8</td><td>11.2</td><td>17.0</td><td>95.1</td><td>94.9</td><td>96.0</td></tr><tr><td rowspan="2">500</td><td>500</td><td>500</td><td>-0.9</td><td>-0.9</td><td>-0.8</td><td>-0.9</td><td>-0.8</td><td>98.5</td><td>3.2</td><td>3.2</td><td>3.2</td><td>5.4</td><td>5.4</td><td>99.1</td><td>95.1</td><td>94.8</td><td>99.5</td></tr><tr><td>10</td><td>5</td><td>-0.3</td><td>-27.8</td><td>27.3</td><td>17.7</td><td>55.9</td><td>56.3</td><td>63.5</td><td>71.0</td><td>69.1</td><td>53.5</td><td>86.2</td><td>85.7</td><td>94.7</td><td>96.9</td><td>96.8</td></tr><tr><td rowspan="2">10</td><td>10</td><td>10</td><td>-0.7</td><td>-16.5</td><td>14.6</td><td>15.7</td><td>31.1</td><td>31.9</td><td>38.5</td><td>42.1</td><td>40.6</td><td>37.9</td><td>49.8</td><td>50.3</td><td>95.6</td><td>96.9</td><td>97.1</td></tr><tr><td>100</td><td>100</td><td>-0.3</td><td>-4.4</td><td>3.0</td><td>-3.1</td><td>6.2</td><td>8.7</td><td>43.5</td><td>43.1</td><td>43.0</td><td>38.7</td><td>44.3</td><td>44.3</td><td>92.3</td><td>93.3</td><td>94.0</td></tr><tr><td rowspan="2">100</td><td>100</td><td>100</td><td>-0.4</td><td>-2.0</td><td>1.7</td><td>2.2</td><td>3.1</td><td>14.8</td><td>10.1</td><td>10.2</td><td>10.3</td><td>11.2</td><td>11.5</td><td>19.4</td><td>95.2</td><td>95.0</td><td>96.2</td></tr><tr><td>500</td><td>500</td><td>-0.1</td><td>-0.1</td><td>0.1</td><td>0.5</td><td>0.3</td><td>100.9</td><td>3.2</td><td>3.0</td><td>2.9</td><td>6.3</td><td>5.1</td><td>102.4</td><td>95.1</td><td>94.7</td><td>99.4</td></tr><tr><td rowspan="2">20</td><td>5</td><td>5</td><td>-0.3</td><td>-60.9</td><td>60.8</td><td>66.8</td><td>121.0</td><td>122.5</td><td>82.5</td><td>109.1</td><td>100.5</td><td>94.6</td><td>146.3</td><td>149.0</td><td>97.5</td><td>98.9</td><td>98.9</td></tr><tr><td>10</td><td>10</td><td>1.0</td><td>-45.9</td><td>44.5</td><td>67.1</td><td>89.8</td><td>92.2</td><td>49.3</td><td>68.1</td><td>65.4</td><td>80.9</td><td>102.6</td><td>105.4</td><td>98.2</td><td>99.0</td><td>98.8</td></tr><tr><td rowspan="2">10</td><td>100</td><td>0.7</td><td>-10.9</td><td>12.3</td><td>13.1</td><td>24.0</td><td>28.7</td><td>43.2</td><td>44.7</td><td>43.8</td><td>41.3</td><td>49.4</td><td>52.3</td><td>95.1</td><td>95.8</td><td>96.0</td></tr><tr><td>100</td><td>100</td><td>-0.1</td><td>-6.2</td><td>6.1</td><td>11.4</td><td>12.2</td><td>25.6</td><td>10.1</td><td>11.8</td><td>11.8</td><td>16.0</td><td>16.7</td><td>28.6</td><td>96.0</td><td>96.3</td><td>97.2</td></tr><tr><td rowspan="2">500</td><td>500</td><td>500</td><td>-0.5</td><td>-1.2</td><td>0.2</td><td>0.9</td><td>1.1</td><td>103.6</td><td>3.2</td><td>3.4</td><td>3.2</td><td>5.6</td><td>5.7</td><td>104.2</td><td>95.1</td><td>95.1</td><td>99.6</td></tr></tbody></table>
<table><caption>Table 1: Relative biases, relative stabilities and coverage probabilities of multiple estimators of the variance of Ŷ$_{π}$ in percentage</caption><thead><tr><th rowspan="2">σ<sub>12</sub></th><th rowspan="2">n<sub>1</sub></th><th rowspan="2">n<sub>2</sub></th><th colspan="8">RB</th><th colspan="4">RS</th><th colspan="4">CI<sub>95%</sub></th></tr><tr><th>Unbiased</th><th>Simpl</th><th>Simpl2</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Unbiased</th><th>Simpl</th><th>Simpl2</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Gross</th><th>RaoWu</th><th>Skinner</th></tr></thead><tbody><tr><td rowspan="2">5</td><td>5</td><td>5</td><td>2.2</td><td>-8.6</td><td>9.4</td><td>-7.2</td><td>17.5</td><td>18.1</td><td>55.3</td><td>56.2</td><td>55.9</td><td>44.8</td><td>57.6</td><td>58.4</td><td>91.1</td><td>94.2</td><td>94.2</td></tr><tr><td>10</td><td>10</td><td>1.2</td><td>-4.3</td><td>4.7</td><td>-1.7</td><td>8.7</td><td>10.2</td><td>34.5</td><td>35.1</td><td>35.0</td><td>31.4</td><td>36.0</td><td>37.1</td><td>93.4</td><td>94.5</td><td>94.9</td></tr><tr><td rowspan="2">10</td><td>100</td><td>0.6</td><td>0.2</td><td>2.2</td><td>-6.6</td><td>1.8</td><td>5.6</td><td>47.8</td><td>47.6</td><td>47.5</td><td>44.2</td><td>48.0</td><td>49.1</td><td>91.6</td><td>93.0</td><td>93.0</td></tr><tr><td>100</td><td>100</td><td>-0.1</td><td>-0.6</td><td>0.3</td><td>-0.2</td><td>0.7</td><td>11.9</td><td>10.1</td><td>10.1</td><td>10.2</td><td>10.8</td><td>11.2</td><td>17.0</td><td>95.1</td><td>94.9</td><td>96.0</td></tr><tr><td rowspan="2">500</td><td>500</td><td>500</td><td>-0.9</td><td>-0.9</td><td>-0.8</td><td>-0.9</td><td>-0.8</td><td>98.5</td><td>3.2</td><td>3.2</td><td>3.2</td><td>5.4</td><td>5.4</td><td>99.1</td><td>95.1</td><td>94.8</td><td>99.5</td></tr><tr><td>10</td><td>5</td><td>-0.3</td><td>-27.8</td><td>27.3</td><td>17.7</td><td>55.9</td><td>56.3</td><td>63.5</td><td>71.0</td><td>69.1</td><td>53.5</td><td>86.2</td><td>85.7</td><td>94.7</td><td>96.9</td><td>96.8</td></tr><tr><td rowspan="2">10</td><td>10</td><td>10</td><td>-0.7</td><td>-16.5</td><td>14.6</td><td>15.7</td><td>31.1</td><td>31.9</td><td>38.5</td><td>42.1</td><td>40.6</td><td>37.9</td><td>49.8</td><td>50.3</td><td>95.6</td><td>96.9</td><td>97.1</td></tr><tr><td>100</td><td>100</td><td>-0.3</td><td>-4.4</td><td>3.0</td><td>-3.1</td><td>6.2</td><td>8.7</td><td>43.5</td><td>43.1</td><td>43.0</td><td>38.7</td><td>44.3</td><td>44.3</td><td>92.3</td><td>93.3</td><td>94.0</td></tr><tr><td rowspan="2">100</td><td>100</td><td>100</td><td>-0.4</td><td>-2.0</td><td>1.7</td><td>2.2</td><td>3.1</td><td>14.8</td><td>10.1</td><td>10.2</td><td>10.3</td><td>11.2</td><td>11.5</td><td>19.4</td><td>95.2</td><td>95.0</td><td>96.2</td></tr><tr><td>500</td><td>500</td><td>-0.1</td><td>-0.1</td><td>0.1</td><td>0.5</td><td>0.3</td><td>100.9</td><td>3.2</td><td>3.0</td><td>2.9</td><td>6.3</td><td>5.1</td><td>102.4</td><td>95.1</td><td>94.7</td><td>99.4</td></tr><tr><td rowspan="2">20</td><td>5</td><td>5</td><td>-0.3</td><td>-60.9</td><td>60.8</td><td>66.8</td><td>121.0</td><td>122.5</td><td>82.5</td><td>109.1</td><td>100.5</td><td>94.6</td><td>146.3</td><td>149.0</td><td>97.5</td><td>98.9</td><td>98.9</td></tr><tr><td>10</td><td>10</td><td>1.0</td><td>-45.9</td><td>44.5</td><td>67.1</td><td>89.8</td><td>92.2</td><td>49.3</td><td>68.1</td><td>65.4</td><td>80.9</td><td>102.6</td><td>105.4</td><td>98.2</td><td>99.0</td><td>98.8</td></tr><tr><td rowspan="2">10</td><td>100</td><td>0.7</td><td>-10.9</td><td>12.3</td><td>13.1</td><td>24.0</td><td>28.7</td><td>43.2</td><td>44.7</td><td>43.8</td><td>41.3</td><td>49.4</td><td>52.3</td><td>95.1</td><td>95.8</td><td>96.0</td></tr><tr><td>100</td><td>100</td><td>-0.1</td><td>-6.2</td><td>6.1</td><td>11.4</td><td>12.2</td><td>25.6</td><td>10.1</td><td>11.8</td><td>11.8</td><td>16.0</td><td>16.7</td><td>28.6</td><td>96.0</td><td>96.3</td><td>97.2</td></tr><tr><td rowspan="2">500</td><td>500</td><td>500</td><td>-0.5</td><td>-1.2</td><td>0.2</td><td>0.9</td><td>1.1</td><td>103.6</td><td>3.2</td><td>3.4</td><td>3.2</td><td>5.6</td><td>5.7</td><td>104.2</td><td>95.1</td><td>95.1</td><td>99.6</td></tr></tbody></table>
maintained when the variance estimators are approximately unbiased, they tend to be conservative when the variance estimators are conservative, as one would expect.

## 7.2 Variance estimation for a ratio

In the second population, the values for three variables of interest are generated according to the model

$$z_{k_1 k_2} = \sigma_1 U_{k_1} + \sigma_2 V_{k_2} + \sigma_{12} W_{k_1 k_2}, \quad (7.4)$$

$$x_{k_1 k_2} = \mu_x + \alpha z_{k_1, k_2} + (1 - \alpha)(\sigma_1' U'_{k_1} + \sigma_2' V'_{k_2} + \sigma_{12}' W'_{k_1 k_2}) \quad (7.5)$$

$$y_{k_1 k_2} = \mu_y + \beta z_{k_1, k_2} + (1 - \beta)(\sigma_1'' U''_{k_1} + \sigma_2'' V''_{k_2} + \sigma_{12}'' W''_{k_1 k_2}). \quad (7.6)$$

The variables $U_{k_1}$, $U'_{k_1}$, $U''_{k_1}$, $V_{k_2}$, $V'_{k_2}$, $V''_{k_2}$, and $W_{k_1 k_2}$, $W'_{k_1 k_2}$, $W''_{k_1 k_2}$, are independent standard normal variables. The parameter $\alpha$ in equation (7.5) is used to control the correlation between variables $x$ and $z$. Similarly, the parameter $\beta$ in equation (7.5) is used to control the correlation between variables $y$ and $z$. Thus, high values of $\alpha$ and $\beta$ lead to high correlations between $x$ and $y$. We used $\mu_x = 100$, $\mu_y = 300$, and $\sigma_1 = \sigma_2 = \sigma_{12} = \sigma'_1 = \sigma''_1 = 5$, $\sigma'_2 = \sigma''_2 = 10$, $\sigma'_{12} = \sigma''_{12} = 15$. Finally we considered $\alpha \in \{0.1, 0.5, 0.8\}$ and $\beta = 0.5$.

A Monte-carlo procedure similar to the one in the previous section is then applied to the substitution estimator $\hat{R} = \hat{Y}/\hat{X}$ of the ratio $R = Y/X$. For each sample $t \in \{1, \dots, T\}$, we produce bootstrap estimations $\hat{R}^{(b*)}(t) = \hat{Y}^{(b*)}(t)/\hat{X}^{(b*)}(t)$ from the bootstrap versions $(\hat{Y}^{(b*)}(t), \hat{X}^{(b*)}(t))$ of $(\hat{Y}^{(t)}, \hat{X}^{(t)})$. We obtain a bootstrap variance estimator, and a confidence interval by means of the reverse percentile bootstrap. For the analytic variance estimators, linearized variance estimators are obtained by replacing the variable $y_{k_1 k_2}$ by the linearized variable $\nu_{k_1 k_2} = (y_{k_1 k_2} - \hat{R} x_{k_1 k_2})/\hat{X}$ in the three possible variance estimators.

The simulation results are presented in Table 2. The results are very similar to those obtained for total estimation. We have compiled in Table 3 the Monte-carlo mean values for the dispersions $S_{1,\eta}^2$, $S_{2,\eta}^2$ and $S_{12,\eta}^2$ associated with the variable $\eta_{k_1 k_2} = y_{k_1 k_2} - Rx_{k_1 k_2}$. Results show that $S_{12,\eta}^2$ is comparatively larger as $\alpha$ decreases, and as for the estimation of a total, the bias of the variance estimators is larger in this situation.
maintained when the variance estimators are approximately unbiased, they tend to be conservative when the variance estimators are conservative, as one would expect.

## 7.2 Variance estimation for a ratio

In the second population, the values for three variables of interest are generated according to the model

$$z_{k_1 k_2} = \sigma_1 U_{k_1} + \sigma_2 V_{k_2} + \sigma_{12} W_{k_1 k_2}, \qquad (7.4)$$

$$x_{k_1 k_2} = \mu_x + \alpha z_{k_1, k_2} + (1 - \alpha)(\sigma_1' U'_{k_1} + \sigma_2' V'_{k_2} + \sigma_{12}' W'_{k_1 k_2}) \qquad (7.5)$$

$$y_{k_1 k_2} = \mu_y + \beta z_{k_1, k_2} + (1 - \beta)(\sigma_1'' U''_{k_1} + \sigma_2'' V''_{k_2} + \sigma_{12}'' W''_{k_1 k_2}). \qquad (7.6)$$

The variables $U_{k_1}$, $U'_{k_1}$, $U''_{k_1}$, $V_{k_2}$, $V'_{k_2}$, $V''_{k_2}$, and $W_{k_1 k_2}$, $W'_{k_1 k_2}$, $W''_{k_1 k_2}$, are independent standard normal variables. The parameter $\alpha$ in equation (7.5) is used to control the correlation between variables $x$ and $z$. Similarly, the parameter $\beta$ in equation (7.5) is used to control the correlation between variables $y$ and $z$. Thus, high values of $\alpha$ and $\beta$ lead to high correlations between $x$ and $y$. We used $\mu_x = 100$, $\mu_y = 300$, and $\sigma_1 = \sigma_2 = \sigma_{12} = \sigma'_1 = \sigma''_1 = 5$, $\sigma'_2 = \sigma''_2 = 10$, $\sigma'_{12} = \sigma''_{12} = 15$. Finally we considered $\alpha \in \{0.1, 0.5, 0.8\}$ and $\beta = 0.5$.

A Monte-carlo procedure similar to the one in the previous section is then applied to the substitution estimator $\hat{R} = \hat{Y}/\hat{X}$ of the ratio $R = Y/X$. For each sample $t \in \{1, \dots, T\}$, we produce bootstrap estimations $\hat{R}^{(b*)}(t) = \hat{Y}^{(b*)}(t)/\hat{X}^{(b*)}(t)$ from the bootstrap versions $(\hat{Y}^{(b*)}(t), \hat{X}^{(b*)}(t))$ of $(\hat{Y}^{(t)}, \hat{X}^{(t)})$. We obtain a bootstrap variance estimator, and a confidence interval by means of the reverse percentile bootstrap. For the analytic variance estimators, linearized variance estimators are obtained by replacing the variable $y_{k_1 k_2}$ by the linearized variable $\nu_{k_1 k_2} = (y_{k_1 k_2} - \hat{R} x_{k_1 k_2})/\hat{X}$ in the three possible variance estimators.

The simulation results are presented in Table 2. The results are very similar to those obtained for total estimation. We have compiled in Table 3 the Monte-carlo mean values for the dispersions $S_{1,\eta}^2$, $S_{2,\eta}^2$ and $S_{12,\eta}^2$ associated with the variable $\eta_{k_1 k_2} = y_{k_1 k_2} - Rx_{k_1 k_2}$. Results show that $S_{12,\eta}^2$ is comparatively larger as $\alpha$ decreases, and as for the estimation of a total, the bias of the variance estimators is larger in this situation.
<table><caption>Table 2: Relative biases, relative stabilities and coverage probabilities of multiple estimators of the variance of $\hat{R}$ in percentage</caption><thead><tr><th rowspan="2">α</th><th rowspan="2">n<sub>1</sub></th><th rowspan="2">n<sub>2</sub></th><th colspan="6">RB</th><th colspan="6">RS</th><th colspan="6">CI<sub>95%</sub></th></tr><tr><th>Unbiased</th><th>Simpl</th><th>Simpl2</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Unbiased</th><th>Simpl</th><th>Simpl2</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Gross</th><th>RaoWu</th><th>Skinner</th></tr></thead><tbody><tr><td rowspan="6">0.1</td><td>5</td><td>5</td><td>-0.1</td><td>-25.6</td><td>26.1</td><td>22.4</td><td>60.6</td><td>61.1</td><td>70.7</td><td>76.6</td><td>77.5</td><td>67.7</td><td>60.6</td><td>61.1</td><td>93.8</td><td>96.8</td><td>96.3</td><td>93.8</td><td>96.8</td><td>96.3</td></tr><tr><td>10</td><td>10</td><td>0.2</td><td>-14.0</td><td>14.2</td><td>16.6</td><td>30.6</td><td>31.6</td><td>45.6</td><td>46.8</td><td>48.7</td><td>47.1</td><td>57.9</td><td>58.7</td><td>94.9</td><td>95.8</td><td>95.8</td><td>94.9</td><td>95.8</td><td>95.8</td></tr><tr><td>10</td><td>100</td><td>0.4</td><td>-6.1</td><td>5.9</td><td>3.4</td><td>12.1</td><td>16.4</td><td>35.4</td><td>35.5</td><td>35.8</td><td>33.8</td><td>39.0</td><td>40.8</td><td>94.1</td><td>94.8</td><td>95.2</td><td>94.1</td><td>94.8</td><td>95.2</td></tr><tr><td>100</td><td>100</td><td>0.3</td><td>-1.1</td><td>2.2</td><td>2.3</td><td>3.3</td><td>15.0</td><td>11.6</td><td>11.7</td><td>11.7</td><td>11.7</td><td>12.5</td><td>13.1</td><td>20.5</td><td>95.1</td><td>95.1</td><td>95.1</td><td>95.1</td><td>95.2</td></tr><tr><td>500</td><td>500</td><td>-0.6</td><td>-0.7</td><td>-0.4</td><td>0.4</td><td>0.4</td><td>101.5</td><td>3.9</td><td>4.0</td><td>3.9</td><td>6.0</td><td>6.0</td><td>102.2</td><td>94.9</td><td>95.2</td><td>95.2</td><td>94.9</td><td>95.2</td><td>99.4</td></tr><tr><td>0.5</td><td>5</td><td>5</td><td>0.4</td><td>-24.2</td><td>24.5</td><td>17.3</td><td>50.8</td><td>51.8</td><td>66.0</td><td>70.8</td><td>70.2</td><td>58.4</td><td>85.2</td><td>87.6</td><td>93.9</td><td>96.3</td><td>96.3</td><td>96.4</td><td>96.4</td><td>96.4</td></tr><tr><td rowspan="6">0.8</td><td>10</td><td>10</td><td>-0.1</td><td>-14.6</td><td>13.7</td><td>15.5</td><td>29.1</td><td>31.2</td><td>42.1</td><td>44.2</td><td>43.9</td><td>42.3</td><td>52.5</td><td>54.1</td><td>95.3</td><td>96.0</td><td>96.0</td><td>96.1</td><td>96.1</td><td>96.1</td><td>96.1</td></tr><tr><td>10</td><td>100</td><td>-0.9</td><td>-5.1</td><td>3.9</td><td>-0.1</td><td>8.2</td><td>12.5</td><td>38.7</td><td>38.9</td><td>39.0</td><td>35.5</td><td>39.4</td><td>41.4</td><td>93.5</td><td>94.6</td><td>94.6</td><td>94.8</td><td>94.7</td><td>96.3</td><td>96.3</td></tr><tr><td>100</td><td>100</td><td>0.1</td><td>-1.2</td><td>1.6</td><td>2.1</td><td>3.3</td><td>14.8</td><td>11.4</td><td>11.5</td><td>11.6</td><td>12.5</td><td>12.9</td><td>20.3</td><td>94.8</td><td>95.4</td><td>95.4</td><td>95.2</td><td>95.2</td><td>99.4</td><td>99.4</td></tr><tr><td>500</td><td>500</td><td>0.1</td><td>0.1</td><td>0.3</td><td>1.3</td><td>1.5</td><td>103.3</td><td>3.9</td><td>4.0</td><td>4.0</td><td>6.1</td><td>6.2</td><td>104.0</td><td>95.4</td><td>95.4</td><td>95.2</td><td>95.2</td><td>95.2</td><td>99.4</td><td>99.4</td></tr><tr><td>5</td><td>5</td><td>-0.7</td><td>-14.6</td><td>14.1</td><td>1.5</td><td>30.5</td><td>31.8</td><td>60.6</td><td>62.4</td><td>62.5</td><td>49.2</td><td>69.4</td><td>70.0</td><td>92.1</td><td>94.9</td><td>94.9</td><td>94.6</td><td>94.6</td><td>94.6</td><td>94.6</td></tr><tr><td>10</td><td>10</td><td>1.3</td><td>-7.6</td><td>7.8</td><td>4.0</td><td>16.1</td><td>17.3</td><td>37.2</td><td>37.2</td><td>37.0</td><td>34.1</td><td>40.8</td><td>41.5</td><td>93.9</td><td>95.3</td><td>95.3</td><td>95.1</td><td>95.1</td><td>95.1</td><td>95.1</td></tr><tr><td rowspan="6">100</td><td>10</td><td>100</td><td>-0.5</td><td>-1.7</td><td>1.4</td><td>-6.1</td><td>3.4</td><td>6.5</td><td>39.2</td><td>39.6</td><td>39.2</td><td>36.1</td><td>39.7</td><td>40.6</td><td>92.0</td><td>93.4</td><td>93.4</td><td>93.5</td><td>93.5</td><td>93.5</td><td>93.5</td></tr><tr><td>100</td><td>100</td><td>-0.4</td><td>-0.1</td><td>0.4</td><td>1.0</td><td>1.9</td><td>13.2</td><td>10.5</td><td>10.5</td><td>10.5</td><td>11.3</td><td>11.6</td><td>18.4</td><td>95.1</td><td>94.7</td><td>94.7</td><td>94.7</td><td>96.1</td><td>96.1</td><td>96.1</td></tr><tr><td>500</td><td>500</td><td>0.6</td><td>0.6</td><td>0.8</td><td>0.2</td><td>0.3</td><td>100.8</td><td>3.2</td><td>3.3</td><td>3.3</td><td>5.5</td><td>5.5</td><td>101.4</td><td>94.7</td><td>94.9</td><td>94.9</td><td>94.9</td><td>99.3</td><td>99.3</td><td>99.3</td></tr></tbody></table>
<table><caption>Table 2: Relative biases, relative stabilities and coverage probabilities of multiple estimators of the variance of $\hat{R}$ in percentage</caption><thead><tr><th rowspan="2">α</th><th rowspan="2">n<sub>1</sub></th><th rowspan="2">n<sub>2</sub></th><th colspan="6">RB</th><th colspan="6">RS</th><th colspan="6">CI<sub>95%</sub></th></tr><tr><th>Unbiased</th><th>Simpl</th><th>Simpl2</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Unbiased</th><th>Simpl</th><th>Simpl2</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Gross</th><th>RaoWu</th><th>Skinner</th><th>Gross</th><th>RaoWu</th><th>Skinner</th></tr></thead><tbody><tr><td rowspan="6">0.1</td><td>5</td><td>5</td><td>-0.1</td><td>-25.6</td><td>26.1</td><td>22.4</td><td>60.6</td><td>61.1</td><td>70.7</td><td>76.6</td><td>77.5</td><td>67.7</td><td>60.6</td><td>61.1</td><td>93.8</td><td>96.8</td><td>96.3</td><td>93.8</td><td>96.8</td><td>96.3</td></tr><tr><td>10</td><td>10</td><td>0.2</td><td>-14.0</td><td>14.2</td><td>16.6</td><td>30.6</td><td>31.6</td><td>45.6</td><td>46.8</td><td>48.7</td><td>47.1</td><td>57.9</td><td>58.7</td><td>94.9</td><td>95.8</td><td>95.8</td><td>94.9</td><td>95.8</td><td>95.8</td></tr><tr><td>10</td><td>100</td><td>0.4</td><td>-6.1</td><td>5.9</td><td>3.4</td><td>12.1</td><td>16.4</td><td>35.4</td><td>35.5</td><td>35.8</td><td>33.8</td><td>39.0</td><td>40.8</td><td>94.1</td><td>94.8</td><td>95.2</td><td>94.1</td><td>94.8</td><td>95.2</td></tr><tr><td>100</td><td>100</td><td>0.3</td><td>-1.1</td><td>2.2</td><td>2.3</td><td>3.3</td><td>15.0</td><td>11.6</td><td>11.7</td><td>11.7</td><td>11.7</td><td>12.5</td><td>13.1</td><td>20.5</td><td>95.1</td><td>95.1</td><td>95.1</td><td>95.1</td><td>95.2</td></tr><tr><td>500</td><td>500</td><td>-0.6</td><td>-0.7</td><td>-0.4</td><td>0.4</td><td>0.4</td><td>101.5</td><td>3.9</td><td>4.0</td><td>3.9</td><td>6.0</td><td>6.0</td><td>102.2</td><td>94.9</td><td>95.2</td><td>95.2</td><td>94.9</td><td>95.2</td><td>99.4</td></tr><tr><td>0.5</td><td>5</td><td>5</td><td>0.4</td><td>-24.2</td><td>24.5</td><td>17.3</td><td>50.8</td><td>51.8</td><td>66.0</td><td>70.8</td><td>70.2</td><td>58.4</td><td>85.2</td><td>87.6</td><td>93.9</td><td>96.3</td><td>96.3</td><td>96.4</td><td>96.4</td><td>96.4</td></tr><tr><td rowspan="6">0.8</td><td>10</td><td>10</td><td>-0.1</td><td>-14.6</td><td>13.7</td><td>15.5</td><td>29.1</td><td>31.2</td><td>42.1</td><td>44.2</td><td>43.9</td><td>42.3</td><td>52.5</td><td>54.1</td><td>95.3</td><td>96.0</td><td>96.0</td><td>96.1</td><td>96.1</td><td>96.1</td><td>96.1</td></tr><tr><td>10</td><td>100</td><td>-0.9</td><td>-5.1</td><td>3.9</td><td>-0.1</td><td>8.2</td><td>12.5</td><td>38.7</td><td>38.9</td><td>39.0</td><td>35.5</td><td>39.4</td><td>41.4</td><td>93.5</td><td>94.6</td><td>94.6</td><td>94.8</td><td>94.7</td><td>96.3</td><td>96.3</td></tr><tr><td>100</td><td>100</td><td>0.1</td><td>-1.2</td><td>1.6</td><td>2.1</td><td>3.3</td><td>14.8</td><td>11.4</td><td>11.5</td><td>11.6</td><td>12.5</td><td>12.9</td><td>20.3</td><td>94.8</td><td>95.4</td><td>95.4</td><td>95.2</td><td>95.2</td><td>99.4</td><td>99.4</td></tr><tr><td>500</td><td>500</td><td>0.1</td><td>0.1</td><td>0.3</td><td>1.3</td><td>1.5</td><td>103.3</td><td>3.9</td><td>4.0</td><td>4.0</td><td>6.1</td><td>6.2</td><td>104.0</td><td>95.4</td><td>95.4</td><td>95.2</td><td>95.2</td><td>95.2</td><td>99.4</td><td>99.4</td></tr><tr><td>5</td><td>5</td><td>-0.7</td><td>-14.6</td><td>14.1</td><td>1.5</td><td>30.5</td><td>31.8</td><td>60.6</td><td>62.4</td><td>62.5</td><td>49.2</td><td>69.4</td><td>70.0</td><td>92.1</td><td>94.9</td><td>94.9</td><td>94.6</td><td>94.6</td><td>94.6</td><td>94.6</td></tr><tr><td>10</td><td>10</td><td>1.3</td><td>-7.6</td><td>7.8</td><td>4.0</td><td>16.1</td><td>17.3</td><td>37.2</td><td>37.2</td><td>37.0</td><td>34.1</td><td>40.8</td><td>41.5</td><td>93.9</td><td>95.3</td><td>95.3</td><td>95.1</td><td>95.1</td><td>95.1</td><td>95.1</td></tr><tr><td rowspan="6">100</td><td>10</td><td>100</td><td>-0.5</td><td>-1.7</td><td>1.4</td><td>-6.1</td><td>3.4</td><td>6.5</td><td>39.2</td><td>39.6</td><td>39.2</td><td>36.1</td><td>39.7</td><td>40.6</td><td>92.0</td><td>93.4</td><td>93.4</td><td>93.5</td><td>93.5</td><td>93.5</td><td>93.5</td></tr><tr><td>100</td><td>100</td><td>-0.4</td><td>-0.1</td><td>0.4</td><td>1.0</td><td>1.9</td><td>13.2</td><td>10.5</td><td>10.5</td><td>10.5</td><td>11.3</td><td>11.6</td><td>18.4</td><td>95.1</td><td>94.7</td><td>94.7</td><td>94.7</td><td>96.1</td><td>96.1</td><td>96.1</td></tr><tr><td>500</td><td>500</td><td>0.6</td><td>0.6</td><td>0.8</td><td>0.2</td><td>0.3</td><td>100.8</td><td>3.2</td><td>3.3</td><td>3.3</td><td>5.5</td><td>5.5</td><td>101.4</td><td>94.7</td><td>94.9</td><td>94.9</td><td>94.9</td><td>99.3</td><td>99.3</td><td>99.3</td></tr></tbody></table>
<table><caption>Table 3: Monte Carlo mean value for the population dispersions $S_{1,\eta}^2$, $S_{2,\eta}^2$ and $S_{12,\eta}^2$</caption><thead><tr><th>α</th><th>S<sup>2</sup><sub>1,η</sub></th><th>S<sup>2</sup><sub>2,η</sub></th><th>S<sup>2</sup><sub>12,η</sub></th></tr></thead><tbody><tr><td>0.1</td><td>192</td><td>760</td><td>1,696</td></tr><tr><td>0.5</td><td>88</td><td>276</td><td>586</td></tr><tr><td>0.8</td><td>105</td><td>151</td><td>227</td></tr></tbody></table>

# 8 Application to the ELFE survey

ELFE is a French longitudinal survey that tracks children from birth to adulthood (Pirus et al., 2010; Charles et al., 2020). It commenced in France in 2011, with a survey of babies born in maternity wards, excluding the most premature infants. The sample was selected using two-dimensional cross-classified sampling, achieved by combining the population of 544 maternity wards with the population of 365 days. The population of maternity wards was divided into 5 strata of equal size based on the number of births in 2008. Similarly, the population of days was divided into 4 strata, each representing a season. The sampling design can be viewed as stratified simple random sampling in each dimension (Juillard et al., 2017). All babies born in the selected maternity wards on the designated days were surveyed, resulting in approximately 18,000 children after adjusting for unit non-response due to parental refusal to participate.

In this application, we are interested in the sub-population of maternity wards that hosted between 1,010 and 1,418 births in 2018. It represents our population $U_1$ in the first dimension, of size $N_1 = 109$. The responding maternity wards in $U_1$ are seen as coming from a simple random sample, corresponding to a sample $S_1$ of size $n_1 = 62$. The population $U_2$ in the second dimension is that of all days, of size $N_2 = 365$. The sample of days selected for ELFE is seen as coming from a simple random sample, corresponding to a sample $S_2$ of size $n_2 = 25$. Thus, we obtain a sample $S = S_1 \times S_2$ of 1,550 maternity wards \times days.

We are interested in estimating the total number of births, the total number and percentage of twin births, the total number and percentage of births according to the marital status of the mother (married or remarried; divorced, single, or widowed; civil solidarity pact), and the average birth weight of the babies. We are also interested in variables related to the mother: we estimate the
<table><caption>Table 3: Monte Carlo mean value for the population dispersions $S_{1,\eta}^2$, $S_{2,\eta}^2$ and $S_{12,\eta}^2$</caption><thead><tr><th>α</th><th>S<sup>2</sup><sub>1,η</sub></th><th>S<sup>2</sup><sub>2,η</sub></th><th>S<sup>2</sup><sub>12,η</sub></th></tr></thead><tbody><tr><td>0.1</td><td>192</td><td>760</td><td>1,696</td></tr><tr><td>0.5</td><td>88</td><td>276</td><td>586</td></tr><tr><td>0.8</td><td>105</td><td>151</td><td>227</td></tr></tbody></table>

# 8 Application to the ELFE survey

ELFE is a French longitudinal survey that tracks children from birth to adulthood (Pirus et al., 2010; Charles et al., 2020). It commenced in France in 2011, with a survey of babies born in maternity wards, excluding the most premature infants. The sample was selected using two-dimensional cross-classified sampling, achieved by combining the population of 544 maternity wards with the population of 365 days. The population of maternity wards was divided into 5 strata of equal size based on the number of births in 2008. Similarly, the population of days was divided into 4 strata, each representing a season. The sampling design can be viewed as stratified simple random sampling in each dimension (Juillard et al., 2017). All babies born in the selected maternity wards on the designated days were surveyed, resulting in approximately 18,000 children after adjusting for unit non-response due to parental refusal to participate.

In this application, we are interested in the sub-population of maternity wards that hosted between 1,010 and 1,418 births in 2018. It represents our population $U_1$ in the first dimension, of size $N_1 = 109$. The responding maternity wards in $U_1$ are seen as coming from a simple random sample, corresponding to a sample $S_1$ of size $n_1 = 62$. The population $U_2$ in the second dimension is that of all days, of size $N_2 = 365$. The sample of days selected for ELFE is seen as coming from a simple random sample, corresponding to a sample $S_2$ of size $n_2 = 25$. Thus, we obtain a sample $S = S_1 \times S_2$ of 1,550 maternity wards \times days.

We are interested in estimating the total number of births, the total number and percentage of twin births, the total number and percentage of births according to the marital status of the mother (married or remarried; divorced, single, or widowed; civil solidarity pact), and the average birth weight of the babies. We are also interested in variables related to the mother: we estimate the
average weight gain between the beginning and the end of pregnancy, and the average BMI at the end of pregnancy. For a total $Y$, we use the Horvitz-Thompson estimator $\hat{Y}_{\pi}$ given in equation (2.8). For a ratio $R = Y^1/Y^2$, we use the substitution estimator $\hat{R}_{\pi} = \hat{Y}_{\pi}^1/\hat{Y}_{\pi}^2$.

We first consider analytic variance estimators. For the estimator of the total $\hat{Y}_{\pi}$, we consider the *Unbiased* variance estimator (see equation 4.11), the simplified variance estimator $\hat{V}^{SIMPL}(\hat{Y}_{\pi})$ (*Simpl*) given in equation (4.13), and the simplified variance estimator $\hat{V}^{SIMPL2}(\hat{Y}_{\pi})$ (*Simpl2*) given in equation (4.14). For the substitution estimator $\hat{R}_{\pi} = \hat{Y}_{\pi}^1/\hat{Y}_{\pi}^2$, we use linearization variance estimators obtained by replacing in the three previous estimators the variable of interest $y_{k_1k_2}$ with the linearized variable $e_{k_1k_2} = \{\hat{Y}_{\pi}^2\}^{-1}\{y_{k_1k_2}^2 - \hat{R}_{\pi}y_{k_1k_2}^1\}$. We also consider two bootstrap variance estimators for both $\hat{Y}_{\pi}$ and $\hat{R}_{\pi}$. The first procedure is the pseudo-population bootstrap (see Section 6.1), with the variant proposed by Booth et al. (1994) to take into account the fact that $N_1/n_1$ and $N_2/n_2$ are not integers. The second procedure is the rescaled bootstrap (see Section 6.2). The results are presented in terms of coefficients of variation. The analytic and the bootstrap variance estimators are also plugged into normality-based confidence intervals, with a nominal one-tailed error rate of 2.5 %. For both bootstrap procedures, we also use a confidence interval based on the reverse percentile method (see equation 7.3).

The findings are presented in Table 4 concerning the estimation of counts for variables. It is observed that, in general, the various methods for estimating variance yield similar results, with the exception being the estimation of twin counts. Once again, this scenario demonstrates a notable disparity in the interaction term, with $\hat{S}_{12}^2 = 4.8 \ 10^{-2}$, $\hat{S}_1^2 = 4.4 \ 10^{-4}$, and $\hat{S}_2^2 = 8.2 \ 10^{-4}$. Notably, the variance estimator *Simpl* exhibits a negative bias in this case. Conversely, the variance estimator *Simpl2* displays a positive bias, albeit less pronounced than the bootstrap variance estimators. These observations extend to the estimation of ratios, as evidenced in Table 5.

# 9 Conclusion

In this paper, we investigated CSS designs across an arbitrary number of dimensions. We established the weak consistency and asymptotic normality of the HT-estimator under mild conditions. Specifically, we highlighted that these asymptotic properties necessitate the sample size to approach infinity in each dimension, presenting a significant limitation of CSS designs compared to multi-
average weight gain between the beginning and the end of pregnancy, and the average BMI at the end of pregnancy. For a total $Y$, we use the Horvitz-Thompson estimator $\hat{Y}_{\pi}$ given in equation (2.8). For a ratio $R = Y^1/Y^2$, we use the substitution estimator $\hat{R}_{\pi} = \hat{Y}_{\pi}^1/\hat{Y}_{\pi}^2$.

We first consider analytic variance estimators. For the estimator of the total $\hat{Y}_{\pi}$, we consider the *Unbiased* variance estimator (see equation 4.11), the simplified variance estimator $\hat{V}^{SIMPL}(\hat{Y}_{\pi})$ (*Simpl*) given in equation (4.13), and the simplified variance estimator $\hat{V}^{SIMPL2}(\hat{Y}_{\pi})$ (*Simpl2*) given in equation (4.14). For the substitution estimator $\hat{R}_{\pi} = \hat{Y}_{\pi}^1/\hat{Y}_{\pi}^2$, we use linearization variance estimators obtained by replacing in the three previous estimators the variable of interest $y_{k_1k_2}$ with the linearized variable $e_{k_1k_2} = \{\hat{Y}_{\pi}^2\}^{-1}\{y_{k_1k_2}^2 - \hat{R}_{\pi}y_{k_1k_2}^1\}$. We also consider two bootstrap variance estimators for both $\hat{Y}_{\pi}$ and $\hat{R}_{\pi}$. The first procedure is the pseudo-population bootstrap (see Section 6.1), with the variant proposed by Booth et al. (1994) to take into account the fact that $N_1/n_1$ and $N_2/n_2$ are not integers. The second procedure is the rescaled bootstrap (see Section 6.2). The results are presented in terms of coefficients of variation. The analytic and the bootstrap variance estimators are also plugged into normality-based confidence intervals, with a nominal one-tailed error rate of 2.5 %. For both bootstrap procedures, we also use a confidence interval based on the reverse percentile method (see equation 7.3).

The findings are presented in Table 4 concerning the estimation of counts for variables. It is observed that, in general, the various methods for estimating variance yield similar results, with the exception being the estimation of twin counts. Once again, this scenario demonstrates a notable disparity in the interaction term, with $\hat{S}_{12}^2 = 4.8 \ 10^{-2}$, $\hat{S}_1^2 = 4.4 \ 10^{-4}$, and $\hat{S}_2^2 = 8.2 \ 10^{-4}$. Notably, the variance estimator *Simpl* exhibits a negative bias in this case. Conversely, the variance estimator *Simpl2* displays a positive bias, albeit less pronounced than the bootstrap variance estimators. These observations extend to the estimation of ratios, as evidenced in Table 5.

# 9 Conclusion

In this paper, we investigated CSS designs across an arbitrary number of dimensions. We established the weak consistency and asymptotic normality of the HT-estimator under mild conditions. Specifically, we highlighted that these asymptotic properties necessitate the sample size to approach infinity in each dimension, presenting a significant limitation of CSS designs compared to multi-
<table><caption>Table 4: Estimations of the numbers for count variables, with estimated coefficient of variations and confidence intervals</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th rowspan="2">Birth</th><th rowspan="2">Twins</th><th colspan="3">Marital status</th></tr><tr><th>Married/remarried</th><th>Divorced/single<br>/widowed</th><th>Civil solidarity<br>pact</th></tr></thead><tbody><tr><td></td><td>ŷ<sub>π</sub></td><td>49 487</td><td>975</td><td>21 766</td><td>513</td><td>27 20</td></tr><tr><td rowspan="2">Unbiased</td><td>CV̂</td><td>10.5%</td><td>27.7%</td><td>11.6%</td><td>10.6%</td><td>26.0%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 291, 59 683]</td><td>[446, 1 504]</td><td>[16 807, 26 725]</td><td>[21 559, 32 857]</td><td>[252, 774]</td></tr><tr><td rowspan="2">Simpl</td><td>CV̂</td><td>10.4%</td><td>23.7%</td><td>11.4%</td><td>10.4%</td><td>21.8%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 366, 59 608]</td><td>[523, 1 427]</td><td>[16 888, 26 644]</td><td>[21 641, 32 775]</td><td>[293, 733]</td></tr><tr><td rowspan="2">Simpl2</td><td>CV̂</td><td>10.6%</td><td>31.2%</td><td>11.8%</td><td>10.7%</td><td>29.5%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 216, 59 758]</td><td>[379, 1 571]</td><td>[16 727, 26 805]</td><td>[21 477, 32 939]</td><td>[216, 810]</td></tr><tr><td rowspan="3">Gross</td><td>CV̂</td><td>10.6%</td><td>35.2%</td><td>11.9%</td><td>10.9%</td><td>32.8%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 203, 59 771]</td><td>[302, 1 648]</td><td>[16 694, 26 838]</td><td>[21 382, 33 034]</td><td>[183, 843]</td></tr><tr><td>IC<sub>rev</sub></td><td>[39 580, 60 166]</td><td>[205, 1 540]</td><td>[16 710, 26 874]</td><td>[21 510, 33 086]</td><td>[154, 821]</td></tr><tr><td rowspan="3">Rao-Wu</td><td>CV̂</td><td>10.6%</td><td>34.1%</td><td>11.9%</td><td>10.8%</td><td>33.2%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 248, 59 726]</td><td>[324, 1 626]</td><td>[16 701, 26 831]</td><td>[21 447, 32 969]</td><td>[179, 847]</td></tr><tr><td>IC<sub>rev</sub></td><td>[39 635, 59 822]</td><td>[249, 1 525]</td><td>[16 816, 26 841]</td><td>[21 594, 33 009]</td><td>[140, 801]</td></tr></tbody></table>
<table><caption>Table 4: Estimations of the numbers for count variables, with estimated coefficient of variations and confidence intervals</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th rowspan="2">Birth</th><th rowspan="2">Twins</th><th colspan="3">Marital status</th></tr><tr><th>Married/remarried</th><th>Divorced/single</th><th>Civil solidarity</th></tr><tr><th></th><th></th><th></th><th></th><th colspan="3">/widowed</th><th>pact</th></tr></thead><tbody><tr><td></td><td>ŷ<sub>π</sub></td><td>49 487</td><td>975</td><td>21 766</td><td>513</td><td>27 20</td></tr><tr><td rowspan="2">Unbiased</td><td>CV̂</td><td>10.5%</td><td>27.7%</td><td>11.6%</td><td>10.6%</td><td>26.0%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 291, 59 683]</td><td>[446, 1 504]</td><td>[16 807, 26 725]</td><td>[21 559, 32 857]</td><td>[252, 774]</td></tr><tr><td rowspan="2">Simpl</td><td>CV̂</td><td>10.4%</td><td>23.7%</td><td>11.4%</td><td>10.4%</td><td>21.8%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 366, 59 608]</td><td>[523, 1 427]</td><td>[16 888, 26 644]</td><td>[21 641, 32 775]</td><td>[293, 733]</td></tr><tr><td rowspan="2">Simpl2</td><td>CV̂</td><td>10.6%</td><td>31.2%</td><td>11.8%</td><td>10.7%</td><td>29.5%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 216, 59 758]</td><td>[379, 1 571]</td><td>[16 727, 26 805]</td><td>[21 477, 32 939]</td><td>[216, 810]</td></tr><tr><td rowspan="3">Gross</td><td>CV̂</td><td>10.6%</td><td>35.2%</td><td>11.9%</td><td>10.9%</td><td>32.8%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 203, 59 771]</td><td>[302, 1 648]</td><td>[16 694, 26 838]</td><td>[21 382, 33 034]</td><td>[183, 843]</td></tr><tr><td>IC<sub>rev</sub></td><td>[39 580, 60 166]</td><td>[205, 1 540]</td><td>[16 710, 26 874]</td><td>[21 510, 33 086]</td><td>[154, 821]</td></tr><tr><td rowspan="3">Rao-Wu</td><td>CV̂</td><td>10.6%</td><td>34.1%</td><td>11.9%</td><td>10.8%</td><td>33.2%</td></tr><tr><td>IC<sub>norm</sub></td><td>[39 248, 59 726]</td><td>[324, 1 626]</td><td>[16 701, 26 831]</td><td>[21 447, 32 969]</td><td>[179, 847]</td></tr><tr><td>IC<sub>rev</sub></td><td>[39 635, 59 822]</td><td>[249, 1 525]</td><td>[16 816, 26 841]</td><td>[21 594, 33 009]</td><td>[140, 801]</td></tr></tbody></table>
<table><caption>Table 5: Estimations of ratios, with estimated coefficient of variations and confidence intervals</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Twins</th><th colspan="2">Marital status</th><th rowspan="2">Baby<br>Birth<br>weight</th><th colspan="2">Mother</th></tr><tr><th>Married/remarried</th><th>Divorced/single<br>/widowed</th><th>Weight<br>gain</th><th>IMC</th></tr></thead><tbody><tr><td>Î<sub>π</sub></td><td>1.97%</td><td>43.98%</td><td>54.98%</td><td>3 234</td><td>13.39</td><td>28.70</td></tr><tr><td>Unbiased<br>IC<sub>norm</sub></td><td>24.3%</td><td>4.0%</td><td>3.2%</td><td>0.3%</td><td>0.9%</td><td>0.5%</td></tr><tr><td>Simpl</td><td>19.7%</td><td>3.6%</td><td>2.9%</td><td>0.2%</td><td>0.6%</td><td>0.4%</td></tr><tr><td>Simpl2</td><td>28.1%</td><td>4.3%</td><td>3.4%</td><td>0.3%</td><td>1.1%</td><td>0.6%</td></tr><tr><td>Gross</td><td>32.6%</td><td>4.7%</td><td>3.7%</td><td>0.4%</td><td>1.3%</td><td>0.6%</td></tr><tr><td>IC<sub>norm</sub></td><td>[0.71%, 3.23%]</td><td>[39.96%, 48.00%]</td><td>[50.96%, 59.00%]</td><td>[3 209, 3 260]</td><td>[13.05, 13.73]</td><td>[28.34, 29.05]</td></tr><tr><td>IC<sub>rev</sub></td><td>[0.60%, 3.07%]</td><td>[39.89%, 47.99%]</td><td>[50.93%, 59.11%]</td><td>[3 209, 3 260]</td><td>[13.04, 13.73]</td><td>[28.33, 29.06]</td></tr><tr><td>Rao-Wu</td><td>31.1%</td><td>4.6%</td><td>3.7%</td><td>0.4%</td><td>1.3%</td><td>0.6%</td></tr><tr><td>IC<sub>norm</sub></td><td>[0.77%, 3.17%]</td><td>[40.04%, 47.92%]</td><td>[51.04%, 58.92%]</td><td>[3 209, 3 259]</td><td>[13.06, 13.72]</td><td>[28.35, 29.04]</td></tr><tr><td>IC<sub>rev</sub></td><td>[0.63%, 3.01%]</td><td>[39.92%, 47.92%]</td><td>[51.07%, 58.96%]</td><td>[3 209, 3 259]</td><td>[13.04, 13.70]</td><td>[28.34, 29.04]</td></tr></tbody></table>
<table><caption>Table 5: Estimations of ratios, with estimated coefficient of variations and confidence intervals</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Twins</th><th colspan="2">Marital status</th><th>Baby</th><th colspan="2">Mother</th></tr><tr><th>Married/remarried</th><th>Divorced/single<br>/widowed</th><th>Birth<br>weight</th><th>Weight<br>gain</th><th>IMC</th></tr></thead><tbody><tr><td>Î<sub>π</sub></td><td>1.97%</td><td>43.98%</td><td>54.98%</td><td>3 234</td><td>13.39</td><td>28.70</td></tr><tr><td>Unbiased<br>IC<sub>norm</sub></td><td>24.3%</td><td>4.0%</td><td>3.2%</td><td>0.3%</td><td>0.9%</td><td>0.5%</td></tr><tr><td>Simpl</td><td>19.7%</td><td>3.6%</td><td>2.9%</td><td>0.2%</td><td>0.6%</td><td>0.4%</td></tr><tr><td>Simpl2</td><td>28.1%</td><td>4.3%</td><td>3.4%</td><td>0.3%</td><td>1.1%</td><td>0.6%</td></tr><tr><td>Gross</td><td>32.6%</td><td>4.7%</td><td>3.7%</td><td>0.4%</td><td>1.3%</td><td>0.6%</td></tr><tr><td>IC<sub>norm</sub></td><td>[0.71%, 3.23%]</td><td>[39.96%, 48.00%]</td><td>[50.96%, 59.00%]</td><td>[3 209, 3 260]</td><td>[13.05, 13.73]</td><td>[28.34, 29.05]</td></tr><tr><td>IC<sub>rev</sub></td><td>[0.60%, 3.07%]</td><td>[39.89%, 47.99%]</td><td>[50.93%, 59.11%]</td><td>[3 209, 3 260]</td><td>[13.04, 13.73]</td><td>[28.33, 29.06]</td></tr><tr><td>Rao-Wu</td><td>31.1%</td><td>4.6%</td><td>3.7%</td><td>0.4%</td><td>1.3%</td><td>0.6%</td></tr><tr><td>IC<sub>norm</sub></td><td>[0.77%, 3.17%]</td><td>[40.04%, 47.92%]</td><td>[51.04%, 58.92%]</td><td>[3 209, 3 259]</td><td>[13.06, 13.72]</td><td>[28.35, 29.04]</td></tr><tr><td>IC<sub>rev</sub></td><td>[0.63%, 3.01%]</td><td>[39.92%, 47.92%]</td><td>[51.07%, 58.96%]</td><td>[3 209, 3 259]</td><td>[13.04, 13.70]</td><td>[28.34, 29.04]</td></tr></tbody></table>
stage sampling designs, where only the number of PSUs needs to approach infinity to ensure these
properties (Chauvet and Vallée, 2020).

Additionally, we delved into the intricate matter of variance estimation for CSS (Juillard et al., 2017). We derived the Hoeffding-Sobol decomposition of the HT-estimator and proposed two straightforward variance estimators by focusing on the leading terms in this decomposition. Theoretical and empirical analyses demonstrate the consistency of these estimators as sample sizes increase. Our empirical findings indicate that the plug-in variance estimator tends to be conservative with smaller sample sizes, exhibiting a moderate positive bias as long as interaction dispersion remains moderate. Therefore, if an analytic variance estimator is preferred, we recommend its use.

Furthermore, we explored bootstrap variance estimation, an area which is little addressed in the literature concerning CSS. By employing a conditional Hoeffding-Sobol decomposition, we established that weighted bootstrap techniques are applicable for CSS, provided a weighted bootstrap technique is available for the HT-estimator in each dimension. This greatly facilitates the implementation of simple bootstrap procedures for CSS. Particularly, we proposed adaptations of the population bootstrap and the rescaled bootstrap for simple random sampling in each dimension. Empirical results validate the asymptotic unbiasedness of both bootstrap techniques as sample sizes increase. Additionally, they suggest that the population bootstrap yields variance estimators that are less biased and more stable with smaller sample sizes compared to the rescaled bootstrap.

## Acknowledgement

The ELFE survey is a joint project between the French Institute for Demographic Studies (INED) and the National Institute of Health and Medical Research (INSERM), in partnership with the French blood transfusion service (Etablissement français du sang, EFS), Santé publique France, the National Institute for Statistics and Economic Studies (INSEE), the Direction générale de la santé (DGS, part of the Ministry of Health and Social Affairs), the Direction générale de la prévention des risques (DGPR, Ministry for the Environment), the Direction de la recherche, des études, de l'évaluation et des statistiques (DREES, Ministry of Health and Social Affairs), the Département des études, de la prospective et des statistiques (DEPS, Ministry of Culture), and the Caisse nationale des allocations familiales (CNAF), with the support of the Ministry of Higher Education
stage sampling designs, where only the number of PSUs needs to approach infinity to ensure these
properties (Chauvet and Vallée, 2020).

Additionally, we delved into the intricate matter of variance estimation for CSS (Juillard et al., 2017). We derived the Hoeffding-Sobol decomposition of the HT-estimator and proposed two straightforward variance estimators by focusing on the leading terms in this decomposition. Theoretical and empirical analyses demonstrate the consistency of these estimators as sample sizes increase. Our empirical findings indicate that the plug-in variance estimator tends to be conservative with smaller sample sizes, exhibiting a moderate positive bias as long as interaction dispersion remains moderate. Therefore, if an analytic variance estimator is preferred, we recommend its use.

Furthermore, we explored bootstrap variance estimation, an area which is little addressed in the literature concerning CSS. By employing a conditional Hoeffding-Sobol decomposition, we established that weighted bootstrap techniques are applicable for CSS, provided a weighted bootstrap technique is available for the HT-estimator in each dimension. This greatly facilitates the implementation of simple bootstrap procedures for CSS. Particularly, we proposed adaptations of the population bootstrap and the rescaled bootstrap for simple random sampling in each dimension. Empirical results validate the asymptotic unbiasedness of both bootstrap techniques as sample sizes increase. Additionally, they suggest that the population bootstrap yields variance estimators that are less biased and more stable with smaller sample sizes compared to the rescaled bootstrap.

## Acknowledgement

The ELFE survey is a joint project between the French Institute for Demographic Studies (INED) and the National Institute of Health and Medical Research (INSERM), in partnership with the French blood transfusion service (Etablissement français du sang, EFS), Santé publique France, the National Institute for Statistics and Economic Studies (INSEE), the Direction générale de la santé (DGS, part of the Ministry of Health and Social Affairs), the Direction générale de la prévention des risques (DGPR, Ministry for the Environment), the Direction de la recherche, des études, de l'évaluation et des statistiques (DREES, Ministry of Health and Social Affairs), the Département des études, de la prospective et des statistiques (DEPS, Ministry of Culture), and the Caisse nationale des allocations familiales (CNAF), with the support of the Ministry of Higher Education
and Research and the Institut national de la jeunesse et de l'éducation populaire (INJEP). Via the RECONAI platform, it receives a government grant managed by the National Research Agency under the "Investissements d'avenir" programme (ANR-11-EQPX-0038 and ANR-19-COHO-0001).

# A Proofs for Section 3.1

## A.1 Proof of equation (3.1)

From equation (2.6), a straightforward computation leads to the result:

$$V_p(\hat{Y}) = \sum_{\mathbf{k},l \in \mathcal{U}} \mathrm{Cov}_p \left( y_{\mathbf{k}} \prod_{d=1}^{D} w_{k_d}^d, y_{l} \prod_{d=1}^{D} w_{l_d}^d \right) \quad (A.1)$$

$$= \sum_{\mathbf{k},l \in \mathcal{U}} y_{\mathbf{k}} y_{l} \left[ E_p \left( \prod_{d=1}^{D} w_{k_d}^d w_{l_d}^d \right) - E_p \left( \prod_{d=1}^{D} w_{k_d}^d \right) E_p \left( \prod_{d=1}^{D} w_{l_d}^d \right) \right] \quad (A.2)$$

$$= \sum_{\mathbf{k},l \in \mathcal{U}} y_{\mathbf{k}} y_{l} \left[ \prod_{d=1}^{D} E_p(w_{k_d}^d w_{l_d}^d) - \prod_{d=1}^{D} E_p(w_{k_d}^d) E_p(w_{l_d}^d) \right] \quad (A.3)$$

where the last line follows from the independence of the weights in each dimension.

## A.2 Proof of Proposition 1

Equations (3.2) as well as (3.5) in the first part of the proposition are general properties of the Hoeffding-Sobol decomposition applied to $\hat{Y}$ seen as a function of the independent variables $(S_d)_{d=1}^D$, the general formula corresponding to (3.3). It remains to check if the equation (3.4) is true in our case by computing for each $I \subseteq \{1, \dots, D\}$ the associated $\hat{Y}^I$ starting from the equation (3.3)

$$\hat{Y}^I = \sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} E_p(\hat{Y}|(S_d)_{d \in I'}) \quad (A.4)$$

$$= \sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d=1}^{D} \left[ (w_{k_d}^d - E_p(w_{k_d}^d)) \mathbf{1}_{(d \in I')} + E_p(w_{k_d}^d) \right] \quad (A.5)$$

$$= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} \sum_{I'' \subseteq \{1, \dots, D\}, d \in I''} \prod_{d \in I''} (w_{k_d}^d - E_p(w_{k_d}^d)) \mathbf{1}_{(d \in I')} \prod_{d \notin I''} E_p(w_{k_d}^d). \quad (A.6)$$

By observing that for any $I'' \subseteq \{1, \dots, D\}$, $\prod_{d \in I''} \mathbf{1}_{(d \in I')} = \mathbf{1}_{(I'' \subseteq I')}$ and that $\sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} \mathbf{1}_{(I'' \subseteq I')} = \mathbf{1}_{(I''=I)}$, one can conclude by rearranging the sum that

$$\hat{Y}^I = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \sum_{I'' \subseteq \{1, \dots, D\}, d \in I''} \prod_{d \in I''} (w_{k_d}^d - E_p(w_{k_d}^d)) \prod_{d \notin I''} E_p(w_{k_d}^d) \sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} \mathbf{1}_{(I'' \subseteq I')} \quad (A.7)$$

$$= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d \in I} (w_{k_d}^d - E_p(w_{k_d}^d)) \prod_{d \notin I} E_p(w_{k_d}^d), \quad (A.8)$$
and Research and the Institut national de la jeunesse et de l'éducation populaire (INJEP). Via the RECONAI platform, it receives a government grant managed by the National Research Agency under the "Investissements d'avenir" programme (ANR-11-EQPX-0038 and ANR-19-COHO-0001).

# A Proofs for Section 3.1

## A.1 Proof of equation (3.1)

From equation (2.6), a straightforward computation leads to the result:

$$V_p(\hat{Y}) = \sum_{\mathbf{k},l \in \mathcal{U}} \mathrm{Cov}_p \left( y_{\mathbf{k}} \prod_{d=1}^{D} w_{k_d}^{d}, y_{l} \prod_{d=1}^{D} w_{l_d}^{d} \right) \quad (A.1)$$

$$= \sum_{\mathbf{k},l \in \mathcal{U}} y_{\mathbf{k}} y_{l} \left[ E_p \left( \prod_{d=1}^{D} w_{k_d}^{d} w_{l_d}^{d} \right) - E_p \left( \prod_{d=1}^{D} w_{k_d}^{d} \right) E_p \left( \prod_{d=1}^{D} w_{l_d}^{d} \right) \right] \quad (A.2)$$

$$= \sum_{\mathbf{k},l \in \mathcal{U}} y_{\mathbf{k}} y_{l} \left[ \prod_{d=1}^{D} E_p(w_{k_d}^d w_{l_d}^d) - \prod_{d=1}^{D} E_p(w_{k_d}^d) E_p(w_{l_d}^d) \right] \quad (A.3)$$

where the last line follows from the independence of the weights in each dimension.

## A.2 Proof of Proposition 1

Equations (3.2) as well as (3.5) in the first part of the proposition are general properties of the Hoeffding-Sobol decomposition applied to $\hat{Y}$ seen as a function of the independent variables $(S_d)_{d=1}^D$, the general formula corresponding to (3.3). It remains to check if the equation (3.4) is true in our case by computing for each $I \subseteq \{1, \dots, D\}$ the associated $\hat{Y}^I$ starting from the equation (3.3)

$$\hat{Y}^I = \sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} E_p(\hat{Y}|(S_d)_{d \in I'}) \quad (A.4)$$

$$= \sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d=1}^{D} \left[ (w_{k_d}^d - E_p(w_{k_d}^d)) \mathbf{1}_{(d \in I')} + E_p(w_{k_d}^d) \right] \quad (A.5)$$

$$= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} \sum_{I'' \subseteq \{1, \dots, D\}, d \in I''} \prod_{d \in I''} (w_{k_d}^d - E_p(w_{k_d}^d)) \mathbf{1}_{(d \in I')} \prod_{d \notin I''} E_p(w_{k_d}^d). \quad (A.6)$$

By observing that for any $I'' \subseteq \{1, \dots, D\}$, $\prod_{d \in I''} \mathbf{1}_{(d \in I')} = \mathbf{1}_{(I'' \subseteq I')}$ and that $\sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} \mathbf{1}_{(I'' \subseteq I')} = \mathbf{1}_{(I''=I)}$, one can conclude by rearranging the sum that

$$\hat{Y}^I = \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \sum_{I'' \subseteq \{1, \dots, D\}, d \in I''} \prod_{d \in I''} (w_{k_d}^d - E_p(w_{k_d}^d)) \prod_{d \notin I''} E_p(w_{k_d}^d) \sum_{I' \in \mathcal{P}(I)} (-1)^{|I|-|I'|} \mathbf{1}_{(I'' \subseteq I')} \quad (A.7)$$

$$= \sum_{\mathbf{k} \in \mathcal{U}} y_{\mathbf{k}} \prod_{d \in I} (w_{k_d}^d - E_p(w_{k_d}^d)) \prod_{d \notin I} E_p(w_{k_d}^d), \quad (A.8)$$
from which we can deduce the value of $V_p(\hat{Y}^I)$ for $I \neq \emptyset$:

$$V_p(\hat{Y}^I) = \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} \mathrm{Cov}_p \left( \prod_{d \in I} (w_{k_d}^d - E_p(w_{k_d}^d)), \prod_{d \in I} (w_{l_d}^d - E_p(w_{l_d}^d)) \right) \prod_{d \notin I} E_p(w_{k_d}^d) E_p(w_{l_d}^d) \quad (A.9)$$

$$= \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} E_p \left( \prod_{d \in I} (w_{k_d}^d - E_p(w_{k_d}^d))(w_{l_d}^d - E_p(w_{l_d}^d)) \right) \prod_{d \notin I} E_p(w_{k_d}^d) E_p(w_{l_d}^d) \quad (A.10)$$

$$= \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} \prod_{d \in I} \mathrm{Cov}_p(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} E_p(w_{k_d}^d) E_p(w_{l_d}^d). \quad (A.11)$$

where the last line follows from the independence of the weights in each dimension.

### A.3 Proof of Corollary 1

It suffices to show that for every non-empty subset $I \subseteq \{1, \dots, D\}$, $V_p(\hat{Y}^I)$ is unbiasedly estimated by $\hat{V}_p(\hat{Y}^I)$. It is possible to rewrite the sum over $\mathcal{U}$ by introducing the sample membership indicators. We then obtain

$$E_p(\hat{V}_p(\hat{Y}^I)) = \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} E_p \left[ \prod_{d \in I} \widehat{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} \frac{\delta_{k_d}^d \delta_{l_d}^d}{\pi_{k_d, l_d}^d} \right] \quad (A.12)$$

$$= \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} \prod_{d \in I} E_p \left[ \widehat{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} \frac{E_p(\delta_{k_d}^d \delta_{l_d}^d)}{\pi_{k_d, l_d}^d} \right] \quad (A.13)$$

$$= \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} \prod_{d \in I} \mathrm{Cov}_p(w_{k_d}^d, w_{l_d}^d) \quad (A.14)$$

$$= V_p(\hat{Y}^I). \quad (A.15)$$

where the second line follows from the fact that for a fixed $\mathbf{k}$ and $\mathbf{l} \in \mathcal{U}$, each $\widehat{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d)$ and $\delta_{k_d}^d \delta_{l_d}^d$ are built from the sample $S_d$ only and are thus all independent from each other as long as they correspond to different dimensions.

## B Proofs for Section 3.2

### B.1 Proof of Corollary 2

Equation (3.12) is a direct consequence of equation (3.4), applied for $w_{k_d}^d = \delta_{k_d}^d / \pi_{k_d}^d$. Similarly, equation (3.13) follows from equation (3.6). Finally, equation (3.14) follows by applying equation (3.9) with

$$\widehat{\mathrm{Cov}}^d(w_k^d, w_l^d) = \frac{\Delta_{k_d l_d}^d \delta_{k_d}^d \delta_{l_d}^d}{\pi_{k_d l_d}^d \pi_{k_d}^d \pi_{l_d}^d}. \quad (B.1)$$
from which we can deduce the value of $V_p(\hat{Y}^I)$ for $I \neq \emptyset$:

$$V_p(\hat{Y}^I) = \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} \mathrm{Cov}_p \left( \prod_{d \in I} (w_{k_d}^d - E_p(w_{k_d}^d)), \prod_{d \in I} (w_{l_d}^d - E_p(w_{l_d}^d)) \right) \prod_{d \notin I} E_p(w_{k_d}^d) E_p(w_{l_d}^d) \quad (A.9)$$

$$= \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} E_p \left( \prod_{d \in I} (w_{k_d}^d - E_p(w_{k_d}^d))(w_{l_d}^d - E_p(w_{l_d}^d)) \right) \prod_{d \notin I} E_p(w_{k_d}^d) E_p(w_{l_d}^d) \quad (A.10)$$

$$= \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} \prod_{d \in I} \mathrm{Cov}_p(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} E_p(w_{k_d}^d) E_p(w_{l_d}^d). \quad (A.11)$$

where the last line follows from the independence of the weights in each dimension.

### A.3 Proof of Corollary 1

It suffices to show that for every non-empty subset $I \subseteq \{1, \dots, D\}$, $V_p(\hat{Y}^I)$ is unbiasedly estimated by $\hat{V}_p(\hat{Y}^I)$. It is possible to rewrite the sum over $\mathcal{U}$ by introducing the sample membership indicators. We then obtain

$$E_p(\hat{V}_p(\hat{Y}^I)) = \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} E_p \left[ \prod_{d \in I} \widehat{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} \frac{\delta_{k_d}^d \delta_{l_d}^d}{\pi_{k_d, l_d}^d} \right] \quad (A.12)$$

$$= \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} \prod_{d \in I} E_p \left[ \widehat{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d) \prod_{d \notin I} \frac{E_p(\delta_{k_d}^d \delta_{l_d}^d)}{\pi_{k_d, l_d}^d} \right] \quad (A.13)$$

$$= \sum_{\mathbf{k},\mathbf{l} \in \mathcal{U}} y_{\mathbf{k}} y_{\mathbf{l}} \prod_{d \in I} \mathrm{Cov}_p(w_{k_d}^d, w_{l_d}^d) \quad (A.14)$$

$$= V_p(\hat{Y}^I). \quad (A.15)$$

where the second line follows from the fact that for a fixed $\mathbf{k}$ and $\mathbf{l} \in \mathcal{U}$, each $\widehat{\mathrm{Cov}}^d(w_{k_d}^d, w_{l_d}^d)$ and $\delta_{k_d}^d \delta_{l_d}^d$ are built from the sample $S_d$ only and are thus all independent from each other as long as they correspond to different dimensions.

## B Proofs for Section 3.2

### B.1 Proof of Corollary 2

Equation (3.12) is a direct consequence of equation (3.4), applied for $w_{k_d}^d = \delta_{k_d}^d / \pi_{k_d}^d$. Similarly, equation (3.13) follows from equation (3.6). Finally, equation (3.14) follows by applying equation (3.9) with

$$\widehat{\mathrm{Cov}}^d(w_k^d, w_l^d) = \frac{\Delta_{k_d l_d}^d \delta_{k_d}^d \delta_{l_d}^d}{\pi_{k_d l_d}^d \pi_{k_d}^d \pi_{l_d}^d}. \quad (B.1)$$
## B.2 Proof of Proposition 2

Let us introduce the partition $\mathcal{U}_I \times \mathcal{U}_I = \cup_{I' \in \mathcal{P}(I)} P_{I'}$ for a given non-empty $I \subseteq \{1, \dots, D\}$, where we defined for each $I' \in \mathcal{P}(I)$, the set $P_{I'}$ by

$$ P_{I'} = \{ \mathbf{k}', \mathbf{l}' \in \mathcal{U}_I, \forall d \in I, (k'_d \neq l'_d \Leftrightarrow d \in I') \} \quad (B.2) $$

furthermore, note that we can identify every element $\mathbf{k} \in \mathcal{U}$ with a couple $(\mathbf{k}_1, \mathbf{k}_2) \in \mathcal{U}_I \times \mathcal{U}_{I^c}$ where we introduced $\mathcal{U}_{I^c} = \prod_{d \notin I} \mathcal{U}_d$. More precisely, the coordinates of $\mathbf{k}_1$ are given by the coordinates of $\mathbf{k}$ associated to the dimensions in $I$ and the coordinates of $\mathbf{k}_2$ are given by the coordinates of $\mathbf{k}$ associated to the dimensions that are not in $I$. We will therefore denote by $y_{\mathbf{k}_1 \mathbf{k}_2}$ the quantity $y_{\mathbf{k}}$. This splitting allows us for example to rewrite the subtotal $Y_{\mathbf{k}'}$ in the form $Y_{\mathbf{k}'} = \sum_{\mathbf{l}' \in \mathcal{U}_{I^c}} y_{\mathbf{k}' \mathbf{l}'}$. Now let us fix some subset $I' \in \mathcal{P}(I)$. We can verify using the Cauchy-Schwarz inequality that

$$ \sum_{\mathbf{k}' \in \mathcal{U}_I} Y_{\mathbf{k}'}^2 = \sum_{\mathbf{k}' \in \mathcal{U}_I} \left( \sum_{\mathbf{l}' \in \mathcal{U}_{I^c}} y_{\mathbf{k}' \mathbf{l}'} \right)^2 \le \left( \prod_{d \notin I} N_d \right) \sum_{\mathbf{k}' \in \mathcal{U}_I} \sum_{\mathbf{l}' \in \mathcal{U}_{I^c}} y_{\mathbf{k}' \mathbf{l}'}^2 \le \alpha \frac{N^2}{\prod_{d \in I} N_d} \quad (B.3) $$

by using the assumption (H1).

Now we can similarly identify an element $\mathbf{k}' \in \mathcal{U}_I$ with a couple $(\mathbf{k}'_1, \mathbf{k}'_2) \in \mathcal{U}_{I'} \times \mathcal{U}_{I \setminus I'}$. We will therefore denote by $Y_{\mathbf{k}'_1 \mathbf{k}'_2}$ the subtotal $Y_{\mathbf{k}'}$. Thus, using again the Cauchy-Scharz inequality, it is now easy to verify that

$$ \sum_{(\mathbf{k}', \mathbf{l}') \in P_{I'}} Y_{\mathbf{k}'} Y_{\mathbf{l}'} \le \sum_{\substack{\mathbf{k}'', \mathbf{l}'' \in \mathcal{U}_{I'} \\ \mathbf{m}'' \in \mathcal{U}_{I \setminus I'}}} Y_{\mathbf{k}'' \mathbf{m}''} Y_{\mathbf{l}'' \mathbf{m}''} \quad (B.4) $$

$$ = \sum_{\mathbf{m}'' \in \mathcal{U}_{I \setminus I'}} \left( \sum_{\mathbf{k}'' \in \mathcal{U}_{I'}} Y_{\mathbf{k}'' \mathbf{m}''} \right)^2 \quad (B.5) $$

$$ \le \left( \prod_{d \in I'} N_d \right) \sum_{\mathbf{m}'' \in \mathcal{U}_{I \setminus I'}} \sum_{\mathbf{k}'' \in \mathcal{U}_{I'}} Y_{\mathbf{k}'' \mathbf{m}''}^2 \quad (B.6) $$

$$ \le \alpha \frac{N^2}{\prod_{d \in I \setminus I'} N_d}. \quad (B.7) $$

Furthermore we know by construction of $P_{I'}$ that for any $(\mathbf{k}', \mathbf{l}') \in P_{I'}$, $k'_d \neq l'_d$ if and only if $d \in I'$ and thus using assumption (H3)-(H4) we get that

$$ \forall (\mathbf{k}', \mathbf{l}') \in P_{I'}, \forall d \in I', \left| \frac{\Delta_{k'_d, l'_d}^d}{\pi_{k'_d}^d \pi_{l'_d}^d} \right| \le \frac{\gamma_d}{\lambda_d^2} \frac{1}{n_d}. \quad (B.8) $$

On the other hand, when $d \notin I'$, we have that $k'_d = l'_d$ and thus

$$ \forall (\mathbf{k}', \mathbf{l}') \in P_{I'}, \forall d \notin I', \left| \frac{\Delta_{k'_d, l'_d}^d}{\pi_{k'_d}^d \pi_{l'_d}^d} \right| = \frac{1 - \pi_{k'_d}^d}{\pi_{k'_d}^d} \le \frac{1}{\lambda_d} \frac{N_d}{n_d}. \quad (B.9) $$
## B.2 Proof of Proposition 2

Let us introduce the partition $\mathcal{U}_I \times \mathcal{U}_I = \cup_{I' \in \mathcal{P}(I)} P_{I'}$ for a given non-empty $I \subseteq \{1, \dots, D\}$, where we defined for each $I' \in \mathcal{P}(I)$, the set $P_{I'}$ by

$$P_{I'} = \{\mathbf{k}', \mathbf{l}' \in \mathcal{U}_I, \forall d \in I, (k'_d \neq l'_d \Leftrightarrow d \in I')\} \quad (B.2)$$

furthermore, note that we can identify every element $\mathbf{k} \in \mathcal{U}$ with a couple $(\mathbf{k}_1, \mathbf{k}_2) \in \mathcal{U}_I \times \mathcal{U}_{I^c}$ where we introduced $\mathcal{U}_{I^c} = \prod_{d \notin I} \mathcal{U}_d$. More precisely, the coordinates of $\mathbf{k}_1$ are given by the coordinates of $\mathbf{k}$ associated to the dimensions in $I$ and the coordinates of $\mathbf{k}_2$ are given by the coordinates of $\mathbf{k}$ associated to the dimensions that are not in $I$. We will therefore denote by $y_{\mathbf{k}_1 \mathbf{k}_2}$ the quantity $y_{\mathbf{k}}$. This splitting allows us for example to rewrite the subtotal $Y_{\mathbf{k}'}$ in the form $Y_{\mathbf{k}'} = \sum_{\mathbf{l}' \in \mathcal{U}_{I^c}} y_{\mathbf{k}'}$. Now let us fix some subset $I' \in \mathcal{P}(I)$. We can verify using the Cauchy-Schwarz inequality that

$$\sum_{\mathbf{k}' \in \mathcal{U}_I} Y_{\mathbf{k}'}^2 = \sum_{\mathbf{k}' \in \mathcal{U}_I} \left( \sum_{\mathbf{l}' \in \mathcal{U}_{I^c}} y_{\mathbf{k}'\mathbf{l}'} \right)^2 \le \left( \prod_{d \notin I} N_d \right) \sum_{\mathbf{k}' \in \mathcal{U}_I} \sum_{\mathbf{l}' \in \mathcal{U}_{I^c}} y_{\mathbf{k}'\mathbf{l}'}^2 \le \alpha \frac{N^2}{\prod_{d \in I} N_d} \quad (B.3)$$

by using the assumption (H1).

Now we can similarly identify an element $\mathbf{k}' \in \mathcal{U}_I$ with a couple $(\mathbf{k}'_1, \mathbf{k}'_2) \in \mathcal{U}_{I'} \times \mathcal{U}_{I \setminus I'}$. We will therefore denote by $Y_{\mathbf{k}'_1 \mathbf{k}'_2}$ the subtotal $Y_{\mathbf{k}'}$. Thus, using again the Cauchy-Scharz inequality, it is now easy to verify that

$$\sum_{(\mathbf{k}', \mathbf{l}') \in P_{I'}} Y_{\mathbf{k}'} Y_{\mathbf{l}'} \le \sum_{\substack{\mathbf{k}'', \mathbf{l}'' \in \mathcal{U}_{I'} \\ \mathbf{m}'' \in \mathcal{U}_{I \setminus I'}}} Y_{\mathbf{k}'' \mathbf{m}''} Y_{\mathbf{l}'' \mathbf{m}''} \quad (B.4)$$

$$= \sum_{\mathbf{m}'' \in \mathcal{U}_{I \setminus I'}} \left( \sum_{\mathbf{k}'' \in \mathcal{U}_{I'}} Y_{\mathbf{k}'' \mathbf{m}''} \right)^2 \quad (B.5)$$

$$\le \left( \prod_{d \in I'} N_d \right) \sum_{\mathbf{m}'' \in \mathcal{U}_{I \setminus I'}} \sum_{\mathbf{k}'' \in \mathcal{U}_{I'}} Y_{\mathbf{k}'' \mathbf{m}''}^2 \quad (B.6)$$

$$\le \alpha \frac{N^2}{\prod_{d \in I \setminus I'} N_d} \quad (B.7)$$

Furthermore we know by construction of $P_{I'}$ that for any $(\mathbf{k}', \mathbf{l}') \in P_{I'}$, $k'_d \neq l'_d$ if and only if $d \in I'$ and thus using assumption (H3)-(H4) we get that

$$\forall (\mathbf{k}', \mathbf{l}') \in P_{I'}, \forall d \in I', \left| \frac{\Delta_{k'_d, l'_d}^d}{\pi_{k'_d}^d \pi_{l'_d}^d} \right| \le \frac{\gamma_d}{\lambda_d^2} \frac{1}{n_d}. \quad (B.8)$$

On the other hand, when $d \notin I'$, we have that $k'_d = l'_d$ and thus

$$\forall (\mathbf{k}', \mathbf{l}') \in P_{I'}, \forall d \notin I', \left| \frac{\Delta_{k'_d, l'_d}^d}{\pi_{k'_d}^d \pi_{l'_d}^d} \right| = \frac{1 - \pi_{k'_d}^d}{\pi_{k'_d}^d} \le \frac{1}{\lambda_d} \frac{N_d}{n_d}. \quad (B.9)$$
From there, we can build an upper bound of the sum

$$ \sum_{(\mathbf{k}', \mathbf{I}') \in P_{I'}} Y_{\mathbf{k}'} Y_{\mathbf{I}'} \prod_{d \in I} \frac{\Delta_{\mathbf{k}'d l'_d}^d}{\pi_{\mathbf{k}'d}^d \pi_{\mathbf{I}'d}^d} = \sum_{(\mathbf{k}', \mathbf{I}) \in P_{I'}} Y_{\mathbf{k}'} Y_{\mathbf{I}} \prod_{d \in I'} \frac{\Delta_{\mathbf{k}'d l'_d}^d}{\pi_{\mathbf{k}'d}^d \pi_{\mathbf{I}'d}^d} \prod_{d \in I \setminus I'} \frac{\Delta_{\mathbf{k}'d l'_d}^d}{\pi_{\mathbf{k}'d}^d \pi_{\mathbf{I}'d}^d} \quad (B.10) $$

$$ = O \left( \frac{N^2}{\prod_{d \in I \setminus I'} N_d} \prod_{d \in I'} \frac{1}{n_d} \prod_{d \in I \setminus I'} \frac{N_d}{n_d} \right) \quad (B.11) $$

$$ = O \left( \frac{N^2}{\prod_{d \in I} n_d} \right). \quad (B.12) $$

Now, we can finally conclude that

$$ V_p(\hat{Y}^I) = \sum_{I' \in \mathcal{P}(I)} \left( \sum_{(\mathbf{k}', \mathbf{I}) \in P_{I'}} Y_{\mathbf{k}'} Y_{\mathbf{I}}' \prod_{d \in I} \frac{\Delta_{\mathbf{k}'d l'_d}^d}{\pi_{\mathbf{k}'d}^d \pi_{\mathbf{I}'d}^d} \right) = O \left( \frac{N^2}{\prod_{d \in I} n_d} \right). \quad (B.13) $$

## B.3 Proof of Proposition 3

We can write

$$ \hat{Y}_{\pi} - Y = \sum_{d=1}^{D} \hat{Y}_{\pi}^{\{d\}} + \Delta, \text{ where } \Delta = \sum_{\substack{I \subset \{1, \dots, D\} \\ \mathrm{Card}(I) \ge 2}} \hat{Y}_{\pi}^{I}. $$

We obtain

$$ V_p(\hat{Y}_\pi) = \sum_{d=1}^{D} V_p(\hat{Y}_\pi^{\{d\}}) + V_p(\Delta). \quad (B.14) $$

It follows from Proposition 2 that

$$ V_p(N^{-1}\hat{Y}_\pi^{\{d\}}) = O(n_m^{-1}) \quad \text{and} \quad V_p(N^{-1}\Delta) = o(n_m^{-1}). \quad (B.15) $$

Therefore, we obtain (3.19). Note that from equations (B.14) and (B.15), the constants in assumption (H7) are such that $\sum_{d=1}^{D}(\gamma_d)^2 = 1$.

From equation (B.14), we also obtain

$$ \frac{\hat{Y}_{\pi} - Y}{\sqrt{V_p(\hat{Y}_{\pi})}} = \underbrace{\sum_{d=1}^{D} \gamma_d \frac{\hat{Y}_{\pi}^{\{d\}}}{\sqrt{V_p(\hat{Y}_{\pi}^{\{d\}})}}}_{\Delta_1} + \underbrace{\sum_{d=1}^{D} \left( \sqrt{\frac{V_p(\hat{Y}_{\pi}^{\{d\}})}{V_p(\hat{Y}_{\pi})}} - \gamma_d \right) \frac{\hat{Y}_{\pi}^{\{d\}}}{\sqrt{V_p(\hat{Y}_{\pi}^{\{d\}})}}}_{\Delta_2} + \underbrace{\frac{\Delta}{\sqrt{V_p(\hat{Y}_{\pi})}}}_{\Delta_3}. \quad (B.16) $$

By using assumption (H5) and the right-hand side of (B.14), we have $\Delta_3 \to_{Pr} 0$, where $\to_{Pr}$ stands for the convergence in probability. Also, by using Assumptions (H6) and (H7) and the Slutsky theorem, we obtain that $\Delta_2 \to_{Pr} 0$. Finally, since the variables $\hat{Y}_{\pi}^{\{d\}}$ are independent, we obtain from assumption (H6) by standard arguments that $\Delta_1 \to \mathcal{N}(0, 1)$, which completes the proof.
From there, we can build an upper bound of the sum

$$ \sum_{(\mathbf{k}', \mathbf{I}') \in P_{I'}} Y_{\mathbf{k}'} Y_{\mathbf{I}'} \prod_{d \in I} \frac{\Delta_{\mathbf{k}'d l'_d}^d}{\pi_{\mathbf{k}'d}^d \pi_{l'_d}^d} = \sum_{(\mathbf{k}', \mathbf{I}) \in P_{I'}} Y_{\mathbf{k}'} Y_{\mathbf{I}} \prod_{d \in I'} \frac{\Delta_{\mathbf{k}'d l'_d}^d}{\pi_{\mathbf{k}'d}^d \pi_{l'_d}^d} \prod_{d \in I \setminus I'} \frac{\Delta_{\mathbf{k}'d l'_d}^d}{\pi_{\mathbf{k}'d}^d \pi_{l'_d}^d} \quad (B.10) $$

$$ = O \left( \frac{N^2}{\prod_{d \in I \setminus I'} N_d} \prod_{d \in I'} \frac{1}{n_d} \prod_{d \in I \setminus I'} \frac{N_d}{n_d} \right) \quad (B.11) $$

$$ = O \left( \frac{N^2}{\prod_{d \in I} n_d} \right). \quad (B.12) $$

Now, we can finally conclude that

$$ V_p(\hat{Y}^I) = \sum_{I' \in \mathcal{P}(I)} \left( \sum_{(\mathbf{k}', \mathbf{I}') \in P_{I'}} Y_{\mathbf{k}'} Y_{\mathbf{I}'} \prod_{d \in I} \frac{\Delta_{\mathbf{k}'d l'_d}^d}{\pi_{\mathbf{k}'d}^d \pi_{l'_d}^d} \right) = O \left( \frac{N^2}{\prod_{d \in I} n_d} \right). \quad (B.13) $$

## B.3 Proof of Proposition 3

We can write

$$ \hat{Y}_{\pi} - Y = \sum_{d=1}^{D} \hat{Y}_{\pi}^{\{d\}} + \Delta, \text{ where } \Delta = \sum_{\substack{I \subset \{1, \dots, D\} \\ \mathrm{Card}(I) \ge 2}} \hat{Y}_{\pi}^{I}. $$

We obtain

$$ V_p(\hat{Y}_\pi) = \sum_{d=1}^{D} V_p(\hat{Y}_\pi^{\{d\}}) + V_p(\Delta). \quad (B.14) $$

It follows from Proposition 2 that

$$ V_p(N^{-1}\hat{Y}_\pi^{\{d\}}) = O(n_m^{-1}) \quad \text{and} \quad V_p(N^{-1}\Delta) = o(n_m^{-1}). \quad (B.15) $$

Therefore, we obtain (3.19). Note that from equations (B.14) and (B.15), the constants in assumption (H7) are such that $\sum_{d=1}^{D}(\gamma_d)^2 = 1$.

From equation (B.14), we also obtain

$$ \frac{\hat{Y}_{\pi} - Y}{\sqrt{V_p(\hat{Y}_{\pi})}} = \underbrace{\sum_{d=1}^{D} \gamma_d \frac{\hat{Y}_{\pi}^{\{d\}}}{\sqrt{V_p(\hat{Y}_{\pi}^{\{d\}})}}}_{\Delta_1} + \underbrace{\sum_{d=1}^{D} \left( \sqrt{\frac{V_p(\hat{Y}_{\pi}^{\{d\}})}{V_p(\hat{Y}_{\pi})}} - \gamma_d \right) \frac{\hat{Y}_{\pi}^{\{d\}}}{\sqrt{V_p(\hat{Y}_{\pi}^{\{d\}})}}}_{\Delta_2} + \underbrace{\frac{\Delta}{\sqrt{V_p(\hat{Y}_{\pi})}}}_{\Delta_3}. \quad (B.16) $$

By using assumption (H5) and the right-hand side of (B.14), we have $\Delta_3 \to_{P_r} 0$, where $\to_{P_r}$ stands for the convergence in probability. Also, by using Assumptions (H6) and (H7) and the Slutsky theorem, we obtain that $\Delta_2 \to_{P_r} 0$. Finally, since the variables $\hat{Y}_{\pi}^{\{d\}}$ are independent, we obtain from assumption (H6) by standard arguments that $\Delta_1 \to \mathcal{N}(0, 1)$, which completes the proof.
## B.4 Proof of Proposition 4

Let us fix a non-empty subset $I \subseteq \{1, \dots, D\}$. We can write using the definitions of $\hat{V}^{PLUG}(\hat{Y}_{\pi}^{I})$

$$E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^{I})] = \sum_{\mathbf{k}', \mathbf{l}' \in \mathcal{U}_I} \frac{E_p(\hat{Y}_{\mathbf{k}'} \hat{Y}_{\mathbf{l}'} \delta_{\mathbf{k}'} \delta_{\mathbf{l}'}) \Delta_{\mathbf{k}'\mathbf{l}'}}{\pi_{\mathbf{k}'} \pi_{\mathbf{l}'}}, \quad (B.17)$$

$$= \sum_{\mathbf{k}', \mathbf{l}' \in \mathcal{U}_I} \frac{E_p(\hat{Y}_{\mathbf{k}'} \hat{Y}_{\mathbf{l}'}) E_p(\delta_{\mathbf{k}'} \delta_{\mathbf{l}'}) \Delta_{\mathbf{k}'\mathbf{l}'}}{\pi_{\mathbf{k}'} \pi_{\mathbf{l}'} \pi_{\mathbf{k}'} \pi_{\mathbf{l}'}} \quad (B.18)$$

$$= \sum_{\mathbf{k}', \mathbf{l}' \in \mathcal{U}_I} \frac{E_p(\hat{Y}_{\mathbf{k}'} \hat{Y}_{\mathbf{l}'}) \Delta_{\mathbf{k}'\mathbf{l}'}}{\pi_{\mathbf{k}'} \pi_{\mathbf{l}'} \pi_{\mathbf{k}'} \pi_{\mathbf{l}'}} \quad (B.19)$$

Where the line (B.18) follows from the fact that $\delta_{\mathbf{k}'}$ and $\delta_{\mathbf{l}'}$ are functions of $(S_d)_{d \in I}$ whereas $\hat{Y}_{\mathbf{k}'}$ and $\hat{Y}_{\mathbf{l}'}$ are functions of $(S_d)_{d \notin I}$, and are therefore independent. We can then similarly expand the sum by replacing $\hat{Y}_{\mathbf{k}'}$ and $\hat{Y}_{\mathbf{l}'}$ by their definitions in (3.25):

$$E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^{I})] = \sum_{\mathbf{k}', \mathbf{l}' \in \mathcal{U}_I} \sum_{\substack{\mathbf{k}, \mathbf{l} \in \mathcal{U} \\ \forall d \in I, k_d=k_{d'}, l_d=l_{d}'}} \frac{y_{\mathbf{k}} y_{\mathbf{l}}}{\pi_{\mathbf{k}} \pi_{\mathbf{l}}} \prod_{d \in I} \Delta_{k_d l_d}^d E_p \left( \prod_{d \notin I} \delta_{k_d}^d \delta_{l_d}^d \right) \quad (B.20)$$

$$= \sum_{\mathbf{k}, \mathbf{l} \in \mathcal{U}} \frac{y_{\mathbf{k}} y_{\mathbf{l}}}{\pi_{\mathbf{k}} \pi_{\mathbf{l}}} \prod_{d \in I} \Delta_{k_d l_d}^d \prod_{d \notin I} \pi_{k_d l_d}^d. \quad (B.21)$$

Now we can furthermore observe that for any $\mathbf{k}, \mathbf{l} \in \mathcal{U}$, we have the following identity

$$\prod_{d \notin I} \pi_{k_d l_d}^d = \prod_{d \notin I} (\Delta_{k_d l_d}^d + \pi_{k_d}^d \pi_{l_d}^d) = \sum_{I' \subseteq I^c} \prod_{d \in I'} \Delta_{k_d l_d}^d \prod_{d \in I^c \setminus I'} \pi_{k_d}^d \pi_{l_d}^d \quad (B.22)$$

which can then be substituted in the previous result to give

$$E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^{I})] = \sum_{\mathbf{k}, \mathbf{l} \in \mathcal{U}} \frac{y_{\mathbf{k}} y_{\mathbf{l}}}{\pi_{\mathbf{k}} \pi_{\mathbf{l}}} \prod_{d \in I} \Delta_{k_d l_d}^d \sum_{I' \subseteq I^c} \prod_{d \in I'} \Delta_{k_d l_d}^d \prod_{d \in I^c \setminus I'} \pi_{k_d}^d \pi_{l_d}^d \quad (B.23)$$

$$= \sum_{I' \subseteq I^c} \sum_{\mathbf{k}, \mathbf{l} \in \mathcal{U}} \frac{y_{\mathbf{k}}}{\prod_{d \in I \cup I'} \pi_{k_d}^d} \frac{y_{\mathbf{l}}}{\prod_{d \in I \cup I'} \pi_{l_d}^d} \prod_{d \in I \cup I'} \Delta_{k_d l_d}^d \quad (B.24)$$

Finally we can note that $\{I \cup I', I' \subseteq I^c\} = \{I', I \subseteq I' \subseteq \{1, \dots, D\}\}$. From this observation, we can make a change of variable and conclude

$$E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^{I})] = \sum_{I' \supseteq I} \sum_{\mathbf{k}, \mathbf{l} \in \mathcal{U}} \frac{y_{\mathbf{k}}}{\prod_{d \in I'} \pi_{k_d}^d} \frac{y_{\mathbf{l}}}{\prod_{d \in I'} \pi_{l_d}^d} \prod_{d \in I'} \Delta_{k_d l_d}^d \quad (B.25)$$

$$= \sum_{I' \supseteq I} \sum_{\substack{\mathbf{k}'' , \mathbf{l}'' \in \mathcal{U}_{I'} \\ \mathbf{k}, \mathbf{l} \in \mathcal{U} \\ \forall d \in I', k_d=k_d'', l_d=l_d''}} \frac{y_{\mathbf{k}}}{\prod_{d \in I'} \pi_{k_d}^d} \frac{y_{\mathbf{l}}}{\prod_{d \in I'} \pi_{l_d}^d} \prod_{d \in I'} \Delta_{k_d l_d}^d \quad (B.26)$$

$$= \sum_{I' \supseteq I} \sum_{\mathbf{k}'' , \mathbf{l}'' \in \mathcal{U}_{I'}} \frac{Y_{\mathbf{k}''} Y_{\mathbf{l}''}}{\pi_{\mathbf{k}''} \pi_{\mathbf{l}''}} \Delta_{\mathbf{k}'' \mathbf{l}''} \quad (B.27)$$

$$= \sum_{I' \supseteq I} V_p(\hat{Y}_{\pi}^{I'}). \quad (B.28)$$
## B.4 Proof of Proposition 4

Let us fix a non-empty subset $I \subseteq \{1, \dots, D\}$. We can write using the definitions of $\hat{V}^{PLUG}(\hat{Y}_{\pi}^I)$

$$E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^I)] = \sum_{\mathbf{k}', \mathbf{l}' \in \mathcal{U}_I} \frac{E_p(\hat{Y}_{\mathbf{k}'} \hat{Y}_{\mathbf{l}'} \delta_{\mathbf{k}'} \delta_{\mathbf{l}'}) \Delta_{\mathbf{k}'\mathbf{l}'}}{\pi_{\mathbf{k}'} \pi_{\mathbf{l}'}}, \quad (B.17)$$

$$= \sum_{\mathbf{k}', \mathbf{l}' \in \mathcal{U}_I} \frac{E_p(\hat{Y}_{\mathbf{k}'} \hat{Y}_{\mathbf{l}'}) E_p(\delta_{\mathbf{k}'} \delta_{\mathbf{l}'}) \Delta_{\mathbf{k}'\mathbf{l}'}}{\pi_{\mathbf{k}'} \pi_{\mathbf{l}'} \pi_{\mathbf{k}'} \pi_{\mathbf{l}'}} \quad (B.18)$$

$$= \sum_{\mathbf{k}', \mathbf{l}' \in \mathcal{U}_I} \frac{E_p(\hat{Y}_{\mathbf{k}'} \hat{Y}_{\mathbf{l}'}) \Delta_{\mathbf{k}'\mathbf{l}'}}{\pi_{\mathbf{k}'} \pi_{\mathbf{l}'}}, \quad (B.19)$$

Where the line (B.18) follows from the fact that $\delta_{\mathbf{k}'}$ and $\delta_{\mathbf{l}'}$ are functions of $(S_d)_{d \in I}$ whereas $\hat{Y}_{\mathbf{k}'}$ and $\hat{Y}_{\mathbf{l}'}$ are functions of $(S_d)_{d \notin I}$, and are therefore independent. We can then similarly expand the sum by replacing $\hat{Y}_{\mathbf{k}'}$ and $\hat{Y}_{\mathbf{l}'}$ by their definitions in (3.25):

$$E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^I)] = \sum_{\mathbf{k}', \mathbf{l}' \in \mathcal{U}_I} \sum_{\substack{\mathbf{k}, \mathbf{l} \in \mathcal{U} \\ \forall d \in I, k_d=k_{d'}, l_d=l_{d}'}} \frac{y_{\mathbf{k}} y_{\mathbf{l}}}{\pi_{\mathbf{k}} \pi_{\mathbf{l}}} \prod_{d \in I} \Delta_{k_d l_d}^d E_p \left( \prod_{d \notin I} \delta_{k_d}^d \delta_{l_d}^d \right) \quad (B.20)$$

$$= \sum_{\mathbf{k}, \mathbf{l} \in \mathcal{U}} \frac{y_{\mathbf{k}} y_{\mathbf{l}}}{\pi_{\mathbf{k}} \pi_{\mathbf{l}}} \prod_{d \in I} \Delta_{k_d l_d}^d \prod_{d \notin I} \pi_{k_d l_d}^d. \quad (B.21)$$

Now we can furthermore observe that for any $\mathbf{k}, \mathbf{l} \in \mathcal{U}$, we have the following identity

$$\prod_{d \notin I} \pi_{k_d l_d}^d = \prod_{d \notin I} (\Delta_{k_d l_d}^d + \pi_{k_d}^d \pi_{l_d}^d) = \sum_{I' \subseteq I^c} \prod_{d \in I'} \Delta_{k_d l_d}^d \prod_{d \in I^c \setminus I'} \pi_{k_d}^d \pi_{l_d}^d \quad (B.22)$$

which can then be substituted in the previous result to give

$$E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^I)] = \sum_{\mathbf{k}, \mathbf{l} \in \mathcal{U}} \frac{y_{\mathbf{k}} y_{\mathbf{l}}}{\pi_{\mathbf{k}} \pi_{\mathbf{l}}} \prod_{d \in I} \Delta_{k_d l_d}^d \sum_{I' \subseteq I^c} \prod_{d \in I'} \Delta_{k_d l_d}^d \prod_{d \in I^c \setminus I'} \pi_{k_d}^d \pi_{l_d}^d \quad (B.23)$$

$$= \sum_{I' \subseteq I^c} \sum_{\mathbf{k}, \mathbf{l} \in \mathcal{U}} \frac{y_{\mathbf{k}}}{\prod_{d \in I \cup I'} \pi_{k_d}^d} \frac{y_{\mathbf{l}}}{\prod_{d \in I \cup I'} \pi_{l_d}^d} \prod_{d \in I \cup I'} \Delta_{k_d l_d}^d \quad (B.24)$$

Finally we can note that $\{I \cup I', I' \subseteq I^c\} = \{I', I \subseteq I' \subseteq \{1, \dots, D\}\}$. From this observation, we can make a change of variable and conclude

$$E_p[\hat{V}^{PLUG}(\hat{Y}_{\pi}^I)] = \sum_{I' \supseteq I} \sum_{\mathbf{k}, \mathbf{l} \in \mathcal{U}} \frac{y_{\mathbf{k}}}{\prod_{d \in I'} \pi_{k_d}^d} \frac{y_{\mathbf{l}}}{\prod_{d \in I'} \pi_{l_d}^d} \prod_{d \in I'} \Delta_{k_d l_d}^d \quad (B.25)$$

$$= \sum_{I' \supseteq I} \sum_{\substack{\mathbf{k}'' , \mathbf{l}'' \in \mathcal{U}_{I'} \\ \mathbf{k}, \mathbf{l} \in \mathcal{U} \\ \forall d \in I', k_d=k_{d}'' , l_d=l_{d}''}} \frac{y_{\mathbf{k}}}{\prod_{d \in I'} \pi_{k_d}^d} \frac{y_{\mathbf{l}}}{\prod_{d \in I'} \pi_{l_d}^d} \prod_{d \in I'} \Delta_{k_d l_d}^d \quad (B.26)$$

$$= \sum_{I' \supseteq I} \sum_{\mathbf{k}'' , \mathbf{l}'' \in \mathcal{U}_{I'}} \frac{Y_{\mathbf{k}''}}{\pi_{\mathbf{k}''}} \frac{Y_{\mathbf{l}''}}{\pi_{\mathbf{l}''}} \Delta_{\mathbf{k}'' \mathbf{l}''} \quad (B.27)$$

$$= \sum_{I' \supseteq I} V_p(\hat{Y}_{\pi}^{I'}). \quad (B.28)$$
References

Beaumont, J.-F. and Patak, Z. (2012). On the generalized bootstrap for sample surveys with special attention to poisson sampling. *International Statistical Review*, 80:127–148.

Booth, J. G., Butler, R. W., and Hall, P. (1994). Bootstrap methods for finite populations. *Journal of the American Statistical Association*, 89(428):1282–1289.

Charles, M. A., Thierry, X., Lanoe, J.-L., Bois, C., Dufourg, M.-N., Popa, R., Cheminat, M., Zaros, C., and Geay, B. (2020). Cohort profile: the french national cohort of children (elfe): birth to 5 years. *International journal of epidemiology*, 49(2):368–369j.

Chauvet, G. and Vallée, A.-A. (2020). Inference for two-stage sampling designs. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 82(3):797–815.

Cochran, W. G. (1977). *Sampling techniques*. John Wiley & Sons, New York-London-Sydney, third edition. Wiley Series in Probability and Mathematical Statistics.

Dalén, J. and Ohlsson, E. (1995). Variance estimation in the swedish consumer price index. *Journal of Business & Economic Statistics*, 13(3):347–356.

Fuller, W. A. (2011). *Sampling statistics*, volume 560. John Wiley & Sons.

Gross, S. (1980). Median estimation in sample surveys. *Proceedings of the Section on Survey Research Methods*, pages 181–184.

Hájek, J. (1964). Asymptotic theory of rejective sampling with varying probabilities from a finite population. *The Annals of Mathematical Statistics*, pages 1491–1523.

Hoeffding, W. (1948). A Class of Statistics with Asymptotically Normal Distribution. *The Annals of Mathematical Statistics*, 19(3):293 – 325.

Juillard, H., Chauvet, G., and Ruiz-Gazen, A. (2017). Estimation under cross-classified sampling with application to a childhood survey. *Journal of the American Statistical Association*, 112(518):850–858.

Ohlsson, E. (1996). Cross-classified sampling. *Journal of Official Statistics*, 12(3):241.
References

Beaumont, J.-F. and Patak, Z. (2012). On the generalized bootstrap for sample surveys with special attention to poisson sampling. *International Statistical Review*, 80:127–148.

Booth, J. G., Butler, R. W., and Hall, P. (1994). Bootstrap methods for finite populations. *Journal of the American Statistical Association*, 89(428):1282–1289.

Charles, M. A., Thierry, X., Lanoe, J.-L., Bois, C., Dufourg, M.-N., Popa, R., Cheminat, M., Zaros, C., and Geay, B. (2020). Cohort profile: the french national cohort of children (elfe): birth to 5 years. *International journal of epidemiology*, 49(2):368–369j.

Chauvet, G. and Vallée, A.-A. (2020). Inference for two-stage sampling designs. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 82(3):797–815.

Cochran, W. G. (1977). *Sampling techniques*. John Wiley & Sons, New York-London-Sydney, third edition. Wiley Series in Probability and Mathematical Statistics.

Dalén, J. and Ohlsson, E. (1995). Variance estimation in the swedish consumer price index. *Journal of Business & Economic Statistics*, 13(3):347–356.

Fuller, W. A. (2011). *Sampling statistics*, volume 560. John Wiley & Sons.

Gross, S. (1980). Median estimation in sample surveys. *Proceedings of the Section on Survey Research Methods*, pages 181–184.

Hájek, J. (1964). Asymptotic theory of rejective sampling with varying probabilities from a finite population. *The Annals of Mathematical Statistics*, pages 1491–1523.

Hoeffding, W. (1948). A Class of Statistics with Asymptotically Normal Distribution. *The Annals of Mathematical Statistics*, 19(3):293 – 325.

Juillard, H., Chauvet, G., and Ruiz-Gazen, A. (2017). Estimation under cross-classified sampling with application to a childhood survey. *Journal of the American Statistical Association*, 112(518):850–858.

Ohlsson, E. (1996). Cross-classified sampling. *Journal of Official Statistics*, 12(3):241.
Pirus, C., Bois, C., Dufourg, M.-N., Lanoë, J.-L., Vandentorren, S., Leridon, H., Depledge, R., and the Elfe team (2010). Constructing a cohort: Experience with the french elfe project. *Population*, 65:637–670.

Rao, J., Wu, C., and Yue, K. (1992). Some recent work on resampling methods for complex surveys. *Survey methodology*, 18(2):209–217.

Särndal, C.-E., Swensson, B., and Wretman, J. (1992). *Model Assisted Survey Sampling*. Springer New York.

Skinner, C. (2015). Cross-classified sampling: some estimation theory. *Statistics & Probability Letters*, 104:163–168.
Pirus, C., Bois, C., Dufourg, M.-N., Lanoë, J.-L., Vandentorren, S., Leridon, H., Depledge, R., and the Elfe team (2010). Constructing a cohort: Experience with the french elfe project. *Population*, 65:637-670.

Rao, J., Wu, C., and Yue, K. (1992). Some recent work on resampling methods for complex surveys. *Survey methodology*, 18(2):209-217.

Särndal, C.-E., Swensson, B., and Wretman, J. (1992). *Model Assisted Survey Sampling*. Springer New York.

Skinner, C. (2015). Cross-classified sampling: some estimation theory. *Statistics & Probability Letters*, 104:163-168.