Let us show that for the optimal weights derived in Proposition 2.2.1, we have :

$$ \mathrm{Var}(\hat{t}_{\tilde{\theta}_k}) - \mathrm{Var}(\hat{t}_{\tilde{\theta}_k^{\mathrm{opt}}}) = \mathrm{Var}(\hat{t}_{\tilde{\theta}_k}) - \hat{t}_{\tilde{\theta}_k^{\mathrm{opt}}}). $$

The optimal set of weights is given by $\tilde{\theta}_k^{\mathrm{opt}} = \Delta_k^{-1}1_k (1_k^t \Delta_k^{-1}1_k)^{-1}$, $k \in U_T$.

$$ \begin{aligned} \mathrm{Var}(\hat{t}_{\tilde{\theta}_k}) - \mathrm{Var}(\hat{t}_{\tilde{\theta}_k^{\mathrm{opt}}}) &= \tilde{\theta}_k^t \Delta_k \tilde{\theta}_k - (\tilde{\theta}_k^{\mathrm{opt}})^t \Delta_k \tilde{\theta}_k^{\mathrm{opt}} \\ &= (\tilde{\theta}_k - \tilde{\theta}_k^{\mathrm{opt}})^t \Delta_k (\tilde{\theta}_k - \tilde{\theta}_k^{\mathrm{opt}}) + \tilde{\theta}_k^t \Delta_k \tilde{\theta}_k^{\mathrm{opt}} + (\tilde{\theta}_k^{\mathrm{opt}})^t \Delta_k \tilde{\theta}_k \\ &\quad - 2(\tilde{\theta}_k^{\mathrm{opt}})^t \Delta_k \tilde{\theta}_k^{\mathrm{opt}}. \end{aligned} $$

Now, $\tilde{\theta}_k^t \Delta_k \tilde{\theta}_k^{\mathrm{opt}} = (\tilde{\theta}_k^{\mathrm{opt}})^t \Delta_k \tilde{\theta}_k$ since they are real quantities. Straightforward calculations using the standardization constraint (2.14), give us that $\tilde{\theta}_k^t \Delta_k \tilde{\theta}_k^{\mathrm{opt}} = (1_k^t \Delta_k^{-1}1_k)^{-1}$. Moreover, $(\tilde{\theta}_k^{\mathrm{opt}})^t \Delta_k \tilde{\theta}_k^{\mathrm{opt}} = (1_k^t \Delta_k^{-1}1_k)^{-1}$, and we finally get that

$$ \begin{aligned} \mathrm{Var}(\hat{t}_{\tilde{\theta}_k}) - \mathrm{Var}(\hat{t}_{\tilde{\theta}_k^{\mathrm{opt}}}) &= (\tilde{\theta}_k - \tilde{\theta}_k^{\mathrm{opt}})^t \Delta_k (\tilde{\theta}_k - \tilde{\theta}_k^{\mathrm{opt}}) \\ &= \mathrm{Var}(\hat{t}_{\tilde{\theta}_k} - \hat{t}_{\tilde{\theta}_k^{\mathrm{opt}}}) = \mathrm{Var}(\hat{t}_{\tilde{\theta}_k} - \tilde{\theta}_k^{\mathrm{opt}}) \end{aligned} $$

which ends the proof of Proposition 2.2.2. □

### Proof of Proposition 2.2.3

For all $i \in U_F$, we recall that we assume $0 < \pi_i < 1$, and $\pi_k = 1 - \Pi_{i \in U_F} (1 - \pi_i)^{l_{ik}}$. For Poisson sampling,

$$ \tilde{\theta}_{ik}^{\mathrm{opt}} = \frac{l_{ik} \frac{\pi_i}{1-\pi_i}}{\sum_{i' \in U_F} l_{i'k} \frac{\pi_{i'}}{1-\pi_{i'}}}. $$

Using the variance expressions for Poisson sampling, we have that proving $\mathrm{Var}(\hat{t}_y) < \mathrm{Var}(\hat{t}_{y1}^{\mathrm{opt}})$ is equivalent to prove

$$ \sum_{k \in U_T} y_k^2 \frac{1 - \pi_k}{\pi_k} < \sum_{k \in U_T} y_k^2 \sum_{i \in U_F} \frac{1 - \pi_i}{\pi_i} (\tilde{\theta}_{ik}^{\mathrm{opt}})^2. $$

This inequality is true if, for all $k \in U_T$, we have :

$$ \frac{1 - \pi_k}{\pi_k} < \sum_{i \in U_F} \frac{1 - \pi_i}{\pi_i} (\tilde{\theta}_{ik}^{\mathrm{opt}})^2 $$
which is true as proved next.

$$ 
\begin{align*}
\frac{1 - \pi_k}{\pi_k} < \sum_{i \in U_F} \frac{1 - \pi_i}{\pi_i} (\tilde{\theta}_{ik}^{opt})^2 &\Leftrightarrow \frac{\prod_{i \in U_F} (1 - \pi_i)^{l_{ik}}}{1 - \prod_{i \in U_F} (1 - \pi_i)^{l_{ik}}} < \sum_{i \in U_F} \frac{1 - \pi_i}{\pi_i} (\tilde{\theta}_{ik}^{opt})^2 \\
&\Leftrightarrow \sum_{i \in U_F} l_{ik} \frac{\pi_i}{1 - \pi_i} < \frac{1 - \prod_{i \in U_F} (1 - \pi_i)^{l_{ik}}}{\prod_{i \in U_F} (1 - \pi_i)^{l_{ik}}} \\
&\Leftrightarrow 1 + \sum_{i \in U_F} l_{ik} \frac{\pi_i}{1 - \pi_i} < \frac{1}{\prod_{i \in U_F} (1 - \pi_i)^{l_{ik}}} \\
&\Leftrightarrow 1 + \sum_{i \in U_F} l_{ik} \frac{\pi_i}{1 - \pi_i} < \prod_{i \in U_F} \left(1 + \frac{\pi_i}{1 - \pi_i}\right)^{l_{ik}}.
\end{align*}
$$

This last inequality is true since $\frac{\pi_i}{1 - \pi_i} > 0$, and this finishes the proof of Proposition 2.2.3. □
# Chapitre 3

## Optimal Weights for double
Many-To-One Generalized Weight
Share Method

### Abstract

In probabilistic surveys, the sampling frame of the target population is not always available. If there is a frame population linked to the target population, a sample can be drawn in the target population through indirect sampling. The sampling weights of the target population units can be determined using the *Generalized Weight Share Method* (*GWSM*). These sampling weights take into account the sampling weights of the frame population units, together with some link weights that are assigned by the statistician to the links between the frame and the target populations. This chapter focuses on “many-to-one” links where it is possible to derive a set of optimal link weights that minimize the variance of the *GWSM* estimator, under constraints on the sampling design which are not very restrictive (Medous et al., 2023a). To get unbiased estimators of finite population totals, the link weights have to be standardized, which implies that the links between the frame population and the sample units in the target population have to be known or observed. If the number of links is large, the link weights are thus difficult to retrieve. A solution is to consider an intermediate population, linked to both the frame and the target populations, and to use a double indirect sampling. The sampling weights can be determined using a double *GWSM*, first between the frame and the intermediate populations, and then between the intermediate and the target populations. When both link weights are standardized (double standardization), the double *GWSM* allows for a reduction in the number of observed links. This method is therefore easier to apply than the *GWSM*, but it may deteriorate the precision of the estimator. The objective of the present paper is to derive doubly standardized link
weights that are optimal in the sense that the variance of the double GWSM is minimized. When such weights cannot be computed, as is the case in the French postal traffic survey, alternative weights can be used to improve the precision of the double GWSM estimator. Different alternative weights are compared through Monte Carlo simulations and an application to the French postal traffic estimation is discussed.

**Keywords**: complex sampling design, finite population, generalized weight share method, indirect sampling, surveys.

## 3.1 Introduction

For economic and marketing purposes, the French postal service (La Poste) is interested in obtaining information not only on the monthly postal traffic but also on the types of letters being sent. As certain letters cannot be processed automatically due to either confidentiality concerns or challenges in identifying certain types of letters, La Poste conducts a probability-based survey of postman rounds to estimate the monthly postal traffic. However, due to frequent reorganizations of postal offices, the population of postman rounds is too unstable to be directly observed and La Poste uses the more stable population of addresses as a proxy to indirectly draw samples from the population of rounds, as detailed in Medous et al. (2023a).

Indirect sampling consists in drawing a sample from a frame population linked to the target population and observing all target units linked to this sample. This sampling method requires that all units in the target population are linked to at least one unit of the frame population. In La Poste example, the target population comprises postman rounds departing in a particular month, while the frame population consists of all pairs of French addresses and days of the month, referred to as “address-days.” A round is considered linked to an address-day if the postman delivers letters to that address on that specific day. In this example, all units in the target population are linked to the frame population since each postman round serves at least one address on a particular day. Furthermore, a single address can only be served by one round at most on a given day, making the links of the Many-to-One (MtO) type.

Indirect sampling raises a challenge in computing sampling weights of the target population units, but this can be addressed by using the Generalized Weight Share Method (GWSM). GWSM was introduced by Deville and Lavallée (2006) and Lavallée (2007) and was inspired by the weight sharing method used by Ernst (1986) and Kalton and Brick (1995) for longitudinal household studies. It consists of allocating weights to the links bet-
between frame and target populations, which are then used to compute the final weights. GWSM has been widely studied in the literature, particularly for its utility in sampling hard-to-reach populations, as demonstrated in Deville and Lavallée (2006) and De Vitiis et al. (2014).

Deville and Lavallée (2006) discuss the existence of link weights minimizing the variance of the GWSM estimator for any variable of interest and demonstrate that such optimal link weights exist for Simple Random Sampling Without Replacement and Poisson sampling, if the links are of type MtO. Medous et al. (2023a) extend the results of Deville and Lavallée (2006) to a more general setting and show that optimal link weights exist for the complex sampling design used by La Poste.

However, in MtO setups, a large number of links has to be observed in order to implement the GWSM. In la Poste case, GWSM requires to identify all addresses per round, which can be relatively high. To reduce this number, La Poste has adopted the strategy of drawing double indirect samples using the intermediate population of sorting boxes and double GWSM estimators of finite population totals. According to Medous et al. (2023a), employing double GWSM for MtO-MtO links (MtO links between frame and intermediate and between intermediate and target populations) can significantly reduce the required number of observed links when compared to a single GWSM. However, this reduction comes at the expense of imposing constraints on the link weights, which requires to reformulate the optimal weights described by Deville and Lavallée (2006) and Medous et al. (2023a).

This work intends to investigate the existence of optimal weights that minimize the variance of double GWSM estimators regardless of the variable of interest while respecting the weights constraints defined by Medous et al. (2023a) to ensure a reduction in the number of observed links compared to simple GWSM estimators. We show that, in the case of La Poste, optimal weights for the double GWSM exist but cannot be computed, as the required information is difficult to obtain. Consequently, the current double GWSM estimator used by La Poste suffers from a significant loss of efficiency compared to an optimal double GWSM estimator. To limit this loss of efficiency, we suggest the use of alternative weights computed using auxiliary information.

Section 3.2 provides a brief overview of simple and double indirect sampling and their corresponding simple and double GWSM estimators, along with their variances. In section 3.3, we derive optimal double GWSM weights for MtO-MtO links. We discuss the conditions for the existence of such weights and provide their values for common sampling designs, such as Simple Random Sampling without Replacement or Poisson sampling. In section 3.4, we address the situation where optimal double GWSM weights cannot be computed and presents alternative weight options, including those currently used at La Poste. We examine the loss of precision between an optimal double GWSM estimator and a double GWSM estimator
<figure><img src="image_43.png" /><figcaption>FIGURE 3.1 – Example of MtO links at La Poste.</figcaption></figure>

using alternative weights, and provide advice on choosing the best alternative weights. In Section 3.5, we report results of a simulation study that mimics La Poste situation, where we test various link weights possibilities. Finally, we provide details on the alternative weights currently being tested at La Poste and those planned for the future. Proofs are gathered in the Appendix.

## 3.2 Double Generalized Weight Share Method

### 3.2.1 Indirect sampling and simple GWSM estimator

La Poste estimates its monthly traffic by observing postman rounds. However, the organization of post offices makes it difficult to know which rounds will depart in advance. Direct sampling methods cannot be used to draw a sample of rounds and La Poste has to use alternative sampling methods such as indirect sampling. Indirect sampling involves selecting a sample from a frame population that is linked to the target population.

Let $U_T$ be the target population of size $N_T$ and $U_F$ be a sampling frame of size $N_F$ linked to $U_T$. Following Deville and Lavallée (2006), we assume that every link between $U_F$ and $U_T$ is observable and that each unit of $U_T$ is linked to at least one unit of $U_F$.

Various types of links between populations exist and have been studied in the literature; see Deville (1999), Deville and Lavallée (2006), Medous et al. (2023a), among others. In this paper, we concentrate on a specific type of links, represented in figure 3.1, where each unit in $U_F$ is connected to at most one unit in $U_T$. This type of links is known as the “Many-to-One” (MtO) link.

We denote by $y_k$ the value of the variable of interest $y$ for unit $k$ from $U_T$ and focus on the estimation of the total $t_y = \sum_{k \in U_T} y_k$ over $U_T$. A sample $s_F$ is drawn from $U_F$ using a sampling design $p_F(\cdot)$. The indirect sample $s_T$ from $U_T$ is composed of units from $U_T$ linked
to at least a unit in $s_F$. We denote by $\pi_i = \mathbb{P}(i \in s_F)$ the first-order inclusion probability of unit $i \in U_F$ and by $\pi_{ii'} = \mathbb{P}(i, i' \in s_F)$ the second-order inclusion probability of units $i$ and $i' \in U_F$. We assume that $\pi_i > 0$ for all units $i$ in $U_F$ and note $d_i = 1/\pi_i > 0$ the sampling weight of unit $i$. We denote $\Delta_{ii'} = \frac{\pi_{ii'} - \pi_i\pi_{i'}}{\pi_i\pi_{i'}}$, $i, i' \in U_F$.

Three sampling designs $p_F(\cdot)$ are considered below: Poisson sampling, Simple Random Sampling Without Replacement (SRSWOR) of size $n_F$ and Stratified Simple Random Sampling Without Replacement (SSRSWOR) with $H$ strata of size $N_h$, $h = 1...H$. Poisson sampling consists in performing an independent Bernoulli trial of parameter $\pi_i$ for each $i \in U_F$. SRSWOR of size $n_F$ consists in assigning equal and positive probabilities to each of the without replacement samples of size $n_F$ and a null probability otherwise. The sampling weights are equal to $d_i = N_F/n_F$ for all $i$ in $U_F$. Finally, SSRSWOR with $H$ strata of size $N_h$, $h = 1...H$, consists in performing an independent SRSWOR of size $n_h$ in each stratum $h$, $h = 1...H$. The sampling weights are equal to $d_i = N_h/n_h$ for all $i$ in strata $h$, $h = 1...H$.

As mentioned in Deville and Lavallée (2006), the main difficulty with indirect sampling is the computation of the sampling weights of units $k \in U_T$. In order to deal with this challenge, we assign a link weight $\theta_{ik}$ to each link between $i \in U_F$ and $k \in U_T$, with $\theta_{ik} > 0$ when $i$ and $k$ are linked and $\theta_{ik} = 0$ otherwise. We define the standardized link weights by $\tilde{\theta}_{ik} = \theta_{ik} / \sum_{i' \in U_F} \theta_{i'k}$ for $i \in U_F$, $k \in U_T$, such that $\sum_{i \in U_F} \tilde{\theta}_{ik} = 1$ for all $k \in U_T$. The simple Generalized Weight Share Method (GWSM) estimator of $t_y$, introduced by Deville and Lavallée (2006), is given by:

$$ \hat{t}_{y1} = \sum_{i \in s_F} d_i \left( \sum_{k \in U_T} \tilde{\theta}_{ik} y_k \right) \quad (3.1) $$

and its variance by:

$$ \mathrm{Var}(\hat{t}_{y1}) = \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \sum_{k \in U_T} \tilde{\theta}_{ik} y_k \sum_{k' \in U_T} \tilde{\theta}_{i'k'} y_{k'} \quad (3.2) $$

It can be shown that the simple GWSM estimator is unbiased as long as the link weights are standardized. However, as we can see in equation (3.2), the variance of the simple GWSM estimator depends on the standardized link weights. Deville and Lavallée (2006) and Lavallée (2007) seek for optimal standardized link weights, namely link weights minimizing the variance given by equation (3.2) for any variable of interest $y$. When the links are of type MtO and the samples are drawn using SRSWOR or Poisson sampling, Deville and Lavallée (2006) and Lavallée (2007) demonstrate the existence and uniqueness of a set of optimal standardized link weights. Medous et al. (2023a) extend this result to a wider range of sam-
<figure><img src="image_44.png" /><figcaption>FIGURE 3.2 – Example of MtO-MtO at La Poste.</figcaption></figure>

pling designs satisfying a specific condition. In the following, we refer to the optimal simple GWSM estimator as the simple GWSM estimator obtained with the optimal standardized link weights.

MtO links allow optimal weights to be derived, but the large number of observed links required to compute the standardized link weights, whether they are optimal or not, can raise implementation challenges, as it is the case of La Poste (see Medous et al. (2023a)). To deal with this issue, one can use double indirect sampling as suggested in Deville and Lavallée (2006) and Medous et al. (2023a).

### 3.2.2 Double GWSM estimator

Double indirect sampling consists of introducing a second intermediate population $U_M$, linked to both the frame population $U_F$ and the target population $U_T$. Using the links between these populations, two indirect samples can be obtained, one from $U_M$ and the other from $U_T$. The estimators of totals on $U_T$ are then obtained by applying the GWSM twice, first between $U_F$ and $U_M$, and then between $U_M$ and $U_T$. This method is known as double GWSM (Medous et al., 2023a).

Let $N_M$ be the size of the intermediate population $U_M$. Following Medous et al. (2023a), we make the assumption that the links between $U_F$ and $U_M$, respectively $U_M$ and $U_T$, are both of type MtO. This type of links is referred to as MtO-MtO, which is a specific case of MtO.

As in section 3.2.1, a sample $s_F$ is drawn from $U_F$ with a sampling design $p_F(\cdot)$. The indirect sample $s_M$ from $U_M$ (resp. $s_T$ from $U_T$) is composed of units from $U_M$ linked to at least a unit $i \in s_F$ (resp. of units from $U_T$ linked to at least a unit $j \in s_M$). This sampling technique is called double indirect sampling and is currently used in La Poste survey, with an intermediate population of sorting boxes. These boxes are used to sort the mail before the beginning of a postman round.

For each $i \in U_F$ and $j \in U_M$ (resp. $j \in U_M$, $k \in U_T$), we assign a link weight $\theta_{ij}^{FM}$ (resp.
$\theta_{jk}^{MT}$) with $\theta_{ij}^{FM} > 0$ when i and j are linked and $\theta_{ij}^{FM} = 0$ otherwise (resp. $\theta_{jk}^{MT} > 0$ when j and k are linked and $\theta_{jk}^{MT} = 0$ otherwise). Let $\tilde{\theta}_{ij}^{FM}$ (resp. $\tilde{\theta}_{jk}^{MT}$) be the standardized link weight between $i \in U_F$ and $j \in U_M$ (resp. $j \in U_M$ and $k \in U_T$) such that $\sum_{i \in U_F} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} = 1$. In what follows, we call FM weights (resp. MT weights) the set of standardized weights between $U_F$ and $U_M$ (resp. between $U_M$ and $U_T$).

The double GWSM estimator of $t_y$ is given by

$$ \hat{t}_{y2} = \sum_{i \in s_F} d_i \left( \sum_{k \in U_T} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} y_k \right). \quad (3.3) $$

The double GWSM estimator is a simple GWSM with standardized link weights $\tilde{\theta}_{ik}$, $i \in U_F$, $k \in U_T$, such that

$$ \tilde{\theta}_{ik} = \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT}. $$

The condition $\sum_{i \in U_F} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} = 1$, $k \in U_T$, produces unbiased double GWSM estimators and is verified for various standardization methods for FM and MT link weights. However, not all standardization methods allow for a reduction of the number of observed links. Medous et al. (2023a) show that a reduction of the number of observed links occurs when we consider FM and MT weights such that $\sum_{i \in U_F} \tilde{\theta}_{ij}^{FM} = 1$, $j \in U_M$, and $\sum_{j \in U_M} \tilde{\theta}_{jk}^{MT} = 1$, $j \in U_T$. This is known as double standardization and we say that the link weights are twice (or double) standardized.

The double GWSM with double standardization offers an advantage when it comes to the number of links to observe. Yet, the optimal link weights defined for simple GWSM have to be rewritten to ensure that they satisfy the double standardization constraint. Without this reformulation of the optimal simple weights, the variance of the double GWSM estimator with double standardization can be higher than the variance of the optimal simple GWSM estimator, as shown by Medous et al. (2023a). In order to get an estimator with both a small number of links to observe and a small variance, we look for twice standardized link weights that minimize the variance of the double GWSM estimator for all variable of interest $y$. Thus we look for a set of FM weights and MT weights such that $\sum_{i \in U_F} \tilde{\theta}_{ij}^{FM} = 1$, $j \in U_M$, and $\sum_{j \in U_M} \tilde{\theta}_{jk}^{MT} = 1$, $j \in U_T$, that minimizes the variance given by equation (3.2) for any variable of interest, with $\tilde{\theta}_{ik} = \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT}$. In what follows, we will refer to those weights as the optimal weights for double GWSM estimators.

In what follows, we only consider twice standardized FM and MT weights.
<figure><img src="image_45.png" /><figcaption>FIGURE 3.3 – $U_F$ ordered according to $U_{Fj}$, $j \in U_M$.</figcaption></figure>

## 3.3 Optimality of double GWSM

Medous et al. (2023a) show that optimal weights for the simple GWSM estimator exist for sampling designs satisfying a condition that depends on the links between $U_F$ and $U_T$. Likewise, we can show that optimal weights for double GWSM estimator with double standardization exist for sampling designs that meet two conditions: one depending on the links between $U_F$ and $U_M$, and another depending on the links between $U_M$ and $U_T$.

Let us consider first the links between $U_F$ and $U_M$. For a given unit $j \in U_M$, we denote by $U_{Fj}$ the set of size $N_{Fj}$ containing all units $i$ in $U_F$ that are linked to $j$. If the links between $U_F$ and $U_M$ are MtO, these subpopulations form a partition of $U_F$. In what follows, we assume that units in $U_F$ are ordered according to the subpopulations $U_{Fj}$, $j \in U_M$, as in figure 3.3. Let $\Delta = (\Delta_{ii'})_{i,i' \in U_F}$ be the covariance matrix of size $N_F \times N_F$ and $\Delta_{jj'}^M = (\Delta_{ii'})_{i \in U_{Fj}, i' \in U_{Fj'}}$ be the submatrix of $\Delta$ corresponding to the elements in position $i$ and $i'$ such that $i$ (resp. $i'$) is linked to $j$ (resp. to $j'$). We denote by $\Delta_j^M$ the square submatrix $\Delta_{jj'}^M$ of size $N_{Fj} \times N_{Fj}$. If the links between $U_F$ and $U_M$ are MtO, then the submatrices $\Delta_{jj'}^M$ form a partition of $\Delta$ as in figure 3.3. Let us denote $\mathbf{1}_j^M$ the $N_{Fj}$-dimensional vector of ones and by $u^t$ the transpose of any vector $u$.

*Definition 3.3.1* ($\Delta$-property for $U_M$). For MtO links between $U_F$ and $U_M$, a sampling design is said to satisfy the $\Delta$-property for $U_M$ if, for all $j \in U_M$, $\Delta_j^M$ is invertible and, for $j \neq j' \in U_M$, we have

$$ \Delta_{j,j' \neq j} = c_{jj'} \mathbf{1}_j^M \mathbf{1}_{j'}^{Mt} \quad \text{with } c_{jj'} \text{ not depending on } i \text{ and } i'. \qquad (3.4) $$

Following Medous et al. (2023a), we can show that the $\Delta$-property for $U_M$ holds for
Poisson sampling and SRSWOR with the following constants $c_{jj'}$, $j, j' \in U_M$, $j \neq j'$. For Poisson sampling from $U_F$ with inclusion probabilities $\pi_i$ for $i \in U_F$,

$$c_{jj'} = 0 \text{ for all } j \neq j' \in U_T. \qquad (3.5)$$

For SRSWOR of size $n_F$ from $U_F$,

$$c_{jj'} = -\frac{1-f}{f} \frac{1}{N_F-1} \qquad (3.6)$$

For stratified SRSWOR with H strata of size $N_h$ in $U_F$, $h = 1, \dots, H$, the $\Delta$-property for $U_M$ holds if the following condition is verified :

(C1) for all $j \in U_M$, all units $i \in U_F$ linked to $j$ belong to the same stratum and for each stratum $h$, there are at least two units in $U_M$ linked to $h$.

Then

$$c_{jj'} = \frac{1 - f_h}{f_h} \frac{1}{N_h - 1}, \qquad (3.7)$$

if all units $i, i' \in U_F$ linked to respectively $j$ and $j'$ belong to the same stratum and 0 otherwise, for all $j \neq j' \in U_M$.

Let us denote $\lambda_j = 1_j^{Mt}(\Delta_j^M)^{-1}1_j^M$, $j \in U_M$. The $\Delta$-property for $U_M$ allows us to derive the following proposition on the FM weights.

**Proposition 3.3.1.** If the links between $U_F$ and $U_M$ are MtO and the sampling design satisfies the $\Delta$-property for $U_M$, then there is a unique set of optimal FM weights $\tilde{\theta}_{ij}^{FM, opt}$, $i \in U_F$, $j \in U_M$ minimizing the variance of $\hat{t}_{y2}$ for any variable of interest and any set of weights $\tilde{\theta}_{jk}^{MT}$, $j \in U_M$, $k \in U_T$. These weights are given by :

$$(\tilde{\theta}_{ij}^{FM, opt})_{i \in U_{Fj}} = \frac{(\Delta_j^M)^{-1}1_j^M}{1_j^{Mt}(\Delta_j^M)^{-1}1_j^M} = \frac{1}{\lambda_j}(\Delta_j^M)^{-1}1_j^M, \text{ if } i \text{ is linked to } j, \text{ and } 0 \text{ otherwise.}$$

Interestingly, the FM weights given by Proposition 3.3.1 do not depend on the MT weights used. Let us now consider the links between $U_M$ and $U_T$. For a given unit $k \in U_T$, let us denote by $U_{Mk}$ the set of size $N_{Mk}$ containing all units $j$ in $U_M$ that are linked to $k$. In the MtO-MtO setup, a unit $j$ from $U_M$ can only be linked to one and only one unit $k$ from $U_T$ and the subpopulations $U_{Mk}$, $k \in U_T$, form a partition of $U_M$ (see figure 3.4). In what follows we assume that the units in $U_M$ are ordered according to the subpopulations $U_{Mk}$, $k \in U_T$, as in figure 3.4.

We can define the following condition.

**Definition 3.3.2 (c_k-condition).** If the links are MtO-MtO and the $\Delta$-property holds for $U_M$,
<figure><img src="image_46.png" /><figcaption>FIGURE 3.4 – $U_F$ ordered according to $U_{Fj}$, $j \in U_M$ and $U_{Mk}$, $k \in U_T$, with $\Delta_{kk}^T$ the submatrix of $\Delta$ corresponding to the element in position $i$ and $i'$ such that $i$ and $i'$ are linked to $k$.</figcaption></figure>

$\Delta$ is said to verify the $c_k$-condition if, for all $j$ (resp. $j') \in U_M$ linked to $k$ (resp. $k') \in U_T$ such that $j \neq j'$, the constants $c_{jj'}$ given by the $\Delta$-property for $U_M$ verify

$$c_{jj'} = c_{kk'}$$

with $c_{kk'}$ not depending on $j$ and $j'$.

Let us denote $c_k$, $k \in U_T$ the value of $c_{kk'}$ for $k' = k$. The $c_k$-condition holds for Poisson sampling and SRSWOR. For Poisson sampling and SRSWOR, the values of $c_{jj'}$, $j \neq j' \in U_M$ given by the $\Delta$-property for $U_M$ are constant and do not depend on $j$ or $j'$. Thus the $c_k$-property holds and the values $c_{kk'}$, $k \neq k' \in U_T$, are equal to the values $c_{jj'}$, $j \neq j' \in U_M$ given respectively in equations (3.5) and (3.6). For stratified SRSWOR with $H$ strata of size $N_h$ in $U_F$, $h = 1, \dots, H$, the $c_k$-condition holds if the following condition is verified :
(C2) for all $k \in U_T$, all units $i \in U_F$ linked to k belong to the same stratum.

Then the values $c_{kk'}$, $k \neq k' \in U_T$, are equal to the values $c_{jj'}$, $j \neq j' \in U_M$ given in equation (3.7).

Let $\tilde{\Delta}$ be the matrix of size $N_M \times N_M$ composed of the values

$$ \tilde{\Delta}_{jj'} = \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{i'j'}^{FM}, j, j' \in U_M. \quad (3.8) $$

Let $\tilde{\Delta}_k$ be the submatrix of $\tilde{\Delta}$ corresponding to the elements in position $j$ and $j'$ such that $j$ and $j'$ are linked to $k$. Let $\tilde{1}_k$ be the $N_{Mk}$ dimensional vector of ones. If the $\Delta$-property for $U_M$ and the $c_k$-condition are verified, we can derive the following proposition on the MT weights.

**Proposition 3.3.2.** If the links between $U_F$ and $U_T$ are MtO-MtO, the sampling design satisfies the $\Delta$-property for $U_M$, the $c_k$-condition is verified, and the FM weights $\tilde{\theta}_{ij}^{FM}$, $i \in U_F$, $j \in U_M$, verify :

$$ (\tilde{\theta}_{ij}^{FM})_{i \in U_{Fj}}^t \Delta_j^M (\tilde{\theta}_{i'j}^{FM})_{i' \in U_{Fj}} - c_k \neq 0 \text{ for all } j \in U_{Mk} \\ \text{and } c_k \tilde{1}_k^t \text{diag} \left( ((\tilde{\theta}_{ij}^{FM})_{i \in U_{Fj}}^t \Delta_j^M (\tilde{\theta}_{i'j}^{FM})_{i' \in U_{Fj}} - c_k)^{-1} \right) \tilde{1}_k \neq -1, $$

then there is a unique set of MT weights $\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM})$, $j \in U_M$, $k \in U_T$, functions of $(\tilde{\theta}_{ij}^{FM})_{i \in U_{Fj}}$ and called $\theta^{FM}$-optimal MT weights, minimizing the variance of $\hat{t}_{y2}$ for any variable of interest. These weights are given by :

$$ (\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM}))_{j \in U_{Mk}} = \tilde{\Delta}_k^{-1} \tilde{1}_k (\tilde{1}_k^t \tilde{\Delta}_k^{-1} \tilde{1}_k)^{-1}, \text{ if } j \text{ is linked to } k, \text{ and } 0 \text{ otherwise.} $$

If we compute the $\theta^{FM}$-optimal MT weights for the optimal FM weights given by Proposition 3.3.1, we get a unique set of FM and MT weights that minimize the variance of the double GWSM estimator for any variable of interest as detailed in Proposition 3.3.3.

**Proposition 3.3.3.** If the links are MtO-MtO, the $\Delta$-property holds for $U_M$, the $c_k$-condition is verified and for all $j \in U_{Mk}$, $k \in U_T$, the following requirement is verified,

$$ (\text{inversibility conditions}) \quad c_k \lambda_j \neq 1 \quad \text{and} \quad c_k \sum_{j \in U_{Mk}} \lambda_j/(1 - c_k \lambda_j) \neq -1, $$

with the value of $c_k$ given by the $c_k$-condition, then the $\theta^{FM}$-optimal MT weights for the optimal FM weights $\tilde{\theta}_{ij}^{FM, opt}$, $i \in U_F$, $j \in U_M$, are given by :
$$ \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}) = \frac{\lambda_j/(1 - c_k\lambda_j)}{\sum_{j' \in U_{Mk}} \lambda_{j'}/(1 - c_k\lambda_{j'})}, \text{ if } j \text{ is linked to } k, \text{ and } 0 \text{ otherwise.} $$

Moreover, we can show that the optimal weights for the simple GWSM estimator $\tilde{\theta}_{ik}^{opt}$ can be rewritten for all $i \in U_F, k \in U_T$ as :

$$ \tilde{\theta}_{ik}^{opt} = \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM, opt} \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}). $$

The proofs of Propositions 3.3.1, 3.3.2 and 3.3.3 are given in the Appendix.

We show in Proposition 3.3.3 the existence of optimal link weights for double GWSM estimators with double standardization. Using these weights allows for a reduction in the number of observed links compared to a simple GWSM while preserving the accuracy of the optimal simple GWSM estimator. In what follows, the double GWSM estimator with double standardization obtained by using optimal FM and $\theta^{FM, opt}$-optimal MT weights is called optimal double GWSM estimator and denoted $\hat{l}_{y2}^{opt}$, and the set of optimal FM weights and $\theta^{FM, opt}$-optimal MT weights is referred to as optimal double GWSM weights.

Let us derive the optimal FM weights and the $\theta^{FM, opt}$-optimal MT weights for Poisson sampling, SRSWOR and stratified SRSWOR. We denote by $l_{ij}^{FM}$ (resp. $l_{jk}^{MT}$) the link between unit $i \in U_F$ and unit $j \in U_M$ (resp. between $j \in U_M$ and unit $k \in U_T$) with $l_{ij}^{FM} = 1$ if $i$ and $j$ are linked and 0 otherwise (resp. $l_{jk}^{MT} = 1$ if $j$ and $k$ are linked and 0 otherwise). For Poisson sampling, we have $c_k = 0$ for all $k \in U_T$, thus the inversibility conditions hold. Let $\Delta_{ii} = (1 - \pi_i)/\pi_i$ for $i \in U_F$. The optimal double GWSM weights are given by

$$ \tilde{\theta}_{ij}^{FM, opt} = \frac{l_{ij}^{FM} / \Delta_{ii}}{\sum_{i' \in U_F} l_{i'j}^{FM} / \Delta_{i'i'}} \quad \text{and} \quad \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}) = \frac{l_{jk}^{MT} \sum_{i \in U_F} l_{ij}^{FM} / \Delta_{ii}}{\sum_{j' \in U_M} l_{j'k}^{MT} \sum_{i' \in U_F} l_{i'j'}^{FM} / \Delta_{i'i'}}, $$

for all $j \in U_M$ and $k \in U_T$. For SRSWOR sampling, it can be shown, following Deville and Lavallée (2006), that $\lambda_j = -(1/c_k) N_{Fj} / (N_F - N_{Fj})$, $j \in U_{Mk}, k \in U_T$. Then $c_k \lambda_j \neq 1$ and

$$ c_k \sum_{j \in U_{Mk}} \lambda_j/(1 - c_k\lambda_j) = - \sum_{j \in U_{Mk}} N_{Fj}/N_F. $$

For any unit $k$ from $U_T$, let us denote $U_{Fk}$ of size $N_{Fk}$ the set of units in $U_F$ linked to $k$. When the links are of type MtO-MtO, $N_{Fk} = \sum_{j \in U_{Mk}} N_{Fj}$ for all $k \in U_T$ and $N_F = \sum_{k \in U_T} N_{Fk}$. If $N_T > 1$, then $-N_{Fk}/N_F \neq -1$ for all $k \in U_T$ and the inversibility conditions hold. The
optimal double GWSM weights are given by

$$ \tilde{\theta}_{ij}^{FM, opt} = l_{ij}^{FM} / N_{Fj} \text{ and } \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}) = N_{Fj} l_{jk}^{MT} / \sum_{j' \in U_F} (N_{Fj'} l_{j'k}^{MT}) $$

for all $j \in U_M$ and $k \in U_T$. Details on the derivation of the value of $(\Delta_j^M)^{-1}1_j^M$ can be found in Deville and Lavallée (2006). The remaining computation is straightforward.

Following the reasoning in Medous et al. (2023a), it can be shown that the optimal link weights $\tilde{\theta}_{ij}^{FM, opt}$ and $\theta^{FM, opt}$-optimal link weights $\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$ for stratified SRSWOR are the same as those for SRSWOR if the previously mentioned conditions (C1) and (C2) are verified.

Following Medous et al. (2023a), we can show that, if the links are MtO-MtO, the $\Delta$-property holds for $U_M$ and the $c_k$-condition is verified, the loss of efficiency between an estimator $\hat{t}_{y2}$ given by (3.3) and $\hat{t}_{y2}^{opt}$ is given by :

$$ \mathrm{Var}(\hat{t}_{y2}) - \mathrm{Var}(\hat{t}_{y2}^{opt}) = \sum_{k \in U_T} y_k^2 \mathrm{Var}(\hat{t}_{\hat{\theta}_k} - \hat{t}_{\hat{\theta}_k}^{opt}) = \sum_{k \in U_T} y_k^2 \mathrm{Var}(\hat{t}_{\hat{\theta}_k - \hat{\theta}_k^{opt}}), $$

where

$$ \hat{t}_{\hat{\theta}_k^{opt}} = \sum_{i \in s_F} d_i \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM, opt} \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}) $$

is the HT estimator of the total

$$ t_{\hat{\theta}_k^{opt}} = \sum_{i \in U_F} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM, opt} \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}) = 1, $$

and

$$ \hat{t}_{\hat{\theta}_k - \hat{\theta}_k^{opt}} = \sum_{i \in s_F} d_i \left( \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} - \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM, opt} \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}) \right) $$

is the HT estimator of the total $t_{\hat{\theta}_k - \hat{\theta}_k^{opt}} = 0$, for all $k \in U_T$.

Consider now the double GWSM estimator with double standardization $\hat{t}_{y2}^{FM, opt}$ computed with the optimal FM weights $\tilde{\theta}_{ij}^{FM, opt}$, $i \in U_F$, $j \in U_M$ and any weights $\tilde{\theta}_{jk}^{MT}$, $k \in U_T$. The loss of efficiency between the optimal double GWSM estimator $\hat{t}_{y2}^{opt}$ and the double GWSM estimator $\hat{t}_{y2}^{FM, opt}$ can be expressed as a function of the differences between MT weights.

**Proposition 3.3.4.** If the links are MtO-MtO, the $\Delta$-property holds for $U_M$ and the $c_k$-condition is verified, then the loss of efficiency between the optimal double GWSM estimator
$\hat{t}_{y2}^{opt}$ and the double GWSM estimator $\hat{t}_{y2}^{FM, opt}$ is given by :

$$Var(\hat{t}_{y2}^{FM, opt}) - Var(\hat{t}_{y2}^{opt}) = \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{MT} - \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}))^2 \left(\frac{1}{\lambda_j} - c_k\right)$$

with the value of $c_k$ given by the $c_k$-condition.

Since Poisson sampling, SRSWOR and, under conditions (C1) and (C2), stratified SR-WOR verify the Δ-property for $U_M$ and the $c_k$-condition, we can compute the loss of efficiency between $\hat{t}_{y2}^{opt}$ and $\hat{t}_{y2}^{FM, opt}$ for these sampling designs as follows.

**Corollary 3.3.3.** If the links are MtO-MtO and the sampling designs are Poisson, SRSWOR, or stratified SRSWOR under conditions (C1) and (C2), then the loss of efficiency between $\hat{t}_{y2}^{opt}$ and $\hat{t}_{y2}^{FM, opt}$ can be expressed using Proposition 3.3.4.

For Poisson sampling, we have :

$$Var(\hat{t}_{y2}^{FM, opt}) - Var(\hat{t}_{y2}^{opt}) = \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{MT} - \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}))^2 \frac{1}{\sum_{i \in U_F} l_{ij} / \Delta_{ii}}$$

For SRSWOR of size $n_F$ from $U_F$, let us denote $f = n_F/N_F$. We have :

$$Var(\hat{t}_{y2}^{FM, opt}) - Var(\hat{t}_{y2}^{opt}) = \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{MT} - \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}))^2 \left( \frac{1}{n_F} \frac{N_F - n_F}{N_F - 1} \frac{N_F - N_{Fj}}{N_{Fj}} + \frac{1 - f}{f} \frac{1}{N_F - 1} \right).$$

For stratified SRSWOR with H strata of size $N_h$ in $U_F$, $h = 1, \dots, H$, let us denote $f_h = n_h/N_h$, $h = 1, \dots, H$. Proposition 3.3.4 holds if the previously mentioned conditions (C1) and (C2) are verified. Then we have :

$$Var(\hat{t}_{y2}^{FM, opt}) - Var(\hat{t}_{y2}^{opt}) = \sum_{h=1}^{H} \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{MT} - \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}))^2 \left( \frac{1}{n_h} \frac{N_h - n_h}{N_h - 1} \frac{N_h - N_{Fj}}{N_{Fj}} + \frac{1 - f_h}{f_h} \frac{1}{N_h - 1} \right).$$

The expressions of the difference of variances for SRSWOR and stratified SRSWOR can be approximated by simpler expressions as follows.

**Corollary 3.3.4.** If Proposition 3.3.4 holds, then, under conditions detailed in the Appendix, we have :
For SRSWOR of size $n_F$ from $U_F$ resulting in an indirect sample of size $n_T$ from $U_T$:

$$n_T \left( \frac{\mathrm{Var}(\hat{t}_{y2}^{FM, opt})}{N_T^2} - \frac{\mathrm{Var}(\hat{t}_{y2}^{opt})}{N_T^2} \right) \\ \approx \frac{n_T}{N_T} \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} \left( \tilde{\theta}_{jk}^{MT} - \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}) \right)^2 \left( \frac{1}{n_F} \frac{N_F - n_F}{N_F - 1} \frac{N_F - N_{Fj}}{N_{Fj}} \right).$$

For stratified SRSWOR with H strata of size $N_h$ in $U_F$, $h = 1, \dots, H$:

$$n_T \left( \frac{\mathrm{Var}(\hat{t}_{y2}^{FM, opt})}{N_T^2} - \frac{\mathrm{Var}(\hat{t}_{y2}^{opt})}{N_T^2} \right) \\ \approx \frac{n_T}{N_T} \sum_{h=1}^{H} \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} \left( \tilde{\theta}_{jk}^{MT} - \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}) \right)^2 \left( \frac{1}{n_h} \frac{N_h - n_h}{N_h - 1} \frac{N_h - N_{Fj}}{N_{Fj}} \right).$$

In La Poste survey, the sampling design is complex but can be approximated by a stratified SRSWOR verifying both the $\Delta$-property and $c_k$-condition. As a result, optimal FM weights and $\theta^{FM, opt}$-optimal MT weights exist. Although the optimal FM weights are currently being used at La Poste, the $\theta^{FM, opt}$-optimal MT weights cannot be computed which results in a loss of precision compared to the optimal double GWSM estimator as discussed in Medous et al. (2023a).

## 3.4 Alternative to double GWSM optimal link weights

In this section, we address the situation where the $\theta^{FM, opt}$-optimal MT weights cannot be computed. We propose alternative weights and discuss their practical computation.

### 3.4.1 Alternative link weights using auxiliary information

In section 3.3, we demonstrate that there exists a set of optimal FM weights and $\theta^{FM, opt}$-optimal MT weights that minimize the variance of double GWSM estimators under certain conditions, irrespective of the variable of interest. To compute the optimal FM weights, the knowledge of all links between $U_F$ and indirectly sampled units in $U_M$, as well as the sampling design on $U_F$, is required. Usually, this information is available in practice as the sampling design is selected by the statistician. Moreover, in cases where the links are MtO-MtO, the number of links between $U_F$ and units in $U_M$ may be reasonably low. For instance, in the case of La Poste, a sorting box typically contains 10 addresses.

However, computing the $\theta^{FM, opt}$-optimal MT weights between the indirectly sampled
units $j \in U_M$ and $k \in U_T$ can be challenging in practice as it requires knowing the value of $\lambda_{j'} = 1_{j'}^{Mt}(\Delta_{j'}^M)^{-1}1_{j'}^M$ for all units $j' \in U_M$ linked to $k$. In La Poste case, for instance, $\lambda_j$ is the number of addresses in box $j'$. Unfortunately, the number of addresses contained in each box is not always explicitly stated, and counting them individually is time-consuming, as mentioned in Medous et al. (2023a). Therefore, it is not possible to compute the $\theta^{FM, opt}$-optimal MT weights for La Poste survey. However, it is not uncommon for auxiliary information on $U_M$ to be available. For example, La Poste has access to datasets containing a list of characteristics for all boxes, including their dimensions and an approximate count of the number of addresses they contain. If such information is available, an alternative MT weight can be used in place of the $\theta^{FM, opt}$-optimal MT weights.

If auxiliary information $x_j$ is available for all $j \in U_M$, then, for all $k \in U_T$, we can compute the following alternative MT weight

$$ \tilde{\theta}_{jk}^{MT,alt} = \frac{x_j}{\sum_{j' \in U_{Mk}} x_{j'}}, \quad \text{if } j \text{ is linked to } k, \text{ and } 0 \text{ otherwise.} \qquad (3.9) $$

In what follows, the **GWSM** estimator computed using $\tilde{\theta}_{ik}^{MT,alt}$ is called alternative double **GWSM** estimator and is denoted by $\hat{t}_{y2}^{alt}$. In practice, in order to compute an alternative double **GWSM** estimate, it is enough to know $x_j, j \in U_{Mk}$ for all $k \in U_T$ indirectly sampled.

It is possible to choose $x_j, j \in U_M$, constant for all $j \in U_M$. In that case the weight given by (3.9) reduces to the usual MT weight $\tilde{\theta}_{jk}^{MT,alt} = 1/N_{Mk}, k \in U_T, j \in U_{Mk}$. The usual MT weights are the easiest MT weights to compute, as they only require to know the size of the set $U_{Mk}$ for a given $k \in U_T$. However, these weights generally have poor performance when compared to the $\theta^{FM, opt}$-optimal MT weights. Further information regarding the double **GWSM** estimator using usual MT weights can be found in Medous et al. (2023a).

### 3.4.2 How to choose the best auxiliary variable

In practice, it is not uncommon to have more than one auxiliary variable available for all $j \in U_M$. To minimize the loss of efficiency of the alternative double **GWSM** estimator, we aim to select the variable $x_j, j \in U_M$, that results in the smallest difference between the alternative MT weights and the $\theta^{FM, opt}$-optimal MT weights. We propose two approaches to identify such a variable.

One approach to select the best auxiliary variable $x_j$ for each $j \in U_M$ is to choose a variable that is highly correlated with the non-standardized $\theta^{FM, opt}$-optimal MT weight $\theta_{jk}^{MT}(\theta^{FM, opt})$ for all $j \in U_{Mk}, k \in U_T$. This method allows to compute the correlations between the non-standardized $\theta^{FM, opt}$-optimal MT weight and the available auxiliary va-
riables. For instance, in La Poste case, the non-standardized $\theta^{FM, opt}$-optimal MT weight is the number of addresses in a box, which is known for all indirectly sampled boxes. The survey provides additional variables, such as the location and dimensions of the boxes, which can be used to determine the relationship between the number of addresses and the auxiliary variables. However, the loss of efficiency between the optimal and alternative double GWSM also depends on the value of the variable of interest, as demonstrated in Proposition 3.3.4. In other words, when the variable of interest takes high values, the difference in variance between the optimal and alternative double GWSM estimators can still be substantial, even when the difference in weights is small. This suggests that the correlation between the chosen auxiliary variable and the non-standardized $\theta^{FM, opt}$-optimal MT weight may not be strong enough to minimize the loss of efficiency.

The second approach involves using a variable $x_j$ that is correlated with the variable $y_j = y_k \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$ for all $j \in U_{Mk}$. This approach takes into account the impact of the variable of interest and typically results in better weights than the first approach for lower correlation, as demonstrated in section 3.5. However, this approach requires to obtain the $\theta^{FM, opt}$-optimal MT weights for many indirectly sampled boxes, which can be challenging. In practice, the $\theta^{FM, opt}$-optimal MT weights are unknown and the choice of the auxiliary variable is based on expert knowledge. Moreover, since the alternative MT weights obtained depend on the values $y_k$, $k \in U_T$, they have to be computed for each variable of interest, unlike the first approach which gives a unique set of weights suitable for all variables of interest.

When the number of auxiliary variables available on $U_M$ is large, it can be unrealistic to expect that there exists a single variable $x_j$ that is strongly correlated with either $\theta_{jk}^{MT, opt}$ or $y_j$ for all $j \in U_{Mk}$ and $k \in U_T$. Instead, it is more practical to assume that $x_j$ is highly correlated with $\theta_{jk}^{MT, opt}$ or $y_j$ for only a subset of units $k \in U_T$. To enable more flexibility in the choice of auxiliary variables, an alternative standardized weight can be computed by using different variables for different subsets of $U_T$.

**Definition 3.4.1** (subset weights). Let the population $U_T$ be divided into $G$ subsets $g$, $g = 1, \dots, G$, such that, for each subset $g$, a variable $x_{gj}$ is available for all $j \in U_M$ linked to the units $k \in g$. Then, for all subset $g$, the following subset weights can be used as alternative MT weights :

$$ \tilde{\theta}_{jk}^{MT} = \frac{x_{gj}}{\sum_{j \in U_{Mk}} x_{gj}}, j \in U_{Mk}, k \in g. $$

Using subset weights allows for greater precision in selecting the alternative MT weights as both solutions can be applied to the same estimator.
## 3.5 Simulations

The simulations presented in the following section aim to demonstrate the performance of the alternative double GWSM estimators based on the choice of auxiliary variables used to compute the alternative MT weights. The conclusion of this section is meant as a practical guide for computing alternative MT weights.

### 3.5.1 Populations

The populations used in the simulations are generated to replicate the setup of La Poste, where $U_F$ represents the population of addresses, $U_M$ the population of boxes, and $U_T$ the population of rounds. To ensure the validity of the results presented in sections 3.3 and 3.4, the links between populations are MtO-MtO. An address $i$ is linked to a box $j$ (resp. $j$ is linked to a round $k$) if $j$ contains letters to $i$ (resp. if $k$ deliver the letters contained in $j$).

La Poste has access to extensive historical data on rounds, including information such as traffic volume and details on sampled boxes within each round. In our simulations, we consider a target population $U_T$ of 6,750 rounds selected from this historical data. The frame population $U_F$ of size 5,816,907 and the intermediate population $U_M$ of size 400,011 are obtained by generating numbers of addresses in a box and numbers of boxes in a round from distributions observed in La Poste's historical data.

Our objective is to estimate the total of three variables of interest on the target population $U_T$. The first variable of interest is the traffic $t$, which is derived from La Poste data history after modification for confidentiality reasons. The second variable of interest, $u$, follows a uniform distribution. The third variable of interest, $q$, is proportional to the inverse of the traffic. The bounds of the uniform distributions of variable $u$ and the coefficient of proportionality of variable $q$ are selected in such a way that all variables have the same mean as the traffic variable.

Samples from the target population $U_T$ are obtained through double indirect sampling, with samples from the frame population $U_F$ drawn using SRSWOR. In that case, we saw in Sections 3.3 and 3.4 that, for a box $j \in U_M$, the $\theta^{FM, opt}$-optimal MT weights before standardization reduces to the number $N_{Fj}$ of addresses in box $j$ and the optimal FM weights reduces to $1/N_{Fj}$. In the simulations, we compare five double GWSM estimators, each with optimal FM weights and different MT weights. The five considered MT weights are computed by using the following variables $x_j$, $j \in U_M$:

1. the number of addresses in a box $N_{Fj}$, which gives the $\theta^{FM, opt}$-optimal MT weights
$$ \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt}) = \frac{N_{Fj}}{N_{Fk}}, $$

2. a constant variable $D$, $D \in \mathbb{R}$, which gives the usual MT weight $D_k = \frac{1}{N_{Mk}}$,
3. a variable alt $N_{Fj}$ correlated with $N_{Fj}$,

4. a variable alt $y_j$ correlated with $y_j = y_k \tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$,

5. and a variable alt $\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$ correlated with the optimal weights.

We aim to observe the behavior of the estimators obtained using alt $N_{Fj}$, alt $y_j$, and alt $\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$ when the correlation to respectively $N_{Fj}$, $y_j$, and $\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$ varies. Five correlation values, namely 0.1, 0.3, 0.5, 0.7, and 0.9, are considered for each variable. To achieve this, for each $j \in U_M$, a different value of alt $N_{Fj}$ (resp. alt $y_j$, alt $\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$) is generated for each correlation. The values of alt $N_{Fj}$ (resp. alt $y_j$, alt $\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$), $j \in U_M$, are obtained by taking a random value between 0.5$N_{Fj}$ and 1.5$N_{Fj}$ (resp. 0.5$y_j$ and 1.5$y_j$, 0.5$\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$ and 1.5$\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$) using a uniform distribution. Then, for each $j \in U_M$, a uniformly distributed random number is added to alter the correlation value.

### 3.5.2 Parameters and results

We draw $R = 10,000$ samples of size 1,000 in the frame population using the SRSWOR sampling design. For each r-th sample, $r = 1, \dots, R$, we compute the value of the five double GWSM estimators considered in subsection 3.5.1. Let us denote $\hat{t}_{y2}^{(r)}$ a double GWSM estimate computed for the r-th sample, $r = 1, \dots, R$. For each double GWSM estimator $\hat{t}_{y2}$ introduced in subsection 3.5.1 we have computed the Monte Carlo relative bias :

$$RB_{MC}(\hat{t}_{y2}) = 100 \times \frac{1}{R} \sum_{r=1}^{R} \frac{\hat{t}_{y2}^{(r)} - t_y}{t_y},$$

and the Relative Root Mean Square Error (RRMSE) with respect to the optimal double GWSM estimator $\hat{t}_{y2}^{opt}$:

$$RRMSE_{MC}(\hat{t}_{y2}) = \sqrt{\frac{MSE_{MC}(\hat{t}_{y2})}{MSE_{MC}(\hat{t}_{y2}^{opt})}}$$

where

$$MSE_{MC}(\hat{t}_{y2}) = R^{-1} \sum_{r=1}^{R} \left( \hat{t}_{y2}^{(r)} - R^{-1} \sum_{r=1}^{R} \hat{t}_{y2}^{(r)} \right)^2 .$$

Figures 3.5, 3.6 and 3.7 display the values of the RRMSE of each double GWSM estimator for different correlations, with respectively variables $t$, $u$ and $q$ as the variable of interest. As the absolute value of the relative bias of all estimators considered remains below 0.3%, it is not reported.

The estimators obtained with alt $\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM, opt})$ are consistently less precise than the estimators obtained with alt $N_{Fj}$ in all figures, regardless of the correlation. Noting that the
<figure><img src="image_47.png" /><figcaption>FIGURE 3.5 – Relative efficiency of the alternative double GWSM estimates compared to optimal double GWSM estimates of the total traffic.</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>FIGURE 3.6 – Relative efficiency of the alternative double GWSM estimates compared to optimal double GWSM estimates of the total of the inverse traffic.</figcaption></figure>
Y: uniform variable

<figure><img src="image_49.png" /><figcaption>FIGURE 3.7 – Relative efficiency of the alternative double GWSM estimates compared to optimal double GWSM estimates of the total of the uniform variable.</figcaption></figure>

computation of *alt θ̃<sub>jk</sub><sup>MT, opt</sup>* (*θ*<sup>FM, opt</sup>) requires knowledge of the *θ*<sup>FM, opt</sup>-optimal MT weights *N*<sub>*Fj*</sub>/*N*<sub>*Fk*</sub> for some cases *j*, making it more costly to use than *alt N*<sub>*Fj*</sub>, we do not recommend using *alt θ̃<sub>jk</sub><sup>MT, opt</sup>* (*θ*<sup>FM, opt</sup>) to compute alternative MT weights.

The precision of the estimators obtained with *alt N*<sub>*Fj*</sub> is consistently the best in all figures when the correlation is higher than 0.8. However, in figures 3.6 and 3.7, when the correlation is low, the usual estimator performs better. When the correlation is reasonably high (above 0.1 for figure 3.6 and 0.4 for figure 3.7) yet below 0.8, the estimator obtained with *alt y*<sub>*j*</sub> is the best. In some cases, as in figure 3.5, the estimator obtained with *alt y*<sub>*j*</sub> consistently performs better than the usual estimator.

The results suggest that while the usual double **GWSM** estimator may not be the most efficient alternative in terms of variance, its loss in efficiency compared to the optimal double **GWSM** estimator is still acceptable, with a multiplier ranging from 1.6 to 2. However, it is important to note that these simulations were designed specifically for situations similar to La Poste, where the loss in efficiency between the usual and optimal double **GWSM** estimators is around 1.6 (as discussed in Medous et al. (2023a)).

To summarize, we recommend the use of alternative MT weights obtained with *alt N*<sub>*Fj*</sub>
due to their low computation cost and close proximity to the *θ*<sup>FM, opt</sup>-optimal MT weights,
provided that the correlation between *alt N*<sub>*Fj*</sub> and *N*<sub>*Fj*</sub> is high. Alternatively, when the
correlation between all potential *alt* $N_{Fj}$ and $N_{Fj}$ is too low, it may be beneficial to choose a variable *alt* $y_j$ with a reasonable correlation to $y_j$ instead. However, this method requires to know the value of $y_j$ for many units $j \in s_M$ and for all variables of interest, which can lead to an increase in observation and computation costs. The usual double **GWSM** estimator is the simplest to compute and show an acceptable loss of precision in our simulations. Medous et al. (2023a) study the impact of the distribution of links between populations on the precision of the usual double **GWSM** and identify a few situations where the usual double **GWSM** estimator performs poorly compared to the optimal simple **GWSM** estimator.

### 3.5.3 Perspectives for La Poste application

Currently, La Poste calculates the total quarterly traffic carried by rounds through double indirect sampling and double **GWSM**, using a frame population of addresses and an intermediate population of boxes. The samples in the frame population are obtained through stratified SRSWOR, satisfying both the Δ-property and $c_k$-condition, such that there exists a set of optimal FM weights and $\theta^{FM,opt}$-optimal MT weights given by Proposition 3.3.3.

In La Poste case, $N_{Fj}$ is the number of addresses in the box $j$ and $N_{Fk}$ is the number of addresses in the round $k$. The optimal FM weight for a box $j$ is $1/n_{Fj}$ and the $\theta^{FM,opt}$-optimal MT weight between $j$ and round $k$ is $N_{Fj}/N_{Fk}$ if $k$ and $j$ are linked, 0 otherwise. Unfortunately, in La Poste case, the number of addresses $N_{Fk}$ in round $k$ is unknown for all $k \in U_T$ and the $\theta^{FM,opt}$-optimal MT weights cannot be computed for La Poste estimator. Thus, La Poste currently uses optimal FM weights and usual MT weights to compute its double **GWSM** estimator. However, in the future, the sample size will be reduced, leading to a loss of precision in the estimators. Therefore, La Poste aims to enhance their MT weights to limit this loss of precision.

La Poste intends to use available auxiliary information on the population of boxes to compute alternative MT weights. This information includes dimensions of the boxes, which are observed for all boxes linked to an indirectly sampled round. However, using this variable to compute alternative MT weights did not result in a significant improvement in precision compared to the usual MT weights.

La Poste is currently exploring the use of another observed variable, namely the traffic within the box, which is available for all boxes of an indirectly sampled round. This variable is assumed to be correlated with the variable $y_j$ defined in section 3.4 for a set of round that is yet to be identified. Work is ongoing to evaluate the potential of using this variable to compute alternative MT weights and improve the precision of the estimators.

As an additional source of information, La Poste has access to the RAO dataset, which
contains a non-updated number of addresses in boxes. While this dataset shows potential for computing alternative MT weights, it currently suffers from inconsistencies and cannot be used at the moment.

## 3.6 Conclusion

The use of indirect sampling coupled with double GWSM can reduce the number of observed links, but this comes at the cost of a potential decrease in accuracy of the double GWSM estimator, as discussed in Medous et al. (2023a). In the present paper, we demonstrate the existence of optimal weights for double GWSM that can be used to reduce costs without deteriorating the precision of the estimators relative to an optimal simple GWSM estimator. We have also proposed alternative weights for situations where the optimal weights are not computable, and have demonstrated their relative efficiency in a case similar to La Poste example. To gain a deeper understanding of how alternative weights work, it may be necessary to develop models for both the variable of interest (Särndal et al. (1992)) and the links (Xu and Lavallée (2009)). However, this topic is not in the scope of our paper, which is primarily focused on optimal double GWSM.

## 3.7 Appendix

### Propositions from Medous et al. (2023a)

For any unit $k$ from $U_T$, let us denote $U_{Fk}$ of size $N_{Fk}$ the set of units in $U_F$ linked to $k$. For MtO links, the sets $U_{Fk}$, $k \in U_T$ form a partition of $U_F$ and we can assume units in $U_F$ are ordered with respect to the sub-populations $U_{Fk}$. We denote by $\Delta = (\Delta_{ii'})_{i,i' \in U_F}$ the covariance matrix of size $N_F \times N_F$. For given $k, k'$ from $U_T$, let $\Delta_{kk'} = (\Delta_{ii'})_{i \in U_{Fk}, i' \in U_{Fk'}}$ be the submatrix of $\Delta$ corresponding to elements in position $i$ (resp. $i'$) if $i$ (resp. $i'$) is linked to $k$ (resp. $k'$) and let $\Delta_k$ be the submatrix $\Delta_{kk}$. For MtO links, these submatrices form a partition of $\Delta$:

$$ \Delta = (\Delta_{kk'})_{k,k' \in U_T}. $$

Let us denote by $1_k$ the vector of ones of size $N_{Fk}$ and by $A^t$ the transpose of a given matrix or vector $A$.

The following definition was introduced by Medous et al. (2023a) to characterize all sampling designs for which a set of optimal standardized link weight exists.

*Definition 3.7.1* (Δ-property for $U_T$). For MtO links, a sampling design is said to satisfy the
Δ-property for $U_T$ if, for any $k$ in $U_T$, $\Delta_k$ is invertible and for any $k \neq k'$ in $U_T$,

$$ \Delta_{k,k' \neq k} = c_{kk'} \mathbf{1}_k \mathbf{1}_{k'}^t \quad \text{with } c_{kk'} \text{ not depending on } i \text{ and } i' $$

The following proposition, introduced by Medous et al. (2023a), gives the expression of the set of optimal link weights for all designs satisfying the Δ-property.

**Proposition 3.7.1 (Medous et al. (2023a)).** If the links are MtO and the sampling design satisfies the Δ-property for $U_T$, then there exists a unique set of standardized link weights minimizing the variance of $\hat{t}_{y1}$ for all variable of interest $y$ given by :

$$ (\tilde{\theta}_{ik}^{opt})_{i \in U_{Fk}} = \Delta_k^{-1} \mathbf{1}_k (\mathbf{1}_k^t \Delta_k^{-1} \mathbf{1}_k)^{-1}, \quad \text{for all } k \in U_T. $$

This proposition holds as long as the matrix Δ verify the Δ-property for $U_T$, but Δ does not need to be a covariance matrix. Let us consider a simple GWSM estimator $\hat{t}_{y1}$, given by equation (3.1), obtained using some standardized link weights $\tilde{\theta}_{ik}$, $i \in U_F$, $k \in U_T$ and the optimal simple GWSM estimator $\hat{t}_{y1}^{opt}$ obtained using optimal link wieghs. The following proposition, introduced by Medous et al. (2023a), gives the loss of efficiency between $\hat{t}_{y1}^{opt}$ and $\hat{t}_{y1}$ when the sampling design satisfy the Δ-property for $U_T$.

**Proposition 3.7.2 (Medous et al. (2023a)).** If the links are MtO and the sampling design satisfies the Δ-property for $U_T$, then the loss of efficiency compared with optimal link weights $\tilde{\theta}_{ik}^{opt}$, $i \in U_F$, $k \in U_T$, is given by :

$$ \mathrm{Var}(\hat{t}_{y1}) - \mathrm{Var}(\hat{t}_{y1}^{opt}) = \sum_{k \in U_T} y_k^2 \mathrm{Var}(\hat{t}_{\tilde{\theta}_k} - \hat{t}_{\tilde{\theta}_k^{opt}}) = \sum_{k \in U_T} y_k^2 \mathrm{Var}(\hat{t}_{\tilde{\theta}_k - \tilde{\theta}_k^{opt}}), $$

where $\hat{t}_{\tilde{\theta}_k^{opt}} = \sum_{i \in s_F} d_i \tilde{\theta}_{ik}^{opt}$ is the Horvitz-Thompson estimator of the total $t_{\tilde{\theta}_k^{opt}} = \sum_{i \in U_F} \tilde{\theta}_{ik}^{opt} = 1$, and $\hat{t}_{\tilde{\theta}_k - \tilde{\theta}_k^{opt}} = \sum_{i \in s_F} d_i (\tilde{\theta}_{ik} - \tilde{\theta}_{ik}^{opt})$ is the Horvitz-Thompson estimator of the total $t_{\tilde{\theta}_k - \tilde{\theta}_k^{opt}} = 0$, for all $k \in U_T$.

## Sherman-Morrison formula

This formula, introduced by Sherman and Morrison (1950), is used in the proof of Propositions 3.3.2 and 3.3.3. Let $A$ be a square matrix of size $n \times n$, $n \in \mathbb{N}$, and $u, v$ two vectors of size $n$. If $A$ is invertible and $v^t A^{-1} u \neq -1$, then $A + uv^t$ is invertible and

$$ (A + uv^t)^{-1} = A^{-1} - \frac{A^{-1}uv^tA^{-1}}{1 + v^tA^{-1}u}. $$
### Proof of Proposition 3.3.1

Let us consider the double GWSM estimator $\hat{t}_{y2}$ computed for any variable of interest $y$ and any set of weights $\tilde{\theta}_{jk}^{FM}$, $j \in U_M$, $k \in U_T$. The double GWSM estimator can be rewritten as a GWSM estimator between $U_F$ and $U_M$ of a variable $z_j = \sum_{k \in U_T} \tilde{\theta}_{jk}^{FM} y_k$:

$$ \hat{t}_{y2} = \sum_{i \in s_F} d_i \left( \sum_{k \in U_T} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} y_k \right) = \sum_{i \in s_F} d_i \left( \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} z_j \right). $$

Since the $\Delta$-property holds for $U_M$, we can apply Proposition 3.7.1 to the GWSM estimator $\sum_{i \in s_F} d_i (\sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} z_j)$ of the total $\sum_{j \in U_M} z_j$. We obtain the unique set of optimal weights between $U_F$ and $U_M$:

$$ (\tilde{\theta}_{ij}^{FM, opt})_{i \in U_{Fj}} = (\Delta_j^M)^{-1} \mathbf{1}_j^M (\mathbf{1}_j^{Mt} (\Delta_j^M)^{-1} \mathbf{1}_j^M)^{-1}, \text{ for all } j \in U_M, $$

which concludes the proof. $\square$

### Proof of Proposition 3.3.2

Let us consider any set of FM weights $\tilde{\theta}_{ij}^{FM}$, $i \in U_F$, $j \in U_M$. The double GWSM estimator computed can be rewritten as a simple GWSM estimator between $U_M$ and $U_T$ with sampling weights $\tilde{d}_j = \sum_{i \in U_{Fj}} d_i \tilde{\theta}_{ij}^{FM}$:

$$ \hat{t}_{y2} = \sum_{i \in s_F} d_i \left( \sum_{k \in U_T} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} y_k \right) = \sum_{j \in s_M} \tilde{d}_j \left( \sum_{k \in U_M} \tilde{\theta}_{ij}^{FM} z_j \right). $$

We can get the matrix $\tilde{\Delta} = (\tilde{\Delta}_{jj'})_{j,j' \in U_M}$ of size $N_M \times N_M$, as defined in equation (3.8), by rewriting the variance of $\hat{t}_{y2}$:

$$ \begin{aligned} \text{Var}(\hat{t}_{y1}) &= \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \sum_{k \in U_T} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} y_k \sum_{k' \in U_T} \sum_{j' \in U_M} \tilde{\theta}_{i'j'}^{FM} \tilde{\theta}_{j'k'}^{MT} y_{k'} \\ &= \sum_{j \in U_M} \sum_{j' \in U_M} \left( \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{i'j'}^{FM} \right) \sum_{k \in U_T} \tilde{\theta}_{jk}^{MT} y_k \sum_{k' \in U_T} \tilde{\theta}_{j'k'}^{MT} y_{k'} \\ &= \sum_{j \in U_M} \sum_{j' \in U_M} \tilde{\Delta}_{jj'} \sum_{k \in U_T} \tilde{\theta}_{jk}^{MT} y_k \sum_{k' \in U_T} \tilde{\theta}_{j'k'}^{MT} y_{k'} \end{aligned} $$

$$ \text{with } \tilde{\Delta}_{jj'} = \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{i'j'}^{FM}. $$

Let us show that the $\Delta$-property for $U_T$ holds for $\tilde{\Delta}$. Let $\tilde{\Delta}_k$ be the submatrix of $\tilde{\Delta}$ corresponding to the elements in position $j$ and $j'$ such that $j$ and $j'$ are linked to $k$. We
have :

$$
\tilde{\Delta}_k = \left( \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{i'j'}^{FM} \right)_{j,j' \in U_{Mk}} = \left( (\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_{jj'}^M (\tilde{\theta}_{i'j'}^{FM})_{i' \in U_{Fj'}} \right)_{j,j' \in U_{Mk}}
$$

Since the Δ-property for $U_M$ and the $c_k$-condition are verified, $\Delta_{jj'}^M = \Delta_j^M$ if $j = j'$ and
$c_k 1_j^M 1_{j'}^{Mt}$ otherwise, and, for all $j, j' \in U_{Mk}$,

$$
(\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_{jj'}^M (\tilde{\theta}_{i'j'}^{FM})_{i' \in U_{Fj'}} = \begin{cases} (\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_j^M (\tilde{\theta}_{i'j'}^{FM})_{i' \in U_{Fj}}, & \text{if } j = j', \\ c_k (\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} 1_j^M 1_{j'}^{Mt} (\tilde{\theta}_{i'j'}^{FM})_{i' \in U_{Fj'}} = c_k & \text{otherwise,} \end{cases}
$$

because, for all $j \in U_M$, $(\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} 1_j^M = 1$ by standardisation.

$\tilde{\Delta}_k$ is the square matrix containing the values of $(\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_{jj'}^M (\tilde{\theta}_{i'j'}^{FM})_{i' \in U_{Fj'}}$ for $j, j' \in U_{Mk}$. Since $(\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_j^M (\tilde{\theta}_{i'j'}^{FM})_{i' \in U_{Fj}}, j \in U_{Mk}$, and $c_k$ are real numbers, this means that $\tilde{\Delta}_k$ is a square matrix with the values $(\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_j^M (\tilde{\theta}_{i'j'}^{FM})_{i' \in U_{Fj}}, j \in U_{Mk}$, on its diagonal and $c_k$ everywhere else. Let $\tilde{1}_k$ be the $N_{Mk}$ dimensional vector of ones and, for any vector $u$, let us denote $diag(u)$ the diagonal matrix containing the elements of $u$. We can rewrite $\tilde{\Delta}_k$, $k \in U_T$, as follows :

$$
\tilde{\Delta}_k = \operatorname{diag} \left( \left( (\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_j^M (\tilde{\theta}_{i'j}^{FM})_{i' \in U_{Fj}} - c_k \right)_{j \in U_{Mk}} \right) + c_k \tilde{1}_k \tilde{1}_k^t.
$$

Using the Sherman-Morrison formula, we can show that $\tilde{\Delta}_k$ is invertible if and only if (iif)
$(\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_j^M (\tilde{\theta}_{i'j}^{FM})_{i' \in U_{Fj}} - c_k \neq 0$ for all $j \in U_{Mk}$ and

$$
c_k \tilde{1}_k^t \operatorname{diag} \left( \left( (\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_j^M (\tilde{\theta}_{i'j}^{FM})_{i' \in U_{Fj}} - c_k \right)_{j \in U_{Mk}} \right)^{-1} \tilde{1}_k \neq -1.
$$

Let $\tilde{\Delta}_{kk'}$, $k \neq k' \in U_T$ be the submatrix of $\tilde{\Delta}$ of size $N_{Mk} \times N_{Mk'}$ corresponding to the elements in position $j$ and $j'$ such that $j$ (resp. $j'$) is linked to $k$ (resp. $k'$). Since the links are of type MtO-MtO, $k \neq k'$ implies that, for all $j \in U_{Mk}$ and $j' \in U_{Mk'}$, $j \neq j'$. Thus, by using the $c_k$-condition, we have :

$$
\begin{align*}
\tilde{\Delta}_{kk'} &= \left( (\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} \Delta_{jj'}^M (\tilde{\theta}_{i'j'}^{FM})_{i' \in U_{Fj'}} \right)_{\substack{j \in U_{Mk} \\ j' \in U_{Mk}'}} \\
&= \left( c_{kk'} (\tilde{\theta}_{ij}^{FM})^t_{i \in U_{Fj}} 1_j^M 1_{j'}^{Mt} (\tilde{\theta}_{i'j'}^{FM})_{i' \in U_{Fj'}} \right)_{\substack{j \in U_{Mk} \\ j' \in U_{Mk}'}} = c_{kk'} \tilde{1}_k \tilde{1}_{k'}^t,
\end{align*}
$$

which proves that $\tilde{\Delta}$ satisfies the $\Delta$-property for $U_T$.

Then we can apply Proposition 3.7.1 to the GWSM estimator $\sum_{j \in s_M} \tilde{d}_j (\sum_{k \in U_T} \tilde{\theta}_{jk}^{MT} y_k)$
and we obtain the unique set of optimal weights between $U_M$ and $U_T$:

$$ (\tilde{\theta}_{jk}^{MT, opt}(\theta^{FM}))_{j \in U_{Mk}} = \tilde{\Delta}_k^{-1} \tilde{1}_k (\tilde{1}_k^t \tilde{\Delta}_k^{-1} \tilde{1}_k)^{-1}, \text{ for all } j \in U_M, $$

which concludes the proof. □

### Proof of Proposition 3.3.3

We assume that the $\Delta$-property for $U_M$ and the $c_k$-condition hold. Let us first compute the value of $\tilde{\Delta}_k$, $k \in U_T$, defined above for the optimal FM weights given by Proposition 3.3.1, $(\tilde{\theta}_{ij}^{FM, opt})_{i \in U_{Fj}} = (\Delta_j^M)^{-1} \tilde{1}_j^M \lambda_j^{-1}$, for all $j \in U_M$, with $\lambda_j = \tilde{1}_j^{Mt} (\Delta_j^M)^{-1} \tilde{1}_j^M$:

$$ \begin{aligned} \tilde{\Delta}_k &= ((\tilde{\theta}_{ij}^{FM, opt})_{i \in U_{Fj}}^t \Delta_{jj'}^M (\tilde{\theta}_{i'j'}^{FM, opt})_{i' \in U_{Fj'}})_{j,j' \in U_{Mk}} \\ &= \left( \frac{\tilde{1}_j^{Mt} (\Delta_j^M)^{-1}}{\lambda_j} \Delta_{jj'}^M \frac{(\Delta_{j'}^M)^{-1} \tilde{1}_{j'}^M}{\lambda_{j'}} \right)_{j,j' \in U_{Mk}} \\ &= \mathrm{diag} \left( (1/\lambda_j - c_k)_{j \in U_{Mk}} \right) + c_k \tilde{1}_k \tilde{1}_k^t \end{aligned} $$

since for $j, j' \in U_{Mk}$, $\tilde{1}_j^{Mt} (\Delta_j^M)^{-1} \Delta_{jj'}^M (\Delta_{j'}^M)^{-1} \tilde{1}_{j'}^M = 1/\lambda_j$ if $j = j'$, and $c_k$ otherwise, with the value of $c_k$, $k \in U_T$, given by the $c_k$-condition.

Since $\Delta_j^M$, $j \in U_M$, is invertible thanks to the $\Delta$-property for $U_M$ and $1/\lambda_j - c_k = (1 - c_k \lambda_j)/\lambda_j$, $j \in U_M$, is a real number, we get by using the Sherman-Morrison formula that $\tilde{\Delta}_k$ is invertible if:

$$ \begin{aligned} & c_k \tilde{1}_k^t \mathrm{diag} \left( (1/\lambda_j - c_k)_{j \in U_{Mk}} \right)^{-1} \tilde{1}_k \neq -1 \\ \iff & c_k \sum_{j \in U_{Mk}} \lambda_j / (1 - c_k \lambda_j) \neq -1 \text{ and } c_k \lambda_j \neq 1 \text{ for all } j \in U_{Mk} \end{aligned} $$

and we have:

$$ \tilde{\Delta}_k^{-1} = \tilde{D}^{-1} - \frac{c_k \tilde{D}^{-1} \tilde{1}_k \tilde{1}_k^t \tilde{D}^{-1}}{1 + c_k \tilde{1}_k^t \tilde{D}^{-1} \tilde{1}_k}, $$

with $\tilde{D}^{-1} = \mathrm{diag} \left( (1/\lambda_j - c_k)_{j \in U_{Mk}} \right)^{-1} = \mathrm{diag} \left( (\lambda_j/(1 - c_k \lambda_j)_{j \in U_{Mk}} \right)$. Then we can compute

$$ \tilde{\Delta}_k^{-1} \tilde{1}_k = \tilde{D}^{-1} \tilde{1}_k - \frac{c_k \tilde{D}^{-1} \tilde{1}_k \tilde{1}_k^t \tilde{D}^{-1} \tilde{1}_k}{1 + c_k \tilde{1}_k^t \tilde{D}^{-1} \tilde{1}_k} = \left( 1 - \frac{c_k \tilde{1}_k^t \tilde{D}^{-1} \tilde{1}_k}{1 + c_k \tilde{1}_k^t \tilde{D}^{-1} \tilde{1}_k} \right) \tilde{D}^{-1} \tilde{1}_k $$

and

$$ \tilde{1}_k^t \tilde{\Delta}_k^{-1} \tilde{1}_k = \left( 1 - \frac{c_k \tilde{1}_k^t \tilde{D}^{-1} \tilde{1}_k}{1 + c_k \tilde{1}_k^t \tilde{D}^{-1} \tilde{1}_k} \right) \tilde{1}_k^t \tilde{D}^{-1} \tilde{1}_k. $$
<figure><img src="image_50.png" /><figcaption>FIGURE 3.8 – Expression of $\Delta_k$ when the $\Delta$-property holds for $U_M$ and the $c_k$-condition is verified.</figcaption></figure>

Finally, we get :

$$ \begin{aligned} \left(\tilde{\theta}_{jk}^{\mathrm{MT}}(\theta^{\mathrm{FM, opt}})\right)_{j \in U_{Mk}} &= \tilde{\Delta}_k^{-1} \tilde{1}_k \left(\tilde{1}_k^t \tilde{\Delta}_k^{-1} \tilde{1}_k\right)^{-1} \\ &= \frac{\tilde{D}^{-1} \tilde{1}_k}{\tilde{1}_k^t \tilde{D}^{-1} \tilde{1}_k} = \left( \frac{\lambda_j / (1 - c_k \lambda_j)}{\sum_{j' \in U_{Mk}} \lambda_{j'} / (1 - c_k \lambda_{j'})} \right)_{j \in U_{Mk}}, \end{aligned} $$

since $\tilde{D}^{-1} = \mathrm{diag}((\lambda_j/(1 - c_k \lambda_j)_{j \in U_{Mk}})$.

We will now show that the optimal weights for simple GWSM estimators given by Proposition 3.7.1 can be rewritten as an expression of the optimal FM weights given by Proposition 3.3.1 and the $\theta^{\mathrm{FM, opt}}$-optimal MT weights if the $\Delta$-property for $U_M$ and the $c_k$-condition hold, and if $c_k \sum_{j \in U_{Mk}} \lambda_j / (1 - c_k \lambda_j) \neq 1$ and $c_k \lambda_j \neq 1$ for all $j \in U_{Mk}$.

Let us compute the value of

$$ (\tilde{\theta}_{ik}^{\mathrm{opt}})_{i \in U_{Fk}} = \Delta_k^{-1} 1_k (1_k^t \Delta_k^{-1} 1_k)^{-1}, \quad \text{for all } k \in U_T $$

when the links are MtO-MtO, the $\Delta$-property holds for $U_M$ and the $c_k$-condition is verified.

For a set of square matrices $(A_j)_{j \in U_{Mk}}$, let $\mathrm{diag}((A_j)_{j \in U_{Mk}})$ be the block diagonal matrix whose main-diagonal blocks are the matrices $(A_j)_{j \in U_{Mk}}$. When the links are MtO-MtO, the $\Delta$-property holds for $U_M$ and the $c_k$-condition is verified, the matrix $\Delta_k$ defined in section 3.2.1 is given in figure 3.8 and can be rewritten for all $k \in U_T$ as

$$ \Delta_k = \mathrm{diag} \left( \left( \Delta_j^M - c_k \mathbf{1}_j^M \mathbf{1}_j^M \right)_{j \in U_{Mk}} \right) + c_k \mathbf{1}_k \mathbf{1}_k^t. $$

We will first show that $\Delta_k$ is invertible if $c_k \sum_{j \in U_{Mk}} \lambda_j / (1 - c_k \lambda_j) \neq 1$ and $c_k \lambda_j \neq 1$
1 for all $j \in U_{Mk}$. Let us denote $D = \text{diag}((\Delta_j^M - c_k \mathbf{1}_j^M \mathbf{1}_j^M)_{j \in U_{Mk}})$ the square matrix of size $N_{Fk}$. Using the Sherman-Morrison formula, we can show that $\Delta_k$ is invertible if $D$ is invertible and $c_k \mathbf{1}_k^t D^{-1} \mathbf{1}_k \neq -1$.

Since $D^{-1} = \text{diag}((\Delta_j^M - c_k \mathbf{1}_j^M \mathbf{1}_j^M)_{j \in U_{Mk}})^{-1} = \text{diag}((\Delta_j^M - c_k \mathbf{1}_j^M \mathbf{1}_j^M)_{j \in U_{Mk}}^{-1})$, $D$ is invertible if $\Delta_j^M - c_k \mathbf{1}_j^M \mathbf{1}_j^M$ is invertible for all $j \in U_{Mk}$.

Since $\Delta_j^M, j \in U_M$, is invertible thanks to the $\Delta$-property for $U_M$, we can show using the Sherman-Morrison formula that $\Delta_j^M - c_k \mathbf{1}_j^M \mathbf{1}_j^M$ is invertible if $c_k \mathbf{1}_j^M (\Delta_j^M)^{-1} \mathbf{1}_j^M = c_k \lambda_j \neq 1$, and we have :

$$ (\Delta_j^M - c_k \mathbf{1}_j^M \mathbf{1}_j^M)^{-1} = (\Delta_j^M)^{-1} + \frac{c_k (\Delta_j^M)^{-1} \mathbf{1}_j^M \mathbf{1}_j^M (\Delta_j^M)^{-1}}{1 - c_k \mathbf{1}_j^M (\Delta_j^M)^{-1} \mathbf{1}_j^M} \\ = (\Delta_j^M)^{-1} + \frac{c_k (\Delta_j^M)^{-1} \mathbf{1}_j^M \mathbf{1}_j^M (\Delta_j^M)^{-1}}{1 - c_k \lambda_j} $$

Then we can compute $D^{-1}$:

$$ D^{-1} = \text{diag} \left( \left( (\Delta_j^M)^{-1} + \frac{c_k (\Delta_j^M)^{-1} \mathbf{1}_j^M \mathbf{1}_j^M (\Delta_j^M)^{-1}}{1 - c_k \lambda_j} \right)_{j \in U_{Mk}} \right). $$

Since $\mathbf{1}_k$ is the vector of size $N_{Fk}$ composed of all vectors of size $N_{Fj}$ $\mathbf{1}_j^M, j \in U_{Mk}$, we can rewrite $D^{-1} \mathbf{1}_k$ as :

$$ D^{-1} \mathbf{1}_k = \left( \left( 1 + \frac{c_k \lambda_j}{1 - c_k \lambda_j} \right) (\Delta_j^M)^{-1} \mathbf{1}_j^M \right)_{j \in U_{Mk}} = \left( \frac{(\Delta_j^M)^{-1} \mathbf{1}_j^M}{1 - c_k \lambda_j} \right)_{j \in U_{Mk}} $$

with $((\Delta_j^M)^{-1} \mathbf{1}_j^M) (1 - c_k \lambda_j)^{-1}, j \in U_{Mk}$, vectors of size $N_{Fj}$. We compute $\mathbf{1}_k^t D^{-1} \mathbf{1}_k$:

$$ \mathbf{1}_k^t D^{-1} \mathbf{1}_k = \sum_{j \in U_{Mk}} \frac{\mathbf{1}_j^M (\Delta_j^M)^{-1} \mathbf{1}_j^M}{1 - c_k \lambda_j} = \sum_{j \in U_{Mk}} \frac{\lambda_j}{1 - c_k \lambda_j} $$

Then $\Delta_k^{-1}$ is invertible if $c_k \mathbf{1}_k^t D^{-1} \mathbf{1}_k = c_k \sum_{j \in U_{Mk}} \frac{\lambda_j}{1 - c_k \lambda_j} \neq -1$ and $c_k \lambda_j \neq 1$ for all $j \in U_{Mk}$ and its inverse is given by :

$$ \Delta_k^{-1} = D^{-1} - \frac{c_k D^{-1} \mathbf{1}_k \mathbf{1}_k^t D^{-1}}{1 + c_k \mathbf{1}_k^t D^{-1} \mathbf{1}_k} $$
We can use Proposition 3.7.1 to compute the optimal weights $\tilde{\theta}_{ik}^{opt}$, $i \in U_F$, $k \in U_T$:

$$ \begin{aligned} (\tilde{\theta}_{ik}^{opt})_{i \in U_{Fk}} &= \Delta_k^{-1} \mathbf{1}_k (\mathbf{1}_k^t \Delta_k^{-1} \mathbf{1}_k)^{-1} \\ &= \left(1 - \frac{c_k \mathbf{1}_k^t D^{-1} \mathbf{1}_k}{1 + c_k \mathbf{1}_k^t D^{-1} \mathbf{1}_k}\right) D^{-1} \mathbf{1}_k \times \left(\left(1 - \frac{c_k \mathbf{1}_k^t D^{-1} \mathbf{1}_k}{1 + c_k \mathbf{1}_k^t D^{-1} \mathbf{1}_k}\right) \mathbf{1}_k^t D^{-1} \mathbf{1}_k\right)^{-1} \\ &= \frac{D^{-1} \mathbf{1}_k}{\mathbf{1}_k^t D^{-1} \mathbf{1}_k}. \end{aligned} $$

By plugging the values of $D^{-1}\mathbf{1}_k$ and $\mathbf{1}_k^t D^{-1}\mathbf{1}_k$ computed above, we get:

$$ \begin{aligned} (\tilde{\theta}_{ik}^{opt})_{i \in U_{Fk}} &= \left( \frac{(\Delta_j^M)^{-1} \mathbf{1}_j^M}{1 - c_k \lambda_j} \right)_{j \in U_{Mk}} \times \left( \sum_{j' \in U_{Mk}} \frac{\lambda_{j'}}{1 - c_k \lambda_{j'}} \right)^{-1} \\ &= \left( \frac{(\Delta_j^M)^{-1} \mathbf{1}_j^M / (1 - c_k \lambda_j)}{\sum_{j' \in U_{Mk}} \lambda_{j'} / (1 - c_k \lambda_{j'})} \right)_{j \in U_{Mk}} \\ &= \left( \frac{((\Delta_j^M)^{-1} \mathbf{1}_j^M / \lambda_j) \times (\lambda_j / (1 - c_k \lambda_j))}{\sum_{j' \in U_{Mk}} \lambda_{j'} / (1 - c_k \lambda_{j'})} \right)_{j \in U_{Mk}} \\ &= \left( \frac{(\Delta_j^M)^{-1} \mathbf{1}_j^M}{\lambda_j} \right)_{j \in U_{Mk}} \times \frac{\lambda_j / (1 - c_k \lambda_j)}{\sum_{j' \in U_{Mk}} \lambda_{j'} / (1 - c_k \lambda_{j'})} \\ &= (\tilde{\theta}_j^{FM})_{j \in U_{Mk}} \tilde{\theta}_{jk}^{MT, opt} (\theta^{FM, opt}) \end{aligned} $$

with $\tilde{\theta}_j^{FM, opt} = (\tilde{\theta}_{ij}^{FM, opt})_{i \in U_{Fj}}$, $j \in U_{Mk}$, vectors of size $N_{Fj}$. Recalling that $\sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT}$ has only one non null element when the links are MtO-MtO, this concludes the proof. $\square$
Proof of Proposition 3.3.4

Using Proposition 3.7.2, the difference of variances between $\hat{t}_{y2}^{MT, opt}$ and $\hat{t}_{y2}^{MT, alt}$ can be rewritten as

$$
\begin{align*}
\mathrm{Var}(\hat{t}_{y2}^{\mathrm{alt}}) - \mathrm{Var}(\hat{t}_{y2}^{\mathrm{opt}}) &= \sum_{k \in U_T} y_k^2 \left( \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} (\tilde{\theta}_{ik}^{\mathrm{alt}} - \tilde{\theta}_{ik}^{\mathrm{opt}}) (\tilde{\theta}_{i'k}^{\mathrm{alt}} - \tilde{\theta}_{i'k}^{\mathrm{opt}}) \right) \\
&= \sum_{k \in U_T} y_k^2 \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \sum_{j \in U_M} \left( (\tilde{\theta}_{jk}^{\mathrm{MT, alt}} - \tilde{\theta}_{jk}^{\mathrm{MT, opt}} (\theta^{\mathrm{FM, opt}})) \tilde{\theta}_{ij}^{\mathrm{FM, opt}} \right) \times \\
&\qquad \sum_{j' \in U_M} \left( (\tilde{\theta}_{j'k}^{\mathrm{MT, alt}} - \tilde{\theta}_{j'k}^{\mathrm{MT, opt}}) \tilde{\theta}_{i'j'}^{\mathrm{FM, opt}} \right) \\
&= \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{\mathrm{MT, alt}} - \tilde{\theta}_{jk}^{\mathrm{MT, opt}} (\theta^{\mathrm{FM, opt}})) \sum_{j' \in U_M} (\tilde{\theta}_{j'k}^{\mathrm{MT, alt}} - \tilde{\theta}_{j'k}^{\mathrm{MT, opt}}) \times \\
&\qquad \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \tilde{\theta}_{ij}^{\mathrm{FM, opt}} \tilde{\theta}_{i'j'}^{\mathrm{FM, opt}}
\end{align*}
$$

For all $j \in U_M$, we denote by $\theta_j^{MT, opt}$ the vector of optimal FM weights defined in Proposition 3.3.3. Then we can rewrite

$$
\begin{align*}
\sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \tilde{\theta}_{ij}^{FM, opt} \tilde{\theta}_{i'j'}^{FM, opt} &= (\theta_j^{MT, opt})^t \Delta_{jj'}^M \theta_j^{MT, opt} \\
&= \frac{\mathbf{1}_j^{Mt} (\Delta_j^M)^{-1}}{\mathbf{1}_j^{Mt} (\Delta_j^M)^{-1} \mathbf{1}_j^M} \Delta_{jj'}^M \frac{(\Delta_{j'}^M)^{-1} \mathbf{1}_{j'}^M}{\mathbf{1}_{j'}^{Mt} (\Delta_{j'}^M)^{-1} \mathbf{1}_{j'}^M}.
\end{align*}
$$

Since the sampling design satisfies the $\Delta$-property for $U_M$ and the $c_k$-condition, $\Delta_{jj'} = c_{jj'} \mathbf{1}_j^{Mt} \mathbf{1}_{j'} = c_k \mathbf{1}_j^{Mt} \mathbf{1}_{j'}$ if $j \neq j'$ both linked to $k$ and

$$
\sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \tilde{\theta}_{ij}^{FM, opt} \tilde{\theta}_{i'j'}^{FM, opt} = \begin{cases} \frac{1}{\mathbf{1}_j^{Mt} (\Delta_j^M)^{-1} \mathbf{1}_j^M}, & \text{if } j = j', \\ c_k, & \text{otherwise.} \end{cases}
$$
Then

$$
\begin{align*}
\mathrm{Var}(\hat{t}_{y2}^{\mathrm{alt}}) - \mathrm{Var}(\hat{t}_{y2}^{\mathrm{opt}}) &= \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{\mathrm{MT,alt}} - \tilde{\theta}_{jk}^{\mathrm{MT,opt}}(\theta^{\mathrm{FM,opt}}))^2 \frac{1}{\mathbf{1}_j^{\mathrm{Mt}} (\mathbf{\Delta}_j^M)^{-1} \mathbf{1}_j^M} + \\
&\qquad \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{\mathrm{MT,alt}} - \tilde{\theta}_{jk}^{\mathrm{MT,opt}}(\theta^{\mathrm{FM,opt}})) \sum_{j' \in U_M, j' \neq j} (\tilde{\theta}_{j'k}^{\mathrm{MT,alt}} - \tilde{\theta}_{j'k}^{\mathrm{MT,opt}}) c_k \\
&= \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{\mathrm{MT,alt}} - \tilde{\theta}_{jk}^{\mathrm{MT,opt}}(\theta^{\mathrm{FM,opt}}))^2 \frac{1}{\mathbf{1}_j^{\mathrm{Mt}} (\mathbf{\Delta}_j^M)^{-1} \mathbf{1}_j^M} + \\
&\qquad \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{\mathrm{MT,alt}} - \tilde{\theta}_{jk}^{\mathrm{MT,opt}}(\theta^{\mathrm{FM,opt}})) \sum_{j' \in U_M} (\tilde{\theta}_{j'k}^{\mathrm{MT,alt}} - \tilde{\theta}_{j'k}^{\mathrm{MT,opt}}) c_k - \\
&\qquad \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{\mathrm{MT,alt}} - \tilde{\theta}_{jk}^{\mathrm{MT,opt}}(\theta^{\mathrm{FM,opt}}))^2 c_k.
\end{align*}
$$

Since for all $k \in U_T$ the link weights $\tilde{\theta}_{jk}^{\mathrm{MT,opt}}(\theta^{\mathrm{FM,opt}})$ and $\tilde{\theta}_{jk}^{\mathrm{MT,alt}}$ between $j \in U_M$ and $k \in U_T$ are standardized, we have $\sum_{j \in U_M} (\tilde{\theta}_{jk}^{\mathrm{MT,alt}} - \tilde{\theta}_{jk}^{\mathrm{MT,opt}}(\theta^{\mathrm{FM,opt}})) = 0$. Then

$$
\mathrm{Var}(\hat{t}_{y2}^{\mathrm{alt}}) - \mathrm{Var}(\hat{t}_{y2}^{\mathrm{opt}}) = \sum_{k \in U_T} y_k^2 \sum_{j \in U_M} (\tilde{\theta}_{jk}^{\mathrm{MT,alt}} - \tilde{\theta}_{jk}^{\mathrm{MT,opt}}(\theta^{\mathrm{FM,opt}}))^2 \left( \frac{1}{\mathbf{1}_j^{\mathrm{Mt}} (\mathbf{\Delta}_j^M)^{-1} \mathbf{1}_j^M} - c_k \right).
$$

□

**Proof of Corollary 3.3.4**

To prove the corollary, we have to consider estimator of the mean over $U_T$, since the total
over $U_T$ may not converge. Let us denote $f_N = \frac{n_F}{N_F}$. We consider the asymptotic framework
considered in Isaki and Fuller (1982) for $U_F$, $U_M$ and $U_T$ simultaneously, with MtO-MtO
links between populations. More specifically, we assume that the population sizes of $N_F$, $N_M$
and $N_T$, and the sample sizes $n_F$, $n_M$ and $n_T$ go to infinity while the size of $N_{Mk}$ for all
$k \in U_T$ and $N_{Fj}$ for all $j \in U_M$ is bounded. In the example of La Poste survey, this entails
that the number of boxes in a round or the number of addresses in a box does not go to
infinity. Every time a new round is added to $U_T$, all boxes and addresses linked to this round
are added to $U_M$ and $U_F$ respectively without impacting the preexisting rounds.

We assume the following :

$$
\mathbf{H1} \quad \lim_{N \to \infty} f_N = \pi \in (0, 1),
$$

$$
\textbf{H2} \quad N_F, N_T, N_M, n_F \text{ and } n_T > 1 \text{ and that } \frac{N_F}{N_T} = O(1) \text{ with the Landau notation } O,
$$
H3 there is a positive constant C such that, for all $N_T$ and all $N_M$, $\frac{1}{N_T} \sum_{k \in U} y_k^2 N_{Mk} < C$.

Corollary 3.3.3 gives the expression of the difference of variances of double GWSM estimator of the mean over $U_T$ for samples in $U_F$ drawn using SRSWOR:

$$ \text{Var}\left(\frac{1}{N_T} \hat{t}_{y2}^{\text{alt}}\right) - \text{Var}\left(\frac{1}{N_T} \hat{t}_{y2}^{\text{opt}}\right) = \sum_{k \in U_T} \frac{1}{N_T^2} y_k^2 \sum_{j \in U_M} \left(\tilde{\theta}_{jk}^{\text{MT,alt}} - \tilde{\theta}_{jk}^{\text{MT,opt}}(\theta^{\text{FM,opt}})\right)^2 \left(\frac{1}{n_F} \frac{N_F - n_F}{N_F - 1} \frac{N_F - N_{Fj}}{N_{Fj}} + \frac{1 - f_N}{f_N} \frac{1}{N_F - 1}\right). $$

To show that we can approximate this expression with

$$ \sum_{k \in U_T} \frac{1}{N_T^2} y_k^2 \sum_{j \in U_M} \left( \tilde{\theta}_{jk}^{\text{MT,alt}} - \tilde{\theta}_{jk}^{\text{MT,opt}}(\theta^{\text{FM,opt}}) \right)^2 \left( \frac{1}{n_F} \frac{N_F - n_F}{N_F - 1} \frac{N_F - N_{Fj}}{N_{Fj}} \right) $$

we have to show that the difference

$$ A_N = \sum_{k \in U_T} \frac{1}{N_T^2} y_k^2 \sum_{j \in U_M} \left( \tilde{\theta}_{jk}^{\text{MT,alt}} - \tilde{\theta}_{jk}^{\text{MT,opt}}(\theta^{\text{FM,opt}}) \right)^2 \left( \frac{1 - f_N}{f_N} \frac{1}{N_F - 1} \right) $$

is such that $n_T|A_N|$ converges to 0.

For all $j \in U_M$ and $k \in U_T$, the standardized weights are smaller than 1 thus $|\tilde{\theta}_{jk}^{\text{MT,alt}} - \tilde{\theta}_{jk}^{\text{MT,opt}}(\theta^{\text{FM,opt}})| < 1$ and $\sum_{j \in U_M} (\tilde{\theta}_{jk}^{\text{MT,alt}} - \tilde{\theta}_{jk}^{\text{MT,opt}}(\theta^{\text{FM,opt}}))^2 < N_{Mk}$. We get that:

$$ \begin{align*} n_T|A_N| &= n_T \sum_{k \in U_T} \frac{1}{N_T^2} y_k^2 \sum_{j \in U_M} \left( \tilde{\theta}_{jk}^{\text{MT,alt}} - \tilde{\theta}_{jk}^{\text{MT,opt}}(\theta^{\text{FM,opt}}) \right)^2 \left( \frac{1 - f_N}{f_N} \frac{1}{N_F - 1} \right) \\ &\le n_T \sum_{k \in U_T} \frac{1}{N_T^2} y_k^2 N_{Mk} \left( \frac{1 - f_N}{f_N} \frac{1}{N_F - 1} \right) \\ &\le \frac{n_T}{N_T} \left( \frac{1}{N_T} \sum_{k \in U_T} y_k^2 N_{Mk} \right) \left( \frac{1 - f_N}{f_N} \frac{1}{N_F - 1} \right) \\ &\le \left( \frac{n_T}{N_T} \frac{1 - f_N}{f_N} \right) \left( \frac{1}{N_T} \sum_{k \in U_T} y_k^2 N_{Mk} \right) \frac{1}{N_F - 1} \\ &\le \left( \frac{1 - f_N}{f_N} \right) \left( \frac{1}{N_T} \sum_{k \in U_T} y_k^2 N_{Mk} \right) \frac{1}{N_F - 1} \end{align*} $$
Thus

$$ \lim_{N \to \infty} n_T |A_N| \le \lim_{N \to \infty} \left( \frac{1 - f_N}{f_N} \right) \frac{N_F}{N_T} \left( \frac{1}{N_T} \sum_{k \in U_T} y_k^2 N_{Mk} \right) \frac{1}{N_F - 1} = 0, $$

using **H1**, **H2** and **H3**, which concludes the proof.

For stratified SRSWOR with H strata of size $N_h$ in $U_F$, $h = 1, \dots, H$, we denote by $h_T$ (resp. $h_M$) the strata of size $N_{Th}$ in $U_T$ (resp. of size $N_{Mh}$ in $U_M$) composed of all units $k \in U_T$ (resp. $j \in U_M$) linked to strata h in $U_F$. Since the links are of type MtO-MtO, the strata in $U_M$ and $U_T$ do not overlap. If **H1**, **H2** and **H3** holds for each stratum in $U_F$, $U_M$ and $U_T$, then the previous reasoning can be done for each stratum, which concludes the proof. □
Intégration statistique de données
(ISD)
<figure><img src="image_51.png" /></figure>
# Chapitre 4

## QR Prediction for Statistical Data
Integration<fnref n="1" />

### Abstract

In this paper, we investigate how a big non-probability database can be used to improve estimates of finite population totals from a small probability sample through data integration techniques. In the situation where the study variable is observed in both data sources, Kim and Tam (2021) proposed two design-consistent estimators that can be justified through dual frame survey theory. First, we provide conditions ensuring that these estimators are more efficient than the Horvitz-Thompson estimator when the probability sample is selected using either Poisson sampling or simple random sampling without replacement. Then, we study the class of QR predictors, introduced by Särndal and Wright (1984), to handle the less common case where the non-probability database contains no study variable but auxiliary variables. We also require that the non-probability database is large and can be linked to the probability sample. We provide conditions ensuring that the QR predictor is asymptotically design-unbiased. We derive its asymptotic design variance and provide a consistent design-based variance estimator. We compare the design properties of different predictors, in the class of QR predictors, through a simulation study. This class includes a model-based predictor, a model-assisted estimator and a cosmetic estimator. In our simulation setups, the cosmetic estimator performed slightly better than the model-assisted estimator. These findings are confirmed by an application to La Poste data, which also illustrates that the properties of the cosmetic estimator are preserved irrespective of the observed non-probability sample.

<footnote n="1">1. Article de E.Medous, C.Goga, A. Ruiz-Gazen, J.F Beaumont, A. Dessertaine et P. Puech, à paraître dans *Techniques d'Enquêtes/Survey Methodology*, 2023</footnote>
**Keywords** : cosmetic estimator, dual frame, GREG estimator, non-probability sample, probability sample, variance estimator.

4.1 Introduction

In the field of economics and social sciences, surveys are usually based on probability
sampling methods. At the French postal service (La Poste) for example, the postal traffic
is estimated through quarterly probability surveys. Controlling the sampling design allows
for design-based inference without resorting to modeling of the study variables; this feature
is attractive to many survey statisticians. Neyman (1934) is usually known as the founding
paper of probability sampling theory. Since then, the literature on this topic has grown
rapidly with an interplay between theory and practice (see Rao (2005) for the most important
contributions).

Recently, survey statisticians have observed a decline in response rates together with an
increase of the survey costs, which make probability sampling more challenging. In addition,
large non-probability samples, such as administrative data or web-based surveys, become
available often at low cost (see, e.g., Beaumont (2020) and Rao (2021) for more details). These
observations are also true at La Poste where, for cost reasons, the size of probability samples
is bound to decrease while a big database containing the automatically processed postal
mail is available. Even if non-probability samples are associated with unknown selection
mechanisms and may suffer from selection bias and measurement errors, they provide timely
information on the population of interest. This context leads survey statisticians to study
the integration or combination of data from probability and non-probability samples.

The literature on data integration in survey sampling has grown rapidly recently, and
the reader may refer to several reviews on the subject (see Beaumont (2020), Yang and Kim
(2020), Rao (2021), Kim (2022) and Wu (2022)). If we focus on the problem of combining
probability and non-probability samples, the different data integration methods can be divi-
ded into three groups depending on whether the study variable is observed in the probability
sample only, in the non-probability sample only, or in both samples (see e.g. Rao (2021)).
Most methods tackle the problem of the study variable observed in the non-probability
sample only, e.g. Kim (2022). In this context, the objective is to address the selection bias by
combining data from the non-probability sample with auxiliary data available in a probability
sample.

At La Poste, the problem is rather that the study variables (such as the different types
of mails sent) are only available in the probability sample whereas auxiliary information is
only available in the non-probability database. Such a context is rather rare in practice and
has therefore not been studied in detail so far. The aim of the present paper is to study this particular context thoroughly. The method we recommend and study in detail is applicable nicely if (i) the overlap between the probability and the non-probability samples is ideally large but at least non-empty, and (ii) it is possible to accurately match observations from both samples. At La Poste, the non-probability sample represents more than 80% of the population; as a result, its intersection with the probability sample is large. The matching, however, is a difficult task that La Poste is still investigating.

In the situation where the study variables are measured in both samples, Kim and Tam (2021) propose a design-based dual frame approach to improve the efficiency of the Horvitz-Thompson estimator (Horvitz and Thompson, 1952), which uses the probability sample only. The total of the study variable over the whole population is estimated by summing the true total over the non-probability sample and an estimator of the total over the complementary of the non-probability sample. Kim and Tam (2021) propose several estimators that can be deduced from a calibration perspective.

In Section 4.2, we revisit the approach of Kim and Tam (2021) and derive general results on the efficiency of their proposed dual frame estimators. In the situation where the study variable is not measured in the non-probability sample, we propose to replace the true unknown total over the non-probability sample by some prediction. In Section 4.3, we adapt the general class of QR predictors, introduced in Wright (1983), to data integration. This class of estimators includes the well-known model-assisted (GREG) and model-based estimators, but also the cosmetic estimator (Särndal and Wright, 1984). We first exhibit a condition under which the QR predictors can be written in a projection form. We then derive a QR condition such that these predictors are identical to model-assisted predictors. In Section 4.4, we look at the asymptotic properties of the QR estimators. We show that they are asymptotically unbiased under the model and the sampling design. We also prove that, under the QR condition, the predictors are asymptotically design-unbiased. We derive their asymptotic design variance and provide a design-consistent variance estimator. In Section 4.5, we use Monte Carlo simulations to compare several QR predictors and show that the cosmetic estimator is a good compromise for several setups. In Section 4.6, we consider an application to La Poste data and illustrate the impact of the non-probability sample on the estimators. Finally we conclude and give perspectives in Section 4.7.

## 4.2 Study variable observed in both samples

We are interested in estimating the population total $T = \sum_{k \in U} y_k$, where $y_k$ is the value of the variable of interest $Y$ for unit $k$ of the population $U$. A probability sample $s_P$ is drawn
from $U$ using a sampling design $p(s_P|Z)$, where the population matrix $Z$ contains design information such as strata identifiers. The sample inclusion indicator, $I_k$, $k \in U$, takes the value 1 if unit $k$ is selected in $s_P$, and 0 otherwise. The probability that a given population unit $k$ is selected in $s_P$ is $\pi_k = E_p(I_k|Z)$. We assume in the present section that the variable of interest $Y$ is observed for each unit of the probability sample but also for each unit in the non-probability sample $s_{NP} \subset U$. The inclusion indicator in $s_{NP}$ for the population unit $k \in U$ is denoted as $\delta_k$ (i.e., $\delta_k = 1$, if $k \in s_{NP}$, and $\delta_k = 0$, otherwise). We assume that $\delta_k$ is available for each unit of the probability sample $s_P$. Denote by $N$ (resp. $N_{NP}$) the size of $U$ (resp. $s_{NP}$) and by $n$ the expected size of $s_P$. Let $\hat{T}_{HT} = \sum_{k \in s_P} d_k y_k$ be the well-known expansion or Horvitz-Thompson estimator of $T$ with the sampling weights $d_k = 1/\pi_k$. If $\pi_k > 0$, for all $k \in U$, $\hat{T}_{HT}$ is a design-unbiased estimator of $T$.

The non-probability sample $s_{NP}$ is usually a cheap and large source of data. Its selection mechanism is unknown, and its selection bias cannot be ignored when making inference. On the other hand, the probability sample $s_P$ is assumed representative (without selection bias), yet often expensive and of relatively small size. By combining information from the two samples, we can expect to find an estimator more precise than the expansion estimator obtained using $s_P$.

Kim and Tam (2021) propose two estimators using combined data from $s_P$ and $s_{NP}$ and we propose to revisit the properties of these estimators. The total can be decomposed as :

$$T = T_{NP} + T_C,$$

where $T_{NP} = \sum_{k \in s_{NP}} y_k = \sum_{k \in U} \delta_k y_k$ and $T_C = \sum_{k \in U - s_{NP}} y_k = \sum_{k \in U} (1 - \delta_k) y_k$. Since $y_k$ is measured for all units of $s_{NP}$, $T_{NP}$ is known, and we only have to estimate $T_C$. Kim and Tam (2021) propose the following estimator :

$$\hat{T}_{DI} = T_{NP} + \sum_{k \in s_P} d_k (1 - \delta_k) y_k, \quad (4.1)$$

where $T_C$ is estimated using the expansion estimator. As pointed out by Beaumont (2020), this can be viewed as a dual frame problem, with frames $U$ and $s_{NP}$, where the sample $s_P$ is randomly selected from $U$ and a census is taken from $s_{NP}$. In this context of two sampling frames, $\hat{T}_{DI}$ is an estimator that results from a direct application of the method proposed by Bankier (1986). One may think that $\hat{T}_{DI}$ is more efficient than $\hat{T}_{HT}$, especially if the size of the non-probability sample is large, but this is not true in general. The following Proposition shows that, while the variance of $\hat{T}_{DI}$ is always smaller than the variance of $\hat{T}_{HT}$ for Poisson sampling, the property is only true under a condition on the study variable for
simple random sampling without replacement.

**Proposition 4.2.1.** (i) For Poisson sampling, the variance of $\hat{T}_{DI}$ is less than or equal to the variance of $\hat{T}_{HT}$.

(ii) For simple random sampling without replacement, the variance of $\hat{T}_{DI}$ is less than or equal to the variance of $\hat{T}_{HT}$ if and only if

$$CV_{NP}^2 \geq -\frac{N_{NP}}{N_{NP}-1} \left(1 + \frac{N_{NP}}{N} - 2\frac{\bar{Y}_U}{\bar{Y}_{NP}}\right),$$

where $\bar{Y}_U = \frac{1}{N} \sum_{k \in U} y_k$ is the mean of $Y$ over $U$, $\bar{Y}_{NP} = \frac{1}{N_{NP}} \sum_{k \in U} \delta_k y_k$ is the mean of $Y$ over $s_{NP}$, and $CV_{NP} = \sqrt{S_{Y,NP}^2/\bar{Y}_{NP}}$ is the coefficient of variation of $Y$ in $s_{NP}$, with $S_{Y,NP}^2 = \frac{1}{N_{NP}-1} \sum_{k \in U} \delta_k (y_k - \bar{Y}_{NP})^2$.

The proof of Proposition 4.2.1 is given in the appendix. Intuitively, the result (ii) of Proposition 4.2.1 can be explained by the fact that the size of $s_P$ is fixed for simple random sampling without replacement in the expression of $\hat{T}_{HT}$ while the size of $s_P \cap U - s_{NP}$ is random for $\hat{T}_{DI}$. In other words, the estimator $\hat{T}_{DI}$ is calibrated on $N_{NP}$ and $T_{NP}$, but not on $N$ while $\hat{T}_{HT}$ is calibrated on $N$.

If the size of the population $U$ is known, Kim and Tam (2021) propose to improve $\hat{T}_{DI}$ by using the following estimator :

$$\hat{T}_{PDI} = T_{NP} + \hat{T}_{C}^{(Ha)},$$

where

$$\hat{T}_{C}^{(Ha)} = (N - N_{NP}) \frac{\sum_{k \in s_P} d_k(1 - \delta_k)y_k}{\sum_{k \in s_P} d_k(1 - \delta_k)}$$

is a Hájek-type estimator of the total $T_C$. Kim and Tam (2021) proved that $\hat{T}_{PDI}$ is a Generalized Regression (GREG) estimator calibrated on $N$, $N_{NP}$ and $T_{NP}$. Its expression can be further generalized by including additional calibration variables.

Following Kim and Tam (2021), it is possible to use the linearization approach and derive the approximate variance of $\hat{T}_{PDI}$, denoted as AVar($\hat{T}_{PDI}$). For Poisson sampling, the independence of the inclusion indicators reduces the comparison of $\hat{T}_{PDI}$ and $\hat{T}_{DI}$ to the comparison of Horvitz-Thompson and Hájek estimators of the total $T_C = \sum_{k \in U} (1 - \delta_k)y_k$. The Hájek estimator can be significantly more efficient than the Horvitz-Thompson estimator but it is not true in general (see, e.g., Särndal et al. (1992))

For Poisson sampling, the estimator $\hat{T}_{PDI}$ can be substantially more efficient than the Horvitz-Thompson estimator $\hat{T}_{HT}$, as illustrated in our simulation study in Section 4.5. For
simple random sampling without replacement, the approximate variance of $\hat{T}_{PDI}$ can be compared to the variance of $\hat{T}_{HT}$ in more general conditions than in Kim and Tam (2021). Proposition 4.2.2 below shows that the approximate variance of $\hat{T}_{PDI}$ is smaller than the variance of $\hat{T}_{HT}$ for simple random sampling without replacement, and gives the expression of the difference between the variances.

**Proposition 4.2.2.** For simple random sampling without replacement,

$$Var(\hat{T}_{HT}) - AVar(\hat{T}_{PDI}) = \frac{N^2(1-f)}{(N-1)n} \left( \sum_{k \in U} \delta_k (y_k - \bar{Y}_U)^2 + \sum_{k \in U} (1 - \delta_k) (\bar{Y}_C - \bar{Y}_U)^2 \right),$$

where $\bar{Y}_U = \frac{1}{N} \sum_{k \in U} y_k$ is the mean of Y over U, and $\bar{Y}_C = \frac{1}{N - N_{NP}} \sum_{k \in U} (1 - \delta_k)y_k$ is the mean of Y over $U - s_{NP}$.

In the present section, the study variable Y is assumed to be measured in both samples, $s_P$ and $s_{NP}$. In the next section, we alleviate this assumption by considering that the study variable is not known in the non-probability sample. This situation is the one encountered at La Poste where not all variables of interest are measured in the automatically processed postal mail. The big non-probability database is based on an image recognition process and covers around 80% of the postal mails. This database contains some relevant auxiliary information such as the departure dates from the sending post office. However, such data are subject to selection bias (e.g., mails with atypical shape are not automatically processed), and measurement errors (e.g., errors in barcode scanning during the image recognition process). In such a situation, we propose to use the intersection between the big database and the probability sample, where the auxiliary variables together with the study variable are available, and predict the unknown $y_k$ for $k \in s_{NP} - s_P$.

## 4.3 Prediction estimators for study variable unobserved in the non-probability sample

Recall that the finite population total of Y can be decomposed as $T = T_{NP} + T_C$. The total $T_C$ is estimated as in Section 4.2 by the Hájek-type estimator $\hat{T}_C^{(Ha)}$. In the present section, $y_k$ is unknown for $k \in s_{NP}$, and contrarily to Section 4.2, the total $T_{NP}$ has to be estimated. In order to do so, we introduce a working model for Y and the general QR class of predictors of $T_{NP}$ that does not require $y_k$ to be known for units in $s_{NP}$. We assume that a vector of auxiliary variables $x_k = (X_{k1}, \dots, X_{kp})^\top$ is available for each unit k of a non-probability sample $s_{NP} \subset U$. We also assume that $\delta_k$ and $\delta_k x_k$ are available for each
unit k of the probability sample $s_P$. Table 4.1 gives a summary of the characteristics of the data we consider in the remainder of this paper.

<table><caption>TABLE 4.1 – Data characteristics in the data integration context of Section 4.3.</caption><thead><tr><th>Sample</th><th>y<sub>k</sub> measured</th><th>δ<sub>k</sub> available</th><th>known selection mechanism</th><th>Auxiliary variables available</th></tr></thead><tbody><tr><td>s<sub>P</sub></td><td>Yes</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>s<sub>NP</sub></td><td>No</td><td>Yes</td><td>No</td><td>Yes</td></tr></tbody></table>

The variable $Y$ is not available in $s_{NP}$ and we cannot use anymore $\hat{T}_{PDI}$ since the total $T_{NP} = \sum_{k \in U} \delta_k y_k$ is unknown. The idea behind the class of estimators introduced in this section is to predict $y_k$ for $k \in s_{NP}$ by using regression modelling between $Y$ and the auxiliary variables, and then predict $T_{NP}$. We assume the following working model between the study variable $Y$ and the vector of auxiliary variables $\mathbf{x}_k$:

$$y_k = \mathbf{x}_k^\top \hat{\beta} + \varepsilon_k, \quad k \in s_{NP}, \qquad (4.2)$$

where the errors $\varepsilon_k$ are independent with expectation $E_m(\varepsilon_k) = 0$ and variance $Var_m(\varepsilon_k)$ proportional to $\nu(\mathbf{x}_k) = v_k$ for some known positive constants $v_k$. The subscript $m$ indicates that the expectation and variance are taken with respect to model (4.2) conditionally on observed auxiliary variables $\mathbf{x}_k$, $k \in s_{NP}$. Note that model (4.2) only needs to hold for units in the non-probability sample. A model for $Y$ does not need to be explicitly specified for units $k \in U - s_{NP}$ as we always make inferences conditional on $y_k$, $k \in U - s_{NP}$.

We define a predictor $\hat{y}_k$ of $y_k$ for $k \in s_{NP}$ by $\hat{y}_k = \mathbf{x}_k^\top \hat{\beta}$ with

$$\hat{\beta} = \left( \sum_{k \in s_P} q_k \delta_k \mathbf{x}_k \mathbf{x}_k^\top \right)^{-1} \left( \sum_{k \in s_P} q_k \delta_k \mathbf{x}_k y_k \right), \qquad (4.3)$$

where $q_k$ are known positive constants for $k \in s_{NP}$. We assume that the $p \times p$ dimensional matrix $\sum_{k \in s_P} q_k \delta_k \mathbf{x}_k \mathbf{x}_k^\top$ is nonsingular for all possible samples $s_P$.

We propose to estimate $T_{NP} = \sum_{k \in U} \delta_k y_k$ by a QR predictor as suggested in Wright (1983):

$$\begin{aligned} \hat{T}_{NP}^{(QR)} &= \sum_{k \in U} \delta_k \hat{y}_k + \sum_{k \in s_P} r_k \delta_k (y_k - \hat{y}_k) \\ &= \sum_{k \in U} \delta_k \mathbf{x}_k^\top \hat{\beta} + \sum_{k \in s_P} r_k \delta_k (y_k - \mathbf{x}_k^\top \hat{\beta}), \end{aligned} \qquad (4.4)$$

where $r_k \ge 0$ are predefined constants. The initials Q and R refer to the constants $q_k$ and
$r_k$. The final estimator of $T$ is then given by

$$ \hat{T}^{(\text{QR})} = \hat{T}_{NP}^{(\text{QR})} + \hat{T}_{C}^{(\text{Ha})}. \qquad (4.5) $$

Various choices of $q_k$ and $r_k$ yield predictors $\hat{T}_{NP}^{(\text{QR})}$ with familiar forms as detailed below.

1. For $q_k = d_k v_k^{-1}$ and $r_k = d_k$, we obtain the model-assisted or GREG-type estimator :

$$ \hat{T}_{NP}^{(\text{MA})} = \sum_{k \in U} \delta_k \hat{y}_k^{(\text{MA})} + \sum_{k \in s_P} \delta_k d_k (y_k - \hat{y}_k^{(\text{MA})}), $$

where $\hat{y}_k^{(\text{MA})} = \mathbf{x}_k^{\top} \hat{\beta}^{(\text{MA})}$ with $\hat{\beta}^{(\text{MA})} = (\sum_{k \in s_P} d_k v_k^{-1} \delta_k \mathbf{x}_k \mathbf{x}_k^{\top})^{-1} (\sum_{k \in s_P} d_k v_k^{-1} \delta_k \mathbf{x}_k y_k)$.

2. For $q_k = v_k^{-1}$ and $r_k = 1$, we obtain the model-based type estimator :

$$ \hat{T}_{NP}^{(\text{MB})} = \sum_{k \in U} \delta_k \hat{y}_k^{(\text{MB})} + \sum_{k \in s_P} \delta_k (y_k - \hat{y}_k^{(\text{MB})}), $$

where $\hat{y}_k^{(\text{MB})} = \mathbf{x}_k^{\top} \hat{\beta}^{(\text{MB})}$ with $\hat{\beta}^{(\text{MB})} = (\sum_{k \in s_P} \delta_k v_k^{-1} \mathbf{x}_k \mathbf{x}_k^{\top})^{-1} (\sum_{k \in s_P} \delta_k v_k^{-1} \mathbf{x}_k y_k)$.

3. For $q_k = (d_k - 1)v_k^{-1}$ and $r_k = 1$, we obtain the cosmetic-type estimator (Särndal and Wright, 1984; Brewer, 1999) :

$$ \hat{T}_{NP}^{(\text{Cos})} = \sum_{k \in U} \delta_k \hat{y}_k^{(\text{Cos})} + \sum_{k \in s_P} \delta_k (y_k - \hat{y}_k^{(\text{Cos})}), $$

where $\hat{y}_k^{(\text{Cos})} = \mathbf{x}_k^{\top} \hat{\beta}^{(\text{Cos})}$ with

$$ \hat{\beta}^{(\text{Cos})} = \left( \sum_{k \in s_P} (d_k - 1) v_k^{-1} \delta_k \mathbf{x}_k \mathbf{x}_k^{\top} \right)^{-1} \left( \sum_{k \in s_P} (d_k - 1) v_k^{-1} \delta_k \mathbf{x}_k y_k \right). $$

Let us derive some properties for this class of QR predictors. Proposition 4.3.1 gives a general condition on the constants $q_k$ and $r_k$ such that the QR predictor can be defined as a sum of predictions over the population. Proposition 4.3.2 gives another general condition on the constants $q_k$ and $r_k$ such that the QR predictor is a model-assisted type estimator. The proofs are given in the Appendix.

**Proposition 4.3.1.** (projection form) Consider the QR predictor $\hat{T}_{NP}^{(\text{QR})}$ given by (4.4). Under the condition that there exists a vector $\mu \in \mathbf{R}^p$ such that

$$ (\text{Proj}) : \quad \mu^{\top} \mathbf{x}_k q_k = r_k \quad \text{for all } k \in s_{NP}, \qquad (4.6) $$
we have $\sum_{k \in s_P} r_k \delta_k (y_k - \hat{y}_k) = 0$. In this case, $\hat{T}_{NP}^{(\text{QR})}$ can be written in the projection form :

$$ \hat{T}_{NP}^{(\text{QR})} = \sum_{k \in U} \delta_k \hat{y}_k. $$

The model-assisted estimator $\hat{T}_{NP}^{(\text{MA})}$ and model-based estimator $\hat{T}_{NP}^{(\text{MB})}$ satisfy Condition (Proj) if there exists a vector $\boldsymbol{\mu} \in \mathbf{R}^p$ such that $\boldsymbol{\mu}^{\top} \boldsymbol{x}_k = v_k$ for all $k \in s_{NP}$. This condition is satisfied when $v_k$ is one of the auxiliary variables in the model. If $v_k = 1$, it is satisfied provided that the intercept is included in the model. Condition (Proj) holds for $\hat{T}_{NP}^{(\text{Cos})}$ if $\boldsymbol{\mu}^{\top} \boldsymbol{x}_k = v_k(d_k - 1)^{-1}$ for all $k \in s_{NP}$. A consequence of Proposition 4.3.1 is that, for equal probability sampling design such as simple random sampling without replacement, the model-assisted, the model-based and the cosmetic estimators are all equal.

Using Theorem 2 from Wright (1983), we derive the following proposition. For $r_k$ satisfying Condition (QR) below and any given $q_k$, the QR predictor of $T_{NP}$ is identical to the model-assisted predictor of $T_{NP}$ with the same $q_k$.

**Proposition 4.3.2.** Suppose that the constants $r_k$ and $q_k$ are such that there exists some vector $\boldsymbol{\lambda} \in \mathbf{R}^p$ such that

$$ (QR): \quad 1 - \pi_k r_k = \pi_k q_k \boldsymbol{x}_k^{\top} \boldsymbol{\lambda} \quad \text{for all } k \in s_{NP}. \qquad (4.7) $$

Then :

$$ \hat{T}_{NP}^{(\text{QR})} = \hat{T}_{NP}^{(\text{Q}\pi)}, $$

where

$$ \hat{T}_{NP}^{(\text{Q}\pi)} = \sum_{k \in U} \delta_k \boldsymbol{x}_k^{\top} \hat{\boldsymbol{\beta}} + \sum_{k \in s_P} d_k \delta_k (y_k - \boldsymbol{x}_k^{\top} \hat{\boldsymbol{\beta}}) \qquad (4.8) $$

is the model-assisted type predictor of $T_{NP}$ with $\hat{\boldsymbol{\beta}}$ given by (4.3).

Following Wright (1983), we note that the (QR) condition always holds for $\hat{T}_{NP}^{(\text{MA})}$. This condition also holds for the model-based estimator $\hat{T}_{NP}^{(\text{MB})}$ if and only if there exists a vector $\boldsymbol{\lambda} \in \mathbf{R}^p$ such that $v_k(d_k - 1) = \boldsymbol{x}_k^{\top} \boldsymbol{\lambda}$, for all $k \in s_{NP}$. This condition is true if we take $v_k(d_k - 1)$ among the auxiliary variables $\boldsymbol{x}_k$. Condition (QR) holds for the cosmetic estimator $\hat{T}_{NP}^{(\text{Cos})}$ if and only if there exists a vector $\boldsymbol{\lambda} \in \mathbf{R}^p$ such that $v_k = \boldsymbol{x}_k^{\top} \boldsymbol{\lambda}$, for all $k \in s_{NP}$. This condition is true if $v_k$ is included in the vector of auxiliary variables.
## 4.4 Asymptotic properties and variance estimation of QR predictors

Let us consider the class of QR predictors $\hat{T}^{(QR)}$ given in (4.5) and start by studying the prediction error $\hat{T}^{(QR)} - T$ under the model (4.2) and the sampling design $p(\cdot)$. We have

$$ \hat{T}^{(QR)} - T = (\hat{T}_{NP}^{(QR)} - T_{NP}) + (\hat{T}_{C}^{(Ha)} - T_{C}). $$

The first right-hand term depends on the model and the sampling design, while the second right-hand term only depends on the sampling design. Assuming that the sampling design is not informative with respect to the model (4.2), we can prove that the expectation of $\hat{T}_{NP}^{(QR)} - T_{NP}$ computed with respect to the model is equal to 0. Indeed, the model bias of $\hat{T}_{NP}^{(QR)}$ is given by :

$$ E_m(\hat{T}_{NP}^{(QR)} - T_{NP}) = \sum_{k \in U} \delta_k E_m(\boldsymbol{x}_k^{\top} \hat{\boldsymbol{\beta}} - y_k) + \sum_{k \in s_P} r_k \delta_k E_m(y_k - \boldsymbol{x}_k^{\top} \hat{\boldsymbol{\beta}}), $$

with $\hat{\boldsymbol{\beta}} = (\sum_{k \in s_P} q_k \delta_k \boldsymbol{x}_k \boldsymbol{x}_k^{\top})^{-1} (\sum_{k \in s_P} q_k \delta_k \boldsymbol{x}_k y_k)$. Under the model (4.2), $E_m(y_k) = \boldsymbol{x}_k^{\top} \hat{\boldsymbol{\beta}}$ for all $k \in s_{NP}$, $E_m(\hat{\boldsymbol{\beta}}) = \hat{\boldsymbol{\beta}}$ and $E_m(\boldsymbol{x}_k^{\top} \hat{\boldsymbol{\beta}} - y_k) = 0$, implying that

$$ E_m(\hat{T}_{NP}^{(QR)} - T_{NP}) = 0. \qquad (4.9) $$

The estimator $\hat{T}_C^{(Ha)}$ is a Hájek-type estimator of $T_C$ and is not design-unbiased for $T_C$. Following Särndal (1980), we rather look at its asymptotic design-unbiasedness. An estimator $\hat{T}$ is said to be asymptotically design-unbiased for the finite population total $T$ if $\lim_{N \to \infty} N^{-1}[E_p(\hat{T}) - T] = 0$, where $E_p(\cdot)$ is the expectation computed with respect to the sampling design. The asymptotic framework from Isaki and Fuller (1982) can be considered to allow for the population and sample sizes to grow to infinity. Assuming that the probability of observing an empty intersection set $s_P \cap s_{NP}$ is negligible, then $\hat{T}_C^{(Ha)}$ is asymptotically design-unbiased for $T_C$. Combining this property with relation (4.9), and considering a non-informative sampling design, we get that $\hat{T}^{(QR)}$ is asymptotically mp-unbiased for $T$ under model (4.2).
### 4.4.1 Bias properties of $\hat{T}^{(Q\pi)}$

Let us consider now the QR class of predictors that satisfy the (QR) condition given by (4.7). For this class of predictors, the final estimator of $T$ is

$$ \hat{T}^{(Q\pi)} = \hat{T}_{NP}^{(Q\pi)} + \hat{T}_{C}^{(Ha)} $$

and the standardized total error is given by:

$$ \frac{1}{N} (\hat{T}^{(Q\pi)} - T) = \frac{1}{N} (\hat{T}_{NP}^{(Q\pi)} - T_{NP}) + \frac{1}{N} (\hat{T}_{C}^{(Ha)} - T_{C}). $$

The estimator $\hat{T}^{(Q\pi)}$ is not exactly design-unbiased because of the nonlinearity of $\hat{\beta}$ and of the Hájek estimator $\hat{T}_{C}^{(Ha)}$. Wright (1983) proved that the (QR) condition given in Proposition 4.3.2 is a sufficient condition for $\hat{T}_{NP}^{(Q\pi)}$ to be asymptotically design-unbiased for $T_{NP}$, provided that

$$ \lim_{N \to \infty} \frac{1}{N} \mathrm{E}_p \left[ \left( \sum_{k \in U} \delta_k \mathbf{x}_k - \sum_{k \in s_P} d_k \delta_k \mathbf{x}_k \right)^\top (\hat{\beta} - \tilde{\beta}) \right] = 0, \quad (4.10) $$

where $\tilde{\beta} = (\sum_{k \in U} \pi_k q_k \delta_k \mathbf{x}_k \mathbf{x}_k^\top)^{-1} \sum_{k \in U} \pi_k q_k \delta_k \mathbf{x}_k y_k$, and assuming that $\sum_{k \in U} \pi_k q_k \delta_k \mathbf{x}_k \mathbf{x}_k^\top$ is nonsingular. Following Breidt and Opsomer (2000), if the sampling fraction $n/N$ converges to a constant different from 0, assuming mild conditions on the first and second-order inclusion probabilities of the sampling design, and on the auxiliary information vectors $\mathbf{x}_k$ for all $k \in s_{NP}$, it can be shown that:

$$ \lim_{N \to \infty} \mathrm{E}_p \left\| N^{-1} \left( \sum_{k \in U} \delta_k \mathbf{x}_k - \sum_{k \in s_P} d_k \delta_k \mathbf{x}_k \right) \right\|^2 = 0, $$

where $\|\cdot\|$ is the usual Euclidian norm. Equation (4.10) follows by assuming that the regression coefficient estimator satisfies $\lim_{N \to \infty} \mathrm{E}_p \|\hat{\beta} - \tilde{\beta}\|^2 = 0$ (see Cardot et al. (2013) for more details). The estimator $\hat{T}_{C}^{(Ha)}$ is a Hájek-type estimator which can be shown to be asymptotically design-unbiased for $T_C$ if the probability to observe the empty set for $s_P \cap s_{NP}$ is negligible. We conclude that the QR predictor $\hat{T}^{(Q\pi)}$ is asymptotically design-unbiased for $T$.

### 4.4.2 Asymptotic variance and variance estimation of $\hat{T}^{(Q\pi)}$

Because the QR estimator $\hat{T}^{(Q\pi)}$ is asymptotically design-unbiased, we estimate its asymptotic design variance rather than its design mean square error. We can write the standardized
total error as

$$ \frac{1}{N}(\hat{T}^{(Q\pi)} - T) = \frac{1}{N}\left(\sum_{k \in s_p} d_k(E_k + e_k) - \sum_{k \in U} (E_k + e_k)\right) + R_1 + R_2, $$

where

$$ E_k = \delta_k(y_k - \mathbf{x}_k^\top \tilde{\beta}), \quad e_k = (1 - \delta_k) \left( y_k - \frac{\sum_{k' \in U} (1 - \delta_{k'}) y_{k'}}{N - N_{NP}} \right), $$

$$ R_1 = -\frac{1}{N} \left( \sum_{k \in s_p} d_k \delta_k \mathbf{x}_k - \sum_{k \in U} \delta_k \mathbf{x}_k \right)^\top (\hat{\beta} - \tilde{\beta}) $$

and

$$ R_2 = \left( \frac{1}{N} \sum_{k \in s_p} d_k (1 - \delta_k) \right)^{-1} \frac{1}{N} \left( N - N_{NP} - \sum_{k \in s_p} d_k (1 - \delta_k) \right) \frac{1}{N} \left( \sum_{k \in s_p} d_k e_k - \sum_{k \in U} e_k \right). $$

As in Section 4.4.1, we assume the usual conditions on the sampling fraction $n/N$, on the first and second-order inclusion probabilities, on the variable of interest $y_k$ and on the auxiliary information vector $\mathbf{x}_k$. We use the Landau notations big $O_p$ and little $o_p$. If $||\hat{\beta} - \tilde{\beta}||^2 = o_p(1)$ and $(\sum_{k \in s_p} d_k(1 - \delta_k)/N)^{-1} = O_p(1)$, then $R_1 = o_p(n^{-1/2})$, $R_2 = O_p(n^{-1})$ and the standardized total error can be approximated by

$$ \frac{1}{N}(\hat{T}^{(Q\pi)} - T) \approx \frac{1}{N} \left( \sum_{k \in s_p} d_k(E_k + e_k) - \sum_{k \in U} (E_k + e_k) \right), $$

where the right-hand side of the above expression is of order $O_p(n^{-1/2})$. Since $\sum_{k \in s_p} d_k(E_k + e_k)$ is the Horvitz-Thompson estimator of the total $\sum_{k \in U} (E_k + e_k)$, the asymptotic variance of the QR estimator $\hat{T}^{(Q\pi)}$ is given by :

$$ AVar(\hat{T}^{(Q\pi)}) = \sum_{k \in U} \sum_{l \in U} \Delta_{kl} d_k d_l (E_k + e_k) (E_l + e_l). $$

Assuming that $\pi_{kl} > 0$ for all $k, l \in U$, an estimator of the asymptotic variance is given by

$$ \hat{V}(\hat{T}^{(Q\pi)}) = \sum_{k \in s_p} \sum_{l \in s_p} \frac{\Delta_{kl}}{\pi_{kl}} d_k d_l (\hat{E}_k + \hat{e}_k) (\hat{E}_l + \hat{e}_l), $$

where $\hat{E}_k = \delta_k(y_k - \mathbf{x}_k^\top \tilde{\beta})$ and $\hat{e}_k = (1 - \delta_k)(y_k - \sum_{k' \in s_p} d_{k'}(1 - \delta_{k'})y_{k'}/(N - N_{NP}))$. We can show, under assumptions detailed in Breidt and Opsomer (2000) and Goga et al. (2009),
that this estimator is design-consistent for AVar($\hat{T}^{(Q\pi)}$) in the sense that $N^{-2}n(\hat{V}(\hat{T}^{(Q\pi)}) - AVar(\hat{T}^{(Q\pi)})) = o_p(1)$.

# 4.5 Simulations

In this section, we show the results of a Monte-Carlo study that compares the efficiency of three special cases of the QR predictor, $\hat{T}^{(QR)} = \hat{T}_{NP}^{(QR)} + \hat{T}_{C}^{(Ha)}$, given in Section 4.3, namely the model-assisted, the model-based and the cosmetic estimators, assuming that $v_k = 1$ in model (4.2). We also compare these estimators with the expansion estimator and the PDI estimator defined in Section 4.2. To illustrate that the relative superiority of estimators depends on the data structure, we define three different setups based on different artificial populations. As mentioned in Section 4.3, if the probability samples are drawn using simple random sampling without replacement, the three QR estimators are all equal. Therefore, we focus on Poisson sampling with inclusion probabilities proportional to an auxiliary variable.

## 4.5.1 Populations and setups

The variables are generated using Gamma distributions to ensure their positiveness. Similar simulation results were obtained with Gaussian distributions but are not reported below. All populations have a size $N = 1,000$. We generate two auxiliary variables $X_1$ and $X_2$, where $X_1$ (resp $X_2$) follows a Gamma distribution with mean $\nu_1 = 20$ (resp $\nu_2 = 30$) and standard deviation (Std) $\sigma_1 = 15$ (resp $\sigma_2 = 20$). We use different models to generate the variable $Y$ for all population units. For each model, $Y|X_1, X_2$ follows a Gamma distribution with mean $\mu_{Y|X_1, X_2}$ and constant variance $\sigma_{Y|X_1, X_2}^2$, which depend on the model.

1. For Model 1, $\mu_{Y|X_1, X_2}$ is a linear function of $X_1$ and $X_2$:

$$ \mu_{Y|X_1, X_2} = a_0 + a_1 X_1 + a_2 X_2. $$

2. For Model 2, $\mu_{Y|X_1, X_2}$ is a quadratic function of $X_1$ and a linear function of $X_2$:

$$ \mu_{Y|X_1, X_2} = b_0 + b_1(X_1 - \bar{X}_1)^2 + b_2 X_2 \text{ with } \bar{X}_1 \text{ the mean of } X_1 \text{ over } U. $$

3. For Model 3, $\mu_{Y|X_1, X_2}$ is a linear function of $X_2$:

$$ \mu_{Y|X_1, X_2} = c_0 + c_1 X_2. $$

To make the results comparable between the three models, we determine the constants $a_0$,
$a_1, a_2, b_0, b_1, b_2, c_0, c_2$, and $\sigma_{Y|X_1,X_2}^2$ in such a way that the following characteristics are the same :

- the unconditional mean $\mu$ and variance $\sigma^2$ of the variable $Y$,

- the coefficient of determination of the model, denoted as $R^2$,

- the ratio of variances for the explanatory variables :

$$ \gamma = \text{Var}(a_1 X_1) / \text{Var}(a_2 X_2) = \text{Var}(b_1 (X_1 - \bar{X}_1)^2) / \text{Var}(b_2 X_2). $$

This ratio is only relevant for models 1 and 2 since $X_1$ is not included in Model 3.

In the following, we set $\mu = 100$, $\sigma^2 = 100$, and $\gamma = 0.5$. In Section 4.5.2, the $R^2$ value is either fixed to 0.8 or varies between 0.1 and 0.96. The main characteristics of the three population models are summarized in Table 4.2.

<table><caption>TABLE 4.2 – Population models with $\mu = 100$, $\sigma^2=100$, and $\gamma = 0.5$.</caption><thead><tr><th>Model</th><th>Mean of (X<sub>1</sub>, X<sub>2</sub>)</th><th>Std of (X<sub>1</sub>, X<sub>2</sub>)</th><th>Mean of Y|X<sub>1</sub>, X<sub>2</sub></th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>1</td><td rowspan="3">(20,30)</td><td rowspan="3">(15,20)</td><td>μ<sub>Y</sub> = a<sub>0</sub> + a<sub>1</sub>X<sub>1</sub> + a<sub>2</sub>X<sub>2</sub></td><td>equal</td></tr><tr><td>2</td><td>μ<sub>Y</sub> = b<sub>0</sub> + b<sub>1</sub>(X<sub>1</sub> - X̄<sub>1</sub>)<sup>2</sup> + b<sub>2</sub>X<sub>2</sub></td><td>between</td></tr><tr><td>3</td><td>μ<sub>Y</sub> = c<sub>0</sub> + c<sub>2</sub>X<sub>2</sub></td><td>populations</td></tr></tbody></table>

A non-probability sample of size 900 is drawn using simple random sampling without replacement and is the same for all populations. The probability samples are drawn using Poisson sampling with expected size 200 or 50 and probabilities proportional to $X_1$. We consider three setups. In each setup, we generate $Y|X_1, X_2$ using one of the three different population models, and we compute $\hat{y}_k, k \in s_{NP}$ for different QR predictors. The variables used as explanatory variables in the prediction models differ between setups as follows :

1. Setup 1 : Informative case. Population Model 1 is used to generate population $Y$ values and only $X_2$ is used as explanatory variable in the prediction model along with the intercept.

2. Setup 2 : Quadratic case. Population Model 2 is used to generate population $Y$ values and both auxiliary variables $X_1$ and $X_2$ are used as explanatory variables in the prediction model along with the intercept.

3. Setup 3 : Non-informative case. Population Model 3 is used to generate population $Y$ values and only $X_2$ is used as explanatory variable in the prediction model along with the intercept.

For the informative and quadratic setups, the prediction model differs from the population model for $Y$, while the correct model is used in the non-informative setup. Table 4.3 gives a summary of the three setups.
<table><caption>TABLE 4.3 – Studied setups.</caption><thead><tr><th>Setup</th><th>Population</th><th>Variables used in prediction</th><th>Model correctly specified</th></tr></thead><tbody><tr><td>Informative</td><td>μ<sub>Y</sub> = a<sub>0</sub> + a<sub>1</sub>X<sub>1</sub> + a<sub>2</sub>X<sub>2</sub></td><td><b>x</b><sub>k</sub><sup>T</sup> = (1, x<sub>2k</sub>)</td><td>No</td></tr><tr><td>Quadratic</td><td>μ<sub>Y</sub> = b<sub>0</sub> + b<sub>1</sub>(X<sub>1</sub> - X&#773;<sub>1</sub>)<sup>2</sup> + b<sub>2</sub>X<sub>2</sub></td><td><b>x</b><sub>k</sub><sup>T</sup> = (1, x<sub>1k</sub>, x<sub>2k</sub>)</td><td>No</td></tr><tr><td>Non-informative</td><td>μ<sub>Y</sub> = c<sub>0</sub> + c<sub>2</sub>X<sub>2</sub></td><td><b>x</b><sub>k</sub><sup>T</sup> = (1, x<sub>2k</sub>)</td><td>Yes</td></tr></tbody></table>

## 4.5.2 Results

Let us consider the three setups defined above and compare the following estimators :

$$— \hat{T}_{HT} = \sum_{k \in s_P} d_k y_k,$$

$$— \hat{T}_{PDI} = T_{NP} + (N - N_{NP}) \frac{\sum_{k \in s_P} d_k (1 - \delta_k) y_k}{\sum_{k \in s_P} d_k (1 - \delta_k)},$$

$$— \hat{T}^{(\text{MB})} = \sum_{k \in U} \delta_k \hat{y}_k^{(\text{MB})} + \sum_{k \in s_P} \delta_k (y_k - \hat{y}_k^{(\text{MB})}) + (N - N_{NP}) \frac{\sum_{k \in s_P} d_k (1 - \delta_k) y_k}{\sum_{k \in s_P} d_k (1 - \delta_k)},$$

$$— \hat{T}^{(\text{MA})} = \sum_{k \in U} \delta_k \hat{y}_k^{(\text{MA})} + \sum_{k \in s_P} \delta_k d_k (y_k - \hat{y}_k^{(\text{MA})}) + (N - N_{NP}) \frac{\sum_{k \in s_P} d_k (1 - \delta_k) y_k}{\sum_{k \in s_P} d_k (1 - \delta_k)},$$

$$— \hat{T}^{(\text{Cos})} = \sum_{k \in U} \delta_k \hat{y}_k^{(\text{Cos})} + \sum_{k \in s_P} \delta_k (y_k - \hat{y}_k^{(\text{Cos})}) + (N - N_{NP}) \frac{\sum_{k \in s_P} d_k (1 - \delta_k) y_k}{\sum_{k \in s_P} d_k (1 - \delta_k)}.$$

For each setup, $L = 10,000$ probability samples $s_P$ are drawn according to Poisson sampling as detailed above and several Monte Carlo measures are computed. We compute the Monte Carlo relative bias of a given estimator $\hat{T}$ (either $\hat{T}_{HT}$, $\hat{T}^{(\text{MB})}$, $\hat{T}^{(\text{MA})}$, $\hat{T}^{(\text{Cos})}$ or $\hat{T}_{PDI}$) as

$$\text{RB}_{MC}(\hat{T}) = 100 \times \frac{1}{L} \sum_{l=1}^{L} \frac{\hat{T}^{(l)} - T}{T},$$

where $\hat{T}^{(l)}$ is an estimate of $T$ computed for the $l$-th sample, $l = 1, \dots, L$.

As a measure of efficiency, we compute the Monte Carlo relative mean square error (RMSE) of an estimator $\hat{T}$ (relative to $\hat{T}^{(\text{Cos})}$):

$$\text{RMSE}_{MC}(\hat{T}) = 100 \times \frac{\text{MSE}_{MC}(\hat{T})}{\text{MSE}_{MC}(\hat{T}^{(\text{Cos})})},$$

where

$$\text{MSE}_{MC}(\hat{T}) = L^{-1} \sum_{l=1}^{L} (\hat{T}^{(l)} - T)^2.$$

We also compute the Monte Carlo relative variance (RVar) of an estimator $\hat{T}$ (relative to
<table><caption>TABLE 4.4 – Relative bias (in % of the true value), percent relative (to $\hat{T}^{(\text{Cos})}$) variance and MSE of the different estimators for the three different setups; the expected size of the probability sample is 200 and the size of the non-probability sample is 900.</caption><thead><tr><th>Population parameters</th><th>Setup</th><th>Monte Carlo measures</th><th>T&#x0302;<sub>HT</sub></th><th>T&#x0302;<sup>(MB)</sup></th><th>T&#x0302;<sup>(MA)</sup></th><th>T&#x0302;<sup>(Cos)</sup></th><th>T&#x0302;<sub>PDI</sub></th></tr></thead><tbody><tr><td rowspan="9">μ = 100<br>σ² = 100<br>R² = 0.8<br>γ = 0.5</td><td rowspan="3">Setup 1</td><td>RB<sub>MC</sub></td><td>-0.13</td><td>3.34</td><td>0.11</td><td>0.11</td><td>0.03</td></tr><tr><td>RVar<sub>MC</sub></td><td>23,566.93</td><td>55.62</td><td>114.06</td><td>100.00</td><td>20.97</td></tr><tr><td>RMSE<sub>MC</sub></td><td>22,897.58</td><td>2,715.21</td><td>113.91</td><td>100.00</td><td>20.65</td></tr><tr><td rowspan="3">Setup 2</td><td>RB<sub>MC</sub></td><td>-0.07</td><td>-1.65</td><td>-0.06</td><td>-0.05</td><td>0.02</td></tr><tr><td>RVar<sub>MC</sub></td><td>36,947.99</td><td>84.94</td><td>118.21</td><td>100.00</td><td>23.17</td></tr><tr><td>RMSE<sub>MC</sub></td><td>36,638.27</td><td>1,056.44</td><td>118.42</td><td>100.00</td><td>23.15</td></tr><tr><td rowspan="3">Setup 3</td><td>RB<sub>MC</sub></td><td>0.03</td><td>-0.01</td><td>0.01</td><td>0.01</td><td>0.01</td></tr><tr><td>RVar<sub>MC</sub></td><td>41,088.93</td><td>58.38</td><td>100.49</td><td>100.00</td><td>33.47</td></tr><tr><td>RMSE<sub>MC</sub></td><td>41,080.51</td><td>58.39</td><td>100.48</td><td>100.00</td><td>33.51</td></tr></tbody></table>

$\hat{T}^{(\text{Cos})}$):

$$RVar_{MC}(\hat{T}) = 100 \times \frac{\text{Var}_{MC}(\hat{T})}{\text{Var}_{MC}(\hat{T}^{(\text{Cos})})},$$

where

$$\text{Var}_{MC}(\hat{T}) = L^{-1} \sum_{l=1}^{L} (\hat{T}^{(l)})^2 - \left( L^{-1} \sum_{l=1}^{L} \hat{T}^{(l)} \right)^2.$$

Table 4.4 contains the simulation results for the three setups when $R^2 = 0.8$. In all setups, we confirm that both $\hat{T}_{PDI}$ and $\hat{T}_{HT}$ have a small Monte Carlo bias, as expected. In terms of MSE, $\hat{T}_{PDI}$ is the most precise estimator, while $\hat{T}_{HT}$ is the least precise estimator among all estimators. This result is expected since the expansion estimator does not make use of any auxiliary information, while $\hat{T}_{PDI}$ takes into account the true values of the study variable $y_k$ for $k \in s_{NP}$; i.e., it takes into account the true values of Y for 900 units out of the 1,000 population units. In our context, where the study variable is not observed in $s_{NP}$, the estimator $\hat{T}_{PDI}$ is however not computable and serves more as a gold standard.

The Monte Carlo bias of $\hat{T}^{(\text{MA})}$ and $\hat{T}^{(\text{Cos})}$ is negligible in the three setups whereas $\hat{T}^{(\text{MB})}$ is biased in the informative and quadratic setups. In these two setups, the prediction model differs from the population model used to generate Y values. In the non-informative setup, where the prediction model is correctly specified, the bias of $\hat{T}^{(\text{MB})}$ is also negligible. The estimator $\hat{T}^{(\text{MA})}$ has the largest variance of the QR predictors in the informative and quadratic setups, while $\hat{T}^{(\text{MB})}$ has the smallest variance in all setups. In the quadratic setup, the variance of $\hat{T}^{(\text{MB})}$ is similar to the variance of $\hat{T}^{(\text{Cos})}$ but $\hat{T}^{(\text{MB})}$ has the highest MSE amongst the QR predictors in both informative and quadratic setups. This means that the bias of $\hat{T}^{(\text{MB})}$ degrades its MSE a lot despite its small variance. In the non-informative setup,
$\hat{T}^{(\text{MB})}$ has the lowest MSE amongst the QR predictors. We can see in Table 4.4 that this comes from the absence of bias for $\hat{T}^{(\text{MB})}$ in this setup together with its small variance.

In the informative and quadratic setups, $\hat{T}^{(\text{Cos})}$ is more precise in term of variance than $\hat{T}^{(\text{MA})}$. The estimators $\hat{T}^{(\text{MA})}$ and $\hat{T}^{(\text{Cos})}$ are similar in the non-informative setup. Both estimators use weighted regression with slightly different weights ($d_k$ for $\hat{T}^{(\text{MA})}$ and $d_k - 1$ for $\hat{T}^{(\text{Cos})}$).

To summarize, when the prediction model is incorrectly specified, as in the informative and quadratic setups, both $\hat{T}^{(\text{MA})}$ and $\hat{T}^{(\text{Cos})}$ are significantly more efficient than $\hat{T}^{(\text{MB})}$ because of the bias of $\hat{T}^{(\text{MB})}$, even though the bias is not large. On the opposite, if the model is correctly specified but the design weights and $Y$ are uncorrelated, as in the non-informative setup, $\hat{T}^{(\text{MB})}$ is better than $\hat{T}^{(\text{MA})}$ and $\hat{T}^{(\text{Cos})}$ in terms of MSE. In all setups, $\hat{T}^{(\text{Cos})}$ is more efficient or similar to $\hat{T}^{(\text{MA})}$.

To better understand the impact of the $R^2$ on the results, we also plot, on the $y$-axis of Figures 4.1, 4.2 and 4.3, the $\text{RMSE}_{MC}$ for 10 different values of $R^2$ on the $x$-axis: 0.1, 0.2, ..., 0.9, 0.96. In order to do that, we generate for each setup ten populations, one for each $R^2$ value. Figure 4.1 (resp. Figure 4.2 and 4.3) gives the results for Setup 1 (resp. 2 and 3) with the sample size equal to 200 (resp. 50) on the left (resp. right) column plots. On all plots, the curves correspond to the different estimators with a red curve at 100 for $\hat{T}^{(\text{Cos})}$ (since the RMSE are relative to $\hat{T}^{(\text{Cos})}$) and different colors for $\hat{T}_{HT}$, $\hat{T}^{(\text{MA})}$, $\hat{T}^{(\text{MB})}$ and $\hat{T}_{PDI}$. The plots on the top row of the figures include all the estimators while for the second row (and third row for Figures 4.1 and 4.2), $\hat{T}_{HT}$ (and $\hat{T}^{(\text{MB})}$ for Figures 4.1 and 4.2) is removed in order to zoom in and ease the comparison between $\hat{T}^{(\text{Cos})}$, $\hat{T}^{(\text{MA})}$, $\hat{T}^{(\text{MB})}$ and $\hat{T}_{PDI}$. The scale on the $y$-axis is kept fixed for the two columns (sample sizes).

As expected, $\hat{T}_{PDI}$ is by far the best estimator with the smallest MSE in all setups. In all figures, $\hat{T}_{HT}$ has a very bad relative MSE, especially when $R^2$ is high. Note that in fact the absolute MSE of $\hat{T}_{HT}$ remains stable when $R^2$ increases (results not reported), while the MSE of the other estimators improves. This result is expected because $\hat{T}_{HT}$ does not depend on the distribution of $Y|X_1, X_2$, but depends on $\mu$ and $\sigma^2$ which are constant across the populations.

Figure 4.1 (resp. 4.2) shows the evolution of $\text{RMSE}_{MC}$ with respect to $R^2$ in the informative setup (resp. quadratic setup) for sample $s_P$ of expected size 200 (left column) and 50 (right column). In these two setups, not only is $\hat{T}^{(\text{Cos})}$ better than $\hat{T}^{(\text{MB})}$ or $\hat{T}^{(\text{MA})}$, as seen in Table 4.4, but its gain compared to its competitors increases the most with $R^2$. The precision of $\hat{T}^{(\text{MA})}$ also increases, but at a slightly slower pace. The MSE of $\hat{T}^{(\text{MB})}$ worsens with $R^2$ because the prediction model differs too much from the population model in these setups. This fact implies a larger bias of $\hat{T}^{(\text{MB})}$ when $R^2$ increases. For informative and quadratic
<figure><img src="image_52.png" /><figcaption>FIGURE 4.1 – Relative MSE (in %), with the MSE of the cosmetic estimator as the baseline, versus $R^2$ in the informative setup</figcaption></figure>
<figure><img src="image_53.png" /><figcaption>FIGURE 4.2 – Relative MSE (in %), with the MSE of the cosmetic estimator as the baseline, versus $R^2$ in the quadratic setup</figcaption></figure>
<figure><img src="image_54.png" /><figcaption>FIGURE 4.3 – Relative MSE (in %), with the MSE of the cosmetic estimator as the baseline, versus $R^2$ in the non-informative setup</figcaption></figure>

setups, a smaller size reduces the difference between $RMSE_{MC}(\hat{T})$ of QR predictors.

Figure 4.3 shows the evolution of the $RMSE_{MC}$ with respect to the $R^2$ in the non-informative setup. This time, $\hat{T}^{(\text{MB})}$ does not lose precision when $R^2$ increases because the prediction model is the same as the population model. All QR predictors show an increase in precision with $R^2$, with $\hat{T}^{(\text{Cos})}$ and $\hat{T}^{(\text{MA})}$ having similar precision for all values of $R^2$. In this setup, the plots are comparable for the two sample sizes, because the model is correctly specified for all prediction models.

To sum up, if the prediction model is misspecified, the cosmetic estimator is the best choice in our setups. It has the smallest MSE amongst all QR predictors, and its precision increases faster with $R^2$ than the other estimators. The advantage of $\hat{T}^{(\text{Cos})}$ over $\hat{T}^{(\text{MA})}$ might disappear in a scenario where the probability sample size would be a smaller fraction of the population size. The estimator $\hat{T}^{(\text{MB})}$ is biased and has the largest MSE, even for smaller values of $R^2$. If the model is correctly specified, and $Y$ is not correlated to $X_1$ while the first-order inclusion probabilities are proportional to $X_1$, $\hat{T}^{(\text{MB})}$ is the best choice in terms
of MSE. However, the efficiency gain achieved by choosing $\hat{T}^{(\text{MB})}$ over $\hat{T}^{(\text{Cos})}$ in this third setup is significantly smaller than the efficiency loss observed when choosing $\hat{T}^{(\text{MB})}$ over $\hat{T}^{(\text{Cos})}$ in the first two setups. We thus recommend the choice of the Cosmetic estimator as a good compromise in all setups, followed closely by the model-assisted estimator. Similar observations can be made from real data, as shown in the next section.

# 4.6 Application to La Poste data

## 4.6.1 Data presentation

In France, more than 90% of letters sent are sorted out using automatic sorting machines. The information collected by the machines is based on pictures of the letters and form the non-probability database. La Poste has also access to a probability sample and wants to use both databases to estimate monthly totals of different types of letters using the data integration methods proposed above. Some letters such as *"lettres vertes"* (letters with a green stamp and ecologically transported) are not recognized by the sorting machines, which only take black and white pictures, meaning that the variable of interest is not available in the non-probability database. The auxiliary information associated with the letters automatically sorted is not easily linked to the probability sample. This issue is currently being investigated at La Poste. Thus, the illustration below is based on data from previous surveys made over the years at La Poste.

Data from these surveys are collected from postmen rounds and contain in particular the number of *lettres vertes*, the number of letters named *"produits 1b"* (which include different types of letters, including *lettres vertes*) and the total number of letters in the round. The idea is to mimic the situation at La Poste where the number of *lettres vertes* is available in the probability sample but not on the non-probability database whereas the number of *produits 1b* and the total number of letters are available in the non-probability database but not in the probability sample.

The goal of this section is on one hand to see if the conclusion drawn with simulated populations holds with real data and on the other hand to see if the selection method for the non-probability sample has an impact on the estimators. The population of interest consists in 11,906 rounds from historical data. In the following, we assume that the variable of interest is the number of *lettres vertes* and the explanatory variable $X_1$ (resp. $X_2$) is the total number of letters (resp. the number of *produits 1b*).

We consider three configurations. For each configuration, a non-probability sample of size 9,524 is drawn from the population; the non-probability sampling fraction is thus 80%. In
<table><caption>TABLE 4.5 – Relative bias (in % of the true value), percent relative (to $T^{(Cos)}$) variance and MSE of the different estimators for the three different non-probability samples; the expected size of the probability sample is 2,000 and the size of the non-probability sample is 9,524.</caption><thead><tr><th>S<sub>NP</sub></th><th>Monte Carlo measures</th><th>T&#x0302;<sub>HT</sub></th><th>T&#x0302;<sup>(MB)</sup></th><th>T&#x0302;<sup>(MA)</sup></th><th>T&#x0302;<sup>(Cos)</sup></th><th>T&#x0302;<sub>PDI</sub></th></tr></thead><tbody><tr><td rowspan="3">SRSWOR</td><td>RB<sub>MC</sub></td><td>0.08</td><td>8.30</td><td>0.08</td><td>0.08</td><td>0.06</td></tr><tr><td>RVar<sub>MC</sub></td><td>253.67</td><td>122.66</td><td>102.21</td><td>100.00</td><td>87.41</td></tr><tr><td>RMSE<sub>MC</sub></td><td>252.95</td><td>5,489.24</td><td>102.18</td><td>100.00</td><td>87.31</td></tr><tr><td rowspan="3">Highest Y values</td><td>RB<sub>MC</sub></td><td>0.03</td><td>5.51</td><td>0.03</td><td>0.03</td><td>0.01</td></tr><tr><td>RVar<sub>MC</sub></td><td>1,528.57</td><td>197.54</td><td>101.82</td><td>100.00</td><td>3.78</td></tr><tr><td>RMSE<sub>MC</sub></td><td>1,528.57</td><td>13,676.65</td><td>101.87</td><td>100.00</td><td>3.82</td></tr><tr><td rowspan="3">Lowest Y values</td><td>RB<sub>MC</sub></td><td>-0.02</td><td>1.23</td><td>0.13</td><td>0.13</td><td>0.41</td></tr><tr><td>RVar<sub>MC</sub></td><td>197.09</td><td>100.67</td><td>100.26</td><td>100.00</td><td>93.17</td></tr><tr><td>RMSE<sub>MC</sub></td><td>195.36</td><td>184.83</td><td>100.24</td><td>100.00</td><td>92.95</td></tr></tbody></table>

the first configuration, the non-probability sample is drawn using SRSWOR. For the second (resp. third) configuration, the non-probability sample contains the 9,524 rounds with highest (resp. lowest) values of Y. Only $X_2$ is used as explanatory variable in the prediction model along with the intercept.

We compare the same estimators as in section 4.5.2. For each of the 3 non-probability samples, we draw $L = 1,000$ probability samples $s_P$ of expected size 2,000 using Poisson sampling with probabilities proportional to $X_1$. The same Monte-Carlo measures as those in Section 4.5.2 are computed.

## 4.6.2 Results

Table 4.5 contains the simulation results for La Poste data for the three non-probability samples. For all configurations, the prediction model appears to be significantly misspecified, which causes $T^{(MB)}$ to be less efficient, both in terms of variance and MSE, than in the simulations of Section 4.5.2. Different diagnostic statistics could be computed beforehand to highlight the misspecification problem and alternative prediction models could be proposed. The assessment of such models is not in the scope of this paper.

For all configurations, the results are similar to those obtained in the informative and quadratic setups of Section 4.5.2 with $T^{(MA)}$ and $T^{(Cos)}$ having similar efficiency and being both more efficient than $T^{(MB)}$, which is the only biased estimator. It can be noted that, although the choice of $s_{NP}$ does not impact the relative efficiency of $T^{(MA)}$, it impacts the relative precision of the HT and PDI estimators.

In all configurations, $s_P$ is drawn using Poisson sampling and the variance of the data
<table><caption>TABLE 4.6 – Percent relative variance of the estimators when $s_{NP}$ contains the highest or lowest Y values, relative to their variance when $s_{NP}$ is drawn with SRSWOR.</caption><thead><tr><th>S<sub>NP</sub></th><th>T&#x0302;<sub>HT</sub></th><th>T&#x0302;<sub>NP</sub><sup>(MA)</sup></th><th>T&#x0302;<sub>NP</sub><sup>(Cos)</sup></th><th>T&#x0302;<sub>NP</sub><sup>(Ha)</sup></th><th>T&#x0302;<sub>C</sub><sup>(Ha)</sup></th></tr></thead><tbody><tr><td>Highest Y values</td><td>101.97</td><td>111.39</td><td>110.38</td><td>101.61</td><td>0.92</td></tr><tr><td>Lowest Y values</td><td>104.61</td><td>73.07</td><td>69.03</td><td>15.43</td><td>181.29</td></tr></tbody></table>

integration estimators can be simplified as follows :

$$ \mathrm{Var}(\hat{T}) = \mathrm{Var}(\hat{T}_{NP}) + \mathrm{Var}(\hat{T}_{C}^{(\mathrm{Ha})}) $$

with $\hat{T}_{NP}$ the predictor of the total $T_{NP}$ and $\hat{T}_{C}^{(\mathrm{Ha})}$ the Hájek estimator of the total $T_C$. To further understand the impact of the selection of $s_{NP}$ on the estimators, we study the variance of $\hat{T}_{NP}^{(\mathrm{MA})}$, $\hat{T}_{NP}^{(\mathrm{Cos})}$ and $\hat{T}_{C}^{(\mathrm{Ha})}$, and the variance of the Hájek estimator $\hat{T}_{NP}^{(\mathrm{Ha})} = N_{NP} \sum_{k \in s_P} d_k \delta_k y_k / \sum_{k \in s_P} d_k \delta_k$ of the total $T_{NP}$.

Table 4.6 contains the relative variance of $\hat{T}_{HT}$, $\hat{T}_{NP}^{(\mathrm{MA})}$, $\hat{T}_{NP}^{(\mathrm{Cos})}$, $\hat{T}_{C}^{(\mathrm{Ha})}$ and $\hat{T}_{NP}^{(\mathrm{Ha})}$ for the second and third configurations, when $s_{NP}$ contains the highest or lowest values of Y, relative to their variance in the first configuration (SRSWOR). As expected, the precision of the HT estimator does not depend on $s_{NP}$ and is fixed for the three configurations. In the second configuration, $\hat{T}_{C}^{(\mathrm{Ha})}$ has a much smaller variance than in the other configurations, since only the smallest values of Y are left in $U - s_{NP}$. Similarly, its variance is greater in the third configuration, when only the largest values of Y are left in $U - s_{NP}$. A similar reasoning explains the ratio of variances between configurations for $\hat{T}_{NP}^{(\mathrm{MA})}$, $\hat{T}_{NP}^{(\mathrm{Cos})}$ and $\hat{T}_{NP}^{(\mathrm{Ha})}$. The important point is that both $\hat{T}_{NP}^{(\mathrm{MA})}$ and $\hat{T}_{NP}^{(\mathrm{Cos})}$ are not very sensitive to $s_{NP}$ in terms of variance.

We could argue that the non-probability samples selected in the second and third configurations are not realistic. Yet, the results observed for more realistic non-probability samples, selected using Poisson sampling proportional to $X_1$ or Y, were not significantly different from the results observed when $s_{NP}$ is drawn using SRSWOR and are not reported. Similarly, using both auxiliary variables in the prediction model only impacts the bias of $\hat{T}^{(\mathrm{MB})}$, and the results are again not reported.

To further improve the precision, other estimators of $T_C$ may be used. As mentioned before, the data integration methods proposed in this paper can be used when the overlap between the probability and the non-probability samples is non-empty and, ideally, large. In this context, we recommend, based on our empirical results, the choice of the cosmetic or the model-assisted estimators.
## 4.7 Conclusion

Most of the literature on data integration in finite population tackles the problem of unobserved study variable in the probability sample. In this paper, we have proposed to fill the gap and considered the problem of unobserved study variable in the non-probability sample, assuming that it is observed in the probability sample and that auxiliary information is available in both samples. We have defined a general class of prediction estimators, based on the already known QR class, which includes the model-assisted, model-based and cosmetic estimators, and studied theoretically their bias and variance properties. We have also derived a variance estimator and compared the three types of estimators with the usual Horvitz-Thompson estimator in different simulation setups, both in terms of bias and MSE, and concluded that the cosmetic estimator is a good compromise in general.

The main conclusion of our experiments is that significant efficiency gains can be achieved by leveraging a big non-probability database that contains auxiliary information associated with the main study variables. For large domains, the efficiency gains obtained from using model-assisted estimators, including the cosmetic estimator, may be sufficient to obtain high-quality estimates of the population parameters of interest. For smaller domains, these estimators may not achieve precision targets. However, they could be used as direct estimates in a small area estimation model, such as the well-known Fay-Herriot area level model. This model requires area level auxiliary information. The big non-probability database would be a natural candidate for providing the auxiliary information required for producing small area estimates. Small area estimation methods often yield significant precision gains over direct estimators at the expense of introducing model assumptions.

# Appendix

## Proof of Proposition 4.2.1

We recall that $\hat{T}_{DI} = \sum_{k \in U} \delta_k y_k + \sum_{k \in s_P} (1-\delta_k) d_k y_k$ and $\hat{T}_{HT} = \sum_{k \in s_P} d_k y_k = \sum_{k \in s_P} \delta_k d_k y_k + \sum_{k \in s_P} (1-\delta_k) d_k y_k$. Thus, we have :

$$ \mathrm{Var}(\hat{T}_{HT}) - \mathrm{Var}(\hat{T}_{DI}) = \mathrm{Var} \left( \sum_{k \in s_P} \delta_k d_k y_k \right) + 2 \mathrm{Cov} \left( \sum_{k \in s_P} \delta_k d_k y_k, \sum_{k \in s_P} (1-\delta_k) d_k y_k \right). $$
(i) For Poisson sampling, we have :

$$ \mathrm{Cov} \left( \sum_{k \in s_P} \delta_k d_k y_k, \sum_{k \in s_P} (1 - \delta_k) d_k y_k \right) = \sum_{k \in U} \delta_k (1 - \delta_k) (d_k - 1) y_k^2 = 0 $$

and

$$ \mathrm{Var}(\hat{T}_{HT}) - \mathrm{Var}(\hat{T}_{DI}) = \mathrm{Var} \left( \sum_{k \in s_P} \delta_k d_k y_k \right) = \sum_{k \in U} \delta_k (d_k - 1) y_k^2 \geq 0, $$

which proves the first part of the proposition.

(ii) For simple random sampling without replacement, let $\bar{Y}_U = \sum_{k \in U} y_k/N$, $\bar{Y}_{NP} = \sum_{k \in U} \delta_k y_k/N_{NP}$, $S_{Y,NP}^2 = \sum_{k \in U} \delta_k (y_k - \bar{Y}_{NP})^2/(N_{NP}-1)$ and $CV_{NP}^2 = S_{Y,NP}^2/\bar{Y}_{NP}^2$. Using some simple calculus, we have :

$$ \mathrm{Var} \left( \sum_{k \in s_P} \delta_k d_k y_k \right) = \frac{N}{n} \frac{N-n}{N(N-1)} \left( N(N_{NP}-1)S_{Y,NP}^2 + N_{NP}\bar{Y}_{NP}^2(N-N_{NP}) \right), $$

$$ \mathrm{Cov} \left( \sum_{k \in s_P} \delta_k d_k y_k, \sum_{k \in s_P} (1 - \delta_k) d_k y_k \right) = - \frac{N}{n} \frac{N-n}{N(N-1)} N_{NP} \bar{Y}_{NP} (N \bar{Y}_U - N_{NP} \bar{Y}_{NP}), $$

and thus

$$ \mathrm{Var}(\hat{T}_{HT}) - \mathrm{Var}(\hat{T}_{DI}) = \frac{N}{n} \frac{N-n}{N(N-1)} \left( N(N_{NP}-1)S_{Y,NP}^2 + N_{NP}\bar{Y}_{NP} \left( (N + N_{NP})\bar{Y}_{NP} - 2N\bar{Y}_U \right) \right). $$

We conclude that $\mathrm{Var}(\hat{T}_{HT})$ is larger than or equal to $\mathrm{Var}(\hat{T}_{DI})$ if and only if

$$ N(N_{NP}-1)S_{Y,NP}^2 + N_{NP}\bar{Y}_{NP} \left( (N + N_{NP})\bar{Y}_{NP} - 2N\bar{Y}_U \right) \geq 0, $$

which is equivalent to :

$$ CV_{NP}^2 \geq - \frac{N_{NP}}{N_{NP}-1} \left( 1 + \frac{N_{NP}}{N} - 2 \frac{\bar{Y}_U}{\bar{Y}_{NP}} \right), $$

and proves the second part of the proposition.

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

### Proof of Proposition 4.2.2

We have :

$$ \mathrm{Var}(\hat{T}_{HT}) = \mathrm{Var} \left( \sum_{k \in s_P} d_k y_k \right) = N^2(1-f) \frac{S_{Y,U}^2}{n}, $$
$$AVar(\hat{T}_{PDI}) = Var \left( \sum_{k \in s_P} (1 - \delta_k) d_k (y_k - \bar{Y}_C) \right) = Var \left( \sum_{k \in s_P} d_k \tilde{y}_k \right) = N^2(1-f) \frac{S_{Y,U}^2}{n}$$

where

$$S_{Y,U}^2 = \frac{1}{N-1} \sum_{k \in U} (y_k - \bar{Y}_U)^2,$$

$$\tilde{y}_k = (1 - \delta_k)(y_k - \bar{Y}_C), k \in U$$

$$S_{\tilde{Y},U}^2 = \frac{1}{N-1} \sum_{k \in U} (\tilde{y}_k - \bar{\tilde{Y}}_C)^2 = \frac{1}{N-1} \sum_{k \in U} \tilde{y}_k^2.$$

Using some basic but tedious calculus, we obtain :

$$\begin{align*}
\mathrm{Var}(\hat{T}_{HT}) - \mathrm{AVar}(\hat{T}_{PDI}) &= N^2(1-f) \frac{S_{Y,U}^2 - S_{\tilde{Y},U}^2}{n} \\
&= N^2(1-f) \frac{1}{n} \frac{1}{N-1} \left( \sum_{k \in U} \delta_k (y_k - \bar{Y}_U)^2 + (N - N_{NP})(\bar{Y}_C - \bar{Y}_U)^2 \right) \\
&= N^2(1-f) \frac{1}{n} \frac{1}{N-1} \left( S_{Y,NP}^2 (N_{NP} - 1) + N_{NP} \frac{N}{N - N_{NP}} (\bar{Y}_{NP} - \bar{Y}_U)^2 \right).
\end{align*}$$

□

### Proof of Proposition 4.3.1

Let $\mathbf{R}_{s_P} = \mathrm{diag}(r_k \delta_k)_{k \in s_P}$, $\mathbf{X}_{s_P} = (\mathbf{x}_k^\top)_{k \in s_P}$, $\mathbf{y}_{s_P} = (y_k)_{k \in s_P}$ and $\mathbf{Q}_{xs_P}^\top = \mathbf{X}_{s_P}^\top \mathrm{diag}(q_k \delta_k)_{k \in s_P}$. Then $\hat{\beta} = (\mathbf{Q}_{xs_P}^\top \mathbf{X}_{s_P})^{-1} \mathbf{Q}_{xs_P}^\top \mathbf{y}_{s_P}$. We can write the sum $\sum_{k \in s_P} r_k \delta_k (y_k - \hat{y}_k)$ in a matrix form as follows :

$$\sum_{k \in s_P} r_k \delta_k (y_k - \hat{y}_k) = \mathbf{1}_{s_P}^\top \mathbf{R}_{s_P} (\mathbf{y}_{s_P} - \mathbf{X}_{s_P} \hat{\beta}),$$

where $\mathbf{1}_{s_P}$ is a vector of ones with dimension the size of $s_P$. If the condition $\boldsymbol{\mu}^\top \mathbf{x}_k q_k - r_k = 0$ is fulfilled for all $k \in s_{NP}$, then $\delta_k (\boldsymbol{\mu}^\top \mathbf{x}_k q_k - r_k) = 0$ for all $k \in s_P$ and as a consequence, $\boldsymbol{\mu}^\top \mathbf{Q}_{xs_P}^\top = \mathbf{1}_{s_P}^\top \mathbf{R}_{s_P}$. We get then $\mathbf{1}_{s_P}^\top \mathbf{R}_{s_P} (\mathbf{y}_{s_P} - \mathbf{X}_{s_P} \hat{\beta}) = 0$.

### Proof of Proposition 4.3.2

We have

$$\begin{align*}
\hat{T}_{NP}^{(QR)} - \hat{T}_{NP}^{(Q\pi)} &= \sum_{k \in s_P} (r_k - d_k) \delta_k (y_k - \mathbf{x}_k^\top \hat{\beta}) \\
&= -\boldsymbol{\lambda}^\top \sum_{k \in s_P} q_k \delta_k \mathbf{x}_k (y_k - \mathbf{x}_k^\top \hat{\beta}) = 0.
\end{align*}$$
# 4.8 Supplementary material

This section details the results presented in section 4.4.2.

## 4.8.1 Assumptions

To the assumptions already mentioned in the paper :

$$||\hat{\beta} - \tilde{\beta}||^2 = o_p(1)$$

and

$$\left( \sum_{k \in s_p} d_k(1 - \delta_k)/N \right)^{-1} = O_p(1),$$

we add the following assumptions to ensure the convergence of the **HT** estimator as discussed in section 1.1.3:

1. We assume that $\lim_{N \to \infty} \frac{n}{N} \in (0, 1)$.

2. We assume the following assumptions on the inclusion probabilities :

$$\min_{k \in U} \pi_k \ge \lambda > 0 \text{ and } \limsup_{N \to \infty} \max_{k \ne l \in U} |\pi_{kl} - \pi_k \pi_l| < C_1 < \infty;$$

3. We assume that there is a positive constant $C_2$ such that, for all $N$, $\frac{1}{N} \sum_{k \in U} y_k^2 < C_2$.

4. We assume that there is a positive constant $C_3$ such that, for all $k \in U$, $||x_k||^2 < C_3$.

## 4.8.2 Development of the results presented in section 4.4.2

In this section, we show that we can write the standardized total error as

$$\frac{1}{N}(\hat{T}^{(Q\pi)} - T) = \frac{1}{N}\left(\sum_{k \in s_p} d_k(E_k + e_k) - \sum_{k \in U} (E_k + e_k)\right) + R_1 + R_2.$$

By definition, the standardized total error is given by :

$$\frac{1}{N}(\hat{T}^{(Q\pi)} - T) = \frac{1}{N}(\hat{T}_{NP}^{(Q\pi)} - T_{NP}) + \frac{1}{N}(\hat{T}_{C}^{(Ha)} - T_{C}).$$

For any variable $z_k$, $k \in U$, let $T_z = \sum_{k \in U} z_k$ be the total of $z$ over $U$ and $\hat{t}_z = \sum_{k \in U} d_k z_k$ be the **HT** estimator of $T_z$.
Let us first look at $\frac{1}{N}(\hat{T}_{NP}^{(Q\pi)} - T_{NP})$. The QR estimator $\hat{T}_{NP}^{(Q\pi)}$ can be rewritten

$$ \hat{T}_{NP}^{(Q\pi)} = \sum_{k \in s_p} d_k \delta_k y_k - \left( \sum_{k \in s_p} d_k \delta_k \boldsymbol{x}_k^t - \sum_{k \in U} \delta_k \boldsymbol{x}_k^t \right) \hat{\beta} = \hat{t}_{\delta y} - (\hat{t}_{\delta x} - T_{\delta x})^t \hat{\beta}. $$

Then

$$ \begin{aligned} \frac{1}{N}(\hat{T}_{NP}^{(Q\pi)} - T_{NP}) &= \frac{1}{N}(\hat{t}_{\delta y} - T_{NP}) - \frac{1}{N}(\hat{t}_{\delta x} - T_{\delta x})^t (\hat{\beta} - \tilde{\beta}) - \frac{1}{N}(\hat{t}_{\delta x} - T_{\delta x})^t \tilde{\beta} \\ &= \frac{1}{N}(\hat{t}_{\delta y} - \hat{t}_{\delta x}^t \tilde{\beta}) - \frac{1}{N}(T_{NP} - T_{\delta x}^t \tilde{\beta}) + R_1 \\ &= \frac{1}{N} \left( \sum_{k \in s_p} d_k \delta_k (y_k - \boldsymbol{x}_k^t \tilde{\beta}) - \sum_{k \in U} \delta_k (y_k - \boldsymbol{x}_k^t \tilde{\beta}) \right) + R_1 \\ &= \frac{1}{N} (\hat{t}_E - T_E) + R_1 \end{aligned} $$

with $R_1 = -\frac{1}{N}(\hat{t}_{\delta x} - T_{\delta x})^t (\hat{\beta} - \tilde{\beta})$ and for all $k \in U$ $E_k = \delta_k(y_k - \boldsymbol{x}_k^t \tilde{\beta})$. Under the assumed conditions 1, 2 and 4, the **HT** estimator $\hat{t}_{\delta x} = \sum_{k \in U} d_k \delta_k \boldsymbol{x}_k$ satisfy $\frac{1}{N}(\hat{t}_{\delta x} - T_{\delta x}) = O_p(1/\sqrt{n})$. Since we assumed that $||\hat{\beta} - \tilde{\beta}||^2 = o_p(1)$, then $R_1 = o_p(1/\sqrt{n})$.

Let us now look at $\frac{1}{N}(\hat{T}_C^{(Ha)} - T_C)$:

$$ \frac{1}{N}(\hat{T}_C^{(Ha)} - T_C) = \frac{1}{N}T_{1-\delta} \left( \frac{\hat{t}_{(1-\delta)y}}{\hat{t}_{1-\delta}} - \frac{T_{(1-\delta)y}}{T_{1-\delta}} \right). $$

The difference of ratio can be rewritten as:

$$ \begin{aligned} \frac{\hat{t}_{(1-\delta)y}}{\hat{t}_{1-\delta}} - \frac{T_{(1-\delta)y}}{T_{1-\delta}} &= \frac{1}{\hat{t}_{1-\delta}} \left( \hat{t}_{(1-\delta)y} - \frac{T_{(1-\delta)y}}{T_{1-\delta}} \hat{t}_{1-\delta} \right) = \frac{1}{\hat{t}_{1-\delta}} \hat{t}_e = \frac{1}{\hat{t}_{1-\delta}} (\hat{t}_e - T_e) \\ &= \left( \frac{1}{\hat{t}_{1-\delta}} - \frac{1}{T_{1-\delta}} \right) (\hat{t}_e - T_e) + \frac{1}{T_{1-\delta}} (\hat{t}_e - T_e) \\ &= \frac{N}{T_{1-\delta}} R_2 + \frac{N}{T_{1-\delta}} \frac{1}{N} (\hat{t}_e - T_e) \end{aligned} $$

with for all $k \in U$, $e_k = (1 - \delta_k) \left( y_k - \frac{T_{(1-\delta)y}}{T_{1-\delta}} \right) = (1 - \delta_k) \left( y_k - \frac{T_{(1-\delta)y}}{N - N_{NP}} \right)$, $T_e = \sum_{k \in U} e_k = 0$ and $R_2 = \left( \frac{(\hat{t}_{1-\delta} - T_{1-\delta})/N}{\hat{t}_{1-\delta}/N} \right) \frac{1}{N} (\hat{t}_e - T_e)$. Since $(\hat{t}_e - T_e)/N$ and $(\hat{t}_{1-\delta} - T_{1-\delta})/N$ are both $O_p(1/\sqrt{n})$ under the assumed conditions, $R_2 = O_p(1/n)$ if $\frac{1}{\hat{t}_{1-\delta}/N} = O_p(1)$.
Finally we get

$$ \begin{aligned} \frac{1}{N}(\hat{T}^{(Q\pi)} - T) &= \frac{1}{N}(\hat{T}_{NP}^{(Q\pi)} - T_{NP}) + \frac{1}{N}(\hat{T}_{C}^{(Ha)} - T_{C}) \\ &= \frac{1}{N}(\hat{t}_{E} - T_{E}) + \frac{1}{N}(\hat{t}_{e} - T_{e}) + R_{1} + R_{2} \\ &= \frac{1}{N}(\hat{t}_{E} - T_{E}) + \frac{1}{N}(\hat{t}_{e} - T_{e}) + o_{p}(1/\sqrt{n}). \end{aligned} $$

□
# Conclusions et perspectives

Dans cette thèse, j’ai cherché à améliorer la précision en terme de variance des estimateurs du trafic postal, d’une part en améliorant la méthode de partage des poids utilisée à La Poste et d’autre part en utilisant les larges sources de données auxiliaires dont dispose La Poste.

Je me suis d’abord penchée sur l’étude des estimateurs du total utilisés par La Poste afin de comprendre la perte de précision observée depuis quelques années et j’ai montré que l’augmentation de la variance était majoritairement due à l’utilisation de poids MGPP non optimaux. Ce travail a mené à des propositions d’amélioration de la pondération des estimateurs actuels via la modélisation du nombre d’adresses par case en fonction des données auxiliaires. Cette modélisation est à l’étude à La Poste et son effet sur la précision des estimateurs reste à confirmer.

Les travaux présentés dans les chapitres 2 et 3 ont mis en valeur deux avantages liés à l’utilisation d’une MGPP quand les liens sont de type TpU : d’une part, l’existence de poids optimaux et d’autre part, dans le cas TpU-TpU, la possibilité de réduire le nombre de liens à observer via une MGPP double. Il est intéressant de se demander si ces avantages peuvent être étendus à d’autres situations que les cas TpU et TpU-TpU.

Dans mes travaux, je ne m’intéresse pas à l’existence de poids optimaux indépendants de la variable d’intérêt pour les cas non TpU car l’existence de poids optimaux dans ce cas est complexe. Deville et Lavallée montrent que, pour les plans de Poisson et SASSR, ces poids n’existent pas. Cependant, il est possible d’obtenir des poids qui minimisent la variance des estimateurs MGPP dans des cas non TpU en modifiant la population d’intérêt de manière à obtenir des liens de type TpU.

Une façon de faire est d’utiliser une population de clusters, comme indiqué dans la figure 4.4. Dans cette figure, on regroupe les individus de *U* reliés à au moins un individu commun de $U_F$. Par exemple, $U_1$ et $U_2$ sont tous deux reliés à $F_3$, ils sont donc regroupés. On crée ensuite une population $U'$ composée des clusters ainsi formés. L’individu $U'_1$ est donc le cluster composé de $U_1$ et $U_2$, et dont la valeur de la variable d’intérêt est $y'_{U'_1} = y_{U_1} + y_{U_2}$.
<figure><img src="image_56.png" /><figcaption>FIGURE 4.4 – Utilisation d'une population $U'$ composée de clusters.</figcaption></figure>

Le total $T_y$ peut être réécrit comme le total $T_{y'}$ de $y'$ sur $U'$. Les liens entre $U_F$ et $U'$ étant de type TpU, comme indiqué dans la figure 4.5, la propriété d'optimalité décrite dans le chapitre 2 de ce manuscrit s'applique si le plan de sondage satisfait la $\Delta$-propriété pour $U'$. Cependant, cette méthode augmente de manière significative le nombre de liens à observer.

<figure><img src="image_57.png" /><figcaption>FIGURE 4.5 – Liens entre $U_F$ et $U'$.</figcaption></figure>

On a vu dans le chapitre 2 que la MGPP double avec double standardisation permettait de réduire le nombre de liens à identifier. Ce nombre peut être encore diminué en considérant plus de deux populations intermédiaires, telles que les liens entre chaque population soit de type TpU. L'estimateur MGPP du total $T_y$ s'obtient alors similairement à l'estimateur MGPP double, en standardisant les liens entre chaque populations. L'existence de poids
optimaux pour ces liens, ainsi que le gain (ou perte) en nombre de liens observés selon le nombre de populations intermédiaires, est une des perspective de la thèse que j’ai commencé à étudier.

Il serait aussi possible de considérer les liens qui n’ont pas pu être identifiés comme non-répondants. Les techniques de traitement de la non-réponse s’appliquent alors, comme détaillé dans Xu and Lavallée (2009). Une idée pourrait être de s’inspirer des graphes aléatoires en attribuant à chaque couple $(i, k) \in U_F \times U$ une probabilité d’être liés.

Dans un deuxième temps, je me suis penchée sur l’intégration de bases de données massives pour améliorer les estimations postales, dans le cas où les variables d’intérêt ne sont pas disponibles dans les bases non probabilistes. J’ai proposé des améliorations théoriques qui dépendent de la capacité à identifier de manière unique un individu dans les différentes bases de données et de la qualité des variables auxiliaires. Cependant, La Poste fait face à des challenges industriels liés à l’accessibilité et à la qualité des bases de données massives, rendant l’application des mesures proposées délicate.

Un des problèmes vient de la construction du modèle. La méthode proposée dans ma thèse consiste à utiliser un modèle linéaire pour prédire la valeur des variables d’intérêt pour les individus contenus dans la base de données massive. Or, l’étude des corrélations des données disponibles à La Poste montre que celles-ci sont trop faibles pour justifier l’utilisation d’un tel modèle.

Ce problème pourrait être contourné en utilisant un modèle non paramétrique qui serait moins sensible à la qualité (et quantité) des données auxiliaires. L’idée est d’obtenir un modèle et de calculer des prédictions, qui puissent être utilisées pour toutes les variables d’intérêt. Dans ma thèse, j’ai testé, sur données simulées, le comportement des estimateurs QR avec des modèles de type forêts aléatoires. Ce travail est prometteur mais à approfondir.

De son côté, La Poste explore différentes méthodes de calage afin d’améliorer la précision de ses estimateurs du total, mais ces méthodes n’ont pas été étudiées dans le cadre de ma thèse.
Bibliographie

Ardilly, P. (2006). *Les techniques de sondage*. Editions TECHNIP.

Bankier, M. D. (1986). Estimators based on several stratified samples with applications to multiple frame surveys. *Journal of the American Statistical Association*, 81(396):1074-1079.

Beaumont, J.-F. (2020). Are probability surveys bound to disappear for the production of official statistics? *Survey Methodology*, 46(1):1-29.

Breidt, F.-J. and Opsomer, J.-D. (2000). Local polynomial regression estimators in survey sampling. *The Annals of Statistics*, 28(4):1023-1053.

Brewer, K. (1999). Cosmetic calibration with unequal probability sampling. *Survey Methodology*, 25(2):205-212.

Cardot, H., Goga, C., and Lardin, P. (2013). Uniform convergence and asymptotic confidence bands for model-assisted estimators of the mean of sampled functional data. *Electronic Journal of Statistics*, 7:562-596.

De Vitiis, C., Falorsi, S., Inglese, F., Masi, A., Pannuzi, N., and Russo, M. (2014). A methodological approach based on indirect sampling to survey the homeless population. *Rivista di statisticaufficiale*, 1(2):9-30.

Dessertaine, A. and Fluteaux, L. (2004). Utilisation de la méthode généralisée du partage des poids dans le cadre des estimations de flux de courrier à la poste. *In : Ardilly, Pascal, (ed.) Echantillonnage et méthodes d'enquêtes, Science Sup.*

Deville, J. and Maumy-Bertrand, M. (2006). Extension of the indirect sampling method and its application to tourism. *Survey Methodology*, 32(2):177.

Deville, J.-C. (1999). Simultaneous calibration of several surveys. In *Proceedings of Statistics Canada Symposium 99 of Statistics Canada*, pages 207-212.
Deville, J.-C. and Lavallée, P. (2006). Indirect sampling : The foundations of the generalized weight share method. *Survey Methodology*, Vol. 32(2) :165–176.

Deville, J.-C. and Lavallée, P. (2006). Sondage indirect : les fondements de la méthode généralisée du partage des poids. *Techniques d'enquête*, 32(2) :185.

Deville, J.-C. and Särndal, C.-E. (1992). Calibration estimators in survey sampling. *Journal of the American Statistical Association*, 87 :376–382.

Ernst, L. R. (1986). *Weighting issues for longitudinal household and family estimates*. US Bureau of the Census.

Falorsi, P. D., Righi, P., and Lavallée, P. (2019). Cost optimal sampling for the integrated observation of different populations. *Survey methodology*, 45(3) :485–511.

Goga, C., Deville, J.-C., and Ruiz-Gazen, A. (2009). Use of functionals in linearization and composite estimation with application to two-sample data. *Biometrika*, 96(3) :691–709.

Hájek, J. (1971). Comment on a paper by D. Basu. In *Foundations of Statistical Inference* (eds. V.P. Godambe and D.A. Sprott), page 236. Toronto : Holt, Rinehart and Winston.

Hájek, J. (1981). *Sampling from a finite population*. Statistics : Textbooks and Monographs. Marcel Dekker, New York.

Haziza, D. and Beaumont, J.-F. (2017). Construction of weights in surveys : A review. *Statistical Science*, Vol. 32(2) :206–226.

Horvitz, D. and Thompson, D. (1952). A generalization of sampling without replacement from a finite universe. *Journal of the American Statistical Association*, 47 :663–685.

Isaki, C.-T. and Fuller, W.-A. (1982). Survey design under the regression superpopulation model. *J. Amer. Statist. Assoc.*, 77 :49–61.

Kalton, G. and Brick, J. M. (1995). Weighting schemes for household panel surveys. *Survey Methodology*, 21(2) :33–34.

Kiesl, H. (2016). Indirect sampling : a review of theory and recent applications. *AStA Wirtschafts-und Sozialstatistisches Archiv*, 10(4) :289–303.

Kim, J. K. (2022). A gentle introduction to data integration in survey sampling. *The survey statistician*, 85 :19–29.
Kim, J. K., Park, S., Chen, Y., and Wu, C. (2021). Combining non-probability and probability survey samples through mass imputation. *Journal of the Royal Statistical Society Series A : Statistics in Society*, 184(3):941–963.

Kim, J.-K. and Tam, S.-M. (2021). Data integration by combining big data and survey sample data for finite population inference. *International Statistical Review*, 89(2):382–401.

Kokic, P. and Bell, P. (1994). Optimal winsorizing cutoffs for a stratified finite population estimator. *JOURNAL OF OFFICIAL STATISTICS-STOCKHOLM*-, 10:419–419.

Lardin-Puech, P. (2014). Estimation du trafic de courrier distribué en france métropolitaine par sondage indirect. In *the proceedings of the 8ème colloque francophone sur les sondages, Dijon*. http://paperssondages14.sfds.asso.fr/submission_100.pdf.

Lavallée, P. (2007). *Indirect sampling*. Springer Science & Business Media.

Medous, E., Goga, C., Ruiz-Gazen, A., Beaumont, J.-F., Dessertaine, A., and Puech, P. (2023a). Many-to-one indirect sampling with application to the french postal traffic estimation. *The Annals of Applied Statistics*, 17(1):838–859.

Medous, E., Goga, C., Ruiz-Gazen, A., Beaumont, J.-F., Dessertaine, A., and Puech, P. (To be published, 2023b). Qr prediction for statistical data integration. *Survey Methodology*.

Narain, R. (1951). On sampling without replacement with varying probabilities. *Journal of the Indian Society of Agricultural Statistics*, 3(2):169–175.

Neyman, J. (1934). On the two different aspects of the representative method : the method of stratified sampling and the method of purposive selection. *Journal of the Royal Statistical Society*, 97:558–606.

Rao, J. (2005). Interplay between sample survey theory and practice : An appraisal. *Survey Methodology*, 31(2):117.

Rao, J. (2021). On making valid inferences by integrating data from surveys and other sources. *Sankhya B*, 83(1):242–272.

Rendtel, U. and Harms, T. (2009). Weighting and calibration for household panels. *Methodology of longitudinal surveys*, pages 265–286.

Särndal, C.-E. (1980). On the π-inverse weighting best linear unbiased weighting in probability sampling. *Biometrika*, 67:639–650.
Särndal, C.-E., Swensson, B., and Wretman, J. (1992). *Model assisted survey sampling*. Springer Series in Statistics. Springer-Verlag, New York.

Särndal, C.-E., Swensson, B., and Wretman, J. H. (1989). The weighted residual technique for estimating the variance of the general regression estimator of the finite population total. *Biometrika*, 76(3):527-537.

Särndal, C.-E. and Wright, R. (1984). Cosmetic form of estimators in survey sampling. *Scandinavian J. of Statistics*, 11:146-156.

Sherman, J. and Morrison, W. J. (1950). Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix. *The Annals of Mathematical Statistics*, 21(1):124-127.

Tillé, Y. (2019). *Théorie des sondages : échantillonnage et estimation en populations finies*. Dunod, Paris.

Wright, R. (1983). Finite population sampling with multivariate auxiliary information. *Journal of the American Statistical Association*, 78(384):879-884.

Wu, C. (2022). Statistical inference with non-probability survey samples. *Survey Methodology*, 48(2):283-311.

Xu, X. and Lavallée, P. (2009). Treatments for link nonresponse in indirect sampling. *Survey Methodology*, 35(2):153-164.

Yang, S. and Kim, J. K. (2020). Integration of survey data and big observational data for finite population inference using mass imputation. *Japanese Journal of Statistics and Data Science*, 3:625-650.
**Résumé** La Poste française a mis en place une étude par sondage pour effectuer un suivi de la distribution des objets postaux en France Métropolitaine. Le but des travaux réalisés lors de cette thèse est d'améliorer la précision des estimateurs du trafic total postal. La Poste souhaite observer des tournées de facteurs, mais ne peut pas les échantillonner directement. Elle a donc mis en place un tirage indirect et utilise une Méthode Généralisée de Partage des Poids (MGPP) double, inspirée de la MGPP simple proposée par Deville and Lavallée (2006), pour calculer les estimations du trafic postal. Dans un premier temps, nous étudions dans le cadre de La Poste l'existence d'une MGPP simple optimale, au sens que la variance des estimateurs est minimale, et la comparons aux estimateurs postaux actuels. La MGPP double, telle qu'utilisée à La Poste, est plus facile à mettre en place que la MGPP simple, mais est moins précise que la MGPP simple optimale. Nous nous penchons donc sur la question d'une MGPP double optimale dans le cas de La Poste et proposons des alternatives basées sur l'utilisation de données auxiliaires. Dans un deuxième temps, nous nous intéressons à l'Intégration Statistique de Données (ISD), qui consiste à utiliser un échantillon non-probabiliste, généralement une base de données massives, et un échantillon probabiliste pour construire des estimateurs de totaux. La littérature sur l'ISD traite majoritairement de cas où la variable d'intérêt est disponible dans la base de données massives. Les données massives de La Poste, intitulées *Traitement Automatisé de l'Enveloppe (TAE)*, sont récoltées sur les objets triés par machine et ne contiennent pas d'information sur certaines variables d'intérêt. Des méthodologies adaptées au cas particulier de La Poste doivent donc être développées. Nous proposons une amélioration des estimateurs présentés dans Kim and Tam (2021) en utilisant une famille de prédicteurs, dits prédicteurs QR, proposée par Wright (1983), pour prédire le total des variables d'intérêt sur TAE. Nous utilisons un estimateur de Hájek pour estimer le total des variables d'intérêt sur le complémentaire de TAE. Nous étudions les propriétés théoriques des estimateurs QR et préconisons l'utilisation de l'estimateur cosmétique qui est un cas particulier d'estimateur QR. Les résultats présentés dans cette thèse sont illustrés par des études Monte-Carlo basées sur des données simulées et sur des données historiques postales, afin d'évaluer le gain de précision potentiel apporté par l'utilisation des différents estimateurs que nous proposons.

**Mots Clés** : données non probabilistes, estimateur QR, Intégration Statistique de Données, Méthode Généralisée de Partage des Poids, minimisation de la variance, poids de sondage optimaux, population finie, propriétés asymptotiques, sondage.

**Abstract** The French Postal Service (La Poste) has set up a sample survey to monitor mail distribution in Metropolitan France. The aim of the work carried out in this thesis is to improve the accuracy of estimators of total postal traffic. La Poste wants to observe postman rounds, but cannot sample them directly. It has therefore set up an indirect sampling and uses a double Generalized Weight Share Method (GWSM), inspired by the simple GWSM proposed by Deville and Lavallée (2006), to compute the estimates of the postal traffic. First, we study in the situation of La Poste the existence of a simple optimal GWSM, in the sense that the variance of the estimators is minimal, and compare it with current postal estimators. The double GWSM, as used at La Poste, is easier to implement than the simple GWSM, but is less accurate than the optimal simple GWSM. We therefore address the question of an optimal double GWSM in La Poste case and propose alternatives based on the use of auxiliary data. Secondly, we focus on Statistical Data Integration (SDI), which consists in using a non-probability sample, mainly a massive database, and a probability sample to compute estimators of totals. The literature on SDI deals mainly with cases where the variable of interest is available in the massive database used. La Poste massive database, entitled *Traitement Automatisé de l'Enveloppe (TAE)*, is collected on mails sorted by machine and does not contain information for several variables of interest. Methodologies adapted to the specific case of La Poste must therefore be developed. We propose an amelioration of the estimators presented in Kim and Tam (2021) by using a family of predictors, known as QR predictors, proposed by Wright (1983), to predict the total of the variables of interest in TAE. We propose a Hájek estimator to estimate the total of the variables of interest in the complementary set of TAE. We study the theoretical properties of QR estimators and recommend the use of the cosmetic estimator, which is a special case of QR estimator. The results presented in this thesis are illustrated by Monte-Carlo studies based on simulated data and historical postal data, in order to assess the potential gain in accuracy brought by the use of the different proposed estimators.

**Key Words** : asymptotic properties, finite population, Generalized Weight Sharing Method, non-probabilistic data, optimal sampling weights, QR estimator, sampling method, Statistical Data Integration, variance minimisation.